CS289_project.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker\n",
    "from sklearn import svm,neighbors\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.model_selection import cross_validate,KFold\n",
    "\n",
    "from DataPrep import DATA\n",
    "from ATLAS import ATLAS\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.decomposition import FastICA\n",
    "from sklearn.cross_decomposition import CCA\n",
    "\n",
    "Atlas = ATLAS()\n",
    "\n",
    "Data = DATA()\n",
    "Data.Train_Test(0.8,12345)\n",
    "Data.Add_MRI(\"brain\")\n",
    "Data.Split_Data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Do not run\n",
    "def Train_Test_SVM_KCross(C_space,train,valid,train_index,valid_index):\n",
    "    train_score=[]\n",
    "    valid_score=[]\n",
    "    for C in C_space:\n",
    "        SVM = svm.SVC(kernel='linear',C=C)\n",
    "        SVM.fit(train,Data.labels_train[train_index])\n",
    "        train_score.append(SVM.score(train,Data.labels_train[train_index]))\n",
    "        valid_score.append(SVM.score(valid,Data.labels_train[valid_index]))\n",
    "        \n",
    "    return(train_score,valid_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Do not run \n",
    "\n",
    "folds =4\n",
    "\n",
    "kf = KFold(n_splits=folds,shuffle=True, random_state = 3)\n",
    "C_space = np.logspace(-7,0,25)\n",
    "comp=50\n",
    "\n",
    "Train_temp = np.zeros((folds,C_space.shape[0]))\n",
    "Valid_temp = np.zeros((folds,C_space.shape[0]))\n",
    "i=0\n",
    "idxs = [5,13,47]\n",
    "for train_index, valid_index in kf.split(Data.features_train):\n",
    "    pca = PCA(n_components=comp, whiten=False)\n",
    "    pca.fit(Data.features_train[train_index,5:])\n",
    "    X_train_pca = pca.transform(Data.features_train[train_index,5:])\n",
    "    X_valid_pca = pca.transform(Data.features_train[valid_index,5:])\n",
    "\n",
    "    train_pca = np.hstack((X_train_pca,Data.features_train[train_index,0:5]))\n",
    "    valid_pca = np.hstack((X_valid_pca,Data.features_train[valid_index,0:5]))\n",
    "\n",
    "    if len(idxs)==1:\n",
    "        train = np.reshape(train_pca[:,idxs],(-1,1))\n",
    "        valid = np.reshape(valid_pca[:,idxs],(-1,1))\n",
    "    else:\n",
    "        train = train_pca[:,idxs]\n",
    "        valid = valid_pca[:,idxs]\n",
    "    SVM_PCA = Train_Test_SVM_KCross(C_space,train,valid,train_index,valid_index)\n",
    "    Train_temp[i,:] = SVM_PCA[0]\n",
    "    Valid_temp[i,:] = SVM_PCA[1]\n",
    "    i+=1\n",
    "Train_temp = \n",
    "    \n",
    "pca = PCA(n_components=comp, whiten=False)\n",
    "pca.fit(Data.features_train[train_index,5:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DO not run\n",
    "folds =4\n",
    "\n",
    "kf = KFold(n_splits=folds,shuffle=True, random_state = 3)\n",
    "C_space = np.logspace(-7,1,25)\n",
    "comp=50\n",
    "\n",
    "idx_list = np.array([])\n",
    "sweep = np.arange(0,train_pca.shape[1])\n",
    "Train_Max = [] \n",
    "Valid_Max = []\n",
    "Valid_STD = []\n",
    "Test_Max = []\n",
    "for i in range(10):   \n",
    "    Train_Acc = []\n",
    "    Valid_Acc = []\n",
    "    Valid_std = []\n",
    "    Test_Acc = []\n",
    "    for idx in sweep:\n",
    "        idxs = np.hstack((idx_list,idx)).astype('int')\n",
    "        i=0\n",
    "        Train_temp = np.zeros((folds,C_space.shape[0]))\n",
    "        Valid_temp = np.zeros((folds,C_space.shape[0]))\n",
    "        for train_index, valid_index in kf.split(Data.features_train):\n",
    "            pca = PCA(n_components=comp, whiten=False)\n",
    "            pca.fit(Data.features_train[train_index,5:])\n",
    "            X_train_pca = pca.transform(Data.features_train[train_index,5:])\n",
    "            X_valid_pca = pca.transform(Data.features_train[valid_index,5:])\n",
    "\n",
    "            train_pca = np.hstack((X_train_pca,Data.features_train[train_index,0:5]))\n",
    "            valid_pca = np.hstack((X_valid_pca,Data.features_train[valid_index,0:5]))\n",
    "\n",
    "            if len(idxs)==1:\n",
    "                train = np.reshape(train_pca[:,idxs],(-1,1))\n",
    "                valid = np.reshape(valid_pca[:,idxs],(-1,1))\n",
    "            else:\n",
    "                train = train_pca[:,idxs]\n",
    "                valid = valid_pca[:,idxs]\n",
    "            SVM_PCA = Train_Test_SVM_KCross(C_space,train,valid,train_index,valid_index)\n",
    "            Train_temp[i,:] = SVM_PCA[0]\n",
    "            Valid_temp[i,:] = SVM_PCA[1]\n",
    "            i+=1\n",
    "        Train_mean = np.mean(Train_temp,axis=0)\n",
    "        Valid_mean = np.mean(Valid_temp,axis=0)\n",
    "        Valid_stds = np.std(Valid_temp,axis=0)\n",
    "        \n",
    "        Train_Acc.append(Train_mean[Valid_mean.argmax()])\n",
    "        Valid_Acc.append(Valid_mean[Valid_mean.argmax()])\n",
    "        Valid_std.append(Valid_stds[Valid_mean.argmax()])\n",
    " #       Test_Acc.append(SVM_PCA[4])\n",
    "        \n",
    "    Train_Acc = np.array(Train_Acc)\n",
    "    Valid_Acc = np.array(Valid_Acc)\n",
    "    idx_list = np.hstack((idx_list,sweep[Valid_Acc.argmax()]))\n",
    "    sweep = np.delete(sweep, Valid_Acc.argmax())\n",
    "    Train_Max.append(Train_Acc[Valid_Acc.argmax()])\n",
    "    Valid_Max.append(max(Valid_Acc))\n",
    "    Valid_STD.append(Valid_std[Valid_Acc.argmax()])\n",
    "#    Test_Max.append(Test_Acc[Valid_Acc.argmax()])\n",
    "    print(idx_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not run\n",
    "\n",
    "xrange = np.arange(1,31)\n",
    "\n",
    "plt.errorbar(x = xrange, y=Valid_Max,yerr=Valid_STD)\n",
    "plt.plot(xrange,Test_Max,'.')\n",
    "plt.plot(xrange,Train_Max,'--')\n",
    "plt.xlabel(\"Components\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"Test\",\"Train (Avg)\",\"Valid\"],loc='lower right')\n",
    "plt.title(\"Forward Step Selection PCA\")\n",
    "# plt.savefig(\"ForwardStep.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nilearn import masking\n",
    "\n",
    "def Plot_EigenBrain(slices,minV,maxV):\n",
    "    fig,axes= plt.subplots(ncols=3)\n",
    "    axes[0].imshow(np.rot90(slices[0]), vmin=minV, vmax=maxV,cmap='seismic')\n",
    "    axes[0].set_title(\"Sagittal\")\n",
    "    \n",
    "    axes[1].imshow(np.rot90(slices[1]), vmin=minV, vmax=maxV,cmap='seismic')\n",
    "    axes[1].set_title(\"Axial\")\n",
    "    \n",
    "    a = axes[2].imshow(slices[2], vmin=minV, vmax=maxV,cmap='seismic')\n",
    "    axes[2].set_title(\"Coronal\")\n",
    "    \n",
    "    cbar = fig.colorbar(a,ax=axes[2], shrink=0.4)\n",
    "    tick_locator = ticker.MaxNLocator(nbins=4)\n",
    "    cbar.locator = tick_locator\n",
    "    cbar.ax.yaxis.set_major_locator(ticker.AutoLocator())\n",
    "    cbar.update_ticks()\n",
    "    \n",
    "    axes[0].axis(\"off\")\n",
    "    axes[1].axis(\"off\")\n",
    "    axes[2].axis(\"off\")\n",
    "    \n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Atlas.Print_ROIs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "comp = 60\n",
    "\n",
    "pca = PCA(n_components=comp, whiten=False)\n",
    "pca.fit(Data.features_train[:,5:])\n",
    "\n",
    "X_train_pca = pca.transform(Data.features_train[:,5:])\n",
    "X_test_pca = pca.transform(Data.features_test[:,5:])\n",
    "\n",
    "train_pca = np.hstack((X_train_pca,Data.features_train[:,0:5]))\n",
    "test_pca = np.hstack((X_test_pca,Data.features_test[:,0:5]))\n",
    "\n",
    "#Most predictive features (5. 13. 47. 31. 45.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "\n",
    "plotting.view_img(masking.unmask(pca.components_[49], Data.masker))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_img(Atlas.Mask([37]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Atlas.Print_ROIs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_slices = [9,5,49] #[5,13,47,31] \n",
    "for sliver in print_slices:\n",
    "    slicer = masking.unmask(pca.components_[sliver], Data.masker).get_fdata()\n",
    "    slices = [0,0,0]\n",
    "    slices[0] = slicer[45,:,:]\n",
    "    slices[1] = slicer[:,54,:]\n",
    "    slices[2] = slicer[:,:,45]\n",
    "    a = Plot_EigenBrain(slices,np.min(slicer),np.max(slicer))\n",
    "    a.savefig(\"Video/PCA\"+str(sliver)+\".pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data.idx_AD.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data.idx_CN.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Train_Test_SVM(train,test):\n",
    "#     print(\"Training SVM...\")\n",
    "    train_score=[]\n",
    "    valid_score=[]\n",
    "    std_valid=[]\n",
    "    C_space = np.logspace(-7,0,25)\n",
    "    for C in C_space:\n",
    "        SVM = svm.SVC(kernel='linear',C=C)\n",
    "        cvs = cross_validate(SVM,train,Data.labels_train, cv=4,return_train_score=True)\n",
    "        train_score.append(np.mean(cvs[\"train_score\"]))\n",
    "        valid_score.append(np.mean(cvs[\"test_score\"]))\n",
    "        std_valid.append(np.std(cvs[\"test_score\"]))\n",
    "\n",
    "    max_idx = valid_score.index(max(valid_score))\n",
    "    maxC = C_space[max_idx]\n",
    "    \n",
    "#     if maxC == C_space[0] or maxC == C_space[-1]:\n",
    "#         print(\"Error: C_range limit hit\")\n",
    "#         print(maxC)\n",
    "\n",
    "    SVM = svm.SVC(kernel='linear',C=maxC)\n",
    "    SVM.fit(train,Data.labels_train)\n",
    "    test_acc = SVM.score(test,Data.labels_test)\n",
    "\n",
    "    return(maxC,train_score[max_idx],valid_score[max_idx],std_valid[max_idx],test_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_list = np.array([])\n",
    "sweep = np.arange(0,train_pca.shape[1])\n",
    "Train_Max = [] \n",
    "Valid_Max = []\n",
    "Valid_STD = []\n",
    "Test_Max = []\n",
    "for i in range(50):   \n",
    "    Train_Acc = []\n",
    "    Valid_Acc = []\n",
    "    Valid_std = []\n",
    "    Test_Acc = []\n",
    "    for idx in sweep:\n",
    "        idxs = np.hstack((idx_list,idx)).astype('int')\n",
    "        if len(idxs)==1:\n",
    "            train = np.reshape(train_pca[:,idxs],(-1,1))\n",
    "            test = np.reshape(test_pca[:,idxs],(-1,1))\n",
    "        else:\n",
    "            train = train_pca[:,idxs]\n",
    "            test = test_pca[:,idxs]\n",
    "        SVM_PCA = Train_Test_SVM(train,test)\n",
    "        Train_Acc.append(SVM_PCA[1])\n",
    "        Valid_Acc.append(SVM_PCA[2])\n",
    "        Valid_std.append(SVM_PCA[3])\n",
    "        Test_Acc.append(SVM_PCA[4])\n",
    "        \n",
    "    Train_Acc = np.array(Train_Acc)\n",
    "    Valid_Acc = np.array(Valid_Acc)\n",
    "    idx_list = np.hstack((idx_list,sweep[Valid_Acc.argmax()]))\n",
    "    sweep = np.delete(sweep, Valid_Acc.argmax())\n",
    "    Train_Max.append(Train_Acc[Valid_Acc.argmax()])\n",
    "    Valid_Max.append(max(Valid_Acc))\n",
    "    Valid_STD.append(Valid_std[Valid_Acc.argmax()])\n",
    "    Test_Max.append(Test_Acc[Valid_Acc.argmax()])\n",
    "    print(idx_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#18 per row\n",
    "xrange = np.arange(1,51)\n",
    "\n",
    "plt.errorbar(x = xrange, y=Valid_Max,yerr=Valid_STD)\n",
    "plt.plot(xrange,Test_Max,'.')\n",
    "plt.plot(xrange,Train_Max,'--')\n",
    "plt.xlabel(\"Components\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"Test\",\"Train (Avg)\",\"Valid\"],loc='lower right')\n",
    "# plt.title(\"Forward Step Selection PCA\")\n",
    "plt.savefig(\"ForwardStep.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(Test_Max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Test_Max[Valid_Max.index(max(Valid_Max))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(Valid_Max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_pca = pca.transform(Data.features_train[:,5:])\n",
    "X_test_pca = pca.transform(Data.features_test[:,5:])\n",
    "\n",
    "train_pca = np.hstack((X_train_pca,Data.features_train[:,0:5]))\n",
    "test_pca = np.hstack((X_test_pca,Data.features_test[:,0:5]))\n",
    "\n",
    "print(\"Training SVM...\")\n",
    "train_score_pca=[]\n",
    "train_score_ica=[]\n",
    "test_score_pca=[]\n",
    "test_score_ica=[]\n",
    "C_space = np.logspace(-4,15,30)\n",
    "for C in C_space:\n",
    "    SVM = svm.LinearSVC(penalty='l2',C=C,max_iter=1000000,dual=False)\n",
    "    cvs = cross_validate(SVM,train_pca,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_pca.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_pca.append(np.mean(cvs[\"test_score\"]))\n",
    "    \n",
    "    cvs = cross_validate(SVM,train_ica,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_ica.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_ica.append(np.mean(cvs[\"test_score\"]))\n",
    "\n",
    "maxC_pca = C_space[test_score_pca.index(max(test_score_pca))]\n",
    "maxC_ica = C_space[test_score_ica.index(max(test_score_ica))]\n",
    "\n",
    "plt.semilogx(C_space,test_score_pca)\n",
    "plt.plot(C_space,test_score_ica)\n",
    "plt.xlabel(\"C\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"PCA\",\"ICA\"])\n",
    "\n",
    "print(\"Training NN...\")\n",
    "train_score_pca=[]\n",
    "train_score_ica=[]\n",
    "test_score_pca=[]\n",
    "test_score_ica=[]\n",
    "neighbor_space = np.arange(1,5)\n",
    "for neigh in neighbor_space:\n",
    "    NN = neighbors.KNeighborsClassifier(neigh)\n",
    "    cvs = cross_validate(NN,train_pca,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_pca.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_pca.append(np.mean(cvs[\"test_score\"]))\n",
    "    \n",
    "    cvs = cross_validate(NN,train_ica,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_ica.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_ica.append(np.mean(cvs[\"test_score\"]))\n",
    "\n",
    "plt.figure()\n",
    "plt.semilogx(neighbor_space,test_score_pca)\n",
    "plt.plot(neighbor_space,test_score_ica)\n",
    "plt.xlabel(\"Neighbors\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"PCA\",\"ICA\"])\n",
    "\n",
    "print(\"Training AdaBoost...\")\n",
    "train_score_pca=[]\n",
    "train_score_ica=[]\n",
    "test_score_pca=[]\n",
    "test_score_ica=[]\n",
    "neighbor_space = np.arange(1,5)\n",
    "for neigh in neighbor_space:\n",
    "    NN = neighbors.KNeighborsClassifier(neigh)\n",
    "    cvs = cross_validate(NN,train_pca,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_pca.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_pca.append(np.mean(cvs[\"test_score\"]))\n",
    "    \n",
    "    cvs = cross_validate(NN,train_ica,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_ica.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_ica.append(np.mean(cvs[\"test_score\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([1,2,3,4,5])\n",
    "y=np.vstack((2.1*x,2*x,1.9*x))\n",
    "yerr = np.std(y,axis=0)\n",
    "plt.figure()\n",
    "plt.errorbar(x,np.mean(y,axis=0),yerr)\n",
    "plt.xscale(\"log\", nonposx='clip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yerr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#L1\n",
    "print(max(test_score_pca),max(test_score_ica))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#L2\n",
    "print(max(test_score_pca),max(test_score_ica))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cca = CCA(n_components=comp)\n",
    "cca.fit(Data.features_train[:,5:], Data.labels_train)\n",
    "X_train_cca = cca.transform(Data.features_train[:,5:])\n",
    "X_test_cca = cca.transform(Data.features_test[:,5:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_score_cca=[]\n",
    "test_score_cca=[]\n",
    "\n",
    "C_space = np.logspace(-4,15,30)\n",
    "for C in C_space:\n",
    "    SVM = svm.LinearSVC(penalty='l2',C=C,max_iter=1000000,dual=False)\n",
    "    cvs = cross_validate(SVM,train_cca,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score_cca.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score_cca.append(np.mean(cvs[\"test_score\"]))\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Data.features_train[:,5:].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "PCA_space = np.arange(1,25)*6\n",
    "Train_Best = []\n",
    "Valid_Best = []\n",
    "Test_Best = []\n",
    "MaxCs = []\n",
    "\n",
    "for comp in PCA_space:\n",
    "    pca = PCA(n_components=comp)\n",
    "    pca.fit(Data.features_train[:,5:])\n",
    "\n",
    "    X_train_pca = pca.transform(Data.features_train[:,5:])\n",
    "    X_test_pca = pca.transform(Data.features_test[:,5:])\n",
    "\n",
    "    train = np.hstack((X_train_pca,Data.features_train[:,0:5]))\n",
    "    test = np.hstack((X_test_pca,Data.features_test[:,0:5]))\n",
    "\n",
    "    print(\"Training....\")\n",
    "    train_score = []\n",
    "    test_score = []\n",
    "    C_space = np.logspace(-7,5,30)\n",
    "    for C in C_space:\n",
    "        SVM = svm.SVC(kernel='poly',degree=1, C=C)\n",
    "        cvs = cross_validate(SVM,train,Data.labels_train, cv=4,return_train_score=True)\n",
    "        train_score.append(np.mean(cvs[\"train_score\"]))\n",
    "        test_score.append(np.mean(cvs[\"test_score\"]))\n",
    "    #     print(\"C=\"+str(C)+\"Train\"+str(np.mean(cvs[\"train_score\"]))+\"Valid\"+str(np.mean(cvs[\"test_score\"])))\n",
    "\n",
    "    maxC = C_space[test_score.index(max(test_score))]\n",
    "\n",
    "    SVM = svm.SVC(kernel='linear', C=maxC)\n",
    "    SVM.fit(train,Data.labels_train)\n",
    "    pred_test=SVM.predict(test)\n",
    "\n",
    "    Train_Best.append(train_score[test_score.index(max(test_score))])\n",
    "    Valid_Best.append(max(test_score))\n",
    "    Test_Best.append(accuracy_score(pred_test,Data.labels_test))\n",
    "    MaxCs.append(maxC)\n",
    "\n",
    "    print(\"PCA = \"+str(comp))\n",
    "    print(\"Max Acc = \"+str(np.around(max(test_score),4))+\", with C = \"+str(maxC)+\" Test Acc = \"+str(accuracy_score(pred_test,Data.labels_test)))\n",
    "\n",
    "plt.plot(PCA_space,Train_Best)\n",
    "plt.plot(PCA_space,Valid_Best)\n",
    "plt.plot(PCA_space,Test_Best)\n",
    "plt.xlabel(\"PCA components\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"Train\",\"Valid\",\"Test\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "PCA_space = np.arange(1,15)*2\n",
    "Train_Best = []\n",
    "Valid_Best = []\n",
    "Test_Best = []\n",
    "MaxCs = []\n",
    "\n",
    "for comp in PCA_space:\n",
    "    pca = PCA(n_components=comp)\n",
    "    pca.fit(Data.features_train[:,5:])\n",
    "\n",
    "    X_train_pca = pca.transform(Data.features_train[:,5:])\n",
    "    X_test_pca = pca.transform(Data.features_test[:,5:])\n",
    "\n",
    "    train = np.hstack((X_train_pca,Data.features_train[:,0:5]))\n",
    "    test = np.hstack((X_test_pca,Data.features_test[:,0:5]))\n",
    "\n",
    "    print(\"Training....\")\n",
    "    train_score = []\n",
    "    test_score = []\n",
    "    C_space = np.logspace(-7,3,30)\n",
    "    for C in C_space:\n",
    "        SVM = svm.SVC(kernel='linear',C=C)\n",
    "        cvs = cross_validate(SVM,train,Data.labels_train, cv=4,return_train_score=True)\n",
    "        train_score.append(np.mean(cvs[\"train_score\"]))\n",
    "        test_score.append(np.mean(cvs[\"test_score\"]))\n",
    "    #     print(\"C=\"+str(C)+\"Train\"+str(np.mean(cvs[\"train_score\"]))+\"Valid\"+str(np.mean(cvs[\"test_score\"])))\n",
    "\n",
    "    maxC = C_space[test_score.index(max(test_score))]\n",
    "\n",
    "    SVM = svm.SVC(kernel='linear', C=maxC)\n",
    "    SVM.fit(train,Data.labels_train)\n",
    "    pred_test=SVM.predict(test)\n",
    "\n",
    "    Train_Best.append(train_score[test_score.index(max(test_score))])\n",
    "    Valid_Best.append(max(test_score))\n",
    "    Test_Best.append(accuracy_score(pred_test,Data.labels_test))\n",
    "    MaxCs.append(maxC)\n",
    "\n",
    "    print(\"PCA = \"+str(comp))\n",
    "    print(\"Max Acc = \"+str(np.around(max(test_score),4))+\", with C = \"+str(maxC)+\" Test Acc = \"+str(accuracy_score(pred_test,Data.labels_test)))\n",
    "\n",
    "plt.plot(PCA_space,Train_Best)\n",
    "plt.plot(PCA_space,Valid_Best)\n",
    "plt.plot(PCA_space,Test_Best)\n",
    "plt.xlabel(\"PCA components\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"Train\",\"Valid\",\"Test\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train = Data.features_train[:,0:5]\n",
    "test = Data.features_test[:,0:5]\n",
    "\n",
    "train_score = []\n",
    "test_score = []\n",
    "C_space = np.logspace(-7,3,30)\n",
    "for C in C_space:\n",
    "    SVM = svm.SVC(kernel='linear',degree=2, C=C)\n",
    "    cvs = cross_validate(SVM,train,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score.append(np.mean(cvs[\"train_score\"]))\n",
    "    test_score.append(np.mean(cvs[\"test_score\"]))\n",
    "#     print(\"C=\"+str(C)+\"Train\"+str(np.mean(cvs[\"train_score\"]))+\"Valid\"+str(np.mean(cvs[\"test_score\"])))\n",
    "\n",
    "maxC = C_space[test_score.index(max(test_score))]\n",
    "\n",
    "SVM = svm.SVC(kernel='linear', C=maxC)\n",
    "SVM.fit(train,Data.labels_train)\n",
    "pred_test=SVM.predict(test)\n",
    "\n",
    "Train_Best.append(train_score[test_score.index(max(test_score))])\n",
    "Valid_Best.append(max(test_score))\n",
    "Test_Best.append(accuracy_score(pred_test,Data.labels_test))\n",
    "MaxCs.append(maxC)\n",
    "\n",
    "print(\"Max Acc = \"+str(np.around(max(test_score),4))+\", with C = \"+str(maxC)+\" Test Acc = \"+str(accuracy_score(pred_test,Data.labels_test)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVM vs Ada comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "FileDir = os.getcwd()+'/SVM_Outputs/'\n",
    "file_object = open(FileDir+'SVM_Log.txt', \"r\")\n",
    "Data_Loader = np.zeros((45,5))\n",
    "for (i,line) in enumerate(file_object):\n",
    "    split = line.split(',')\n",
    "    Data_Loader[i,0]=split[0] #ATLAS\n",
    "    Data_Loader[i,1]=split[1] #Features\n",
    "    Data_Loader[i,2]=split[2] #C\n",
    "    Data_Loader[i,3]=split[3] #Valid\n",
    "    Data_Loader[i,4]=split[4] #Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "FileDir = os.getcwd()+'/Ada_Outputs/'\n",
    "file_object = open(FileDir+'Ada_Log.txt', \"r\")\n",
    "Data_Loader_Ada = np.zeros((45,5))\n",
    "for (i,line) in enumerate(file_object):\n",
    "    split = line.split(',')\n",
    "    Data_Loader_Ada[i,0]=split[0] #ATLAS\n",
    "    Data_Loader_Ada[i,1]=split[1] #Features\n",
    "    Data_Loader_Ada[i,2]=split[2] #C\n",
    "    Data_Loader_Ada[i,3]=split[3] #Valid\n",
    "    Data_Loader_Ada[i,4]=split[4] #Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(Data_Loader[:,0],Data_Loader[:,3],'b')\n",
    "plt.plot(Data_Loader[:,0],Data_Loader[:,4],'--b')\n",
    "plt.plot(Data_Loader_Ada[:,0],Data_Loader_Ada[:,3],'r')\n",
    "plt.plot(Data_Loader_Ada[:,0],Data_Loader_Ada[:,4],'--r')\n",
    "plt.xlabel(\"Atlas ROI\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"SVM Valid\",\"SVM Test\",\"Ada Valid\",\"Ada Test\"],loc='center left', bbox_to_anchor=(1, 0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idxs = np.argsort(Data_Loader[:,3]+Data_Loader_Ada[:,3])\n",
    "print(Data_Loader[idxs,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idxs = np.argsort(Data_Loader[:,3])\n",
    "print(Data_Loader[idxs,0])\n",
    "print(Data_Loader[idxs,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idxs = np.argsort(Data_Loader_Ada[:,3])\n",
    "print(Data_Loader_Ada[idxs,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVM Ensemble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn import svm\n",
    "from sklearn.model_selection import cross_validate\n",
    "\n",
    "from DataPrep import DATA\n",
    "\n",
    "FileDir = os.getcwd()+'/SVM_Outputs/'\n",
    "file_object = open(FileDir+'SVM_Log.txt', \"r\")\n",
    "Data_Loader = np.zeros((45,5))\n",
    "for (i,line) in enumerate(file_object):\n",
    "    split = line.split(',')\n",
    "    Data_Loader[i,0]=split[0] #ATLAS\n",
    "    Data_Loader[i,1]=split[1] #Features\n",
    "    Data_Loader[i,2]=split[2] #C\n",
    "    Data_Loader[i,3]=split[3] #Valid\n",
    "    Data_Loader[i,4]=split[4] #Test\n",
    "\n",
    "Data = DATA()\n",
    "Data.Train_Test(0.8,12345)\n",
    "#test on 3,8, 26 - DL 2,7,25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ensemble = [] \n",
    "for (i,DL) in enumerate(Data_Loader):\n",
    "    if DL[3]>=0.58:\n",
    "        Ensemble.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Predictions = np.zeros((len(Ensemble),Data.idx_test.shape[0]))\n",
    "for (i,Ens) in enumerate(Ensemble):\n",
    "    Selector = Data_Loader[Ens][0]\n",
    "    maxC = Data_Loader[Ens][2]\n",
    "    \n",
    "    Data.Add_MRI([Selector])\n",
    "    Data.Split_Data()\n",
    "    \n",
    "    SVM = svm.SVC(kernel='linear', C=maxC)\n",
    "    SVM.fit(Data.features_train,Data.labels_train)\n",
    "    \n",
    "    pred_test = SVM.predict(Data.features_test)\n",
    "    acc = accuracy_score(pred_test,Data.labels_test)\n",
    "    \n",
    "    Predictions[i]=pred_test\n",
    "    \n",
    "    print(Selector,acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "accuracy_score(scipy.stats.mode(Predictions[:,:],axis=0)[0][0],Data.labels_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# View Brain (Video Prep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn import svm,neighbors\n",
    "from sklearn.model_selection import cross_validate\n",
    "\n",
    "from DataPrep import DATA\n",
    "from ATLAS import ATLAS\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.decomposition import FastICA\n",
    "\n",
    "from sklearn.ensemble import AdaBoostClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data = DATA()\n",
    "Data.Train_Test(0.8,12345)\n",
    "Data.Add_MRI(\"brain\") #\"brain\" or [34]\n",
    "Data.Split_Data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = Data.features_train\n",
    "test = Data.features_test\n",
    "\n",
    "train_score=[]\n",
    "valid_score=[]\n",
    "std_valid=[]\n",
    "#0.0004281332398719391 whole brain\n",
    "C_space = np.logspace(-2,1,30)\n",
    "C_space = [0.0004281332398719391]\n",
    "for C in C_space:\n",
    "    SVM = svm.SVC(kernel='linear',C=C)\n",
    "    cvs = cross_validate(SVM,train,Data.labels_train, cv=4,return_train_score=True)\n",
    "    train_score.append(np.mean(cvs[\"train_score\"]))\n",
    "    valid_score.append(np.mean(cvs[\"test_score\"]))\n",
    "    std_valid.append(np.std(cvs[\"test_score\"]))\n",
    "    print(C,np.mean(cvs[\"train_score\"]),np.mean(cvs[\"test_score\"]),np.std(cvs[\"test_score\"]))\n",
    "\n",
    "max_idx = valid_score.index(max(valid_score))\n",
    "maxC = C_space[max_idx]\n",
    "\n",
    "SVM = svm.SVC(kernel='linear',C=maxC)\n",
    "SVM.fit(train,Data.labels_train)\n",
    "test_acc = SVM.score(test,Data.labels_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(x = C_space, y=valid_score,yerr=std_valid)\n",
    "# plt.plot(C_space,Test_Max,'.')\n",
    "plt.semilogx(C_space,train_score,'-')\n",
    "plt.xlabel(\"C\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend([\"Train (Avg)\",\"Valid\"],loc='lower right')\n",
    "plt.title(\"Whole Brain\") #SVM_34_Parahippocampal Gyrus, anterior division\n",
    "plt.ylim([0.4,1.02])\n",
    "print(test_acc)\n",
    "# plt.savefig(\"Video/SVM_WholeBrain.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(valid_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import masking\n",
    "from nilearn import plotting\n",
    "\n",
    "a = masking.unmask(SVM.coef_[0,5:], Data.masker)\n",
    "plotting.view_img(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie = a.get_fdata()\n",
    "plt.imshow(np.rot90(movie[32,:,:]), vmin=-np.max(np.abs(movie)), vmax=np.max(np.abs(movie)),cmap='seismic')\n",
    "plt.colorbar()\n",
    "plt.axis(\"off\")\n",
    "plt.savefig(\"Video/SVM_Hippo34_Weights.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie.shape[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from celluloid import Camera\n",
    "\n",
    "movie = a.get_fdata()\n",
    "\n",
    "fig,axes= plt.subplots(ncols=2)\n",
    "camera = Camera(fig)\n",
    "for i in range(movie.shape[2]):\n",
    "    axes[0].imshow(movie[:,:,i], vmin=np.min(movie), vmax=np.max(movie),cmap='seismic')\n",
    "    axs = axes[1].imshow(np.rot90(movie[i,:,:]), vmin=np.min(movie), vmax=np.max(movie),cmap='seismic')\n",
    "    axes[0].axis(\"off\")\n",
    "    axes[1].axis(\"off\")\n",
    "    camera.snap()\n",
    "animation = camera.animate()\n",
    "animation.save('Video/SVM_BrainCoeff.mp4', fps=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(movie[:,:,i], vmin=np.min(movie), vmax=np.max(movie),cmap='seismic')\n",
    "cbar = plt.colorbar()\n",
    "tick_locator = ticker.MaxNLocator(nbins=5)\n",
    "cbar.locator = tick_locator\n",
    "cbar.ax.yaxis.set_major_locator(ticker.AutoLocator())\n",
    "cbar.update_ticks()\n",
    "plt.savefig(\"Video/Colobar_for_BrainCoeff.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from celluloid import Camera\n",
    "\n",
    "movie = Atlas.mask_img.get_fdata()\n",
    "\n",
    "fig,axes= plt.subplots(ncols=2)\n",
    "camera = Camera(fig)\n",
    "for i in range(shape[2]):\n",
    "    axes[0].imshow(movie[:,:,i], vmin=np.min(movie), vmax=np.max(movie),cmap='seismic')\n",
    "    axes[1].imshow(np.rot90(movie[i,:,:]), vmin=np.min(movie), vmax=np.max(movie),cmap='seismic')\n",
    "#     plt.colorbar()\n",
    "    axes[0].axis(\"off\")\n",
    "    axes[1].axis(\"off\")\n",
    "    camera.snap()\n",
    "animation = camera.animate()\n",
    "animation.save('Brain_Atlas2.mp4', fps=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Atlas=ATLAS()\n",
    "print(mask_img.shape)\n",
    "plotting.view_img(Atlas.mask_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slicer = Atlas.Mask([35]).get_fdata()\n",
    "slices = [0,0,0]\n",
    "slices[0] = slicer[45,:,:]\n",
    "slices[1] = slicer[:,54,:]\n",
    "slices[2] = slicer[:,:,45]\n",
    "a = Plot_EigenBrain(slices,np.min(slicer),np.max(slicer))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from celluloid import Camera\n",
    "import nibabel as nib\n",
    "gm = nib.load(Data.gm_imgs[0])\n",
    "wm = nib.load(Data.wm_imgs[0])\n",
    "movie = gm.get_fdata()#+wm.get_fdata()\n",
    "\n",
    "fig,axes= plt.subplots(ncols=2)\n",
    "camera = Camera(fig)\n",
    "for i in range(shape[2]):\n",
    "    axes[0].imshow(movie[:,:,i], vmin=np.min(movie), vmax=np.max(movie),cmap='gray')\n",
    "    axes[1].imshow(np.rot90(movie[i,:,:]), vmin=np.min(movie), vmax=np.max(movie),cmap='gray')\n",
    "#     plt.colorbar()\n",
    "    axes[0].axis(\"off\")\n",
    "    axes[1].axis(\"off\")\n",
    "    camera.snap()\n",
    "animation = camera.animate()\n",
    "animation.save('GMScan.mp4', fps=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data.idx_CN.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data.labels_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "27/39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}