{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# BirdCLEF23 EDA with 3D embeddings\n\nIn this notebook, I will demonstrate how explore the BirdCLEF 2023 dataset using 3D projections of the embeddings produced by [Google's Bird Vocalization Classifier](https://www.kaggle.com/models/google/bird-vocalization-classifier), directly in IPython! I'll cover:\n1. How to project embeddings into 3D with [UMAP](https://umap-learn.readthedocs.io/en/latest/) and [t-SNE](https://en.wikipedia.org/wiki/T-distributed_stochastic_neighbor_embedding)\n2. How to make 3D scatterplots using plotly\n3. Discovering noisy recordings in the projection space\n4. Visualizing metadata breakdowns in 3D with color\n5. Some musings about these embeddings\n\nI was greatly inspired by [SPS444's notebook](https://www.kaggle.com/code/sps444/birdclef-2023-interactive-eda-windowed-data), which demonstrated how to use their Spotlight tool to explore datasets with 2D embeddings and other visualizations.","metadata":{}},{"cell_type":"markdown","source":"# But first: what are embeddings?\n\nEmbeddings are a rich feature representation of high dimensional data, often created by taking the output of an intermediate layer of a neural network that is trained on the data. The idea is that this intermediate layer has extracted higher level semantics and thrown away noise or other unuseful parts of the data, in order for later layers to process it more easily. What's mind blowing about this is that these represetations are often useful for doing things that aren't directly related to the task that the network was trained to do. For example, a common usecase in image processing is to train a model as a classifier on the imagenet dataset, and then use one of its intermediate layer outputs to do [image comparison](https://rom1504.medium.com/image-embeddings-ed1b194d113e). Even though the network wasn't directly trained to determine image similarity, its internal representations turn out to be very useful in doing that.\n\nIn our case, the input data is a 5 second audio recording at 32,000 sample/sec, or 160,000 samples (or 160,000 dimensions). Google's Bird Vocalization Classifier takes in this data and produces a vector with just 1280 dimensions--quite a bit of compression! Like the imagenet models, this model was also trained as a classifier, but the embeddings that it outputs should also be useful more generally for understanding the content of the data. If we can get them down even lower to just 2 or 3 dimensions, then we might be able to directly visualize the semantic information that they contain.","metadata":{}},{"cell_type":"markdown","source":"# Data preparation\n\nFirst, let's download all the metadata and embeddings for the BirdCLEF 2023 dataset, which I precomputed [here](https://www.kaggle.com/code/robbynevels/gbvc-embeddings) using V2 of the [Google's Bird Vocalization Classifier](https://www.kaggle.com/models/google/bird-vocalization-classifier).","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path\nimport torch\n\ndf = pd.read_csv('/kaggle/input/birdclef-2023/train_metadata.csv')\ndf['index'] = df.index\nembeddings = torch.load('/kaggle/input/gbvc-embeddings/embeddings.pt')","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:00:35.679613Z","iopub.execute_input":"2023-04-08T22:00:35.680055Z","iopub.status.idle":"2023-04-08T22:00:56.993092Z","shell.execute_reply.started":"2023-04-08T22:00:35.680017Z","shell.execute_reply":"2023-04-08T22:00:56.991666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These embeddings are computed on non-overlapping 5 second intervals of audio in the dataset. For example, the first recording is just under 50 seconds long, so we have 10 embeddings for it, 1 for each 5 second interval of the recording:","metadata":{}},{"cell_type":"code","source":"filename = df.filename[0]\nfilename, embeddings[filename].shape","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:00:56.995610Z","iopub.execute_input":"2023-04-08T22:00:56.996089Z","iopub.status.idle":"2023-04-08T22:00:57.005411Z","shell.execute_reply.started":"2023-04-08T22:00:56.996047Z","shell.execute_reply":"2023-04-08T22:00:57.004120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before we start exploring the embeddings directly, here's a helper function for getting more info about a 5 second recording, including the audio, mel spectrogram, and all the metadata. We can use this to dive into a particular point in the projection space.","metadata":{}},{"cell_type":"code","source":"from IPython.display import Audio\nimport torchaudio\nimport matplotlib.pyplot as plt\n\nSAMPLE_RATE = 32_000\nTRAIN_PATH = Path('/kaggle/input/birdclef-2023/train_audio')\n\ncompute_melspec = torchaudio.transforms.MelSpectrogram(\n    sample_rate=SAMPLE_RATE,\n    n_mels=128,\n    n_fft=2048, \n    hop_length=512,\n    f_min=0,\n    f_max=SAMPLE_RATE // 2,\n)\n\npower_to_db = torchaudio.transforms.AmplitudeToDB(\n    stype=\"power\",\n    top_db=80.0,\n)\n\ndef show_bird(index, start=0):\n    audio = torchaudio.load(TRAIN_PATH / df.filename[index], start, start+32_000*5)[0][0]\n    display(Audio(audio, rate=SAMPLE_RATE))\n    plt.figure(figsize=(12, 2.5))\n    plt.subplot(121)\n    plt.plot(audio)\n    plt.gca().get_xaxis().set_visible(False)\n    plt.subplot(122)\n    plt.imshow(power_to_db(compute_melspec(audio)))\n    plt.show()\n    return df.iloc[index]\n\n\nshow_bird(0)","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:00:57.006856Z","iopub.execute_input":"2023-04-08T22:00:57.007294Z","iopub.status.idle":"2023-04-08T22:01:00.610283Z","shell.execute_reply.started":"2023-04-08T22:00:57.007258Z","shell.execute_reply":"2023-04-08T22:01:00.609038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analyzing the first 5 seconds of all recordings\n\nTo start with, let's take just the first embedding of every recording, which just contains information about only the first 5 seconds of audio.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfirst5sec_embeddings = np.stack([embeddings[filename][0] for filename in df.filename])\nfirst5sec_embeddings.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:01:23.884313Z","iopub.execute_input":"2023-04-08T22:01:23.885754Z","iopub.status.idle":"2023-04-08T22:01:24.002079Z","shell.execute_reply.started":"2023-04-08T22:01:23.885697Z","shell.execute_reply":"2023-04-08T22:01:24.000251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To make it easier to see these embeddings, we can \"project\" each 1280-dimensional vector into a 3 dimensional point using dimensionality reduction techniques. There are [many ways to do this](https://en.wikipedia.org/wiki/Nonlinear_dimensionality_reduction), here I will just explore UMAP since it's fast, but I also included the code below for using t-SNE. Note that this may take a minute or two!","metadata":{}},{"cell_type":"code","source":"def compute_umap(df, embeddings):\n    import umap\n    \n    reducer = umap.UMAP(\n        random_state=42,\n        n_components=3\n    )\n    umap_embedding = reducer.fit_transform(embeddings)\n    df['umap_x'] = umap_embedding[:, 0]\n    df['umap_y'] = umap_embedding[:, 1]\n    df['umap_z'] = umap_embedding[:, 2]\n    \n\ndef compute_tsne(df, embeddings):\n    from sklearn.manifold import TSNE\n    \n    tsne_embedding = TSNE(\n        random_state=42,\n        n_components=3,\n        verbose=0,\n        perplexity=40,\n        n_iter=300\n    ).fit_transform(embeddings)\n    df['tsne_x'] = tsne_embedding[:, 0]\n    df['tsne_y'] = tsne_embedding[:, 1]\n    df['tsne_z'] = tsne_embedding[:, 2]\n    \n\nprint('computing umap...')\ncompute_umap(df, first5sec_embeddings)\n#print('computing tsne...')\n#compute_tsne(df, first5sec_embeddings)\nprint('done!')","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:01:25.938173Z","iopub.execute_input":"2023-04-08T22:01:25.938688Z","iopub.status.idle":"2023-04-08T22:02:42.502655Z","shell.execute_reply.started":"2023-04-08T22:01:25.938647Z","shell.execute_reply":"2023-04-08T22:02:42.500987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's throw these onto a 3D scatter plot! `plotly.express` makes this really easy with its integration with pandas dataframes. We can set the color with primary label of the bird, and also show the index on hover so it's easy to look up the recording later.","metadata":{}},{"cell_type":"code","source":"import plotly.express as px\n\nfig = px.scatter_3d(\n    title=\"Google Bird Classifier embeddings on the first 5 seconds of audio of every recording in BirdCLEF 2023, projected to 3D with UMAP\",\n    data_frame=df, \n    x='umap_x', \n    y='umap_y', \n    z='umap_z', \n    color='primary_label',\n    hover_data=['index'],\n    width=None, \n    height=800,\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:08:32.900865Z","iopub.execute_input":"2023-04-08T23:08:32.901305Z","iopub.status.idle":"2023-04-08T23:08:34.151672Z","shell.execute_reply.started":"2023-04-08T23:08:32.901267Z","shell.execute_reply":"2023-04-08T23:08:34.150390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Pan, zoom, and explore! If you find a particular point of interest, hover over it and copy the index into the `show_bird(<index>)` function. Here's a random outlier that I found particularly interesting :)","metadata":{}},{"cell_type":"code","source":"show_bird(15440)","metadata":{"execution":{"iopub.status.busy":"2023-04-08T22:02:47.033200Z","iopub.execute_input":"2023-04-08T22:02:47.033801Z","iopub.status.idle":"2023-04-08T22:02:49.800483Z","shell.execute_reply.started":"2023-04-08T22:02:47.033762Z","shell.execute_reply":"2023-04-08T22:02:49.799240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Discovering noise\n\nOne of the first things I noticed in the plot above was a couple of weird line/string-like structures containing a variety of different birds all in a row. What could these be? Let's zoom in on one and find out. We can create boxes around them by copying the coordinates shown when hovering over the endpoints of the lines:","metadata":{}},{"cell_type":"code","source":"def add_cube_mesh(fig, p0, p1):\n    x0, y0, z0 = p0\n    x1, y1, z1 = p1\n    fig.add_mesh3d(\n        x=[x0, x0, x0, x0, x1, x1, x1, x1], \n        y=[y0, y0, y1, y1, y0, y0, y1, y1], \n        z=[z0, z1, z0, z1, z0, z1, z0, z1],\n        alphahull=0,\n        opacity=0.1,\n        color='blue'\n    )\n    \ndef is_inside_box(x, y, z, p0, p1):\n    return (\n        (x > p0[0]) & (x < p1[0]) &\n        (y > p0[1]) & (y < p1[1]) &\n        (z > p0[2]) & (z < p1[2])\n    )\n\nmargin = 0.2\ncorner0 = np.array([6.798, -1.25, 6.083]) - margin\ncorner1 = np.array([7.77, -0.8935, 6.62]) + margin\non_line_1 = is_inside_box(df.umap_x, df.umap_y, df.umap_z, corner0, corner1)\n\ncorner2 = np.array([12.383, -1.007, 2.989]) - margin\ncorner3 = np.array([12.721, -1.239, 3.301]) + margin\non_line_2 = is_inside_box(df.umap_x, df.umap_y, df.umap_z, corner2, corner3)\n\nfig = px.scatter_3d(\n    data_frame=df[on_line_1 | on_line_2], \n    x='umap_x',\n    y='umap_y',\n    z='umap_z',\n    color='primary_label',\n    hover_data=['index'],\n)\nadd_cube_mesh(fig, corner0, corner1)\nadd_cube_mesh(fig, corner2, corner3)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-08T23:10:37.728231Z","iopub.execute_input":"2023-04-08T23:10:37.728743Z","iopub.status.idle":"2023-04-08T23:10:38.064901Z","shell.execute_reply.started":"2023-04-08T23:10:37.728696Z","shell.execute_reply":"2023-04-08T23:10:38.063947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, let's sample some audio data from each line and compare them to the rest of the dataset:","metadata":{}},{"cell_type":"code","source":"import torchvision\n\ndef get_specs(df):\n    specs = []\n    audios = []\n    for filename in df.filename.sample(16):\n        audio = torchaudio.load(TRAIN_PATH / filename, 0, 32_000*5)[0][0]\n        if len(audio) < 32_000*5:\n            audio = torch.concat([audio, torch.zeros(32_000*5 - len(audio))])\n        audios.append(audio)\n        spec = power_to_db(compute_melspec(audio))[None, :]\n        specs.append(spec)\n    return torchvision.utils.make_grid(specs, nrow=4)[0], audios\n\n\nnormal_specs, normal_audios = get_specs(df[~(on_line_1 | on_line_2)])\nline1_specs, line1_audios = get_specs(df[on_line_1])\nline2_specs, line2_audios = get_specs(df[on_line_2])\n\nplt.figure(figsize=(15, 5))\n\nfor i, (specs, title) in enumerate([\n    (normal_specs, 'Recordings outside of weird lines'),\n    (line1_specs, 'Recordings on line 1'),\n    (line2_specs, 'Recordings on line 2'),\n]):\n    plt.subplot(231 + i)\n    plt.axis('off')\n    plt.title(title)\n    plt.imshow(specs)\n    plt.colorbar()\n\nfor i, audios in enumerate([normal_audios, line1_audios, line2_audios]):\n    plt.subplot(234 + i)\n    plt.axis('off')\n    plt.ylim(-.2, .2)\n    for audio in audios:\n        plt.plot(audio[32000:33000], alpha=0.5)","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:11:10.034316Z","iopub.execute_input":"2023-04-08T23:11:10.034747Z","iopub.status.idle":"2023-04-08T23:11:12.127592Z","shell.execute_reply.started":"2023-04-08T23:11:10.034707Z","shell.execute_reply":"2023-04-08T23:11:12.126607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, these structures contain mostly silence and/or noise! It may be worth filtering out this part of the audio when we're training or doing validation. I'll save this list of recordings to \"outliers.csv\" if you'd like to look into them more (see the output section of this notebook).","metadata":{}},{"cell_type":"code","source":"df['outlier_type'] = on_line_1 * 1 + on_line_2 * 2\noutliers = df[df.outlier_type != 0]\noutliers.to_csv('outliers.csv')\nprint('Saved list of', len(outliers), 'recordings to outliers.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:12:52.944879Z","iopub.execute_input":"2023-04-08T23:12:52.946307Z","iopub.status.idle":"2023-04-08T23:12:52.960799Z","shell.execute_reply.started":"2023-04-08T23:12:52.946249Z","shell.execute_reply":"2023-04-08T23:12:52.959031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Breaking down metadata (song vs call example)\n\nOne fun way to explore metadata in the dataset is to map additional metadata to colors, shapes, or sizes of points when plotting embeddings. If the metadata is somehow contained within the embeddings, then we should expect to recognize larger structures based on the characteristics of the points. Let's test this idea with an example.\n\nMany of the descriptions of each recording (inside the \"type\" column) mention either the words \"song\" or \"call\". It turns out that there's [quite a distiction](https://en.wikipedia.org/wiki/Bird_vocalization) between these two kinds of vocalizations. But is this distinction clearly represented in the embeddings too? Let's focus in on just the top 5 most common birds in the dataset to find out! First, we can add a column to our dataframe to indicate whether the description contains \"song\" or \"call\":","metadata":{}},{"cell_type":"code","source":"def get_song_call(t):\n    has_call = 'call' in t\n    has_song = 'song' in t\n    if has_call and has_song:\n        return 'both'\n    elif has_call:\n        return 'call'\n    elif has_song:\n        return 'song'\n    return 'other'\n\ndf['song_call'] = df.type.map(get_song_call)\ndf.song_call.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:14:29.038444Z","iopub.execute_input":"2023-04-08T23:14:29.039861Z","iopub.status.idle":"2023-04-08T23:14:29.057879Z","shell.execute_reply.started":"2023-04-08T23:14:29.039809Z","shell.execute_reply":"2023-04-08T23:14:29.056752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we can get just the top 5 most common birds, and then rerun UMAP on only their embeddings in order to refocus our 3D projection space just for them:","metadata":{}},{"cell_type":"code","source":"top5_labels = set(df.primary_label.value_counts()[:5].index)\ntop5_df = df[df.primary_label.isin(top5_labels)].reset_index()\ncompute_umap(top5_df, first5sec_embeddings[np.array(top5_df['index'])])","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:14:31.798472Z","iopub.execute_input":"2023-04-08T23:14:31.799377Z","iopub.status.idle":"2023-04-08T23:14:43.750862Z","shell.execute_reply.started":"2023-04-08T23:14:31.799325Z","shell.execute_reply":"2023-04-08T23:14:43.749744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And now we can indeed see clearly defined clusters for each common bird:","metadata":{}},{"cell_type":"code","source":"px.scatter_3d(\n    data_frame=top5_df, \n    x='umap_x', \n    y='umap_y', \n    z='umap_z',\n    color='primary_label',\n    hover_data=['index'],\n    width=None, \n    height=800,\n).show()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:14:43.753058Z","iopub.execute_input":"2023-04-08T23:14:43.753423Z","iopub.status.idle":"2023-04-08T23:14:43.847297Z","shell.execute_reply.started":"2023-04-08T23:14:43.753389Z","shell.execute_reply":"2023-04-08T23:14:43.846114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Coloring the points by `song_call` instead reveals very separable structures *within* each bird cluster! In fact, \"thrnig1\" in particular has completely different clusters between songs and calls.","metadata":{}},{"cell_type":"code","source":"px.scatter_3d(\n    data_frame=top5_df, \n    x='umap_x', \n    y='umap_y', \n    z='umap_z',\n    color='song_call',\n    hover_data=['primary_label', 'index'],\n    width=None, \n    height=800,\n).show()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:14:43.848846Z","iopub.execute_input":"2023-04-08T23:14:43.849265Z","iopub.status.idle":"2023-04-08T23:14:43.960309Z","shell.execute_reply.started":"2023-04-08T23:14:43.849233Z","shell.execute_reply":"2023-04-08T23:14:43.958923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also see those noisy recordings here too, still showing up in a weirdly-straight line:","metadata":{}},{"cell_type":"code","source":"px.scatter_3d(\n    data_frame=top5_df, \n    x='umap_x', \n    y='umap_y', \n    z='umap_z',\n    color='outlier_type',\n    hover_data=['primary_label', 'index'],\n).show()","metadata":{"execution":{"iopub.status.busy":"2023-04-08T23:18:46.889096Z","iopub.execute_input":"2023-04-08T23:18:46.889591Z","iopub.status.idle":"2023-04-08T23:18:46.987085Z","shell.execute_reply.started":"2023-04-08T23:18:46.889548Z","shell.execute_reply":"2023-04-08T23:18:46.986070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Musing\n\nOne thing I find particularly interesting about the 3D projection of the embedding space is that it feels like it has very strict hierarchical organization: **bird class** is the primary way that clusters are formed, with other semantic information is organized *within* each of those clusters. Although we're only looking at 3D projections of a much higher dimensional space, the projection techniques try to preserve the structure of the high dimensions, so I think it's likely that the embeddings themselves also contain this emphasis on bird type, maybe through the magnitudes of certain dimensions. \n\nThe feature extractor was trained as part of a classifier model, so it's not particularly surprsing that the features would primarily emphasize distinctions between classes. But I wonder whether unsupervised models would do the same; and if not, maybe they'd produce more generalizable semantic representations. Alternatively, maybe there's a way erase or dampen the class information in the embeddings somehow, in order to reveal other clusters of inforation, such as features of bird vocalizations that are invariant to bird class, or background noises.\n\n\n# Conclusion\n\nI hope this exploration was as interesting to you as it was for me! This really only scratches the surface of what we can find here--there's so many directions to go from here. What other metadata would be interesting to visualize? How about the recording quality, or model score? What other ways can we use these embeddings? Feel free to copy this code and it try out!","metadata":{}}]}