roc_curves.ipynb

{
 "metadata": {
  "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.6.7"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.6.7 64-bit ('pytorch': conda)"
  },
  "interpreter": {
   "hash": "62eb3f0dc0d696e2dd85d645fd9023f02c025c88b4bd5b0bd938b54b93374391"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Plot ROC curves\r\n",
    "\r\n",
    "Maintain by Xiaotong.Chen"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import torch\r\n",
    "import torch.nn as nn\r\n",
    "import torch.nn.functional as F\r\n",
    "import torch.optim as optim\r\n",
    "import torch.optim.lr_scheduler as lr_scheduler\r\n",
    "from torch.optim.lr_scheduler import _LRScheduler\r\n",
    "import torch.utils.data as data\r\n",
    "import torchvision.transforms as transforms\r\n",
    "import torchvision.datasets as datasets\r\n",
    "import torchvision.models as models\r\n",
    "from torch.utils.data import Dataset, DataLoader\r\n",
    "from torchvision import  utils\r\n",
    "\r\n",
    "from sklearn import decomposition\r\n",
    "from sklearn import manifold\r\n",
    "from sklearn import metrics\r\n",
    "from sklearn import model_selection\r\n",
    "from sklearn.metrics import confusion_matrix\r\n",
    "\r\n",
    "import matplotlib.pyplot as plt\r\n",
    "%matplotlib inline\r\n",
    "\r\n",
    "\r\n",
    "import copy\r\n",
    "from collections import namedtuple\r\n",
    "import os, re, time\r\n",
    "import random\r\n",
    "import shutil\r\n",
    "import math\r\n",
    "\r\n",
    "import pandas as pd\r\n",
    "import numpy as np\r\n",
    "from sklearn import svm\r\n",
    "\r\n",
    "# set seed to make sure the results are reproducible\r\n",
    "SEED = 99\r\n",
    "random.seed(SEED)\r\n",
    "np.random.seed(SEED)\r\n",
    "torch.manual_seed(SEED)\r\n",
    "torch.cuda.manual_seed(SEED)\r\n",
    "torch.backends.cudnn.deterministic = True\r\n",
    "\r\n",
    "# change directory of data\r\n",
    "CV_ID = 2\r\n",
    "CV_ID = str(CV_ID)\r\n",
    "datapath = os.path.join(\"D:\\\\0jiazhuangxian2020\\\\jiazhuangxian\",\".data\",\"hcb\",\"cv\" + CV_ID)\r\n",
    "load_name = 'basic-model-cv'+CV_ID+'.pt'\r\n",
    "train_dir = os.path.join(datapath, 'train')\r\n",
    "val_dir = os.path.join(datapath, 'val')\r\n",
    "test_dir = os.path.join(datapath, 'test')\r\n",
    "\r\n",
    "pretrained_size = (224, 224)\r\n",
    "pretrained_means = [0.485, 0.456, 0.406]\r\n",
    "pretrained_stds= [0.229, 0.224, 0.225]\r\n",
    "test_transforms = transforms.Compose([\r\n",
    "                           transforms.Resize(pretrained_size),\r\n",
    "                           transforms.ToTensor(),\r\n",
    "                           transforms.Normalize(mean = pretrained_means, \r\n",
    "                                                std = pretrained_stds)\r\n",
    "                       ])\r\n",
    "# inorder to load img with it's label\r\n",
    "class MyImageFolder(datasets.ImageFolder):\r\n",
    "    def __getitem__(self, index):\r\n",
    "        path, _ = self.imgs[index] #img path, label\r\n",
    "        return super(MyImageFolder, self).__getitem__(index), path # return image path\r\n",
    "    \r\n",
    "\r\n",
    "BATCH_SIZE = 1000 #4\r\n",
    "train_data = MyImageFolder(root = train_dir,\r\n",
    "                                transform = test_transforms)\r\n",
    "test_data = MyImageFolder(root = test_dir,\r\n",
    "                                transform = test_transforms)\r\n",
    "valid_data = MyImageFolder(root = val_dir,\r\n",
    "                                transform = test_transforms)\r\n",
    "train_iterator = data.DataLoader(train_data, \r\n",
    "                                shuffle = True,\r\n",
    "                                drop_last = False,\r\n",
    "                                batch_size = BATCH_SIZE) \r\n",
    "valid_iterator = data.DataLoader(valid_data, drop_last = False,\r\n",
    "                                batch_size = BATCH_SIZE)\r\n",
    "test_iterator = data.DataLoader(test_data, drop_last = False,\r\n",
    "                                batch_size = BATCH_SIZE)     \r\n",
    "\r\n",
    "\r\n",
    "# read csv data\r\n",
    "path_to_text = '.data\\\\text.csv'\r\n",
    "text_table = pd.read_csv(path_to_text, header=0, index_col=0)\r\n",
    "Dimension_Text = len(text_table.columns)\r\n",
    "\r\n",
    "def imgid2index(id):\r\n",
    "    return 1000*int(id[0])+int(id[2:])\r\n",
    "\r\n",
    "def imgid2textinfo(imgid):\r\n",
    "    # convert img path into text info in .csv\r\n",
    "    # input: ['.data\\\\train\\\\zy\\\\1_152.bmp','2_8.bmp']\r\n",
    "    # output: tensor([[ 1,  4,  7, 10],\r\n",
    "    #    [ 3,  6,  9, 12]])\r\n",
    "    return torch.tensor(text_table.loc[[imgid2index(re.search('\\d_\\d+',i.split('\\\\')[-1]).group()) for i in imgid],:].values, dtype=torch.float32, device=device)\r\n",
    "\r\n",
    "basic_model = torch.hub.load('pytorch/vision:v0.6.0', 'resnet50')\r\n",
    "# change output dimension to what we need\r\n",
    "IN_FEATURES = basic_model.fc.in_features \r\n",
    "OUTPUT_DIM = 2\r\n",
    "basic_model.fc = nn.Linear(IN_FEATURES, OUTPUT_DIM)\r\n",
    "fc1=nn.Linear(IN_FEATURES, 32)\r\n",
    "fc2=nn.Linear(32,OUTPUT_DIM)\r\n",
    "basic_model.fc = nn.Sequential(fc1, fc2)\r\n",
    "basic_model.load_state_dict(torch.load(load_name))\r\n",
    "basic_model.fc[1] = nn.Identity()\r\n",
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\r\n",
    "# text_table.head()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# load train data\r\n",
    "num_train = len(train_data)\r\n",
    "train_data_combine_model = data.DataLoader(train_data, \r\n",
    "                                shuffle = True, \r\n",
    "                                batch_size = num_train) \r\n",
    "xy, train_id = next(iter(train_data_combine_model))\r\n",
    "x, y = xy\r\n",
    "img_output = basic_model(x)\r\n",
    "train_x = torch.cat((img_output, imgid2textinfo(train_id).cpu()), dim=-1).detach().numpy()\r\n",
    "train_y = y\r\n",
    "# load val data and combine with train\r\n",
    "num_val = len(valid_data)\r\n",
    "val_data_combine_model = data.DataLoader(valid_data, \r\n",
    "                                shuffle = True, \r\n",
    "                                batch_size = num_val) \r\n",
    "xy, val_id = next(iter(val_data_combine_model))\r\n",
    "x, y = xy\r\n",
    "img_output = basic_model(x)\r\n",
    "val_x = torch.cat((img_output, imgid2textinfo(val_id).cpu()), dim=-1).detach().numpy()\r\n",
    "val_y = y\r\n",
    "\r\n",
    "train_x = np.concatenate((train_x, val_x), axis=0)\r\n",
    "train_y = np.concatenate((train_y, val_y), axis=0)\r\n",
    "# load test data\r\n",
    "num_test = len(test_data)\r\n",
    "test_data_combine_model = data.DataLoader(test_data, \r\n",
    "                                batch_size = num_test) \r\n",
    "xy, test_id = next(iter(test_data_combine_model))\r\n",
    "x, y = xy\r\n",
    "img_output = basic_model(x)\r\n",
    "test_x = torch.cat((img_output, imgid2textinfo(test_id).cpu()), dim=-1).detach().numpy()\r\n",
    "test_y = y.numpy()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "tmpt=[imgid2index(re.search('\\d_\\d+',i.split('\\\\')[-1]).group()) for i in train_id]\r\n",
    "tmpv=[imgid2index(re.search('\\d_\\d+',i.split('\\\\')[-1]).group()) for i in val_id]\r\n",
    "train_ID=tmpt+tmpv\r\n",
    "test_ID=[imgid2index(re.search('\\d_\\d+',i.split('\\\\')[-1]).group()) for i in test_id]\r\n",
    "\r\n",
    "x_train=text_table.loc[train_ID,:].values\r\n",
    "y_train=[2-i//1000 for i in train_ID]\r\n",
    "# y_train=[i//1000-1 for i in train_ID]\r\n",
    "# print(x_train,y_train)\r\n",
    "\r\n",
    "x_test=text_table.loc[test_ID,:].values\r\n",
    "y_test=[2-i//1000 for i in test_ID]\r\n",
    "# y_test=[i//1000-1 for i in test_ID]\r\n",
    "# print(x_test,y_test)\r\n",
    "\r\n",
    "all_sen={}\r\n",
    "all_spe={}"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# SVM TEXT+IMG"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "#search best params for SVM\r\n",
    "svm_para = svm.SVC()    \r\n",
    "param_grid = {'C': range(5,30), 'gamma': [1e-2, 7e-3, 5e-3, 3e-3, 1e-3, 7e-4, 5e-4, 3e-4, 1e-4, 7e-5], \r\n",
    "              'kernel': ['rbf', 'linear', 'poly', 'sigmoid']}    \r\n",
    "grid_search = model_selection.GridSearchCV(svm_para, param_grid)    \r\n",
    "grid_search.fit(train_x, train_y)    \r\n",
    "best_parameters = grid_search.best_estimator_.get_params()\r\n",
    "for para, val in list(best_parameters.items()):    \r\n",
    "    print(para, val) "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "def roc_results(y_true, y_pred, path):\r\n",
    "    r\"\"\"\r\n",
    "    Return: sensitivity, specificity, ppv, npv, accuracy, AUC, threshold\r\n",
    "    this function is used to calculate AUC,sensitivity,specificity,ppv,npv,accuracy\r\n",
    "    baesed on each threshold and plot the ROC curve. It will also choose the best threshold, which means that this threshold lead to the highest accuracy.\r\n",
    "    input: true, pred, path\r\n",
    "    \"\"\"\r\n",
    "# calculate the target vector from roc\r\n",
    "    fpr, sen, threshold = metrics.roc_curve(y_true, y_pred)\r\n",
    "    spe = 1-fpr\r\n",
    "    pos_num = np.sum(y_true == 1)\r\n",
    "    neg_num = len(y_true)-pos_num\r\n",
    "    tp = sen*pos_num\r\n",
    "    tn = spe*neg_num\r\n",
    "    fn = pos_num-tp\r\n",
    "    fp = neg_num-tn\r\n",
    "    ppv = tp/(tp+fp+1e-16)\r\n",
    "    npv = tn/(tn+fn+1e-16)\r\n",
    "    acc = (tp+tn)/len(y_true)\r\n",
    "# the best point data\r\n",
    "    j_statistic=sen+spe-1#(based on J_statistic)\r\n",
    "    ind=np.argmax(j_statistic)\r\n",
    "    # ind = np.max(np.where(acc == np.max(acc)))#(based on accuracy)\r\n",
    "    pred_p = tp[ind]+fp[ind]\r\n",
    "    pred_n = tn[ind]+fn[ind]\r\n",
    "    best_acc = acc[ind]\r\n",
    "    best_sen = sen[ind]\r\n",
    "    best_spe = spe[ind]\r\n",
    "    best_ppv = tp[ind]/pred_p\r\n",
    "    best_npv = tn[ind]/pred_n\r\n",
    "    AUC = metrics.auc(fpr, sen)\r\n",
    "# plot\r\n",
    "    plt.plot(fpr, sen)\r\n",
    "    plt.plot(np.linspace(0, 1, 10), np.linspace(0, 1, 10))\r\n",
    "    plt.plot(fpr[ind], sen[ind], 'r.', markersize=9)\r\n",
    "    text = \"AUC={:.2f}\\nSEN={:.2f}\\nSPE={:.2f}\\nPPV={:.2f}\\nNPV={:.2f}\\nACC={:.2f}\".format(\r\n",
    "        AUC, sen[ind], spe[ind], best_ppv, best_npv, best_acc)\r\n",
    "    plt.text(0.82, 0.1, text, fontsize=10,\r\n",
    "             style=\"italic\", horizontalalignment=\"center\")\r\n",
    "    plt.xlabel('1-Specificity')\r\n",
    "    plt.ylabel('Sensitivity')\r\n",
    "    plt.title('ROC Curve')\r\n",
    "    path = path+\"/\"+\"cv\"+str(CV_ID)+\".jpg\"\r\n",
    "    plt.savefig(path)\r\n",
    "    plt.show()\r\n",
    "    return sen, spe, ppv, npv, acc, AUC, threshold"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# combine text info and output of img model with Support Vector Machine\r\n",
    "rbf_svc = svm.SVC(C=best_parameters['C'], gamma=best_parameters['gamma'], kernel=best_parameters['kernel'])\r\n",
    "# rbf_svc = svm.SVC(kernel='rbf',C = 18, gamma = 0.001)\r\n",
    "rbf_svc.fit(train_x, train_y)\r\n",
    "path='D:/0jiazhuangxian2020/jiazhuangxian/fig/svmcombine/train'\r\n",
    "svm_pred = rbf_svc.predict(train_x)\r\n",
    "svm_pred_prob=rbf_svc.decision_function(train_x)\r\n",
    "sen, spe, ppv, npv, acc, auc, _=roc_results(1-train_y,1-svm_pred_prob,path)\r\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "svm_pred = rbf_svc.predict(test_x)\r\n",
    "path='D:/0jiazhuangxian2020/jiazhuangxian/fig/svmcombine/test'\r\n",
    "sen, spe, ppv, npv, acc, auc, _=roc_results(1-test_y,1-rbf_svc.decision_function(test_x),path)\r\n",
    "all_sen['svm_c']=sen\r\n",
    "all_spe['svm_c']=spe"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# SVM TEXT"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "#search best params for SVM\r\n",
    "svm_para = svm.SVC()    \r\n",
    "param_grid = {'C': range(5,8), 'gamma': [1e-2, 1e-1,1e-3], \r\n",
    "              'kernel': ['rbf', 'sigmoid']}    \r\n",
    "grid_search = model_selection.GridSearchCV(svm_para, param_grid)    \r\n",
    "grid_search.fit(x_train, y_train)    \r\n",
    "best_parameters = grid_search.best_estimator_.get_params()\r\n",
    "for para, val in list(best_parameters.items()):    \r\n",
    "    print(para, val) "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# combine text info and output of img model with Support Vector Machine\r\n",
    "rbf_svc = svm.SVC(C=best_parameters['C'], gamma=best_parameters['gamma'], kernel=best_parameters['kernel'])\r\n",
    "# rbf_svc = svm.SVC(kernel='sigmoid',C = 5, gamma = 0.1)\r\n",
    "rbf_svc.fit(x_train, y_train)\r\n",
    "\r\n",
    "path='D:/0jiazhuangxian2020/jiazhuangxian/fig/svmtext/train'\r\n",
    "svm_pred = rbf_svc.predict(x_train)\r\n",
    "y_train=np.array(y_train)\r\n",
    "\r\n",
    "sen, spe, ppv, npv, acc, auc, _=roc_results(1-np.array(y_train),1-np.array(rbf_svc.decision_function(x_train)),path)\r\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "svm_pred = rbf_svc.predict(x_test)\r\n",
    "path='D:/0jiazhuangxian2020/jiazhuangxian/fig/svmtext/test'\r\n",
    "sen, spe, ppv, npv, acc, auc, _=roc_results(1-np.array(y_test),1-np.array(rbf_svc.decision_function(x_test)),path)\r\n",
    "all_sen['svm_t']=sen\r\n",
    "all_spe['svm_t']=spe"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# basic model(only image)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "def get_predictions(model, iterator):\r\n",
    "    model.eval()\r\n",
    "    images = []\r\n",
    "    labels = []\r\n",
    "    probs = []\r\n",
    "    top_pred = []\r\n",
    "    with torch.no_grad():\r\n",
    "        for (xy, _) in iterator:\r\n",
    "            x, y = xy\r\n",
    "            y_pred = model(x)\r\n",
    "            y_prob = F.softmax(y_pred, dim=-1)\r\n",
    "            _, _top_pred = torch.max(y_pred, 1)\r\n",
    "            images.append(x.cpu())\r\n",
    "            labels.append(y.cpu())\r\n",
    "            probs.append(y_prob.cpu())\r\n",
    "            top_pred.append(_top_pred.cpu())\r\n",
    "    images = torch.cat(images, dim=0)\r\n",
    "    labels = torch.cat(labels, dim=0)\r\n",
    "    probs = torch.cat(probs, dim=0)\r\n",
    "    top_pred = torch.cat(top_pred, dim=0)\r\n",
    "    return images, labels, probs, top_pred"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "basic_model.to(\"cpu\")\r\n",
    "\r\n",
    "_, labels, probs, _ =get_predictions(basic_model, test_iterator)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "pred=probs[:,1]\r\n",
    "# sen, spe, ppv, npv, acc, auc, _ = roc_results(1-labels.numpy(),pred.numpy(),path)\r\n",
    "sen, spe, ppv, npv, acc, auc, _ = roc_results(labels.numpy(),pred.numpy(),path)\r\n",
    "all_sen['img']=sen\r\n",
    "all_spe['img']=spe"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# all roc"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "plt.style.use('ggplot')\r\n",
    "plt.plot(1-all_spe['svm_c'],all_sen['svm_c'],linewidth=1.5)\r\n",
    "plt.plot(1-all_spe['svm_t'],all_sen['svm_t'],linewidth=1.5)\r\n",
    "plt.plot(1-all_spe['img'],all_sen['img'],color='goldenrod',linewidth=1.5)\r\n",
    "plt.plot(np.linspace(0, 1, 10), np.linspace(0, 1, 10),'--',label=None,color='forestgreen',linewidth=1.5)\r\n",
    "plt.xlim([-0.02, 1.0])\r\n",
    "plt.ylim([0.0, 1.02])\r\n",
    "plt.xlabel('1-Specificity')\r\n",
    "plt.ylabel('Sensitivity')\r\n",
    "plt.title('ROC Curves')\r\n",
    "plt.legend((r\"$Ours$\",r\"$SVM_C$\",r\"$DNN_I$\"),loc='lower right')\r\n",
    "# plt.axis('equal')\r\n",
    "plt.gca().set_aspect('equal', adjustable='box')\r\n",
    "\r\n",
    "plt.savefig(\"three_curve_roc11.svg\")\r\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  }
 ]
}