{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install /kaggle/input/pyfftwwheel/pyFFTW-0.12.0-cp37-cp37m-manylinux1_x86_64.whl","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nfrom pathlib import Path\nimport pandas as pd\nimport librosa\nimport matplotlib.pyplot as plt\nimport soundfile as sf\nfrom torch.utils.data import Dataset, DataLoader\nimport numpy as np\nimport torch\nimport torchvision\nfrom torch import nn\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\n\nfrom pyfftw.builders import rfft as rfft_builder\nfrom pyfftw import empty_aligned\n\nMEL_BANDS=80\nMEL_MIN=27.5\nMEL_MAX=10000\nSAMPLE_RATE=32_000","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"classes = pd.read_pickle('../input/espstarterpackv3res34minmax/classes.pkl')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def lme_pool(x, alpha=1.0): # log-mean-exp pool\n    '''alpha -> approximates maxpool, alpha -> 0 approximates mean pool'''\n    T = x.shape[1]\n    mult_log = torch.log(torch.tensor(1/T))\n    return 1/alpha * (mult_log + torch.logsumexp((alpha * x), dim=1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Model(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.cnn = nn.Sequential(*list(torchvision.models.resnet34(False).children())[:-2])\n        self.classifier = nn.Sequential(*[\n            nn.Linear(512, 512), nn.ReLU(), nn.Dropout(p=0.5), nn.BatchNorm1d(512),\n            nn.Linear(512, 512), nn.ReLU(), nn.Dropout(p=0.5), nn.BatchNorm1d(512),\n            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] # scaling to between 0 and 1\n        x /= max_per_example[:, None, None, None, None]     # per example!\n        bs, im_num = x.shape[:2]\n        x = x.view(-1, x.shape[2], x.shape[3], x.shape[4])\n        x = self.cnn(x)\n        x = x.mean((2,3))\n        x = self.classifier(x)\n        x = x.view(bs, im_num, -1)\n        x = lme_pool(x)\n        return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = Model()\nmodel.load_state_dict(torch.load('../input/espstarterpackv3res34minmax/135_lmepool_simple_minmax_0.73.pth'))\nmodel.cuda()\nmodel.eval();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The below adapts [code](https://github.com/f0k/birdclef2018/blob/master/experiments/audio.py) by Jan Schlüter.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"filterbank = create_mel_filterbank(SAMPLE_RATE, 256, MEL_BANDS, MEL_MIN, MEL_MAX)\n\ndef audio_to_melspec(audio):\n    spec = spectrogram(audio, SAMPLE_RATE, 256, 128)\n    return (spec @ filterbank).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"TEST_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')\ntest_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class AudioDataset(Dataset):\n    def __init__(self, items, classes, rec, example_dur_mult=1):\n        self.items = items\n        self.vocab = classes\n        self.rec = rec\n        self.example_dur_mult = example_dur_mult\n    \n    def __getitem__(self, idx):\n        _, rec_fn, start = self.items[idx]\n        x = self.rec[start*SAMPLE_RATE:(start+(5 * self.example_dur_mult))*SAMPLE_RATE]\n        example = self.get_specs(x)\n        imgs = example.reshape(-1, 3, 80, 212)\n        return imgs.astype(np.float32)\n    \n    def get_specs(self, x):\n        xs = []\n        for i in range(self.example_dur_mult * 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)\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":{"trusted":true},"cell_type":"code","source":"%%time\n\nrow_ids = []\nall_preds = []\n\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 = 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\nEXAMPLE_DUR_MULT = 10 # example duration = 5 sec x EXAMPLE_DUR_MULT\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    mod_duration = duration - (duration % (5 * EXAMPLE_DUR_MULT))\n    items = [(current_row.row_id, current_row.audio_id, start_sec) for start_sec in [0 + i * 5 * EXAMPLE_DUR_MULT for i in range(mod_duration // (5 * EXAMPLE_DUR_MULT))]]\n    if duration > EXAMPLE_DUR_MULT * 5: items.append((current_row.row_id, current_row.audio_id, duration-(EXAMPLE_DUR_MULT * 5)))\n    test_ds = AudioDataset(items, classes, rec, EXAMPLE_DUR_MULT)\n    dl = 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.08\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}