{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Trying \"acoustic\" features\n[avloss](https://www.kaggle.com/avloss/audio-analysis-with-animation) and [eigrad](https://www.kaggle.com/eigrad/wip-some-audio-digging) has nicely shown us that our \"acoustic_data\" is literally acoustic audio data.  \nSo, I feel like trying methods from Audio/Music Information Retrieval (AIR/MIR).  \nIn this notebook, I try a famous and successful audio feature called [Mel-frequency cepstral coefficients (MFCC)](https://en.wikipedia.org/wiki/Mel-frequency_cepstrum). "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":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 in \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 \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\n\n# Any results you write to the current directory are saved as output.\n\nfrom tqdm import tqdm_notebook as tqdm\nimport matplotlib.pyplot as plt\n\n# train.csv is huge, so I implement csv_fragments() function\n# which yields DataFrame of the specified length while scaning a csv file from start to end.\n\nimport builtins\n\nrandom_seed = 4126\n\ncast = {\n    'acoustic_data': 'int',\n    'time_to_failure': 'float'\n}\n\ndef df_fragments(path, length, skip=1):\n    with open(path, 'r') as f:\n        m = {}\n        cols = []\n        count = 0\n        index = 0\n        for line in f:\n            if len(cols) == 0:\n                for col in line.strip(\"\\n\\r \").split(','):\n                    cols.append(col)\n                continue\n            if count == 0:\n                for col in cols:\n                    m[col] = []\n            if index % skip == 0:\n                for j, cell in enumerate(line.strip(\"\\n\\r \").split(',')):\n                    col = cols[j]\n                    m[col].append(getattr(builtins, cast[col])(cell))\n            count += 1\n            if count == length:\n                if index % skip == 0:\n                    yield pd.DataFrame(m)\n                index += 1\n                count = 0\n\ndef count_rows(path):\n    with open(path, 'r') as f:\n        i = -1\n        for _ in f:\n            i += 1\n        return i","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"[**LibROSA**](https://librosa.github.io/librosa/index.html) is an easy-to-use library to calculate audio features."},{"metadata":{"trusted":true},"cell_type":"code","source":"import librosa, librosa.display","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see MFCC of train data (first 150,000 records)"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"for df in df_fragments('../input/train.csv', 150000):\n    mfcc = librosa.feature.mfcc(df['acoustic_data'].values.astype('float32'))\n    plt.figure(figsize=(25, 5))\n    librosa.display.specshow(mfcc, x_axis='time')\n    plt.colorbar()\n    break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The shape of MFCC is (\\[No. of features (20 by default)\\], \\[time\\]).  \nI tentatively create train data by calculating mean values along time axis for each 150000 train records (same size as test data fragments)."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('counting total...')\ntotal = count_rows('../input/train.csv')\nprint('total: {}'.format(total))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print('generating train data...')\nfragment_size = 150000\nskip = 1\n# you can reduce train data to process for some quick experiments\n# skip = 10\n\nmfcc_ttf_map = {}\nfor df in tqdm(df_fragments('../input/train.csv', length=fragment_size, skip=skip), total=(total//fragment_size)//skip):\n    mfcc = librosa.feature.mfcc(df['acoustic_data'].values.astype('float32'))\n    mfcc_mean = mfcc.mean(axis=1)\n    for i, each_mfcc_mean in enumerate(mfcc_mean):\n        key = 'mfcc_{}'.format(i)\n        if key not in mfcc_ttf_map:\n            mfcc_ttf_map[key] = []\n        mfcc_ttf_map[key].append(each_mfcc_mean)\n    key = 'time_to_failure'\n    if key not in mfcc_ttf_map:\n        mfcc_ttf_map[key] = []\n    mfcc_ttf_map[key].append(df.iloc[-1][key])\n\nmfcc_ttf_df = pd.DataFrame(mfcc_ttf_map)\nfname = 'mfcc_train.csv'\nmfcc_ttf_df.to_csv(fname, index=False)\nprint('saved {}.'.format(fname))\n        ","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false},"cell_type":"markdown","source":"Let's visualize train data."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig = plt.figure(figsize=[20,5])\nax1 = fig.add_subplot(111)\nax2 = ax1.twinx()\n\nmfcc_ttf_df['time_to_failure'].plot(ax=ax1, y='time_to_failure', legend=True, color='black')\nax1.legend(loc='upper left')\nmfcc_ttf_df.drop(['time_to_failure'], axis=1).plot(ax=ax2, legend=True)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some of mean MFCC feature values seem to have linear relationship with time_to_failure.  \nLet's try linear regression (cross validation fold=10)."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import cross_val_score\n\ndef report_cv(model):\n    X = mfcc_ttf_df.drop(['time_to_failure'], axis=1).values\n    y = mfcc_ttf_df['time_to_failure'].values\n    scores = cross_val_score(model, X, y, scoring='neg_mean_absolute_error', cv=10)\n    print('Cross Validation scores: {}'.format(abs(scores)))\n    print('Average score: {}'.format(abs(scores.mean())))\n\nreport_cv(LinearRegression())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Also try XGBoost"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from xgboost import XGBRegressor\n\nreport_cv(XGBRegressor(random_state=random_seed))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Following is a graph with cross-validation fold boundary (red vertical line).  \nIt seems both models work poorly when time_to_failure of test data is abnormally large (5th or 9th fold) or abnormally small (4th fold).  \n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"from matplotlib import patches \n\nfig = plt.figure(figsize=[20,5])\nax1 = fig.add_subplot(111)\nax2 = ax1.twinx()\n\nmfcc_ttf_df['time_to_failure'].plot(ax=ax1, y='time_to_failure', legend=True, color='black')\nax1.legend(loc='upper left')\nmfcc_ttf_df.drop(['time_to_failure'], axis=1).plot(ax=ax2, legend=True)\n\nfold_len = len(mfcc_ttf_df)//10\nfor i in range(1, 10):\n    plt.axvline(x=fold_len * i,color='red')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Following is just my imagination...  \nWe can see small spikes between the times of failure, and after such spikes mfcc values go down in some degree. Aren't they small failures which were not recorded as actual failures? Such small failures seem to postpone following actual failures (corresponding to foreshock in real world?). If we could detect such small failures, we may be able to improve our score, but it seems impossible for me to detect small failures from fragmented and shuffled test data..."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=[20,5])\nax1 = fig.add_subplot(111)\nax2 = ax1.twinx()\n\nmfcc_ttf_df['time_to_failure'].plot(ax=ax1, y='time_to_failure', legend=True, color='black')\nax1.legend(loc='upper left')\nmfcc_ttf_df.drop(['time_to_failure'], axis=1).plot(ax=ax2, legend=True)\n\nax2.add_patch(patches.Rectangle((170,320),80, 80,linewidth=3,edgecolor='r',facecolor='none'))\nax2.add_patch(patches.Rectangle((490,320),80,80,linewidth=3,edgecolor='r',facecolor='none'))\nax2.add_patch(patches.Rectangle((1780,320),80,80,linewidth=3,edgecolor='r',facecolor='none'))\nax2.text(1780, 280, 'small failure?', fontsize=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Anyway, let's create files for submission using all the train data."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import re\n\nprint('generating test features...')\ntest_dir = '../input/test'\ntest_map = {}\nfor fname in tqdm(os.listdir(test_dir)):\n    path = test_dir + '/' + fname\n    df = pd.read_csv(path)\n    mfcc = librosa.feature.mfcc(df['acoustic_data'].values.astype('float32'))\n    mfcc_mean = mfcc.mean(axis=1)\n    for i, each_mfcc_mean in enumerate(mfcc_mean):\n        key = 'mfcc_{}'.format(i)\n        if key not in test_map:\n            test_map[key] = []\n        test_map[key].append(each_mfcc_mean)\n    key = 'seg_id'\n    if key not in test_map:\n        test_map[key] = []\n    test_map[key].append(re.sub('.csv$', '', fname))\ntest_df = pd.DataFrame(test_map)\ntest_csv = 'mfcc_test.csv'\ntest_df.to_csv(test_csv, index=False)\nprint('saved {}'.format(test_csv))\n\n\ndef submit(model, file_path):\n    X = mfcc_ttf_df.drop(['time_to_failure'], axis=1).values\n    y = mfcc_ttf_df['time_to_failure'].values\n    model.fit(X, y)\n    \n    X_submit = test_df.drop(['seg_id'], axis=1).values\n    y_submit = model.predict(X_submit)\n    submit_df = pd.DataFrame({\n        'seg_id': test_df['seg_id'].values,\n        'time_to_failure': y_submit\n    })\n    submit_df.to_csv(file_path, index=False)\n    print('saved {}'.format(file_path))\n    ","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"submit(LinearRegression(), 'submit_linear.csv')\nsubmit(XGBRegressor(random_state=random_seed), 'submit_xgb.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I hope specialists of Audio/Music Information Retrieval go into detail of acoustic features. Thanks!"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}