{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Based on a kernel by Bojan Tunguz: https://www.kaggle.com/tunguz/melanoma-tsne-and-umap-embeddings-with-rapids","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"[Rapids](https://rapids.ai) is an open-source GPU accelerated Data Science and Machine Learning library, developed and mainatained by [Nvidia](https://www.nvidia.com). It is designed to be compatible with many existing CPU tools, such as Pandas, scikit-learn, numpy, etc. It enables **massive** acceleration of many data-science and machine learning tasks, oftentimes by a factor fo 100X, or more. \n\nRapids is still undergoing developemnt, and as of right now it's not availabel in the Kaggle Docker environment. If you are interested in installing and running Rapids locally on your own machine, then you should [refer to the followong instructions](https://rapids.ai/start.html).","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Overview\n\nThe idea is to visualize the space of RGB **histograms** of the tumor images. The distance/dissimilarity between the histograms is measured with the Fisher metric (more info inside the notebook). This space is embedded in 3D using UMAP, and displayed with plotly. \n\n# View of data\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-user-content/o/inbox%2F466630%2F0e09c8550d8320bca4a362ba80ce82dc%2Fumap_with_histograms.png?generation=1592226072209541&alt=media)\n\n# More information\nhttps://www.kaggle.com/c/siim-isic-melanoma-classification/discussion/158829\n\n# Interaction \nClicking on a point on the plot below shows 5 closest images (by the UMAP coords) and their histograms.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-output":true,"scrolled":false},"cell_type":"code","source":"%%time\ntry:\n    import cudf, cuml\n    print('rapids already installed')\nexcept:\n    # INSTALL RAPIDS OFFLINE (FROM KAGGLE DATASET). TAKES 1 MINUTE :-)\n    print('installing rapids (should take ~80sec)')\n    import sys\n    !cp ../input/rapids/rapids.0.13.0 /opt/conda/envs/rapids.tar.gz  2>/dev/null\n    \n    !cd /opt/conda/envs/ && tar -xzf rapids.tar.gz \n    sys.path = [\"/opt/conda/envs/rapids/lib\"] + [\"/opt/conda/envs/rapids/lib/python3.6\"] + [\"/opt/conda/envs/rapids/lib/python3.6/site-packages\"] + sys.path\n    !cp /opt/conda/envs/rapids/lib/libxgboost.so /opt/conda/lib/  2>/dev/null\n    print('done installing rapids')\n    import cudf, cuml","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"import cupy as cp\nimport numpy as np\nimport pandas as pd\nimport os\nfrom cuml.manifold import TSNE, UMAP\nimport matplotlib.pyplot as plt\nfrom matplotlib.pyplot import ylim, xlim\n\nimport plotly.express as px\n\nimport plotly.graph_objs as go\nfrom ipywidgets import Output, VBox\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"from PIL import Image\n\ndef central_crop(im):\n    w = im.shape[0]\n    nw = int(np.floor(w/np.sqrt(2)))\n    d = (w-nw)//2\n    return im[d:-d, d:-d, :]\n\ndef to_grayscale(IM):\n    return np.asarray(Image.fromarray(IM).convert('L'))    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(\"../input/siim-isic-melanoma-classification/train.csv\")\nIMS = np.load('../input/siimisic-melanoma-resized-images/x_train_64.npy')\n\nMERGE = True\n\nif MERGE:\n    IMS_test = np.load('../input/siimisic-melanoma-resized-images/x_test_64.npy')\n    IMS = np.concatenate([IMS, IMS_test])\n    test_df = pd.read_csv(\"../input/siim-isic-melanoma-classification/test.csv\")\n    \n    train_df['is_train'] = True\n    train_df = train_df.append(test_df)\n    train_df.is_train.fillna(False, inplace = True)\n    train_df.target.fillna(-1, inplace = True)\n    train_df.reset_index(inplace = True, drop = True)\n\nIMS = list(map(central_crop, IMS))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"%%time\ndef rgb_hists(IMS):\n    HS = []\n    k = 3 if len(IMS[0].shape)==3 else 1\n    for im in IMS:            \n        comps = [np.bincount(im.reshape(-1,k)[:,i], minlength=256) for i in range(k)]    \n        HS.append(np.concatenate(comps))\n        \n    return np.asarray(HS)\n    \nhists = rgb_hists(IMS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fisher metric\n\nIdeally we'd use something like the KL divergence to compare the histograms, but this version of UMAP works only with the Eclidean metric it seems. We use a sqrt-transform which maps the Fisher manifold isometrically onto the Euclidean space. In other words, computing the Euclidean distance in the transformed space (train) corresponds to computing the Fisher metric in the original space (hists). \n\nThe fisher Metric is an approximation (of the sqrt) of the KL divergence. (BTW: This sqrt is not related to the other sqrt.) More precisely, the geodesics (between two histograms/discrete probability distributions) in the Fisher space minimize the integrals of KL divergences along the curve.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"train = np.sqrt(hists)\n#train = hists","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"MEANS = np.median( np.asarray(IMS).reshape(len(IMS), -1, 3), axis = 1)\n\ntrain_df['mean_col'] = [(r/255,g/255,b/255) for r,g,b in MEANS]\ntrain_df['mean_lum'] = [r/255 + g/255 + b/255 for r,g,b in MEANS]\ntrain_df['mean_r'] = [r/255 for r,g,b in MEANS]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def plot_rgb_hist(ax, hist):\n    hs = np.split(hist,3)\n    ax.plot(hs[0], 'r')\n    ax.plot(hs[1], 'g')\n    ax.plot(hs[2], 'b')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"st = 59","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"print(f'using random_state = {st}')\numap = UMAP(n_components=3, random_state=st, n_neighbors = 12, n_epochs = 1_000)\n\nxyz = umap.fit_transform(train)\n\ntrain_df['emb_x'] = xyz[:, 0]\ntrain_df['emb_y'] = xyz[:, 1]\ntrain_df['emb_z'] = xyz[:, 2]\n\ntrain_df.to_csv('tabular_with_umap_coords', index = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"import plotly.graph_objs as go\nimport plotly.express as px\nfrom ipywidgets import Output, VBox\nfrom scipy.spatial import KDTree\n\ndef show_interactive_embedding(train_df, colors = None, sizes = None):\n    X = train_df[['emb_x', 'emb_y', 'emb_z']].values\n    KD = KDTree(X)\n    \n    train_df['ind'] = train_df.index    \n\n    sc = px.scatter_3d(train_df, x = 'emb_x', y = 'emb_y', z = 'emb_z', \n                  size = sizes,               \n                  #size = (PRED != q.target)*10 + 0.5,\n                  color = colors,\n                  #symbol = 'is_train',\n                  hover_data = train_df.columns,\n                  width = 1200, height = 1200,\n                  )\n\n    fig = go.FigureWidget(data=sc)\n\n    out = Output()\n\n    def same_patients(sel):\n        id = train_df[train_df.index == sel].patient_id.values[0]\n        ALL = train_df[train_df.patient_id == id].index.values\n        ALL[list(ALL).index(sel)], ALL[0] = ALL[0], ALL[list(ALL).index(sel)]\n        return ALL\n\n    def similar_in_umap(ind, k = 20):    \n        ns = KD.query(X[ind], k = k)[1]\n        return ns\n\n    @out.capture(clear_output=True)\n    def handle_click(trace, points, state):\n        if not points.point_inds:\n            print('handle_click received empty selection, probably a bug in plotly...')\n            return\n\n        sel = points.point_inds[0]    \n        #ALL = same_patients(sel)    \n        ALL = similar_in_umap(sel, 50)\n\n        _, axs = plt.subplots(len(ALL), 2, figsize = (10,len(ALL)*2))\n        #axs = axs.ravel()\n\n        new_sizes = sizes.copy()\n        for i,x in enumerate(ALL):\n            axs[i,1].imshow(IMS[x])\n            plot_rgb_hist(axs[i,0], hists[x])        \n            axs[i,0].set_title(f\"{train_df.at[x, 'target']}, {train_df.at[x,'ind']} {train_df.at[x, 'diagnosis']} \")\n            axs[i,0].axes.get_xaxis().set_visible(False)        \n            new_sizes[x] = 15    \n\n        fig.update_traces(marker=dict(size=new_sizes))    \n\n        plt.show()        \n\n        \n    # bug in plotly -- if colors are specified something is wrong with selections...\n    if colors is None:\n        fig.data[0].on_click(handle_click)    \n\n    return VBox([fig, out])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"q = train_df[(train_df.emb_x > 5) & (train_df.emb_y > -2) & (train_df.emb_z < -1)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"sizes_by_target = (train_df.target>0).values*10 + 0.5","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.is_train.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"q.is_train.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"#show_interactive_embedding(train_df, sizes = sizes_by_target, colors = train_df.is_train.values)\nshow_interactive_embedding(train_df, sizes = sizes_by_target, colors = None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[(train_df.emb_y >= 3) & (train_df.target>=1)].target.count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"#show_interactive_embedding(train_df, sizes = None, colors = train_df.is_train.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"#show_interactive_embedding(train_df, sizes = sizes_by_target, colors = train_df.sex.fillna('n/a').values)","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}