End2End_PipeLine.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "End2End_PipeLine.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "2mZF2ZLWu_kC"
      },
      "source": [
        "path = \"\"\n",
        "from google.colab import drive, files\n",
        "drive.mount('/content/drive/')\n",
        "path = '/content/drive/My Drive/App/HPA' #\"/content/drive/My Drive/HPA_Validation/\"\n",
        "\n",
        "%tensorflow_version 1.x"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ShLwOReUUMsd"
      },
      "source": [
        "#!pip install mahotas"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rFvlL1Xvv2uL"
      },
      "source": [
        "import mahotas\n",
        "from skimage.feature import hog\n",
        "from skimage.io import imread\n",
        "from skimage.transform import rescale\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import cv2\n",
        "import os\n",
        "from sklearn.decomposition import PCA\n",
        "from skimage.transform import resize\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "\n",
        "def fd_hu_moments(image):\n",
        "    image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)\n",
        "    feature = cv2.HuMoments(cv2.moments(image)).flatten()\n",
        "    return feature\n",
        "    \n",
        "def fd_haralick(image): \n",
        "    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)\n",
        "    haralick = mahotas.features.haralick(gray).mean(axis=0)\n",
        "    return haralick\n",
        "\n",
        "def fd_hog(image):    \n",
        "    img = resize(image, (1024, 1024))\n",
        "    hog_feature = hog(img, orientations=9, pixels_per_cell=(16, 16), cells_per_block=(2, 2), \n",
        "                                 block_norm='L2-Hys', visualize=False, transform_sqrt=False, \n",
        "                                 feature_vector=True, multichannel=True)\n",
        "    \n",
        "    \n",
        "    hog_feature = hog_feature.reshape(hog_feature.shape[0],1)\n",
        "    scaler = MinMaxScaler()\n",
        "    data_rescaled = scaler.fit_transform(hog_feature)\n",
        "    pca = PCA(n_components = 0.95)\n",
        "    pca.fit(data_rescaled)\n",
        "    hog_feature_pca = pca.transform(data_rescaled)\n",
        "    hog_feature_pca = hog_feature_pca.flatten()\n",
        "    return hog_feature\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J905bMSIMagQ"
      },
      "source": [
        "from keras.applications.vgg16 import VGG16\n",
        "from keras.models import Model, Sequential\n",
        "from keras.layers import Dense,Conv2D,Flatten,MaxPooling2D,GlobalAveragePooling2D,Lambda\n",
        "from keras.utils import plot_model\n",
        "from keras import regularizers\n",
        "from keras.optimizers import Adam,SGD,Nadam\n",
        "\n",
        "#from scipy.misc import logsumexp\n",
        "from scipy.special import logsumexp\n",
        "\n",
        "\n",
        "def Vgg16_Net():\n",
        "    # load model without classifier layers\n",
        "    model = VGG16(include_top=False, input_shape=(1024, 1024, 3))\n",
        "    # add new classifier layers\n",
        "    flat1 = Flatten()(model.layers[-1].output)\n",
        "    class1 = Dense(256, activation='relu')(flat1)\n",
        "    output = Dense(8, activation='sigmoid')(class1)\n",
        "    model = Model(inputs=model.inputs, outputs=output)\n",
        "    model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])\n",
        "    return model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kw8YJ08URo4f"
      },
      "source": [
        "layer_name = 'dense'\n",
        "FC_layer_model = Model(inputs=Vgg16_Net.input, outputs=Vgg16_Net.get_layer(layer_name).output)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QoDJRC51O-Sd"
      },
      "source": [
        "### HU, Haralick, HOG and Vgg16 CNN Feature Extraction ###\n",
        "import os\n",
        "import cv2\n",
        "from PIL import Image\n",
        "from numpy import array\n",
        "from scipy.sparse import csr_matrix\n",
        "import scipy\n",
        "from skimage.transform import resize\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "IMG = []\n",
        "hog_images = []\n",
        "hog_features = []\n",
        "global_features = []\n",
        "features_hog = []\n",
        "features_haralick = []\n",
        "features_hu_moments = []\n",
        "all_features_hand_cnn = []\n",
        "all_features = []\n",
        "features_cnn = []\n",
        "features_hand = []\n",
        "p_hat = []\n",
        "batch_features_hand = []\n",
        "read = lambda imname: np.asarray(Image.open(imname).convert(\"RGB\"))\n",
        "fixed_size = (3000,3000)\n",
        "\n",
        "PATH = \"Image file location\" # example: '/content/drive/My Drive/App/HPA/FrenchLab/Image01.jpg'\n",
        "\n",
        "image = cv2.imread(PATH)\n",
        "image = cv2.resize(image, fixed_size)\n",
        "crop_img = image[300:2700, 300:2700]   \n",
        "       \n",
        "##############################################\n",
        "# DCNN Features\n",
        "###############################################\n",
        "img = resize(image, (1024,1024))\n",
        "images = np.array(img) \n",
        "FC_output = FC_layer_model.predict(np.expand_dims(images,0))\n",
        "FC_output = FC_output.reshape(FC_output.shape[1],FC_output.shape[0])\n",
        "\n",
        "#Data = pd.DataFrame(FC_output)\n",
        "#Data.to_csv(\"csv file name\") #'/content/drive/My Drive/App/HPA/FrenchLab/Image01_CNN_Features_X.csv') \n",
        "\n",
        "####################################\n",
        "# Handcrafted Feature extraction\n",
        "####################################\n",
        "fv_hu_moments = fd_hu_moments(image) \n",
        "fv_haralick   = fd_haralick(image)\n",
        "fv_hog        = fd_hog(image) \n",
        "\n",
        "fv_hu_moments= fv_hu_moments.reshape(fv_hu_moments.shape[1],fv_hu_moments.shape[0])\n",
        "fv_haralick= fv_haralick.reshape(fv_haralick.shape[1],fv_haralick.shape[0])\n",
        "\n",
        "features_hand = np.concatenate ((fv_hu_moments,fv_haralick,fv_hog), axis=1)\n",
        "features = np.hstack((features_hand, FC_output)) \n",
        "\n",
        "file = PATH + '.csv'\n",
        "Data = pd.DataFrame(features)\n",
        "Data.to_csv(\"file\") #'/content/drive/My Drive/App/HPA/FrenchLab/Image01_Features_X.csv') #,index=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3EpnsHVhRbF0"
      },
      "source": [
        "def categorical_entropy(proba):\r\n",
        "    \"\"\"Compute the entropy along the last dimension.\"\"\"\r\n",
        "    return - np.sum(proba * np.log(proba + 1), axis=-1)\r\n",
        "\r\n",
        "def hamming_score(y_true, y_pred, normalize=True, sample_weight=None):\r\n",
        "    acc_list = []\r\n",
        "    for i in range(y_true.shape[0]):\r\n",
        "        set_true = set( np.where(y_true[i])[0] )\r\n",
        "        set_pred = set( np.where(y_pred[i])[0] )\r\n",
        "        tmp_a = None\r\n",
        "        if len(set_true) == 0 and len(set_pred) == 0:\r\n",
        "            tmp_a = 1\r\n",
        "        else:\r\n",
        "            tmp_a = len(set_true.intersection(set_pred))/\\\r\n",
        "                    float( len(set_true.union(set_pred)))\r\n",
        "        acc_list.append(tmp_a)\r\n",
        "    return np.mean(acc_list)\r\n",
        "\r\n",
        "def evaluate(classes, y_gt, y_pred):\r\n",
        "   \r\n",
        "    print(\"******* Evaluation Score******************\")\r\n",
        "\r\n",
        "    y_pred_bin = y_pred #>= threshold_value\r\n",
        "\r\n",
        "    score_f1_macro = f1_score(y_gt, y_pred_bin, average=\"macro\")\r\n",
        "    print(\"Macro f1_socre = {:.6f}\".format(score_f1_macro))\r\n",
        "\r\n",
        "    score_f1_micro = f1_score(y_gt, y_pred_bin, average=\"micro\")\r\n",
        "    print(\"Micro f1_socre = {:.6f}\".format(score_f1_micro))\r\n",
        "\r\n",
        "    # hamming loss\r\n",
        "    h_loss = hamming_loss(y_gt, y_pred_bin)\r\n",
        "    print(\"Hamming Loss = {:.6f}\".format(h_loss))\r\n",
        "\r\n",
        "    print('Hamming score: {0}'.format(hamming_score(y_gt, y_pred_bin)))\r\n",
        "\r\n",
        "    # Exact Match ratio:\r\n",
        "    print(\"Exact Match ratio:\",np.all(y_pred_bin == y_gt, axis=1).mean())\r\n",
        "\r\n",
        "    mAP = average_precision_score(y_gt, y_pred)\r\n",
        "    print(\"mean Average Precision(mAP) = {:.2f}%\".format(mAP * 100))\r\n",
        "    ap_classes = []\r\n",
        "    for i, cls in enumerate(classes):\r\n",
        "        ap_cls = average_precision_score(y_gt[:, i], y_pred[:, i])\r\n",
        "        ap_classes.append(ap_cls)\r\n",
        "        print(\"AP({}) = {:.2f}%\".format(cls, ap_cls * 100))\r\n",
        "    print(\"mAP = {:.2f}%\".format(np.mean(ap_classes) * 100))\r\n",
        "\r\n",
        "    print(\"Test Sample Size:\", y_gt.shape)\r\n",
        "\r\n",
        "    print(multilabel_confusion_matrix(y_gt, y_pred))\r\n",
        "\r\n",
        "def mcc(out_proba,y_test):\r\n",
        "    threshold = np.arange(0.01,0.9,0.01)\r\n",
        "    acc = []\r\n",
        "    accuracies = []\r\n",
        "    best_threshold = np.zeros(out_proba.shape[1])\r\n",
        "    for i in range(out_proba.shape[1]):\r\n",
        "        y_prob = np.array(out_proba[:,i])\r\n",
        "        for j in threshold:\r\n",
        "            y_pred = [1 if prob>=j else 0 for prob in y_prob]\r\n",
        "            acc.append( matthews_corrcoef(y_test[:,i],y_pred))\r\n",
        "        acc   = np.array(acc)\r\n",
        "        index = np.where(acc==acc.max()) \r\n",
        "        accuracies.append(acc.max()) \r\n",
        "        best_threshold[i] = threshold[index[0][0]]\r\n",
        "        acc = []\r\n",
        "\r\n",
        "    print(\"Best Threshold\",best_threshold[:])\r\n",
        "    \r\n",
        "    y_pred = np.array([[1 if out_proba[i,j]>=best_threshold[j] else 0 for j in range(y_test.shape[1])] for i in range(len(y_test))])\r\n",
        "    \r\n",
        "    total_correctly_predicted = len([i for i in range(len(y_test)) if (y_test[i]==y_pred[i]).sum() == 8])\r\n",
        "    print(\"total_correctly_predicted\",total_correctly_predicted)\r\n",
        "    \r\n",
        "    return y_pred, best_threshold\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pltMJi2CX0L4"
      },
      "source": [
        "\"\"\"Multilabel cross validators based on an implementation of the algorithm described in the following paper:\r\n",
        "Sechidis K., Tsoumakas G., Vlahavas I. (2011) On the Stratification of Multi-Label Data. \r\n",
        "Machine Learning and Knowledge Discovery in Databases. ECML PKDD 2011. Lecture Notes in Computer Science, vol 6913. Springer, Berlin, Heidelberg.\r\n",
        "https://github.com/trent-b/iterative-stratification/blob/master/iterstrat/ml_stratifiers.py\r\n",
        "\"\"\"\r\n",
        "\r\n",
        "import numpy as np\r\n",
        "\r\n",
        "from sklearn.utils import check_random_state\r\n",
        "from sklearn.utils.validation import _num_samples, check_array\r\n",
        "from sklearn.utils.multiclass import type_of_target\r\n",
        "from sklearn.model_selection._split import _BaseKFold, _RepeatedSplits, BaseShuffleSplit, _validate_shuffle_split\r\n",
        "\r\n",
        "class MultilabelStratifiedShuffleSplit(BaseShuffleSplit):\r\n",
        "    \"\"\"Multilabel Stratified ShuffleSplit cross-validator\r\n",
        "    Examples\r\n",
        "    --------\r\n",
        "    >>> X = np.array([[1,2], [3,4], [1,2], [3,4], [1,2], [3,4], [1,2], [3,4]])\r\n",
        "    >>> y = np.array([[0,0], [0,0], [0,1], [0,1], [1,1], [1,1], [1,0], [1,0]])\r\n",
        "    >>> msss = MultilabelStratifiedShuffleSplit(n_splits=3, test_size=0.2, random_state=0)\r\n",
        "    >>> for train_index, test_index in msss.split(X, y):\r\n",
        "            print(\"TRAIN:\", train_index, \"TEST:\", test_index)\r\n",
        "            X_train, X_test = X[train_index], X[test_index]\r\n",
        "            y_train, y_test = y[train_index], y[test_index]\r\n",
        "    ---        \r\n",
        "    TRAIN: [1 2 5 7] TEST: [0 3 4 6]\r\n",
        "    TRAIN: [2 3 6 7] TEST: [0 1 4 5]\r\n",
        "    TRAIN: [1 2 5 6] TEST: [0 3 4 7]\r\n",
        "    \"\"\"\r\n",
        "\r\n",
        "    def __init__(self, n_splits=10, test_size=\"default\", train_size=None,\r\n",
        "                 random_state=None):\r\n",
        "        super(MultilabelStratifiedShuffleSplit, self).__init__(\r\n",
        "            n_splits, test_size, train_size, random_state)\r\n",
        "\r\n",
        "    def _iter_indices(self, X, y, groups=None):\r\n",
        "        n_samples = _num_samples(X)\r\n",
        "        y = check_array(y, ensure_2d=False, dtype=None)\r\n",
        "        y = np.asarray(y, dtype=bool)\r\n",
        "        type_of_target_y = type_of_target(y)\r\n",
        "\r\n",
        "        if type_of_target_y != 'multilabel-indicator':\r\n",
        "            raise ValueError(\r\n",
        "                'Supported target type is: multilabel-indicator. Got {!r} instead.'.format(\r\n",
        "                    type_of_target_y))\r\n",
        "\r\n",
        "        n_train, n_test = _validate_shuffle_split(n_samples, self.test_size,\r\n",
        "                                                  self.train_size)\r\n",
        "\r\n",
        "        n_samples = y.shape[0]\r\n",
        "        rng = check_random_state(self.random_state)\r\n",
        "        y_orig = y.copy()\r\n",
        "\r\n",
        "        r = np.array([n_train, n_test]) / (n_train + n_test)\r\n",
        "\r\n",
        "        for _ in range(self.n_splits):\r\n",
        "            indices = np.arange(n_samples)\r\n",
        "            rng.shuffle(indices)\r\n",
        "            y = y_orig[indices]\r\n",
        "\r\n",
        "            test_folds = IterativeStratification(labels=y, r=r, random_state=rng)\r\n",
        "\r\n",
        "            test_idx = test_folds[np.argsort(indices)] == 1\r\n",
        "            test = np.where(test_idx)[0]\r\n",
        "            train = np.where(~test_idx)[0]\r\n",
        "\r\n",
        "            yield train, test\r\n",
        "\r\n",
        "    def split(self, X, y, groups=None):\r\n",
        "        y = check_array(y, ensure_2d=False, dtype=None)\r\n",
        "        return super(MultilabelStratifiedShuffleSplit, self).split(X, y, groups)\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "cellView": "code",
        "id": "KnwpftYZZWon"
      },
      "source": [
        "def DropWeights(x):\r\n",
        "    import tensorflow as tf\r\n",
        "    p = 0.3\r\n",
        "    return tf.nn.dropout(x, rate = p) * (1 - p)\r\n",
        "\r\n",
        "es = EarlyStopping(monitor='val_acc', mode='max', verbose=1, patience=100)\r\n",
        "mc = ModelCheckpoint('/content/drive/My Drive/App/HPA/Best_HNetModel.h5', monitor='val_acc', mode='max', verbose=1, save_best_only=True)\r\n",
        "reduce_lr = ReduceLROnPlateau(monitor='val_acc', factor=0.2, patience=2, min_lr=0.000001, verbose=1, mode=\"auto\", cooldown=0)\r\n",
        "\r\n",
        "def main():\r\n",
        "  \r\n",
        "    kf = MultilabelStratifiedShuffleSplit(n_splits=10, test_size=0.2, random_state=0)\r\n",
        "\r\n",
        "    for train_index, test_index in kf.split(X, Y):    \r\n",
        "        X_train, X_test = X[train_index], X[test_index]\r\n",
        "        Y_train, Y_test = Y[train_index], Y[test_index]\r\n",
        "\r\n",
        "        model = keras.Sequential()\r\n",
        "        model.add(Dense(2048, input_shape=(3996,),kernel_initializer=\"he_uniform\", activation='relu'))\r\n",
        "        model.add(Dense(2048, kernel_initializer=\"he_uniform\", activation='relu'))\r\n",
        "        model.add(Dense(1024, kernel_initializer=\"he_uniform\", activation='relu'))\r\n",
        "        model.add(Dense(512,  kernel_initializer=\"he_uniform\", activation='relu'))\r\n",
        "        model.add(Lambda(DropWeights))\r\n",
        "        model.add(Dense(8, activation='sigmoid'))        \r\n",
        "        adam = Adam(lr= 0.000001, amsgrad=True)\r\n",
        "\r\n",
        "        model.compile(adam, 'binary_crossentropy', metrics=['accuracy'])\r\n",
        "        \r\n",
        "        model.fit(X_train, Y_train, epochs= 250, batch_size=32, validation_split = 0.2, verbose=0, callbacks=[es, mc, reduce_lr])     \r\n",
        "\r\n",
        "        from keras.models import load_model\r\n",
        "        model = load_model (\"/content/drive/My Drive/App/HPA/Best_HNetModel.h5\")  \r\n",
        "\r\n",
        "        num_samples = 1000\r\n",
        "        ens_preds = []\r\n",
        "        for seed in range(num_samples):\r\n",
        "            ens_preds.append(model.predict(X_test))\r\n",
        "        Yt_hat = np.array(ens_preds)\r\n",
        "\r\n",
        "        prediction = np.mean(Yt_hat, axis = 0)\r\n",
        "        Y_predict = mcc(prediction,Y_test)\r\n",
        "\r\n",
        "        evaluate(classes, Y_test, Y_predict)\r\n",
        "\r\n",
        "main()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "17aZ9Gq9brvM"
      },
      "source": [
        "class_names = ['Leydig Cells', 'Elongated/late Spermatids Cells', 'Pachytene Spermatocytes Cells', 'Peritubular Cells', \r\n",
        "           'Preleptotene spermatocytes Cells', 'Round/early Spermatids Cells', 'Sertoli Cells', 'Spermatogonia Cells']\r\n",
        "classes = class_names\r\n",
        "cell_type = 0 # 0: Leydig, .. 7: Spermatogonia\r\n",
        "T = 1000\r\n",
        "\r\n",
        "for t in range(T):\r\n",
        "    p_hat.append(model.predict(X_test,verbose=0))\r\n",
        "p_hat = np.array(p_hat)\r\n",
        "\r\n",
        "\r\n",
        "y_gt = Y_test[:,cell_type:cell_type+1]\r\n",
        "y_pred = y_predict[:,cell_type:cell_type+1]\r\n",
        "\r\n",
        "tn, fp, fn, tp = confusion_matrix(y_gt.flatten(),y_pred.flatten()).ravel()\r\n",
        " \r\n",
        "# Plug-in estimate of Entropy    \r\n",
        "proba_avg = p_hat.mean(axis=0)\r\n",
        "entropy_expected = categorical_entropy(proba_avg)\r\n",
        "\r\n",
        "# Jackknife correction for Bias\r\n",
        "proba_loo = proba_avg + (proba_avg - p_hat) / (len(p_hat) - 1)\r\n",
        "expected_entropy_p = categorical_entropy(proba_loo).mean(axis=0)\r\n",
        "entropy_expected += (len(p_hat) - 1) * (entropy_expected - expected_entropy_p)\r\n",
        "\r\n",
        "MCDW_samples = entropy_expected - categorical_entropy(p_hat).mean(axis=0)\r\n",
        "MCDW_samples = np.array(MCDW_samples)\r\n",
        "\r\n",
        "#uncertainty_norm = (MCDW_samples - np.min(MCDW_samples)) / (np.max(MCDW_samples) - np.min(MCDW_samples))\r\n",
        "\r\n",
        "# For Multi-Class Classification\r\n",
        "#proba_samples_sorted = np.sort(np.mean(MCDW_samples, axis = 0))[::-1]\r\n",
        "#CPPD = proba_samples_sorted[:,-1] - proba_samples_sorted[:, -2]\r\n",
        "\r\n",
        "# For Multi-Label Classification\r\n",
        "CPPD = (proba_avg[:,0:1])\r\n",
        "CPPD = 2*CPPD.flatten() -1\r\n",
        "GTL = CPPD.flatten()/MCDW_samples.flatten()\r\n",
        "\r\n",
        "GTL_samples =  (GTL - np.min(GTL)) / (np.max(GTL) - np.min(GTL))\r\n",
        "GTL_samples = np.array(GTL_samples)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}