Custom CNN for Population Screening.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code is meant to train a custom CNN on the Chest X-ray data and measure its performance to aid in population screening. However, you can use these codes as the skeleton to make use of your own custom CNN for your task of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin with, let us define a few functions to load the data and convert them to Keras compatible targets. Lets begin loading the libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load libraries\n",
    "import cv2\n",
    "import numpy as np\n",
    "import os\n",
    "from keras.utils import np_utils\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "import time\n",
    "from keras.models import Sequential\n",
    "from keras.callbacks import ModelCheckpoint, TensorBoard\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "from keras.layers import Conv2D, Activation, Dense, MaxPooling2D, Flatten, Dropout\n",
    "from sklearn.metrics import log_loss\n",
    "from sklearn.utils import class_weight\n",
    "from keras.optimizers import SGD\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.metrics import classification_report,confusion_matrix, accuracy_score\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import average_precision_score\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We performed 5-fold cross validation at the patient level. we had train and test splits for each fold to ensure that none of the patienet information in the training data leaks into the test data. We randomly split 10% of the training data for validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define data directories\n",
    "train_data_dir = 'f1_cxr/train' #path to your data\n",
    "valid_data_dir = 'f1_cxr/valid'\n",
    "test_data_dir = 'f1_cxr/test'\n",
    "\n",
    "# declare the number of samples in each category\n",
    "nb_train_samples = 475 #  modify for your dataset\n",
    "nb_valid_samples = 53 #  modify for your dataset\n",
    "nb_test_samples = 134 # modify for your dataset\n",
    "num_classes = 2 # binary classification \n",
    "\n",
    "img_rows_orig = 224 # modify these values for your requirements \n",
    "img_cols_orig = 224 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets define functions to load and resize the training, validation and test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_training_data():\n",
    "    labels = os.listdir(train_data_dir)\n",
    "    total = len(labels)\n",
    "    X_train = np.ndarray((nb_train_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_train = np.zeros((nb_train_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating training images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_train = os.listdir(os.path.join(train_data_dir, label))\n",
    "        total = len(image_names_train)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_train:\n",
    "            img = cv2.imread(os.path.join(train_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_train[i] = img\n",
    "            Y_train[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1    \n",
    "    print(i)                \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.')\n",
    "    Y_train = np_utils.to_categorical(Y_train[:nb_train_samples], num_classes)\n",
    "    np.save('imgs_train.npy', X_train, Y_train) #save as numpy files\n",
    "    return X_train, Y_train\n",
    "    \n",
    "def load_validation_data():\n",
    "    # Load validation images\n",
    "    labels = os.listdir(valid_data_dir)\n",
    "    X_valid = np.ndarray((nb_valid_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_valid = np.zeros((nb_valid_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating validation images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_valid = os.listdir(os.path.join(valid_data_dir, label))\n",
    "        total = len(image_names_valid)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_valid:\n",
    "            img = cv2.imread(os.path.join(valid_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_valid[i] = img\n",
    "            Y_valid[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1\n",
    "    print(i)            \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.');\n",
    "    Y_valid = np_utils.to_categorical(Y_valid[:nb_valid_samples], num_classes)\n",
    "    np.save('imgs_valid.npy', X_valid, Y_valid) #save as numpy files\n",
    "    return X_valid, Y_valid\n",
    "\n",
    "def load_test_data():\n",
    "    labels = os.listdir(test_data_dir)\n",
    "    X_test = np.ndarray((nb_test_samples, img_rows_orig, img_cols_orig, 3), dtype=np.uint8)\n",
    "    Y_test = np.zeros((nb_test_samples,), dtype='uint8')\n",
    "    i = 0\n",
    "    print('-'*30)\n",
    "    print('Creating test images...')\n",
    "    print('-'*30)\n",
    "    j = 0\n",
    "    for label in labels:\n",
    "        image_names_test = os.listdir(os.path.join(test_data_dir, label))\n",
    "        total = len(image_names_test)\n",
    "        print(label, total)\n",
    "        for image_name in image_names_test:\n",
    "            img = cv2.imread(os.path.join(test_data_dir, label, image_name), cv2.IMREAD_COLOR)\n",
    "            img = np.array([img])\n",
    "            X_test[i] = img\n",
    "            Y_test[i] = j\n",
    "            if i % 100 == 0:\n",
    "                print('Done: {0}/{1} images'.format(i, total))\n",
    "            i += 1\n",
    "        j += 1\n",
    "    print(i)            \n",
    "    print('Loading done.')\n",
    "    print('Transform targets to keras compatible format.');\n",
    "    Y_test = np_utils.to_categorical(Y_test[:nb_test_samples], num_classes)\n",
    "    np.save('imgs_test.npy', X_test, Y_test) #save as numpy files\n",
    "    return X_test, Y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will define functions to resize the original images to that dimensions required for the custom CNN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_resized_training_data(img_rows, img_cols):\n",
    "\n",
    "    X_train, Y_train = load_training_data()\n",
    "    X_train = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_train[:nb_train_samples,:,:,:]])\n",
    "    \n",
    "    return X_train, Y_train\n",
    "    \n",
    "def load_resized_validation_data(img_rows, img_cols):\n",
    "\n",
    "    X_valid, Y_valid = load_validation_data()\n",
    "    X_valid = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_valid[:nb_valid_samples,:,:,:]])\n",
    "        \n",
    "    return X_valid, Y_valid   \n",
    "\n",
    "def load_resized_test_data(img_rows, img_cols):\n",
    "\n",
    "    X_test, Y_test = load_test_data()\n",
    "    X_test = np.array([cv2.resize(img, (img_rows,img_cols)) for img in X_test[:nb_test_samples,:,:,:]])\n",
    "    \n",
    "    return X_test, Y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An evaluation script has been written to compute the confusion matrix for the performance of the trained model. This function prints and plots the confusion matrix. Normalization can be applied by setting 'normalize=True'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes,\n",
    "                          normalize=False, #if true all values in confusion matrix is between 0 and 1\n",
    "                          title='Confusion matrix',\n",
    "                          cmap=plt.cm.Blues):\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=45)\n",
    "    plt.yticks(tick_marks, classes)\n",
    "\n",
    "    if normalize:\n",
    "        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "        print(\"Normalized confusion matrix\")\n",
    "    else:\n",
    "        print('Confusion matrix, without normalization')\n",
    "\n",
    "    print(cm)\n",
    "    thresh = cm.max() / 2.\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        plt.text(j, i, cm[i, j],\n",
    "                 horizontalalignment=\"center\",\n",
    "                 color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "    plt.tight_layout()\n",
    "    plt.ylabel('True label')\n",
    "    plt.xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now proceed to train a sequential CNN model on our dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_rows=224 #dimensions of image for the custom CNN\n",
    "img_cols=224\n",
    "channel = 3 #RGB\n",
    "num_classes = 2 #binary classification\n",
    "batch_size = 16 # modify based on the GPUs in your system\n",
    "num_epoch = 300 # modify depending on the model convergence with your data\n",
    "\n",
    "#load data\n",
    "X_train, Y_train = load_resized_training_data(img_rows, img_cols)\n",
    "X_valid, Y_valid = load_resized_validation_data(img_rows, img_cols)\n",
    "X_test, Y_test = load_resized_test_data(img_rows, img_cols)\n",
    "\n",
    "#print the shape of the data\n",
    "print(X_train.shape, Y_train.shape)\n",
    "print(X_valid.shape, Y_valid.shape)\n",
    "print(X_test.shape, Y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#model architecture\n",
    "model = Sequential()\n",
    "model.add(Conv2D(64, (3, 3), padding='same', input_shape=(img_rows, img_cols, channel)))\n",
    "model.add(Activation('relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "model.add(Conv2D(64, (3, 3), padding='same'))\n",
    "model.add(Activation('relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "model.add(Conv2D(128, (3, 3),padding='same'))\n",
    "model.add(Activation('relu'))\n",
    "model.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "model.add(Flatten())\n",
    "model.add(Dense(128, activation='relu'))\n",
    "model.add(Dropout(0.5))\n",
    "model.add(Dense(num_classes, activation='softmax'))\n",
    "model.summary()\n",
    "\n",
    "#fix the optimizer\n",
    "sgd = SGD(lr=1e-3, decay=1e-6, momentum=0.9, nesterov=True) #try varying this for your task and see the best fit\n",
    "\n",
    "#compile the model\n",
    "model.compile(optimizer=sgd,\n",
    "              loss='categorical_crossentropy', metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Its time to train the model. We will store only the best model weights by initializing callbacks. Also we can view the performance of our model during run-time by visualizing the performance graphs with Tensorboard. Create a log directory named 'logs' to store the training logs and a separate folder named 'weights' to store the model weights. You can visualize tensorboard graphs simply by navigating to your working directory and do:\n",
    "\n",
    "$tensorboard --logdir=path/to/log-directory/ --port 6006\n",
    "\n",
    "Then open localhost:6006 in your browser to view the performance graphs, model architecture and other parameters of your interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = 'weights/' + model.name + '.{epoch:02d}-{val_acc:.4f}.hdf5'\n",
    "checkpoint = ModelCheckpoint(filepath, monitor='val_acc', verbose=1, \n",
    "                             save_weights_only=True, save_best_only=True, mode='max', period=1)\n",
    "tensor_board = TensorBoard(log_dir='logs/', histogram_freq=0, batch_size=batch_size)\n",
    "callbacks_list = [checkpoint, tensor_board]\n",
    "\n",
    "#compute training time\n",
    "t=time.time()\n",
    "hist = model.fit(X_train, Y_train, batch_size=batch_size, \n",
    "                 callbacks=callbacks_list,\n",
    "                 epochs=num_epoch, verbose=1, \n",
    "                 shuffle=True, validation_data=[X_valid, Y_valid])\n",
    "\n",
    "#compute the training time\n",
    "print('Training time: %s' % (time.time()-t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to visualize the performance of the model in the console other than with Tensorboard, you can use the following script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_loss=hist.history['loss']\n",
    "val_loss=hist.history['val_loss']\n",
    "train_acc=hist.history['acc']\n",
    "val_acc=hist.history['val_acc']\n",
    "xc=range(num_epoch)\n",
    "\n",
    "plt.figure(1,figsize=(20,10), dpi=100)\n",
    "plt.plot(xc,train_loss)\n",
    "plt.plot(xc,val_loss)\n",
    "plt.xlabel('num of Epochs')\n",
    "plt.ylabel('loss')\n",
    "plt.title('train_loss vs val_loss')\n",
    "plt.grid(True)\n",
    "plt.legend(['train','val'])\n",
    "plt.style.use(['classic'])\n",
    "\n",
    "plt.figure(2,figsize=(20,10), dpi=100)\n",
    "plt.plot(xc,train_acc)\n",
    "plt.plot(xc,val_acc)\n",
    "plt.xlabel('num of Epochs')\n",
    "plt.ylabel('accuracy')\n",
    "plt.title('train_acc vs val_acc')\n",
    "plt.grid(True)\n",
    "plt.legend(['train','val'])\n",
    "plt.style.use(['classic'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the model is trained, load the best model weights to predict on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.load_weights('weights/sequential_7.10-0.5311.hdf5') #path to your model weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#predict on the test data\n",
    "X_test, Y_test = load_resized_test_data(img_rows, img_cols)\n",
    "print(X_test.shape, Y_test.shape)\n",
    "print('-'*30)\n",
    "print('Predicting on the test data...')\n",
    "print('-'*30)\n",
    "y_pred = model.predict(X_test, batch_size=batch_size, verbose=1)\n",
    "\n",
    "# compute the accuracy\n",
    "Test_accuracy = accuracy_score(Y_test.argmax(axis=-1),y_pred.argmax(axis=-1))\n",
    "print(\"Test_Accuracy = \",Test_accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now compute the performance metrics for the custom model with the test data. The performance metrics involve computing the ROC curves, cross-entropy loss score, average precision score, prediction probabilities and storing these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute the ROC-AUC values\n",
    "fpr = dict()\n",
    "tpr = dict()\n",
    "roc_auc = dict()\n",
    "for i in range(num_classes):\n",
    "    fpr[i], tpr[i], _ = roc_curve(Y_test[:, i], y_pred[:, i])\n",
    "    roc_auc[i] = auc(fpr[i], tpr[i])\n",
    "    \n",
    "# Compute micro-average ROC curve and ROC area\n",
    "fpr[\"micro\"], tpr[\"micro\"], _ = roc_curve(Y_test.ravel(), y_pred.ravel())\n",
    "roc_auc[\"micro\"] = auc(fpr[\"micro\"], tpr[\"micro\"])\n",
    "\n",
    "#Plot ROC curves\n",
    "plt.figure(figsize=(20,10), dpi=100)\n",
    "lw = 1\n",
    "plt.plot(fpr[1], tpr[1], color='red',\n",
    "         lw=lw, label='ROC curve (area = %0.4f)' % roc_auc[1])\n",
    "plt.plot([0, 1], [0, 1], color='black', lw=lw, linestyle='--')\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('Receiver operating characteristics')\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.show()\n",
    "\n",
    "# computhe the cross-entropy loss score\n",
    "score = log_loss(Y_test,y_pred)\n",
    "print(score)\n",
    "\n",
    "# compute the average precision score\n",
    "prec_score = average_precision_score(Y_test,y_pred)  \n",
    "print(prec_score)\n",
    "\n",
    "# transfer it back\n",
    "y_pred = np.argmax(y_pred, axis=1)\n",
    "Y_test = np.argmax(Y_test, axis=1)\n",
    "print(y_pred)\n",
    "print(Y_test)\n",
    "\n",
    "#save the predictions as a CSV file for further analysis\n",
    "np.savetxt('custom_model_y_pred.csv',y_pred,fmt='%i',delimiter = \",\")\n",
    "np.savetxt('custom_model_Y_test.csv',Y_test,fmt='%i',delimiter = \",\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us plot the confusion matrix to analyze the custom model's performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_names = ['class 0(abnormal)', 'class 1(normal)'] #decide the labels for your own data\n",
    "print(classification_report(Y_test,y_pred,target_names=target_names))\n",
    "print(confusion_matrix(Y_test,y_pred))\n",
    "cnf_matrix = (confusion_matrix(Y_test,y_pred))\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "# Plot non-normalized confusion matrix\n",
    "plt.figure(figsize=(20,10), dpi=100)\n",
    "plot_confusion_matrix(cnf_matrix, classes=target_names,\n",
    "                  title='Confusion matrix')\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}