{"cells":[{"metadata":{},"cell_type":"markdown","source":"We are provided with a large set of retina images taken using fundus photography under a variety of imaging conditions. A clinician has rated each image for the severity of diabetic retinopathy on a scale of 0 to 4:\n\n0 - No DR\n\n1 - Mild\n\n2 - Moderate\n\n3 - Severe\n\n4 - Proliferative DR\n\nWe are asked to build a machine learning model to speed up disease detection. There are thousands of images collected in rural areas to help identify diabetic retinopathy automatically.\n\nSubmissions are scored based on the quadratic weighted kappa, which measures the agreement between two ratings.\n"},{"metadata":{},"cell_type":"markdown","source":"## Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport cv2\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom matplotlib.cm import get_cmap\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nfrom sklearn.model_selection import train_test_split\n\nfrom scipy.stats import norm\n\nimport keras\nfrom keras import layers\nfrom keras.models import Model\nfrom keras import metrics\nfrom keras import backend as K   # 'generic' backend so ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load data"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# load train and test.csv\ntrain = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# How does the train file look like?\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(8, 6))\nsns.countplot(train['diagnosis'], ax=ax)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Almost half of the training data (1805 / 3662) has the score of 0. So this is a unbalanced dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"# How does the test file look like?\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train images\ntrain_image_list = os.listdir(\"../input/train_images\")\nprint(\"{} train images\".format(len(train_image_list)))\nprint(train_image_list[:10])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test images\ntest_image_list = os.listdir(\"../input/test_images\")\nprint(\"{} test images\".format(len(test_image_list)))\nprint(test_image_list[:10])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualization\n## Score: 0"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize random N images whose score is 0 (No DR)\ndef visualize_images(train, s, N):\n    # select image with the score of s\n    train_image_list = train.loc[train['diagnosis'] == s, 'id_code'].tolist()    \n    print(\"There are {} / 3662 images with the score of {}.\".format(len(train_image_list), s))\n\n    # random N choices of images\n    np.random.seed(42)\n    random_choice = np.random.choice(len(train_image_list), N)\n    \n    # visualize\n    nrow = int(np.floor(np.sqrt(N)))\n    fig, ax = plt.subplots(nrow, int(N/nrow), figsize=(20, 12))\n    ax = ax.flatten()\n    for n in range(N):\n        img = cv2.imread('../input/train_images/' + train_image_list[random_choice[n]] + '.png')\n        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # BGR to RGB\n        img = cv2.resize(img, (512, 512)) # resize to 512 x 512 (for now)\n        ax[n].imshow(img)\n        ax[n].set_title(train_image_list[random_choice[n]])    \n    plt.tight_layout()\n    \n# examples of images (score = 0) \nvisualize_images(train, 0, 30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"They are healthy guys. Again, almost half of the training data (1805 / 3662) has the score of 0."},{"metadata":{},"cell_type":"markdown","source":"## Score: 1"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize random N images whose score is 1 (mild)\nvisualize_images(train, 1, 30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I am not a medical professional, so I don't see the difference between score 0 and 1;("},{"metadata":{},"cell_type":"markdown","source":"## Score: 2"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize random N images whose score is 1 (mild)\nvisualize_images(train, 2, 30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK, maybe it is just a variability of the quality of the images but they do not always look similar in spite of the same score."},{"metadata":{},"cell_type":"markdown","source":"## Score: 3"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize random N images whose score is 1 (mild)\nvisualize_images(train, 3, 30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Score: 4"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize random N images whose score is 4 (Proliferative DR)\nvisualize_images(train, 4, 30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ah...some of them may look indeed not healthy..."},{"metadata":{},"cell_type":"markdown","source":"# Organize matrix"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train image data\nnpix = 224 # resize to npix x npix (for now)\nX_train = np.zeros((train.shape[0], npix, npix))\nfor i in range(train.shape[0]):\n    # load an image\n    img = cv2.imread('../input/train_images/' + train.loc[i, 'id_code'] + '.png')\n    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) \n    X_train[i, :, :] = cv2.resize(img, (npix, npix)) \n    \nprint(\"X_train shape: \" + str(np.shape(X_train)))   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test image data\nX_test = np.zeros((test.shape[0], npix, npix))\nfor i in range(test.shape[0]):\n    # load an image\n    img = cv2.imread('../input/test_images/' + test.loc[i, 'id_code'] + '.png')\n    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) \n    X_test[i, :, :] = cv2.resize(img, (npix, npix)) \n    \nprint(\"X_test shape: \" + str(np.shape(X_test))) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check blurriness of images\nAs these seems to be a variety of the quality of images, it may be worthwhile computing **laplacian variance** to quantify the blur."},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplacian_variance(images):\n    return [cv2.Laplacian(image, cv2.CV_32F).var() for image in images]\n\nlaplacian_variances = laplacian_variance(X_train.astype('uint8'))\n# print(len(laplacian_variances))\n# print(len(train['diagnosis']))\nblur_df = pd.DataFrame({'diagnosis': train['diagnosis'], 'laplacian_variance': laplacian_variances})\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\nsns.boxplot(x='diagnosis', y='laplacian_variance', data=blur_df, ax=ax)\nax.set_title('blurriness')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The quality of images vary a lot within the same class, but fortunately it is not systematically different across diagnosis."},{"metadata":{},"cell_type":"markdown","source":"# PCA"},{"metadata":{},"cell_type":"markdown","source":"Before modeling, let's check if we can find clusters of the data in a PCA subspace."},{"metadata":{"trusted":true},"cell_type":"code","source":"# normalize\nX = X_train / 255\n\n# reshape\nX = X.reshape(X.shape[0], -1)\n\n# use a subset of the data\ntrainX, valX, trainy, valy = train_test_split(X, train['diagnosis'], test_size=0.5, random_state=1220)\n\nX_decomposed = PCA(n_components=2).fit_transform(trainX)\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\ncmap = get_cmap(\"tab10\")\n\nfor i in range(5):\n    marker = \"$\" + str(i) + \"$\"\n    idx = trainy == i\n    ax.scatter(X_decomposed[idx, 0], X_decomposed[idx, 1],\n              marker=marker, color=cmap(i))\nax.set_title(\"PCA\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# TSNE"},{"metadata":{},"cell_type":"markdown","source":"How about t-SNE?"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_decomposed = TSNE(n_components=2).fit_transform(trainX)\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\n\nfor i in range(5):\n    marker = \"$\" + str(i) + \"$\"\n    idx = trainy == i\n    ax.scatter(X_decomposed[idx, 0], X_decomposed[idx, 1],\n              marker=marker, color=cmap(i))\nax.set_title(\"t-SNE\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Good news and bad news. Good news is that we can find distinct clusters separating 0 from others. Bad news is that telling a difference between 1 - 4 is very hard."},{"metadata":{},"cell_type":"markdown","source":"# VAE"},{"metadata":{},"cell_type":"markdown","source":"Highly inspired by [\"Visualizing MNIST using a Variational Autoencoder\" by Rebecca Vislay Wade](https://www.kaggle.com/rvislaywade/visualizing-mnist-using-a-variational-autoencoder)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# hyperparameters\nepochs = 10\nbatch_size = 200\nlatent_dim = 2 # dimensionality of the latent space\n\n# encoder\ninput_img = keras.Input(shape=(npix, npix, 1))\nx = layers.Conv2D(32, 3,\n                  padding='same', \n                  activation='relu')(input_img)\nx = layers.Conv2D(64, 3,\n                  padding='same', \n                  activation='relu',\n                  strides=(2, 2))(x)\nx = layers.Conv2D(64, 3,\n                  padding='same', \n                  activation='relu')(x)\nx = layers.Conv2D(64, 3,\n                  padding='same', \n                  activation='relu')(x)\n# need to know the shape of the network here for the decoder\nshape_before_flattening = K.int_shape(x)\n\nx = layers.Flatten()(x)\nx = layers.Dense(32, activation='relu')(x)\n\n# Two outputs, latent mean and (log)variance\nz_mu = layers.Dense(latent_dim)(x)\nz_log_sigma = layers.Dense(latent_dim)(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sampling function\ndef sampling(args):\n    z_mu, z_log_sigma = args\n    epsilon = K.random_normal(shape=(K.shape(z_mu)[0], latent_dim),\n                              mean=0., stddev=1.)\n    return z_mu + K.exp(z_log_sigma) * epsilon\n\n# sample vector from the latent distribution\nz = layers.Lambda(sampling)([z_mu, z_log_sigma])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# decoder\ndecoder_input = layers.Input(K.int_shape(z)[1:])\n\n# expand to npix x npix pixels\nx = layers.Dense(np.prod(shape_before_flattening[1:]),\n                 activation='relu')(decoder_input)\n\n# reshape\nx = layers.Reshape(shape_before_flattening[1:])(x)\n\n# use Conv2DTranspose to reverse the conv layers from the encoder\nx = layers.Conv2DTranspose(32, 3, padding='same', \n                           activation='relu', strides=(2, 2))(x)\nx = layers.Conv2D(1, 3, padding='same', activation='sigmoid')(x)\n\n# declare decoder\ndecoder = Model(decoder_input, x)\n\n# apply the decoder to the sample from the latent distribution\nz_decoded = decoder(z)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# loss\nclass CustomVariationalLayer(keras.layers.Layer):\n    def vae_loss(self, x, z_decoded):\n        x = K.flatten(x)\n        z_decoded = K.flatten(z_decoded)\n        \n        # reconstruction loss\n        xent_loss = keras.metrics.binary_crossentropy(x, z_decoded)\n        \n        # KL divergence\n        kl_loss = -5e-4 * K.mean(1 + z_log_sigma - K.square(z_mu) - K.exp(z_log_sigma), axis=-1)\n        \n        return K.mean(xent_loss + kl_loss)\n    \n    # adds the custom loss to the class\n    def call(self, inputs):\n        x = inputs[0]\n        z_decoded = inputs[1]\n        loss = self.vae_loss(x, z_decoded)\n        self.add_loss(loss, inputs=inputs)\n        return x\n    \n# apply the custom loss to the input images and the decoded latent distribution sample\ny = CustomVariationalLayer()([input_img, z_decoded])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# VAE compile\nvae = Model(input_img, y)\nvae.compile(optimizer='rmsprop', loss=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train VAE\ntrainX = np.reshape(trainX, (-1, npix, npix, 1))\nvalX = np.reshape(valX, (-1, npix, npix, 1))\n# print(np.shape(valX))\n# print(np.shape(trainX))\nvae.fit(x=trainX, y=None, \n        shuffle=True, epochs=epochs,\n        batch_size=batch_size, \n        validation_data=(valX, None))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# translate into the latent space\nencoder = Model(input_img, z_mu)\nX_encoded = encoder.predict(valX, batch_size=batch_size)\n\n# plot\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\n\nfor i in range(5):\n    marker = \"$\" + str(i) + \"$\"\n    idx = valy == i\n    ax.scatter(X_encoded[idx, 0], X_encoded[idx, 1],\n              marker=marker, color=cmap(i))\nax.set_title(\"VAE's latent space\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Again, it seems that it is relatively easy to separate 0 from others but not between 1 - 4."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display a 2D manifold of the digits\nn = 20  # figure with 20x20 digits\ndigit_size = npix\nfigure = np.zeros((digit_size * n, digit_size * n))\n\n# Construct grid of latent variable values\ngrid_x = norm.ppf(np.linspace(0.05, 0.95, n))\ngrid_y = norm.ppf(np.linspace(0.05, 0.95, n))\n\n# decode for each square in the grid\nfor i, yi in enumerate(grid_x):\n    for j, xi in enumerate(grid_y):\n        z_sample = np.array([[xi, yi]])\n        z_sample = np.tile(z_sample, batch_size).reshape(batch_size, 2)\n        x_decoded = decoder.predict(z_sample, batch_size=batch_size)\n        digit = x_decoded[0].reshape(digit_size, digit_size)\n        figure[i * digit_size: (i + 1) * digit_size,\n               j * digit_size: (j + 1) * digit_size] = digit\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\nax.imshow(figure, cmap='gnuplot2')\nax.set_title('generated images by VAE') ","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}