{"cells":[{"metadata":{"papermill":{"duration":0.045074,"end_time":"2021-01-13T22:39:03.980489","exception":false,"start_time":"2021-01-13T22:39:03.935415","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

Human Protein Atlas - Single Cell Classification

\n

Similarity & Duplicate Detections Using RAPIDS and KNN

\n
CREATED BY: DARIEN SCHETTLER
\n"},{"metadata":{"papermill":{"duration":0.042275,"end_time":"2021-01-13T22:39:04.065992","exception":false,"start_time":"2021-01-13T22:39:04.023717","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

TABLE OF CONTENTS

\n\n---\n\n

0    IMPORTS

\n\n---\n\n

1    BACKGROUND INFORMATION

\n\n---\n\n

2    SETUP

\n\n---\n\n

3    HELPER FUNCTIONS

\n\n---\n\n

4    FEATURE EMBEDDING

\n\n---"},{"metadata":{"papermill":{"duration":0.043617,"end_time":"2021-01-13T22:39:04.153416","exception":false,"start_time":"2021-01-13T22:39:04.109799","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

0  IMPORTS

"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.execute_input":"2021-01-13T22:39:04.243237Z","iopub.status.busy":"2021-01-13T22:39:04.242538Z","iopub.status.idle":"2021-01-13T22:39:39.252127Z","shell.execute_reply":"2021-01-13T22:39:39.252932Z"},"papermill":{"duration":35.05654,"end_time":"2021-01-13T22:39:39.253146","exception":false,"start_time":"2021-01-13T22:39:04.196606","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"print(\"\\n... IMPORTS STARTING ...\\n\")\nprint(\"\\n\\tVERSION INFORMATION\")\n\n# Machine Learning and Data Science Imports\nimport tensorflow as tf; print(f\"\\t\\t– TENSORFLOW VERSION: {tf.__version__}\");\nimport pandas as pd; pd.options.mode.chained_assignment = None;\nimport numpy as np; print(f\"\\t\\t– NUMPY VERSION: {np.__version__}\");\nimport torch\n\n# Built In Imports\nfrom collections import Counter\nfrom datetime import datetime\nimport multiprocessing\nfrom glob import glob\nimport warnings\nimport requests\nimport imageio\nimport IPython\nimport urllib\nimport zipfile\nimport pickle\nimport random\nimport shutil\nimport string\nimport math\nimport tqdm\nimport time\nimport gzip\nimport ast\nimport io\nimport os\nimport gc\nimport re\n\n# Visualization Imports\nfrom matplotlib.colors import ListedColormap\nimport matplotlib.patches as patches\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nimport plotly.express as px\nimport seaborn as sns\nfrom PIL import Image\nimport matplotlib; print(f\"\\t\\t– MATPLOTLIB VERSION: {matplotlib.__version__}\");\nimport plotly\nimport PIL\nimport cv2\n\n# Submission Imports\nimport typing as t\nimport base64\nimport zlib\n\n# PRESETS\nLBL_NAMES = [\"Nucleoplasm\", \"Nuclear Membrane\", \"Nucleoli\", \"Nucleoli Fibrillar Center\", \"Nuclear Speckles\", \"Nuclear Bodies\", \"Endoplasmic Reticulum\", \"Golgi Apparatus\", \"Intermediate Filaments\", \"Actin Filaments\", \"Microtubules\", \"Mitotic Spindle\", \"Centrosome\", \"Plasma Membrane\", \"Mitochondria\", \"Aggresome\", \"Cytosol\", \"Vesicles\", \"Negative\"]\nINT_2_STR = {x:LBL_NAMES[x] for x in np.arange(19)}\nINT_2_STR_LOWER = {k:v.lower().replace(\" \", \"_\") for k,v in INT_2_STR.items()}\nSTR_2_INT_LOWER = {v:k for k,v in INT_2_STR_LOWER.items()}\nSTR_2_INT = {v:k for k,v in INT_2_STR.items()}\nFIG_FONT = dict(family=\"Helvetica, Arial\", size=14, color=\"#7f7f7f\")\nLABEL_COLORS = [px.colors.label_rgb(px.colors.convert_to_RGB_255(x)) for x in sns.color_palette(\"Spectral\", len(LBL_NAMES))]\nLABEL_COL_MAP = {str(i):x for i,x in enumerate(LABEL_COLORS)}\n\nprint(\"\\n\\n... IMPORTS COMPLETE ...\\n\")\n\n\nLIMIT = 12\ngpus = tf.config.experimental.list_physical_devices('GPU')\nif gpus:\n try:\n tf.config.experimental.set_virtual_device_configuration(\n gpus[0],\n [tf.config.experimental.VirtualDeviceConfiguration(memory_limit=1024*LIMIT)])\n logical_gpus = tf.config.experimental.list_logical_devices('GPU')\n print(len(gpus), \"Physical GPUs,\", len(logical_gpus), \"Logical GPUs\")\n except RuntimeError as e:\n print(e)\n\nprint('Restrict TensorFlow to max %iGB GPU RAM'%LIMIT)\nprint('so RAPIDS can use %iGB GPU RAM'%(16-LIMIT))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import cudf, cuml, cupy\nfrom cuml.neighbors import NearestNeighbors\n\nprint('RAPIDS', cuml.__version__)","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.042806,"end_time":"2021-01-13T22:39:39.338543","exception":false,"start_time":"2021-01-13T22:39:39.295737","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

1  BACKGROUND INFORMATION

"},{"metadata":{"papermill":{"duration":0.042139,"end_time":"2021-01-13T22:39:39.423376","exception":false,"start_time":"2021-01-13T22:39:39.381237","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

1.1 THE DATA

\n\n---\n\nBACKGROUND INFORMATION\n\nWIP"},{"metadata":{"papermill":{"duration":0.042063,"end_time":"2021-01-13T22:39:39.507915","exception":false,"start_time":"2021-01-13T22:39:39.465852","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

1.2 THE GOAL

\n\n---\n\nWIP"},{"metadata":{"papermill":{"duration":0.042459,"end_time":"2021-01-13T22:39:39.678624","exception":false,"start_time":"2021-01-13T22:39:39.636165","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

2  NOTEBOOK SETUP

"},{"metadata":{"execution":{"iopub.execute_input":"2021-01-13T22:39:39.779248Z","iopub.status.busy":"2021-01-13T22:39:39.778548Z","iopub.status.idle":"2021-01-13T22:39:40.76653Z","shell.execute_reply":"2021-01-13T22:39:40.767021Z"},"papermill":{"duration":1.046203,"end_time":"2021-01-13T22:39:40.767178","exception":false,"start_time":"2021-01-13T22:39:39.720975","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"# Define the root data directory\nROOT_DIR = \"/kaggle/input\"\nDATA_DIR = os.path.join(ROOT_DIR, \"hpa-tfrecords-512x512-rgb-only\")\n\n# Define the paths to the training and testing tfrecord and image folders respectively\nTRAIN_IMG_DIR = os.path.join(ROOT_DIR, \"hpa-single-cell-image-classification\", \"train\")\nTRAIN_TFREC_DIR = os.path.join(DATA_DIR, \"train_slide_records\")\n\n# Capture all the relevant full tfrec paths\nTRAIN_TFREC_PATHS = tf.io.gfile.glob(os.path.join(TRAIN_TFREC_DIR, '*.tfrec'))\nprint(f\"\\n... The number of training tfrecord files is {len(TRAIN_TFREC_PATHS)} ...\\n\")\n\nTRAIN_CSV_PATH = \"../input/hpa-single-cell-image-classification/train.csv\"\ntrain_df = pd.read_csv(TRAIN_CSV_PATH)\n\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.043842,"end_time":"2021-01-13T22:39:40.85577","exception":false,"start_time":"2021-01-13T22:39:40.811928","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

3  HELPER FUNCTIONS

"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.execute_input":"2021-01-13T22:39:40.980411Z","iopub.status.busy":"2021-01-13T22:39:40.979421Z","iopub.status.idle":"2021-01-13T22:39:40.982404Z","shell.execute_reply":"2021-01-13T22:39:40.981881Z"},"papermill":{"duration":0.080032,"end_time":"2021-01-13T22:39:40.982523","exception":false,"start_time":"2021-01-13T22:39:40.902491","status":"completed"},"tags":[],"trusted":true},"cell_type":"code","source":"def load_image_from_channel_paths(channel_paths):\n rgby = [np.asarray(Image.open(path), np.uint8) for path in channel_paths]\n return np.stack(rgby, axis=-1)\n\ndef load_image(img_id, img_dir):\n \"\"\" Load An Image Using ID and Directory Path - Composes 4 Individual Images \"\"\"\n rgby = [\n np.asarray(Image.open(os.path.join(img_dir, img_id+f\"_{c}.png\")), np.uint8) \\\n for c in [\"red\", \"green\", \"blue\", \"yellow\"]\n ]\n return np.stack(rgby, axis=-1)\n\ndef convert_rgby_to_rgb(arr):\n \"\"\" Convert a 4 channel (RGBY) image to a 3 channel RGB image.\n \n Advice From Competition Host/User: lnhtrang\n\n For annotation (by experts) and for the model, I guess we agree that individual \n channels with full range px values are better. \n In annotation, we toggled the channels. \n For visualization purpose only, you can try blending the channels. \n For example, \n - red = red + yellow\n - green = green + yellow/2\n - blue=blue.\n \n Args:\n arr (numpy array): The RGBY, 4 channel numpy array for a given image\n \n Returns:\n RGB Image\n \"\"\"\n \n rgb_arr = np.zeros_like(arr[..., :-1])\n rgb_arr[..., 0] = arr[..., 0]\n rgb_arr[..., 1] = arr[..., 1]+arr[..., 3]/2\n rgb_arr[..., 2] = arr[..., 2]\n \n return rgb_arr\n\ndef plot_ex(arr, figsize=(20,6), title=None, plot_merged=True, rgb_only=False):\n \"\"\" Plot 4 Channels Side by Side \"\"\"\n if plot_merged and not rgb_only:\n n_images=5 \n elif plot_merged and rgb_only:\n n_images=4\n elif not plot_merged and rgb_only:\n n_images=4\n else:\n n_images=3\n plt.figure(figsize=figsize)\n if type(title) == str:\n plt.suptitle(title, fontsize=20, fontweight=\"bold\")\n\n for i, c in enumerate([\"Red Channel – Microtubles\", \"Green Channel – Protein of Interest\", \"Blue - Nucleus\", \"Yellow – Endoplasmic Reticulum\"]):\n if not rgb_only:\n ch_arr = np.zeros_like(arr[..., :-1]) \n else:\n ch_arr = np.zeros_like(arr)\n if c in [\"Red Channel – Microtubles\", \"Green Channel – Protein of Interest\", \"Blue - Nucleus\"]:\n ch_arr[..., i] = arr[..., i]\n else:\n if rgb_only:\n continue\n ch_arr[..., 0] = arr[..., i]\n ch_arr[..., 1] = arr[..., i]\n plt.subplot(1,n_images,i+1)\n plt.title(f\"{c.title()}\", fontweight=\"bold\")\n plt.imshow(ch_arr)\n plt.axis(False)\n \n if plot_merged:\n plt.subplot(1,n_images,n_images)\n \n if rgb_only:\n plt.title(f\"Merged RGB\", fontweight=\"bold\")\n plt.imshow(arr)\n else:\n plt.title(f\"Merged RGBY into RGB\", fontweight=\"bold\")\n plt.imshow(convert_rgby_to_rgb(arr))\n plt.axis(False)\n \n plt.tight_layout(rect=[0, 0.2, 1, 0.97])\n plt.show()\n \n \n\ndef flatten_list_of_lists(l_o_l):\n return [item for sublist in l_o_l for item in sublist]\n\n\ndef get_class_wts(df, low_idx=4):\n label_counts = Counter([c for sublist in df.Label.str.split(\"|\").to_list() for c in sublist])\n low_val = sorted(label_counts.values())[low_idx-1] # Not the lowest as it is very underrepresented\n class_wts ={int(k):min(1.0, low_val/v) for k,v in label_counts.items()}\n return {i:class_wts[i] for i in sorted(class_wts)}\n\n\ndef decode_image(image_data, n_channels=1, resize_to=(512,512), cast_to=tf.uint8):\n image = tf.image.decode_png(image_data, channels=n_channels) \n image = tf.image.resize(image, resize_to) \n return tf.cast(image, cast_to)\n\n\ndef str_2_multi_hot_encoding(tfstring, n_classes=19):\n ragged_indices = tf.strings.to_number(tf.strings.split(tfstring, sep=\"|\"), out_type=tf.int32)\n one_hot_stack = tf.one_hot(ragged_indices, depth=n_classes)\n return tf.reduce_max(one_hot_stack, axis=-2)\n\n\ndef decode(serialized_example, multihot=False, n_channels=1, resize_to=(512,512)):\n \"\"\" Parses a set of features and label from the given `serialized_example`.\n \n It is used as a map function for `dataset.map`\n\n Args:\n serialized_example (tf.Example): A serialized example\n is_test (bool, optional): Whether to allow for the label feature\n \n Returns:\n A decoded tf.data.Dataset object representing the tfrecord dataset\n \"\"\"\n # Defaults are not specified since both keys are required.\n feature_dict = {\n 'image': tf.io.FixedLenFeature(shape=(), dtype=tf.string),\n 'image_name': tf.io.FixedLenFeature(shape=(), dtype=tf.string),\n 'target': tf.io.FixedLenFeature(shape=(), dtype=tf.string),\n }\n # Define a parser\n features = tf.io.parse_single_example(serialized_example, features=feature_dict) \n image = decode_image(features['image'], n_channels, resize_to)\n image_name = features[\"image_name\"]\n if multihot:\n label = str_2_multi_hot_encoding(features[\"target\"])\n else:\n label = features[\"target\"]\n return image, image_name, label\n\n\ndef preprocess_tfrec_ds(red, green, blue, yellow=None, drop_yellow=True, return_id=True):\n if (yellow is None) or (drop_yellow):\n (ri, rn, rl), (gi, gn, gl), (bi, bn, bl) = red, green, blue\n yi, yn, yl = None, None, None\n else:\n (ri, rn, rl), (gi, gn, gl), (bi, bn, bl), (yi, yn, yl) = red, green, blue, yellow\n \n if yi is None:\n combo_img = tf.stack([ri[..., 0], gi[..., 0], bi[..., 0]], axis=-1)\n else:\n combo_img = tf.stack([ri[..., 0], gi[..., 0], bi[..., 0], yi[..., 0]], axis=-1)\n \n if return_id:\n img_id = tf.strings.substr(rn, pos=0, len=36) # 36 is length of id (always)\n return combo_img, img_id, rl\n else:\n return combo_img, rl\n \ndef augment(img_batch, lbl_batch):\n # SEEDING & KERNEL INIT\n K = tf.random.uniform((1,), minval=0, maxval=4, dtype=tf.dtypes.int32)[0]\n\n img_batch = tf.image.random_flip_left_right(img_batch)\n img_batch = tf.image.random_flip_up_down(img_batch)\n img_batch = tf.image.rot90(img_batch, K)\n \n img_batch = tf.image.random_saturation(img_batch, 0.875, 1.125)\n img_batch = tf.image.random_brightness(img_batch, 0.1125)\n img_batch = tf.image.random_contrast(img_batch, 0.825, 1.175)\n\n return img_batch, lbl_batch","execution_count":null,"outputs":[]},{"metadata":{"papermill":{"duration":0.143736,"end_time":"2021-01-13T22:50:29.721487","exception":false,"start_time":"2021-01-13T22:50:29.577751","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

4  FEATURE EMBEDDING

\n\nWIP"},{"metadata":{"papermill":{"duration":0.144463,"end_time":"2021-01-13T22:50:30.008724","exception":false,"start_time":"2021-01-13T22:50:29.864261","status":"completed"},"tags":[]},"cell_type":"markdown","source":"

4.1 DATASET

\n\n---\n\nWIP"},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 64\n\ntrain_ds = tf.data.TFRecordDataset(TRAIN_TFREC_PATHS, num_parallel_reads=None)\ntrain_ds = train_ds.map(lambda x: decode(x, n_channels=3, multihot=True, resize_to=(256,256)))\n\n# See examples\nfor i, (img, image_name, lbl) in enumerate(train_ds.take(3)):\n print(f\"IMAGE SHAPE : {img.shape}\")\n print(f\"IMG #{i+1} -- IMAGE NAME : {image_name.numpy().decode()}\")\n print(f\"IMAGE LABEL : {lbl}\\n\")\n plt.imshow(img.numpy().astype(np.uint8))\n plt.show()\n \n# get ordered image ids \nimg_ids = list(train_ds.map(lambda x,y,z: (y)).as_numpy_iterator())\nimg_ids = [x.decode() for x in img_ids]\n\ntrain_ds = train_ds.map(lambda x,y,z: (x)).batch(BATCH_SIZE).prefetch(tf.data.AUTOTUNE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

4.2 MODEL

\n\n---\n\nWIP"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_full_model(backbone, preprocessing_fn, input_shape=(256,256,3)):\n inputs = tf.keras.layers.Input(input_shape, dtype=tf.uint8)\n prep_inputs = preprocessing_fn(tf.cast(inputs, tf.float32))\n backbone_outputs = backbone(prep_inputs)\n return tf.keras.Model(inputs=[inputs], outputs=[backbone_outputs])\n\nbb = tf.keras.applications.EfficientNetB0(weights='imagenet', include_top=False, pooling=\"avg\")\npp = tf.keras.applications.efficientnet.preprocess_input\nmodel = get_full_model(bb, pp, input_shape=(256,256,3))\nmodel.summary()\n\nmodel.save(\"/kaggle/working/eb0_embedding_model_at_256.h5\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

4.3 FIND SIMILAR IMAGES WITH RAPIDS KNN AND IMAGE EMBEDDINGS

\n\n---\n\nWIP"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_embeddings = model.predict(train_ds,verbose=1)\nprint('image embeddings shape is', image_embeddings.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"IMG_KNN_N = 75\nmodel = NearestNeighbors(n_neighbors=IMG_KNN_N)\nmodel.fit(image_embeddings)\ndistances, indices = model.kneighbors(image_embeddings)\n\nprint(distances.shape, indices.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = []\nfor k in tqdm(range(len(distances)), total=len(distances)):\n IDX = np.where(distances[k,]<2.5)[0]\n IDS = indices[k,IDX]\n preds.append(list(np.array(img_ids)[IDS]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds_map = {p[0]:p[1:] for p in preds if len(p)>1}\nkeys_to_remove = []\nfor pred_key, pred_vals in tqdm(preds_map.copy().items(), total=len(preds_map)):\n if pred_key not in keys_to_remove:\n keys_to_remove += [v for v in pred_vals if v in preds_map.keys()]\n \nfor k in keys_to_remove:\n try:\n _ = preds_map.pop(k)\n except:\n pass\n \nlen(preds_map)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"MAX_DUPS = 2\n\nfor root_id, duplicate_ids in tqdm(preds_map.items(), total=len(preds_map)):\n print(\"\\n------------------------------------------------------------------------\\n\")\n plot_ex(load_image(root_id, TRAIN_IMG_DIR), title=f\"\\nROOT\\nID : {train_df[train_df.ID==root_id].ID.values[0]}\\nLBL : {train_df[train_df.ID==root_id].Label.values[0]}\")\n print(\"\\n------------------------------------------------------------------------\\n\")\n for i, _id in enumerate(duplicate_ids):\n plot_ex(load_image(_id, TRAIN_IMG_DIR), title=f\"\\nDUPLICATE\\nID : {train_df[train_df.ID==_id].ID.values[0]}\\nLBL : {train_df[train_df.ID==_id].Label.values[0]}\")\n if i==2:\n print(f\"\\n... ONLY DISPLAYING {MAX_DUPS} DUPLICATES FOR BREVITY ...\\n\")\n break","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}