{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Parallelized Audio Feature Extraction with Dask\n\n## Preface\n\nThis notebook is the part of the comparative study on the best parallel computation tool to use to parallelize audio feature extraction from the bird call audio files provided as a part of this competition.\n\nIn the study, three different technologies have been evaluated\n* classical *multiprocessing* library\n* *Dask* (illustrated by this notebook)\n* *Ray*\n\n**Note**: This notebook is a supplementary demo to the discussion thread per https://www.kaggle.com/c/birdsong-recognition/discussion/179662\n\n## Dask Overview\n\nFrom the outside, Dask looks a lot like Ray. It, too, is a library for distributed parallel computing in Python, with its own task scheduling system, awareness of Python data frameworks like NumPy, and the ability to scale from one machine to many.\n\nIf looking for a musical metaphor for Dask, it resembles swing. Both in its style and its attitude. You do not have to study a lot to write your first Dask-based mulitiprocessing or distributed application – much like in composing, playing or dancing swing. As we remember, only a few of the pioneers and world starts of swing were ever formally trained in music.\n\nDask is composed of two parts\n* Dynamic task scheduling optimized for computation. This is similar to Airflow, Luigi, Celery, or Make, but optimized for interactive computational workloads.\n* “Big Data” collections like parallel arrays, dataframes, and lists that extend common interfaces like NumPy, Pandas, or Python iterators to larger-than-memory or distributed environments. These parallel collections run on top of dynamic task schedulers.\n\nDask works in two basic ways. The first is by way of parallelized data structures — essentially, Dask’s own versions of NumPy arrays, lists, or Pandas DataFrames. Whenever you swap in the Dask versions of those constructions for their defaults, and Dask will automatically spread their execution across your cluster. \nTypically, it barely involves little more than changing the name of an import. However, in some case you may find yourself rewriting your code completely to make it work.\n\nThe second way is through Dask’s low-level parallelization mechanisms, including function decorators, that parcel out jobs across nodes and return results synchronously (“immediate” mode) or asynchronously (“lazy”). Both modes can be mixed as needed, too.\n\nOne key difference between Dask and Ray is the scheduling mechanism. Dask uses a centralized scheduler that handles all tasks for a cluster. Ray is, in turn, totally decentralized. It means each machine runs its own scheduler, so any issues with a scheduled task are handled at the level of the individual machine, not the whole cluster.\n\n\nWhen applying the capabilities of Dask to the test drive with the bird call recognition, we will use the second approach. The essence of the implementation was to use Dask’s low-level delayed interface.\n","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from dask import delayed\nimport datetime as dt\nimport pandas as pd\nimport numpy as np\nimport librosa\nfrom typing import List, Dict, Tuple\nimport warnings\n\nwarnings.filterwarnings('ignore')\n\n# config settings\nin_kaggle = True\nNUMBER_OF_MFCC = 20\n\nNUMBER_OF_CPU_IN_POOL = 6\n\n# Fourier Transform Settings\n# Default FFT window size\nN_FFT = 2048  # FFT window size\nHOP_LENGTH = 512  # number audio of frames between STFT columns (looks like a good default)\n\n# data output settings\nTRANSFORMED_DATA_PATH = \"/kaggle/working\"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def get_data_file_path(is_in_kaggle: bool) -> Tuple[str, str]:\n    train_path = ''\n    test_path = ''\n\n    if is_in_kaggle:\n        # running in Kaggle, inside the competition\n        train_path = '../input/birdsong-recognition/train.csv'\n        test_path = '../input/birdsong-recognition/test.csv'\n    else:\n        # running locally\n        train_path = 'data/train.csv'\n        test_path = 'data/test.csv'\n\n    return train_path, test_path\n\ndef get_base_train_audio_folder_path(is_in_kaggle: bool) -> str:\n    folder_path = ''\n    if is_in_kaggle:\n        folder_path = '../input/birdsong-recognition/train_audio/'\n    else:\n        folder_path = 'data/train_audio/'\n    return folder_path\n\ndef extract_feautres(trial_audio_file_path: str):\n    # process data frame\n    function_start_time = dt.datetime.now()\n    print(\"Started a file processing at \", function_start_time)\n\n    df0 = extract_feature_means(trial_audio_file_path)\n\n    function_finish_time = dt.datetime.now()\n    print(\"Fininished the file processing at \", function_finish_time)\n\n    processing = function_finish_time - function_start_time\n    print(\"Processed the file: \", trial_audio_file_path, \"; processing time: \", processing)\n\n    return df0\n\n# inspirations: https://musicinformationretrieval.com/basic_feature_extraction.html\ndef extract_feature_means(audio_file_path: str) -> pd.DataFrame:\n    # config settings\n    number_of_mfcc = NUMBER_OF_MFCC\n\n    # 1. Importing 1 file\n    y, sr = librosa.load(audio_file_path)\n\n    # Trim leading and trailing silence from an audio signal (silence before and after the actual audio)\n    signal, _ = librosa.effects.trim(y)\n\n    # 2. Fourier Transform\n    # Default FFT window size\n    n_fft = N_FFT  # FFT window size\n    hop_length = HOP_LENGTH  # number audio of frames between STFT columns (looks like a good default)\n\n    # Short-time Fourier transform (STFT)\n    d_audio = np.abs(librosa.stft(signal, n_fft=n_fft, hop_length=hop_length))\n\n    # 3. Spectrogram\n    # Convert an amplitude spectrogram to Decibels-scaled spectrogram.\n    db_audio = librosa.amplitude_to_db(d_audio, ref=np.max)\n\n    # 4. Create the Mel Spectrograms\n    s_audio = librosa.feature.melspectrogram(signal, sr=sr)\n    s_db_audio = librosa.amplitude_to_db(s_audio, ref=np.max)\n\n    # 5 Zero crossings\n\n    # #6. Harmonics and Perceptrual\n    # Note:\n    #\n    # Harmonics are characteristichs that represent the sound color\n    # Perceptrual shock wave represents the sound rhythm and emotion\n    y_harm, y_perc = librosa.effects.hpss(signal)\n\n    # 7. Spectral Centroid\n    # Note: Indicates where the ”centre of mass” for a sound is located and is calculated\n    # as the weighted mean of the frequencies present in the sound.\n\n    # Calculate the Spectral Centroids\n    spectral_centroids = librosa.feature.spectral_centroid(signal, sr=sr)[0]\n    spectral_centroids_delta = librosa.feature.delta(spectral_centroids)\n    spectral_centroids_accelerate = librosa.feature.delta(spectral_centroids, order=2)\n\n    # spectral_centroid_feats = np.stack((spectral_centroids, delta, accelerate))  # (3, 64, xx)\n\n    # 8. Chroma Frequencies¶\n    # Note: Chroma features are an interesting and powerful representation\n    # for music audio in which the entire spectrum is projected onto 12 bins\n    # representing the 12 distinct semitones ( or chromas) of the musical octave.\n\n    # Increase or decrease hop_length to change how granular you want your data to be\n    hop_length = HOP_LENGTH\n\n    # Chromogram\n    chromagram = librosa.feature.chroma_stft(signal, sr=sr, hop_length=hop_length)\n\n    # 9. Tempo BPM (beats per minute)¶\n    # Note: Dynamic programming beat tracker.\n\n    # Create Tempo BPM variable\n    tempo_y, _ = librosa.beat.beat_track(signal, sr=sr)\n\n    # 10. Spectral Rolloff\n    # Note: Is a measure of the shape of the signal. It represents the frequency below which a specified\n    #  percentage of the total spectral energy(e.g. 85 %) lies.\n\n    # Spectral RollOff Vector\n    spectral_rolloff = librosa.feature.spectral_rolloff(signal, sr=sr)[0]\n\n    # spectral flux\n    onset_env = librosa.onset.onset_strength(y=signal, sr=sr)\n\n    # Spectral Bandwidth¶\n    # The spectral bandwidth is defined as the width of the band of light at one-half the peak\n    # maximum (or full width at half maximum [FWHM]) and is represented by the two vertical\n    # red lines and λSB on the wavelength axis.\n    spectral_bandwidth_2 = librosa.feature.spectral_bandwidth(signal, sr=sr)[0]\n    spectral_bandwidth_3 = librosa.feature.spectral_bandwidth(signal, sr=sr, p=3)[0]\n    spectral_bandwidth_4 = librosa.feature.spectral_bandwidth(signal, sr=sr, p=4)[0]\n\n    audio_features = {\n        \"file_name\": audio_file_path,\n        \"zero_crossing_rate\": np.mean(librosa.feature.zero_crossing_rate(signal)[0]),\n        \"zero_crossings\": np.sum(librosa.zero_crossings(signal, pad=False)),\n        \"spectrogram\": np.mean(db_audio[0]),\n        \"mel_spectrogram\": np.mean(s_db_audio[0]),\n        \"harmonics\": np.mean(y_harm),\n        \"perceptual_shock_wave\": np.mean(y_perc),\n        \"spectral_centroids\": np.mean(spectral_centroids),\n        \"spectral_centroids_delta\": np.mean(spectral_centroids_delta),\n        \"spectral_centroids_accelerate\": np.mean(spectral_centroids_accelerate),\n        \"chroma1\": np.mean(chromagram[0]),\n        \"chroma2\": np.mean(chromagram[1]),\n        \"chroma3\": np.mean(chromagram[2]),\n        \"chroma4\": np.mean(chromagram[3]),\n        \"chroma5\": np.mean(chromagram[4]),\n        \"chroma6\": np.mean(chromagram[5]),\n        \"chroma7\": np.mean(chromagram[6]),\n        \"chroma8\": np.mean(chromagram[7]),\n        \"chroma9\": np.mean(chromagram[8]),\n        \"chroma10\": np.mean(chromagram[9]),\n        \"chroma11\": np.mean(chromagram[10]),\n        \"chroma12\": np.mean(chromagram[11]),\n        \"tempo_bpm\": tempo_y,\n        \"spectral_rolloff\": np.mean(spectral_rolloff),\n        \"spectral_flux\": np.mean(onset_env),\n        \"spectral_bandwidth_2\": np.mean(spectral_bandwidth_2),\n        \"spectral_bandwidth_3\": np.mean(spectral_bandwidth_3),\n        \"spectral_bandwidth_4\": np.mean(spectral_bandwidth_4),\n    }\n\n    # extract mfcc feature\n    mfcc_df = extract_mfcc_feature_means(audio_file_path,\n                                    signal,\n                                    sample_rate=sr,\n                                    number_of_mfcc=number_of_mfcc)\n\n    df = pd.DataFrame.from_records(data=[audio_features])\n\n    df = pd.merge(df, mfcc_df, on='file_name')\n\n    return df\n\n\ndef extract_mfcc_feature_means(audio_file_name: str,\n                          signal: np.ndarray,\n                          sample_rate: int,\n                          number_of_mfcc: int) -> pd.DataFrame:\n    # another MFCC approach\n    # as suggested by https://github.com/Cocoxili/DCASE2018Task2/blob/master/data_transform.py,\n    # https://arxiv.org/abs/1810.12832, and https://www.kaggle.com/c/freesound-audio-tagging\n    mfcc_alt = librosa.feature.mfcc(y=signal, sr=sample_rate,\n                                    n_mfcc=number_of_mfcc)\n    delta = librosa.feature.delta(mfcc_alt)\n    accelerate = librosa.feature.delta(mfcc_alt, order=2)\n\n    mfcc_features = {\n        \"file_name\": audio_file_name,\n    }\n\n    for i in range(0, number_of_mfcc):\n        # dict.update({'key3': 'geeks'})\n\n        # mfcc coefficient\n        key_name = \"\".join(['mfcc', str(i)])\n        mfcc_value = np.mean(mfcc_alt[i])\n        mfcc_features.update({key_name: mfcc_value})\n\n        # mfcc delta coefficient\n        key_name = \"\".join(['mfcc_delta_', str(i)])\n        mfcc_value = np.mean(delta[i])\n        mfcc_features.update({key_name: mfcc_value})\n\n        # mfcc accelerate coefficient\n        key_name = \"\".join(['mfcc_accelerate_', str(i)])\n        mfcc_value = np.mean(accelerate[i])\n        mfcc_features.update({key_name: mfcc_value})\n\n    df = pd.DataFrame.from_records(data=[mfcc_features])\n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Audio Feature Extraction Flow\n\nAs a part of the effort to classify bird calls (songs) in Cornell Birdcall Identification competition, there is a need to extract audio features from the digital audio records of the bird songs. In this project, librosa library is used for audio feature extraction. Since librosa-based audio feature computations are CPU-intensive, parallelizing such computations is the key architectural pattern to build a high-performance data transformation toolset to process the audio files.\n\nFrom the conceptual standpoint, the audio feature extraction flow was implemented as follows\n\n1.\tThe training set was split by the species\n2.\tEvery per-species batch of audio files with the bird songs have been processed with librosa as follows\n* Features extracted/calculated per the list below\n* Every feature that librosa extracts as np.array is further processed to calculate mean across the numeric values in np.array (actually it is everything except BPM variable)\n* The final results saved as a pandas dataframe\n* The dataframes serialized as CSV files for future use by ML pipeline tools down the road\n\nActivities per step 2 have been parallelized using three technologies compared in this research. When working with Dask and Ray, the output in a Pandas dataframe was consciously preferred over the native Dask and Ray counterparts of Pandsas to maximize the reuse of audio feature extraction code across all of the solutions.\nData files for just three species ('American Avocet', 'American Bittern', and 'American Crow') have been used in order to run the calculation within the reasonable time frame (as well as not to hit the Kaggle timeout issues when running the experiment on Kaggle).\n\nThe list of audio features extracted with librosa is below\n1.\tSpectrogram (decibel-scaled)\n2.\tMel Spectrograms (decibel-scaled)\n3.\tZero-crossing rate\n4.\tHarmonics (harmonics are characteristics that represent the sound color)\n5.\tPerceptual shock wave (it represents the sound rhythm and emotion of a soundtrack)\n6.\tSpectral Centroid (it indicates where the ”center of mass” for a sound is located and is calculated as the weighted mean of the frequencies present in the sound)\n7.\tChroma Frequencies (Chroma features). These are an interesting and powerful representation for music audio in which the entire spectrum is projected onto 12 bins representing the 12 distinct semitones (or chromas) of the musical octave.\n8.\tTempo BPM (beats per minute). Dynamic programming beat tracker.\n9.\tSpectral Rolloff. This is a measure of the shape of the signal. It represents the frequency below which a specified percentage of the total spectral energy (e.g. 85 %) lies.\n10.\tSpectral flux.\n11.\tSpectral Bandwidth. The spectral bandwidth is defined as the width of the band of light at one-half the peak maximum (or full width at half maximum (FWHM) and is represented by the two verticalred lines and λSB on the wavelength axis (in fact, we measure three spectral bandwidth features in this experiment)\n12.\tMFCC features (20 MFCC coefficients, 20 MFCC delta coefficients, and 20 MFCC accelerate coefficients)\n\nTotal of 87 audio features has been extracted.\n\n**Note**: the audio feature extraction flow has been designed to calculate tabular features for each audio file provided in the training set. Such an approach will be helpful to try both DL and classical ML algorithms during the feature importance analysis and ML predictions. However, you can find alternative approaches to audio feature extraction where the entire np.arrays are saved as files for futher processing (you can see examples of such a flow in https://github.com/Cocoxili/DCASE2018Task2/blob/master/data_transform.py etc., for instance)\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"start_time = dt.datetime.now()\nprint(\"Started at \", start_time)\n\n# Import data\ntrain_set_path, test_set_path = get_data_file_path(in_kaggle)\ntrain_csv = pd.read_csv(train_set_path)\ntest_csv = pd.read_csv(test_set_path)\n\n# Create some time features\ntrain_csv['year'] = train_csv['date'].apply(lambda x: x.split('-')[0])\ntrain_csv['month'] = train_csv['date'].apply(lambda x: x.split('-')[1])\ntrain_csv['day_of_month'] = train_csv['date'].apply(lambda x: x.split('-')[2])\n\n\n# Create Full Path so we can access data more easily\nbase_dir = get_base_train_audio_folder_path(in_kaggle)\ntrain_csv['full_path'] = base_dir + train_csv['ebird_code'] + '/' + train_csv['filename']\n\nfinal_data = ['American Avocet', 'American Bittern', 'American Crow',]\n\nfor ebird in final_data:\n    print(\"Starting to process a new species: \", ebird)\n    ebird_data = train_csv[train_csv['species'] == ebird]\n\n    short_file_name = ebird_data['ebird_code'].unique()[0]\n    print(\"Short file name: \", short_file_name)\n\n    result = []\n\n    for index, row in ebird_data.iterrows():\n        # process each audio file\n        f = delayed(extract_feautres)(row['full_path'])\n        result.append(f)\n\n    # combine chunks with transformed data into a single training set\n    extracted_features = delayed(pd.concat)(result)\n\n    df = extracted_features.compute()\n\n    # save extracted features to CSV\n    output_path = \"\".join([TRANSFORMED_DATA_PATH, short_file_name, \".csv\"])\n    df.to_csv(output_path, index=False)\n\n    print(\"Finished processing: \", ebird)\n\nprint('We are done. That is all, folks!')\nfinish_time = dt.datetime.now()\nprint(\"Finished at \", finish_time)\nelapsed = finish_time - start_time\nprint(\"Elapsed time: \", elapsed)","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}