{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os\nimport sklearn\n\nimport IPython.display as ipd\nimport librosa\nimport librosa.display\nimport soundfile as sf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# De la siguiente manera cargamos el archivo de audio en el Notebook, que será una serie temporal a modo de 'array'\n# con una tasa de muestreo 'sr' de 22 kHZ mono.\n\naudio_test_1 = \"../input/rfcx-species-audio-detection/train/0d25045a9.flac\"\n\n# Seleccioné como ejemplo ese audio por tener 1 especie presente\n\nx , sr = librosa.load(audio_test_1)\n\nprint(type(x), type(sr))\nprint(x.shape, sr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Si nos interesa, podemos cambiar la tasa de muestreo a 44.1 kHZ\n\nlibrosa.load(audio_test_1, sr = 44100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Para poder escuchar el audio aquí en el 'Notebook' hacemos lo siguiente:\n\nipd.Audio(audio_test_1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Mostramos el Espectrograma como una gráfica de líneas.\n# Tiempo vs Frecuencia\n\n%matplotlib inline\n\nplt.figure(figsize=(14, 5))\nlibrosa.display.waveplot(x, sr=sr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Mostramos el Espectrograma como un'Heatmap'.\n# Tiempo vs Frecuencia\n\nX = librosa.stft(x)\nXdb = librosa.amplitude_to_db(abs(X))\nplt.figure(figsize=(14, 5))\nlibrosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='hz')\nplt.colorbar()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Aplicamos la escala logarítimica para centrar la atención el los sonidos a estudiar\n\nlibrosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='log')\nplt.colorbar()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sr = 22050 # Tiempo de muestreo\nT = 5.0    # Segundos\nt = np.linspace(0, T, int(T*sr), endpoint=False) # Variable tiempo\nx = 0.5*np.sin(2*np.pi*220*t) # Onda senoidal pura a 220 Hz\n\n# Reproduciendo el audio\nipd.Audio(x, rate=sr) # Generando un NumPy 'array'\n\n# Guardando el audio\nsf.write('tono_0d25045a9_220.wav', x, sr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# El 'spectral_centroid' devolverá un 'array' con el mismo número de \n# columnas que número de 'fotogramas' presentes en la muestra de audio\n\nspectral_centroids = librosa.feature.spectral_centroid(x, sr=sr)[0]\nspectral_centroids.shape\n\n# Calculando la variable de tiempo para la visualización\n\nplt.figure(figsize=(12, 4))\nframes = range(len(spectral_centroids))\nt = librosa.frames_to_time(frames)\n\n# Normalizando el 'centroide espectral' para la visualización\n\ndef normalize(x, axis=0):\n    return sklearn.preprocessing.minmax_scale(x, axis=axis)\n\n# Trazando el 'centroide espectral' a lo largo de la forma de onda\n\nlibrosa.display.waveplot(x, sr=sr, alpha=0.4)\nplt.plot(t, normalize(spectral_centroids), color='b')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spectral_rolloff = librosa.feature.spectral_rolloff(x+0.01, sr=sr)[0]\nplt.figure(figsize=(12, 4))\nlibrosa.display.waveplot(x, sr=sr, alpha=0.4)\nplt.plot(t, normalize(spectral_rolloff), color='black')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spectral_bandwidth_2 = librosa.feature.spectral_bandwidth(x+0.01, sr=sr)[0]\nspectral_bandwidth_3 = librosa.feature.spectral_bandwidth(x+0.01, sr=sr, p=3)[0]\nspectral_bandwidth_4 = librosa.feature.spectral_bandwidth(x+0.01, sr=sr, p=4)[0]\nplt.figure(figsize=(15, 9))\nlibrosa.display.waveplot(x, sr=sr, alpha=0.4)\nplt.plot(t, normalize(spectral_bandwidth_2), color='r')\nplt.plot(t, normalize(spectral_bandwidth_3), color='g')\nplt.plot(t, normalize(spectral_bandwidth_4), color='y')\nplt.legend(('p = 2', 'p = 3', 'p = 4'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x, sr = librosa.load(audio_test_1)\n\n# Para centrarnos en un momento específico y ver el comportamiento de las frecuencias, podemos\n# mostrar la gráfica del audio en cuestión para hacerle un zoom en el momento concreto:\n\nplt.figure(figsize=(14, 5))\nlibrosa.display.waveplot(x, sr=sr)\n\n# Aplicamos el zoom, siendo 9000 y 9100 datos registrados en 'x'.\n\n# Dado que 'x' es un \"array\" de 1.323.000 valores, acotando a los valores \n# mencionados, conseguimos hacer zoom a una zona en concreta del audio.\n\n# Haciendo algunos cálculos rápidamente, si la duración total es de 60 segundos\n# y en esos 60 segundos tenemos 1.323.000 valores, el zoom de 9000 a 9100 equivale\n# al periodo de tiempo entre el segundo 0'4081 y el segundo 0'4126.\n\nn0 = 9000\nn1 = 9100\nplt.figure(figsize=(14, 5))\nplt.plot(x[n0:n1])\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Librosa también nos permite calcular fácilmente el número de veces\n# que la onda pasa por cero en el periodo de tiempo determinado\n\nzero_crossings = librosa.zero_crossings(x[n0:n1], pad=False)\nprint(sum(zero_crossings))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Para obtener una gráfica MFCC, se definen los siguientes parámetros:\n\nfs=10  # Siendo 'fs' un escalar > 0 y que define la tasa de muestreo para el eje \"y\" \nmfccs = librosa.feature.mfcc(x, sr = fs)  # Se define la variable mfccs para poder sacar el gráfico, es una matriz\nprint(mfccs.shape)\n(20, 97)\n\n# Utilizando la libreria 'librosa' sacamos el gráfico correspondiente\n\nplt.figure(figsize=(15, 7))\nlibrosa.display.specshow(mfccs, sr = sr, x_axis = 'time')  # Siendo 'sr' la tasa de muestreo utilizada para determinar la escala de tiempo en el eje \"x\".","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Otra herramienta interesante es 'chromagram', pues permite resaltar las características principales del croma\n\nhop_length=12  # La longitud del 'salto', que también se utiliza para determinar la escala de tiempo en el eje \"x\"\nchromagram = librosa.feature.chroma_stft(x, sr = sr, hop_length = hop_length)  # Se define la variable 'chromagram' para poder sacar el gráfico, es una matriz\n\n# Y se obtiene el gráfico con las variables ya definidas\n\nplt.figure(figsize=(15, 5))\nlibrosa.display.specshow(chromagram,\n                         x_axis = 'time', \n                         y_axis = 'chroma',\n                         hop_length = hop_length, \n                         cmap = 'coolwarm'  #'magma', 'gray_r', 'coolwarm'\n                         )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Trabajando en Linux desde mi ordenador, la manera de instalar el paquete 'resnest' fue utilizando:\n\n#pip install resnest\n\n# Pero para poder subir el Notebook a kaggle, la forma de implementar resnest es la siguiente\n\n!pip install resnest > /dev/null","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from pathlib import Path\nimport librosa as lb\n\nimport torch\nfrom  torch.utils.data import Dataset, DataLoader\n\nfrom tqdm.notebook import tqdm\n\nfrom resnest.torch import resnest50","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Se definen en primer lugar algunas variables y se cargan los datos de la competición para \n# implementarles transformaciones mediante unas funciones que se analizarán más adelante\n\nNUM_CLASSES = 24\nSR = 16_000\nDURATION =  60\n\nDATA_ROOT = Path(\"../input/rfcx-species-audio-detection\")\nTRAIN_AUDIO_ROOT = Path(\"../input/rfcx-species-audio-detection/train\")\nTEST_AUDIO_ROOT = Path(\"../input/rfcx-species-audio-detection/test\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Esta primera función crea 'MelSpecComputer', para generar 'melspec',\n# añadiéndole como parámetros para poder utilizar dicha función:\n\nclass MelSpecComputer:\n    def __init__(self, sr, n_mels, fmin, fmax):\n        \n        self.sr = sr  # sr: Tasa de muestreo del eje \"y\"\n        self.n_mels = n_mels  # n_mels: Número de bandas 'Mel' a generar\n        self.fmin = fmin  # fmin: Frecuencia máxima\n        self.fmax = fmax  # fmax: Frecuencia mínima\n\n    def __call__(self, y):\n\n        melspec = lb.feature.melspectrogram(y,\n                                            sr = self.sr,\n                                            n_mels = self.n_mels,\n                                            fmin = self.fmin,\n                                            fmax = self.fmax,\n                                            )\n\n        melspec = lb.power_to_db(melspec).astype(np.float32)\n        \n        return melspec","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Esta primera función se encarga de transformar 'X' en 'V', que contendrá\n# la información necesaria dentro del 'array' para pasar de \"mono\" a \"color\"\n\ndef mono_to_color(X, eps=1e-6, mean=None, std=None):\n    \n    X = np.stack([X, X, X], axis=-1)\n\n    # Estandarización\n    mean = mean or X.mean()\n    std = std or X.std()\n    X = (X - mean) / (std + eps)\n\n    # Normalización a [0, 255]\n    _min, _max = X.min(), X.max()\n\n    if (_max - _min) > eps:\n        V = np.clip(X, _min, _max)\n        V = 255 * (V - _min) / (_max - _min)\n        V = V.astype(np.uint8)\n    \n    else:\n        \n        V = np.zeros_like(X, dtype=np.uint8)\n\n    return V\n\n\n# Esta función se encarga de normalizar la imagen\n\ndef normalize(image, mean=None, std=None):\n    \n    image = image / 255.0\n    \n    if mean is not None and std is not None:\n        \n        image = (image - mean) / std\n        \n    return np.moveaxis(image, 2, 0).astype(np.float32)\n\n# Y esta determinará la longitud de los audios en función de los \n# datos dados por la competición (tmax - tmin)\n\ndef crop_or_pad(y, length, sr, is_train=True):\n    \n    if len(y) < length:\n        \n        y = np.concatenate([y, np.zeros(length - len(y))])\n        \n    elif len(y) > length:\n        \n        if not is_train:\n            \n            start = 0\n            \n        else:\n            \n            start = np.random.randint(len(y) - length)\n\n        y = y[start:start + length]\n\n    y = y.astype(np.float32, copy=False)\n\n    return y\n# Una vez definidas las 3 funciones en la celda anterior, el siguiente paso es \n# completar un nuevo dataset 'RFCXDataset' utilizando las funciones anteriores\n# a modo de preparación de los datos que se recogerán en el nuevo dataset\n\nclass RFCXDataset(Dataset):\n\n    def __init__(self, \n                 data, \n                 sr, \n                 n_mels = 128, \n                 fmin = 0, \n                 fmax = None,  \n                 is_train = False,\n                 num_classes = NUM_CLASSES, \n                 root = None, \n                 duration = DURATION) :\n\n        self.data = data\n        \n        self.sr = sr\n        self.n_mels = n_mels\n        self.fmin = fmin\n        self.fmax = fmax or self.sr//2\n\n        self.is_train = is_train\n\n        self.num_classes = num_classes\n        self.duration = duration\n        self.audio_length = self.duration*self.sr\n        \n        self.root =  root or (TRAIN_AUDIO_ROOT if self.is_train else TEST_AUDIO_ROOT)\n\n        self.wav_transfos = get_wav_transforms() if self.is_train else None\n\n        self.mel_spec_computer = MelSpecComputer(sr=self.sr, n_mels=self.n_mels, fmin=self.fmin, fmax=self.fmax)  # Función definida en la celda anterior\n\n\n    def __len__(self):\n        return len(self.data)\n    \n    def read_index(self, idx, fill_val=1.0, offset=None, use_offset=True):\n        d = self.data.iloc[idx]\n        record, species = d[\"recording_id\"], d[\"species_id\"]\n        try:\n            if use_offset and (self.duration < d[\"duration\"]+1):\n                offset = offset or np.random.uniform(1, int(d[\"duration\"]-self.duration))\n\n            y, _ = lb.load(self.root.joinpath(record).with_suffix(\".flac\").as_posix(),\n                           sr=self.sr, duration=self.duration, offset=offset)\n            \n            if self.wav_transfos is not None:\n                y = self.wav_transfos(y, self.sr)\n            y = crop_or_pad(y, self.audio_length, sr=self.sr)  # Función definida en la celda anterior\n            t = np.zeros(self.num_classes)\n            t[species] = fill_val\n            \n        except Exception as e:\n#             print(e)\n            raise ValueError()  from  e\n            y = np.zeros(self.audio_length)\n            t = np.zeros(self.num_classes)\n        \n        return y,t\n            \n        \n\n    def __getitem__(self, idx):\n\n        y, t = self.read_index(idx)\n        \n        \n        melspec = self.mel_spec_computer(y)  # Función definida en la celda anterior\n        image = mono_to_color(melspec)  # Función definida en la celda anterior\n        image = normalize(image, mean=None, std=None)  # Función definida en la celda anterior\n\n        return image, t\n# Función para obtener la duración de los audios\n\ndef get_duration(audio_name, root=TEST_AUDIO_ROOT):\n    return lb.get_duration(filename=root.joinpath(audio_name).with_suffix(\".flac\"))\n# Definición del nuevo dataset, para aplicarle posteriormente la variable definida\n# en las celdas anteriores 'RFCXDataset'\n\ndata = pd.DataFrame({\n    \"recording_id\": [path.stem for path in Path(TEST_AUDIO_ROOT).glob(\"*.flac\")],\n})\ndata[\"species_id\"] = [[] for _ in range(len(data))]\n\nprint(data.shape)\ndata[\"duration\"] = data[\"recording_id\"].apply(get_duration)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# En este último paso, utilizando el paquete 'resnest' y aplicando 'RFCXDataset' a los datos se obtiene,\n# para cada 'recording_id', la probabilidad de cada una de las 23 especies de aparecer en el audio en\n# concreto, que viene a ser el 'dataset' requerido por la competición\n\n# Se definen los parámetros necesarios\n\nTEST_BATCH_SIZE = 40\nTEST_NUM_WORKERS = 2\n# Aplicamos las funciones al 'dataset' data\n\ntest_data = RFCXDataset(data=data, sr=SR)\n\n# Utilizamos PyTorch para cargar los datos, iterando sobre 'test_data'\n\ntest_loader = DataLoader(test_data, batch_size=TEST_BATCH_SIZE, num_workers=TEST_NUM_WORKERS)\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\nnet = resnest50(pretrained=True).to(device)\nn_features = net.fc.in_features\nnet.fc = torch.nn.Linear(n_features, NUM_CLASSES)\nnet = net.to(device)\nnet.load_state_dict(torch.load(\"../input/kkiller-rfcx-species-detection-public-checkpoints/rfcx_resnest50.pth\", map_location=device))\nnet = net.eval()\nnet","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = []\nnet.eval()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with torch.no_grad():\n    for (xb, yb) in  tqdm(test_loader):\n        xb, yb = xb.to(device), yb.to(device)\n        o = net(xb)\n        o = torch.sigmoid(o) \n        preds.append(o.detach().cpu().numpy())\npreds = np.vstack(preds)\npreds.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.DataFrame(preds, columns=[f\"s{i}\" for i in range(24)])\nsub[\"recording_id\"] = data[\"recording_id\"].values[:len(sub)]\nsub = sub[[\"recording_id\"] + [f\"s{i}\" for i in range(24)]]\nprint(sub.shape)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}