{"cells":[{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"In this notebook we will start with training and go all the way to submission 🙂 This will put some constraint on our code - the plan is for it to be fairly easy to read and to lend itself well to modifications. Let's get started!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom pathlib import Path","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For training, we will use audio files I preprocessed to spectrograms and saved as numpy arrays. They are not compressed and should be very fast to read. Each spectrogram captures under 30 initial seconds of one of the files in train.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_DIR = Path('../input/birdcalldatasetnpy/train_resampled_npy')\nSAMPLE_RATE = 32_000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is what ~2 seconds of audio converted to a spectrogram looks like.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(np.load('../input/birdcalldatasetnpy/train_resampled_npy/aldfly/XC135454.npy')[:, :256]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To generate the spectrograms, I used [code](https://github.com/f0k/birdclef2018/blob/master/experiments/audio.py) shared by Jan Schlüter.\n\n\nLet me copy it over here, we will need it down the road when we predict on the test set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install /kaggle/input/pyfftwwheel/pyFFTW-0.12.0-cp37-cp37m-manylinux1_x86_64.whl\n\nfrom pyfftw.builders import rfft as rfft_builder\nfrom pyfftw import empty_aligned\n\ndef spectrogram(samples, sample_rate, frame_len, fps, batch=48, dtype=None,\n                bins=None, plans=None):\n    \"\"\"\n    Computes a magnitude spectrogram for a given vector of samples at a given\n    sample rate (in Hz), frame length (in samples) and frame rate (in Hz).\n    Allows to transform multiple frames at once for improved performance (with\n    a default value of 48, more is not always better). Returns a numpy array.\n    Allows to return a limited number of bins only, with improved performance\n    over discarding them afterwards. Optionally accepts a set of precomputed\n    plans created with spectrogram_plans(), required when multi-threading.\n    \"\"\"\n    if dtype is None:\n        dtype = samples.dtype\n    if bins is None:\n        bins = frame_len // 2 + 1\n    if len(samples) < frame_len:\n        return np.empty((0, bins), dtype=dtype)\n    if plans is None:\n        plans = spectrogram_plans(frame_len, batch, dtype)\n    rfft1, rfft, win = plans\n    hopsize = int(sample_rate // fps)\n    num_frames = (len(samples) - frame_len) // hopsize + 1\n    nabs = np.abs\n    naa = np.asanyarray\n    if batch > 1 and num_frames >= batch and samples.flags.c_contiguous:\n        frames = np.lib.stride_tricks.as_strided(\n                samples, shape=(num_frames, frame_len),\n                strides=(samples.strides[0] * hopsize, samples.strides[0]))\n        spect = [nabs(rfft(naa(frames[pos:pos + batch:], dtype) * win)[:, :bins])\n                 for pos in range(0, num_frames - batch + 1, batch)]\n        samples = samples[(num_frames // batch * batch) * hopsize::]\n        num_frames = num_frames % batch\n    else:\n        spect = []\n    if num_frames:\n        spect.append(np.vstack(\n                [nabs(rfft1(naa(samples[pos:pos + frame_len:],\n                                dtype) * win)[:bins:])\n                 for pos in range(0, len(samples) - frame_len + 1, hopsize)]))\n    return np.vstack(spect) if len(spect) > 1 else spect[0]\n\n\ndef create_mel_filterbank(sample_rate, frame_len, num_bands, min_freq,\n                          max_freq):\n    \"\"\"\n    Creates a mel filterbank of `num_bands` triangular filters, with the first\n    filter starting at `min_freq` and the last one stopping at `max_freq`.\n    Returns the filterbank as a matrix suitable for a dot product against\n    magnitude spectra created from samples at a sample rate of `sample_rate`\n    with a window length of `frame_len` samples.\n    \"\"\"\n    # prepare output matrix\n    input_bins = (frame_len // 2) + 1\n    filterbank = np.zeros((input_bins, num_bands))\n\n    # mel-spaced peak frequencies\n    min_mel = 1127 * np.log1p(min_freq / 700.0)\n    max_mel = 1127 * np.log1p(max_freq / 700.0)\n    spacing = (max_mel - min_mel) / (num_bands + 1)\n    peaks_mel = min_mel + np.arange(num_bands + 2) * spacing\n    peaks_hz = 700 * (np.exp(peaks_mel / 1127) - 1)\n    fft_freqs = np.linspace(0, sample_rate / 2., input_bins)\n    peaks_bin = np.searchsorted(fft_freqs, peaks_hz)\n\n    # fill output matrix with triangular filters\n    for b, filt in enumerate(filterbank.T):\n        # The triangle starts at the previous filter's peak (peaks_freq[b]),\n        # has its maximum at peaks_freq[b+1] and ends at peaks_freq[b+2].\n        left_hz, top_hz, right_hz = peaks_hz[b:b + 3]  # b, b+1, b+2\n        left_bin, top_bin, right_bin = peaks_bin[b:b + 3]\n        # Create triangular filter compatible to yaafe\n        filt[left_bin:top_bin] = ((fft_freqs[left_bin:top_bin] - left_hz) /\n                                  (top_bin - left_bin))\n        filt[top_bin:right_bin] = ((right_hz - fft_freqs[top_bin:right_bin]) /\n                                   (right_bin - top_bin))\n        filt[left_bin:right_bin] /= filt[left_bin:right_bin].sum()\n\n    return filterbank\n\ndef spectrogram_plans(frame_len, batch=48, dtype=np.float32):\n    \"\"\"\n    Precompute plans for spectrogram(), for a given frame length, batch size\n    and dtype. Returns two plans (single spectrum and batch), and a window.\n    \"\"\"\n    input_array = empty_aligned((batch, frame_len), dtype=dtype)\n    win = np.hanning(frame_len).astype(dtype)\n    return (rfft_builder(input_array[0]), rfft_builder(input_array), win)\n\nfilterbank = create_mel_filterbank(SAMPLE_RATE, 256, 80, 27.5, 10000)\n\ndef audio_to_melspec(audio):\n    spec = spectrogram(audio, SAMPLE_RATE, 256, 128)\n    return (spec @ filterbank).T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To train, we will use PyTorch. Let's put together a dataset we will be able to use to train our model.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"BS = 100\nMAX_LR = 1e-3\n\nclasses = [directory.name for directory in TRAIN_DIR.iterdir()]\ntrain_items = []\n\nfor directory in TRAIN_DIR.iterdir():\n    ebird_code = directory.name\n    for recording in directory.iterdir():\n        train_items.append((ebird_code, recording))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport torchvision\n\nclass TrainDataset(torch.utils.data.Dataset):    \n    def __getitem__(self, idx):\n        cls, path = train_items[idx]\n        example = self.get_spec(path)\n        return example, self.one_hot_encode(cls)\n    \n    def get_spec(self, path):\n        frames_per_spec = 212  \n        x = np.load(path)\n\n        specs = []\n        for _ in range(3):\n            if x.shape[1] < frames_per_spec:\n                spec = np.zeros((80, frames_per_spec))\n                start_frame = np.random.randint(frames_per_spec-x.shape[1])\n                spec[:, start_frame:start_frame+x.shape[1]] = x\n            else:\n                start_frame = int(np.random.rand() * (x.shape[1] - frames_per_spec))\n                spec = x[:, start_frame:start_frame+frames_per_spec]\n            specs.append(spec)\n        \n        return np.stack(specs).reshape(3, 80, frames_per_spec).astype(np.float32)\n    \n    def show(self, idx):\n        x = self[idx][0]\n        return plt.imshow(x.transpose(1,2,0)[:, :, 0])\n        \n    def one_hot_encode(self, cls):\n        y = classes.index(cls)\n        one_hot = np.zeros((len(classes)))\n        one_hot[y] = 1\n        return one_hot\n    def __len__(self):\n        return len(train_items)\n    \ntrain_ds = TrainDataset()\n\nimport multiprocessing\ntrain_dl = torch.utils.data.DataLoader(train_ds, batch_size=BS, num_workers=multiprocessing.cpu_count(), pin_memory=True, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are nearly ready to train! We still need a model and a training loop.\n\nFor the model, let us use an architecture based on resnet34 with pretrained weights.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import torchvision\n\npretrained_res34 = torchvision.models.resnet34(False)\npretrained_res34.load_state_dict(torch.load('../input/pretrained-pytorch/resnet34-333f7ec4.pth'))\n\nclass Model(torch.nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.bn = torch.nn.BatchNorm2d(3)\n        self.cnn = torch.nn.Sequential(*list(pretrained_res34.children())[:-2], torch.nn.AdaptiveMaxPool2d(1))\n        self.classifier = torch.nn.Sequential(*[\n            torch.nn.Linear(512, 512), torch.nn.ReLU(), torch.nn.Dropout(p=0.2), torch.nn.BatchNorm1d(512),\n            torch.nn.Linear(512, 512), torch.nn.ReLU(), torch.nn.Dropout(p=0.2), torch.nn.BatchNorm1d(512),\n            torch.nn.Linear(512, len(classes))\n        ])\n    \n    def forward(self, x):\n        max_per_example = x.view(x.shape[0], -1).max(1)[0]\n        x[max_per_example != 0] /= max_per_example[max_per_example != 0][:, None, None, None]\n        x = self.cnn(x)[:, :, 0, 0]\n        x = self.classifier(x)\n        return x\n\nmodel = Model().cuda()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the model out of the way, time to implement the training loop and start training!\n\nTo speed up the training and improve our results, let's first train just the new classifier we have created, keeping the convolutional part of our model frozen.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"for param in model.cnn.parameters(): param.requires_grad = False","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here is our training loop","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def train(num_epochs):\n    model.train()\n    criterion = torch.nn.BCEWithLogitsLoss()\n    optimizer = torch.optim.Adam(model.parameters(), 1e-3)\n    scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, num_epochs)\n    for epoch in range(num_epochs):\n        for data in train_dl:\n            inputs, labels = data[0].cuda(), data[1].cuda()\n            optimizer.zero_grad()\n\n            outputs = model(inputs)\n            loss = criterion(outputs, labels)\n\n            if np.isnan(loss.item()): raise Exception(f'!!! nan encountered in loss !!! epoch: {epoch}\\n')\n            loss.backward()\n            optimizer.step()\n            scheduler.step()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And let's train!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\ntrain(30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our classifier head is now trained! This means that when we start training the entire model, gradients update will not initally mess lower layers of our architecture too much. This technique called progressive unfreezing is extremely valuable. You can read more about it in [Universal Language Model Fine-tuning for Text Classification](https://arxiv.org/abs/1801.06146) by Jeremy Howard and Sebastian Ruder.\n\nLet's now unfreeze our model and train the entire arch.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfor param in model.cnn.parameters(): param.requires_grad = True\n\ntrain(60)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With training out of the way, its time to predict on the test set. To help us structure our work, let's use the extremely helpful [custom check phase](https://www.kaggle.com/c/birdsong-recognition/discussion/159993) shared by [Alex Shonenkov](https://www.kaggle.com/shonenkov).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\n\nTEST_PATH = Path('../input/birdsong-recognition') if os.path.exists('../input/birdsong-recognition/test_audio') else Path('../input/birdcall-check')\n\nTEST_AUDIO_PATH = TEST_PATH/'test_audio'\ntest_df = pd.read_csv(TEST_PATH/'test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our dataset that we will use for inference will need to be able to work with audio files.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import librosa\n\nclass AudioDataset(torch.utils.data.Dataset):\n    def __init__(self, items, classes, rec):\n        self.items = items\n        self.vocab = classes\n        self.rec = rec\n    \n    def __getitem__(self, idx):\n        _, rec_fn, start = self.items[idx]\n        x = self.rec[start*SAMPLE_RATE:(start+5)*SAMPLE_RATE]\n        example = self.get_specs(x)\n        return example.astype(np.float32)\n    \n    def get_specs(self, x):\n        xs = []\n        for i in range(3):\n            start_frame = int(i * 1.66 * SAMPLE_RATE)\n            xs.append(x[start_frame:start_frame+int(1.66*SAMPLE_RATE)])\n\n        specs = []\n        for x in xs:\n            specs.append(audio_to_melspec(x))\n        return np.stack(specs).reshape(3, 80, 212)\n    \n    def show(self, idx):\n        x = self[idx][0]\n        return plt.imshow(x.transpose(1,2,0)[:, :, 0])\n    \n    def __len__(self):\n        return len(self.items)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's predict on the test set and output predictions!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\n\nrow_ids = []\nall_preds = []\n\nmodel.eval()\nfor audio_id in test_df[test_df.site.isin(['site_1', 'site_2'])].audio_id.unique():\n    items = [(row.row_id, row.audio_id, int(row.seconds)-5) for idx, row in test_df[test_df.audio_id == audio_id].iterrows()]\n    rec = librosa.load(TEST_AUDIO_PATH/f'{audio_id}.mp3', sr=SAMPLE_RATE, res_type='kaiser_fast')[0]\n    test_ds = AudioDataset(items, classes, rec)\n    dl = torch.utils.data.DataLoader(test_ds, batch_size=64)\n    for batch in dl:\n        with torch.no_grad():\n            preds = model(batch.cuda()).sigmoid().cpu().detach()\n            all_preds.append(preds)\n    row_ids += [item[0] for item in items]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfor audio_id in test_df[test_df.site=='site_3'].audio_id.unique():\n    rec = librosa.load(TEST_AUDIO_PATH/f'{audio_id}.mp3', sr=SAMPLE_RATE, res_type='kaiser_fast')[0]\n    current_row = test_df[test_df.audio_id == audio_id].iloc[0] # assuming only one row per recording for site_3\n    duration = rec.shape[0] // SAMPLE_RATE\n    items = [(current_row.row_id, current_row.audio_id, start_sec) for start_sec in [0 + i * 5 for i in range(duration // 5)]]\n    test_ds = AudioDataset(items, classes, rec)\n    dl = torch.utils.data.DataLoader(test_ds, batch_size=64)\n    \n    preds_for_site = []\n    for batch in dl:\n        with torch.no_grad():\n            preds = model(batch.cuda()).sigmoid().cpu().detach()\n            preds_for_site.append(preds)\n    \n    if preds_for_site:\n        row_ids.append(current_row.row_id)\n        all_preds.append(torch.cat(preds_for_site).max(0)[0].unsqueeze(0))     ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_preds = torch.cat(all_preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nthresh = 1\nminimum_prediction_rate = 0.04\n\nwhile (all_preds > thresh).any(1).float().mean() < minimum_prediction_rate:\n    thresh -= 0.001","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results = []\n\nfor row in all_preds:\n    birds = []\n    for idx in np.where(row > thresh)[0]:\n        birds.append(classes[idx])\n    if not birds: birds = ['nocall']\n    results.append(' '.join(birds))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predicted = pd.DataFrame(data={'row_id': row_ids, 'birds': results})\npredicted.to_csv('submission.csv', index=False)","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}