{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"markdown","source":"# Overview\n\nAs I am new to deep learning and my background is more \"classic\" machine learning, I decided to start with random forest (RF) / xgboost (xgb) / logistic regression and then learn how to use neural nets. I started with skimming through articles on musical instrument detection and making a list of important features, coding these features and using these features as an input into RF and xgb. Training both classifiers on the whole dataset and averaging outputs using geometric mean gave 0.844 on the leaderboard, which was already great.\n\nNext I started to learn how to use CNNs. \"Learn from what is there\" they say, so I checked what was already done and came across kernel by Zafar - [Beginner's Guide to Audio Data\n](https://www.kaggle.com/fizzbuzz/beginner-s-guide-to-audio-data), which was a great starting point. Next I did some studying by listening to course by Andrew Ng [links from here](https://www.deeplearning.ai/) and played with hyperparameter optimization.\n\nAs I had dataset for training RF/xgb, I decided to use that dataset along with dataset for CNN as two distinct inputs into one NN which gave good improvement to mapk. \n\nCurrently I have around 0.92 on the leaderboard. Next steps are to study amazing input done by [daisukelab](https://www.kaggle.com/daisukelab) which can be found [here](https://www.kaggle.com/c/freesound-audio-tagging/discussion/57051).\n\nThings I will try:\n* augmentaions\n* oversampling\n* other NN architectures\n* re-sampling audio at 16k or 24k\n* using RNN or 1d CNN as suggested by Zafar (or time-dependent approach by daisukelab)\n* sequential learning / pseudo-labelling (e.g. definitely it is possible to label all test data in 2-3 days: add to train predictions with highest probability, train, etc.)\n\nWhat I have noticed\n* Cross-validated mapk is way more lower then the one on the leaderboard\n* Geometric mean works, arithmetic doesn't (for ensembling)\n* Training classifers on all data can help (cross-validate to tune parameters, then re-train on all data)\n* Ensembling helps\n* TWO SUBMITS IS NOT ENOUGH :D\n\nI will try to share more insights and ideas with the time.\n"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Part 1: let's listen to the data (training on the wrong data?)\n## Training set\n\nFirst I was curious if I should use only manually verified data or all of the data. So I just listened to around first 100 audio files in the training set and found some interesting examples (which do not sound as their labels).\n\n**Examples are below. compare names of the files with what you will hear.**"},{"metadata":{"trusted":true,"_uuid":"db9a581a07ca9294f03db77a2a9c3cbcb134599f","collapsed":true},"cell_type":"code","source":"import IPython.display as ipd  # To play sound in the notebook\ntelephone = '../input/00d3bba3wav/00d3bba3.wav'   # Telephone\nflute = '../input/00d9fa61wav/00d9fa61.wav'   # Flute\nsqueak = '../input/013264d3wav/013264d3.wav'   # Squeak\ncello = '../input/0184c390wav/0184c390.wav'   # Cello\nshatter = '../input/01a39e95wav/01a39e95.wav'   # Shatter\n","execution_count":5,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7e4d530c0cf22723030be319960c1f67bab1cc96","collapsed":true},"cell_type":"code","source":"ipd.Audio(telephone)","execution_count":6,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"00729b8014ed9f4ebd05cdb0043bddac918ec65b","collapsed":true},"cell_type":"code","source":"ipd.Audio(flute)","execution_count":7,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ada28fc5dd94c58e4d9396cfeb7464307a1f0d26","collapsed":true},"cell_type":"code","source":"ipd.Audio(squeak)","execution_count":8,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dbff669db821a3408e4329b44d2b0c26240a365a","collapsed":true},"cell_type":"code","source":"ipd.Audio(cello)","execution_count":9,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"51bc8ec7bf3dee049977d5023bada7cb673680c4","collapsed":true},"cell_type":"code","source":"ipd.Audio(shatter)","execution_count":10,"outputs":[]},{"metadata":{"_uuid":"d900f0388a7bdd3a764513efd5a915a24e47c780"},"cell_type":"markdown","source":"First file (telephone) is definitely not a telephone, but a telegraph. Second one (flute) doesn't sound like flute (just some whistlening sound, flute sounds like [this](https://www.youtube.com/watch?v=com5gPoZ8sI&feature=youtu.be)). Third one is a combination of sound or a car engine, smt like doors opening/closing and squek in the end. Fourth one (cello) is most likely electric piano, and the last one (shatter) is keys jangling.\n\nThere are definitely more examples like this, as I just checked top 100 files (and not all of them). For example, file 02267a1a.wav is just a sound of a person walking. \n\nOther examples from test set can be found below. As mentioned by organizers, not all sounds from the test set are used to calculate scores, but anyway it is worth discussing what are we trying to achieve here. \n\nIt seems that part of the training set is labelled in the wrong way and some of the sounds from test set (examples below) belong to none of the classes found in the training set.\n"},{"metadata":{"_uuid":"8422e561462c10d0e1b81a2685982a1722b6b636"},"cell_type":"markdown","source":"## Test set\nExamples of some of the most weird sounds I found are below. It was easy to find them - I had two classifiers (RF and xgb), I compared predictions produced by them and if most likely class was different, I checked what is that sound. Also I checked around 10-20 sound files for which rf/xgb were less sure to make prediction."},{"metadata":{"trusted":true,"_uuid":"13236cee2879172d0aa642f4673ee3cca8222da8","collapsed":true},"cell_type":"code","source":"ipd.Audio('../input/0ce127f9/0ce127f9.wav')","execution_count":11,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"add672a4bc1f058dc97bc9cd76082fd6c7847860","collapsed":true},"cell_type":"code","source":"ipd.Audio('../input/01e6e112wav/01e6e112.wav')","execution_count":12,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9737872f093c01dfce1e19ebbe6954f4f6d0ffd1","collapsed":true},"cell_type":"code","source":"ipd.Audio('../input/026820e6wav/026820e6.wav')","execution_count":13,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ad2fd619488431f28855d6f6289beb6b72bb6607","collapsed":true},"cell_type":"code","source":"ipd.Audio('../input/013264d3wav/013264d3.wav')","execution_count":14,"outputs":[]},{"metadata":{"_uuid":"c3c9bf1e52aa80d6557ac1c161ef8a127ea33c25"},"cell_type":"markdown","source":"Some of the sound from test set is just music. Other examples are 0300b76b.wav, 02a0eb3c.wav, 01f9883e.wav,  03e1e393.wav, 0539bb41.wav, 054eeab6.wav.\n\nWhat is more important is that I checked only around first 100 files (out of around 2000) for which RF and xgb gave different predictions"},{"metadata":{"_uuid":"fc0a9ae545c41419f74a5f37bb05f0167915958f"},"cell_type":"markdown","source":"## Part of the data is wrong?\nClearly, as most of the data is labelled automatically, it will be partly wrong. So we are training classifiers to classify data labelled by another classifier.  Of course, it is still woth doing so, but I suggest that smt like sequential learning would help in solving such tasks - we start we the data labelled correctly, then we classify next N observations, somehow check predictions, re-label f needed, etc. It would give us better dataset.   "},{"metadata":{"_uuid":"5dddd1f704895f5f168a9081a6f016a42bc973a6"},"cell_type":"markdown","source":"# Part 3: RF and xgb (0.844 on LB)"},{"metadata":{"_uuid":"21669514092dbc817963a6b08639936c01f358d2"},"cell_type":"markdown","source":"## Reading data"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"55e9e9f2f41d5982baf7bf398fd79839e8f13888"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nimport librosa\nfrom scipy.io import wavfile\n\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import precision_score\nfrom sklearn.metrics import make_scorer\nfrom sklearn import preprocessing\nfrom collections import Counter\nfrom sklearn.preprocessing import PolynomialFeatures\nfrom sklearn import linear_model\nimport xgboost as xgb\nfrom sklearn.utils import shuffle\n\nfrom sklearn import manifold, datasets\nfrom sklearn.preprocessing import scale\n\nimport os","execution_count":22,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"c7b80cc5147223132691a2911377c25de44c3693"},"cell_type":"code","source":"def apk(actual, predicted, k=10):\n    \"\"\"\n    Computes the average precision at k.\n    This function computes the average prescision at k between two lists of\n    items.\n    Parameters\n    ----------\n    actual : list\n             A list of elements that are to be predicted (order doesn't matter)\n    predicted : list\n                A list of predicted elements (order does matter)\n    k : int, optional\n        The maximum number of predicted elements\n    Returns\n    -------\n    score : double\n            The average precision at k over the input lists\n    \"\"\"\n    if len(predicted)>k:\n        predicted = predicted[:k]\n\n    score = 0.0\n    num_hits = 0.0\n\n    for i,p in enumerate(predicted):\n        if p in actual and p not in predicted[:i]:\n            num_hits += 1.0\n            score += num_hits / (i+1.0)\n\n    if not actual:\n        return 0.0\n\n    return score / min(len(actual), k)\n\ndef mapk(actual, predicted, k=10):\n    \"\"\"\n    Computes the mean average precision at k.\n    This function computes the mean average prescision at k between two lists\n    of lists of items.\n    Parameters\n    ----------\n    actual : list\n             A list of lists of elements that are to be predicted \n             (order doesn't matter in the lists)\n    predicted : list\n                A list of lists of predicted elements\n                (order matters in the lists)\n    k : int, optional\n        The maximum number of predicted elements\n    Returns\n    -------\n    score : double\n            The mean average precision at k over the input lists\n    \"\"\"\n    return np.mean([apk(a,p,k) for a,p in zip(actual, predicted)])\n","execution_count":16,"outputs":[]},{"metadata":{"_uuid":"e9c6a71ff03edbe14cc5f362aea81450f32a7881"},"cell_type":"markdown","source":"Next I created some auxiliary functions\n* to get spectrum (get_spectra)\n* to get width of the spike in spectrum (get_width)\n* to get running mean of the spectra (running_mean)\n* to get mfcc features\n* to get training/test set (get_training_dataset)\n\nBasic idea was to crate as many different predictors as possible, given some knowledge of what should work. For example, shape and width of spikes spectrum, number of zero-crossings, number of peaks in spectrum, etc. - all these are expected to make difference.\n\nEventually I ended up with quite wide dataset which takes a lot of time to calculate (around 4 hours for each train and test). It takes so much time mostly because of massive amount of tricky predictors. It gives only 0.844 on the LB, but keep in mind that it is non-NN solution and we can feed it to tSNE or MDS and have fun.\n"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"4c29f7fe0fb4eea18bd76aa69068e1074af08f14"},"cell_type":"code","source":"def get_spectra_win(y, L, N):\n    dft = np.fft.fft(y)\n    fl = np.abs(dft)\n    xf = np.arange(0.0, N/L, 1/L)\n    return (xf,fl)\n\ndef get_spectra(signal, fs, M = 1000, sM = 500):\n\n    N = signal.shape[0]\n    ind = np.arange(100, N, M)\n\n    spectra = []\n    meanspectrum = np.repeat(0,M)\n\n    for k in range(1,len(ind)):\n        n1 = ind[k-1]\n        n2 = ind[k]\n        y = signal[n1:n2]\n        L = (n2-n1)/fs\n        N = n2-n1\n        (xq, fq) = get_spectra_win(y, L, N)\n        spectra.append(fq)\n\n    spectra = pd.DataFrame(spectra)\n    meanspectrum = spectra.apply(lambda x: np.log(1+np.mean(x)), axis=0)\n    stdspectrum = spectra.apply(lambda x: np.log(1+np.std(x)), axis=0)\n    \n    meanspectrum = meanspectrum[0:sM]\n    stdspectrum = stdspectrum[0:sM]\n    \n    return (meanspectrum, stdspectrum) \n\ndef get_width(w):\n    if np.sum(w) == 0:\n        return [0,0,0]\n    else:\n        z = np.diff(np.where(np.insert(np.append(w,0),0,0)==0))-1\n        z = z[z>0]\n    return [np.log(1+np.mean(z)),np.log(1+np.std(z)),np.log(1+np.max(z)),len(z)]\n\ndef running_mean(x, N):\n    cumsum = np.cumsum(np.insert(x, 0, 0)) \n    return (cumsum[N:] - cumsum[:-N]) / float(N)\n\ndef clear():\n    os.system( 'cls' )","execution_count":17,"outputs":[]},{"metadata":{"_uuid":"ab5d089d74d4f6e69aa7814e94b1bd4e2b92c4b5"},"cell_type":"markdown","source":"Next we created many-many features (if you are interested in more details, let me know in the comments, code is little bit like spaghetti :D )."},{"metadata":{"trusted":true,"_uuid":"6e7623351683a0708c8de435abe33c7d8c1b7055","collapsed":true},"cell_type":"code","source":"def get_training_dataset(training=1, dir_path='D:/python/audio_train/audio_train/'):\n    \n    if training==1:\n        trainnames = pd.read_csv(dir_path + 'train.csv')\n        labelnames = list(trainnames['label'].unique())\n        le = preprocessing.LabelEncoder()\n        le.fit(labelnames)\n        files_labels = zip(trainnames['fname'].values, trainnames['label'].values)\n    elif training==0:\n        score_filelist = [str(x) for x in os.listdir(dir_path)]\n        # 'D:/python/audio_test/audio_test'\n        labelnames = np.repeat('unlabeled',len(score_filelist))\n        le = preprocessing.LabelEncoder()\n        le.fit(labelnames)\n        files_labels = zip(score_filelist, labelnames)\n    else:\n        return []\n\n    df_m = []\n    df_sd = []\n    df_sig = []\n    df_mfcc = []\n    df_fbank = []\n    df_ssc = []    \n    labels_processed = []\n    filenames_processed = []\n\n    i = 0\n    \n    for filename, labelname in files_labels:\n        \n        label = le.transform([labelname])[0]\n        fname  = dir_path + filename\n        fs, rawsignal = wavfile.read(fname)\n        if rawsignal.size == 0:\n            rawsignal = np.random.randint(0,2,44000)\n            \n        if rawsignal.dtype == 'int16':\n            nb_bits = 16 # -> 16-bit wav files\n        elif rawsignal.dtype == 'int32':\n            nb_bits = 32 # -> 32-bit wav files\n        max_nb_bit = float(2 ** (nb_bits - 1))\n        rawsignal = rawsignal/max_nb_bit\n        \n        # signal features\n        rawsignal_sq = rawsignal*rawsignal\n        silenced = []\n        sound = []\n        attack = []\n        for wd in [2000,10000]:\n            rawsignal_sq_rm = running_mean(rawsignal_sq, wd)            \n            w1 = 1*(rawsignal_sq_rm<0.01*np.max(rawsignal_sq_rm))\n            silenced = silenced + get_width(w1)\n            w2 = 1*(rawsignal_sq_rm<0.05*np.max(rawsignal_sq_rm))\n            silenced = silenced + get_width(w2)            \n            w3 = 1*(rawsignal_sq_rm>0.05*np.max(rawsignal_sq_rm))\n            sound = sound + get_width(w3)\n            w4 = 1*(rawsignal_sq_rm>0.25*np.max(rawsignal_sq_rm))\n            sound = sound + get_width(w4)\n            time_to_attack = np.min(np.where(rawsignal_sq_rm>0.99*np.max(rawsignal_sq_rm)))\n            time_rel = np.where(rawsignal_sq_rm<0.2*np.max(rawsignal_sq_rm))[0]\n            if (time_rel.size == 0):\n                time_to_relax = len(rawsignal_sq_rm)\n            elif (time_rel[time_rel>time_to_attack].size==0):\n                time_to_relax = len(rawsignal_sq_rm)\n            else:\n                time_to_relax = np.min(time_rel[time_rel>time_to_attack])\n            attack.append(np.log(1+time_to_attack))\n            attack.append(np.log(1+time_to_relax))\n\n        lr = len(rawsignal)\n        zerocross_tot = np.log(1+np.sum(np.array(rawsignal[0:(lr-1)])*np.array(rawsignal[1:lr])<=0))\n        zerocross_prop = np.sum(np.array(rawsignal[0:(lr-1)])*np.array(rawsignal[1:lr])<=0)/lr\n        df_sig.append(sound+silenced+attack+[zerocross_tot,zerocross_prop])\n\n        (m, sd) = get_spectra(rawsignal, fs, 2000, 1000 )\n        df_m.append(m)\n        df_sd.append(sd)\n\n        labels_processed.append(label)     \n        filenames_processed.append(filename)\n        \n        # mfcc\n        \n        mfcc_feat = librosa.feature.mfcc(rawsignal, sr = fs, n_mfcc=40)\n        mfcc_feat = pd.DataFrame(np.transpose(mfcc_feat))       \n        \n        # mfcc_feat = mfcc(rawsignal, fs, nfft = 1103, numcep = 30 )\n        # mfcc_feat = pd.DataFrame(mfcc_feat)\n        mfcc_mean = mfcc_feat.apply(lambda x: np.mean(x), axis=0)\n        mfcc_sd = mfcc_feat.apply(lambda x: np.std(x), axis=0)\n        mfcc_max = mfcc_feat.apply(lambda x: np.max(x), axis=0)\n        mfcc_med = mfcc_feat.apply(lambda x: np.median(x), axis=0)        \n        mfcc_res = np.array(list(mfcc_mean)+list(mfcc_sd)+list(mfcc_max)+list(mfcc_med)+[np.log(1+len(rawsignal))])\n        df_mfcc.append(mfcc_res)\n             \n        i = i+1\n        labelname = labelname + ' '*(20-len(labelname))\n        label_string = str(label) + ' '*(3-len(str(label)))\n        i_str = str(i) + ' '*(5-len(str(i)))\n            \n        print('\\r', i_str, filename, ' - ',labelname,' - ',label_string, end='', flush=True)\n        \n    # to data frames\n    df_sig = pd.DataFrame(df_sig)\n    df_sig.fillna(0, inplace = True)\n\n    df_sd = pd.DataFrame(df_sd)\n    df_m = pd.DataFrame(df_m)\n    df_mfcc = pd.DataFrame(df_mfcc) \n    \n    # predictors related to peaks \n    def num_peaks(x):\n        x = np.array(x[0:len(x)])\n        n10 = np.sum(x>0.10*np.max(x))\n        n20 = np.sum(x>0.20*np.max(x))\n        n50 = np.sum(x>0.50*np.max(x))\n        n90 = np.sum(x>0.90*np.max(x))\n        n99 = np.sum(x>0.99*np.max(x))\n        lead_min = np.min(np.where(x==np.max(x)))\n        cnt = 0\n        w10 = get_width(1*(x>0.10*np.max(x)))\n        w20 = get_width(1*(x>0.20*np.max(x)))\n        w50 = get_width(1*(x>0.50*np.max(x)))\n        w90 = get_width(1*(x>0.90*np.max(x)))\n        w99 = get_width(1*(x>0.99*np.max(x)))  \n        W = w10+w20+w50+w90+w99\n\n        f_sc = np.sum(np.arange(0,len(x))*(x*x)/np.sum(x*x))\n\n\n        i1 = np.where(x<0.10*np.max(x))[0]\n        if i1.size == 0:\n            lincoef_w = [0,0,0]\n        else:\n            a1 = i1[i1<lead_min]\n            a2 = i1[i1>lead_min]\n\n            if a1.size == 0:\n                i1_left = 0\n            else:\n                i1_left = np.max(i1[i1<lead_min])\n            if a2.size == 0:\n                i1_right = 0\n            else:\n                i1_right = np.min(i1[i1>lead_min])\n\n            lead_min_width = i1_right - i1_left  \n            if (lead_min_width>2):\n                poly_w = PolynomialFeatures(degree=2, include_bias = False)\n                f_ind_w = poly_w.fit_transform(np.arange(i1_left,i1_right,1).reshape(-1, 1))\n                clf_w = linear_model.LinearRegression()\n                linmodel_w = clf_w.fit(f_ind_w, np.array(x[i1_left:i1_right]))\n                lincoef_w = list(linmodel_w.coef_)+[linmodel_w.intercept_]\n            else:\n                lincoef_w = [0,0,0]\n\n        S = np.sum(x)\n        S_n = np.sum(x)/len(x)\n        S2 = np.sqrt(np.sum(x*x))    \n        S2_n = np.sqrt(np.sum(x*x))/len(x)\n        integrals = [S,S_n,S2,S2_n]\n        \n        poly = PolynomialFeatures(degree=2, include_bias = False)\n        f_ind = poly.fit_transform(np.arange(0,len(x)).reshape(-1, 1))\n        clf = linear_model.LinearRegression()\n        linmodel = clf.fit(f_ind, x)\n        lincoef_spectrum = list(linmodel.coef_)+[linmodel.intercept_]\n\n        high_freq_sum_50 = np.sum(x[0:50]>=0.5*np.max(x))\n        high_freq_sum_90 = np.sum(x[0:50]>=0.9*np.max(x))\n\n        r = [f_sc,n10,n20,n50,n90,n99,lead_min,high_freq_sum_50,high_freq_sum_90]+W+lincoef_spectrum+integrals+lincoef_w\n        return r\n\n    def runningMeanFast(x, N=20):\n        return np.convolve(x, np.ones((N,))/N)[(N-1):]\n\n    df_rm = df_m.apply(runningMeanFast, axis=1)\n    df_sc = df_rm.apply(lambda x: x[np.arange(0,len(x),40)],axis=1)\n    df_m_filt = df_m.apply(lambda x: x[np.arange(0,2,1)],axis=1)\n    df_peaks = df_m.apply(num_peaks,axis=1)\n    df_peaks = pd.DataFrame(list(df_peaks)) \n    df_rm = pd.DataFrame(df_rm)\n    df_sc = pd.DataFrame(df_sc)\n    df_m_filt = pd.DataFrame(df_m_filt)\n    \n    df_fbank = pd.DataFrame(df_fbank)\n    df_ssc = pd.DataFrame(df_ssc)    \n    #\n    \n    df_sd.columns = ['fft_sd'+str(i) for i in range(0,len(df_sd.columns))]\n    df_m.columns = ['fft_mean'+str(i) for i in range(0,len(df_m.columns))]\n    df_rm.columns = ['fft_rmean'+str(i) for i in range(0,len(df_rm.columns))]\n    df_sc.columns = ['fft_scaled'+str(i) for i in range(0,len(df_sc.columns))]\n    df_mfcc.columns = ['mfcc'+str(i) for i in range(0,len(df_mfcc.columns))]\n    df_fbank.columns = ['fbank'+str(i) for i in range(0,len(df_fbank.columns))]\n    df_ssc.columns = ['ssc'+str(i) for i in range(0,len(df_ssc.columns))]\n    \n    df_sig.columns = ['snd_wd_2000_mean_th005','snd_wd_2000_sd_th005','snd_wd_2000_max_th005'\n                      ,'snd_wd_2000_len_th005'\n                     ,'snd_wd_2000_mean_th025','snd_wd_2000_sd_th025','snd_wd_2000_max_th025'\n                      ,'snd_wd_2000_len_th025'\n                     ,'snd_wd_10000_mean_th005','snd_wd_10000_sd_th005','snd_wd_10000_max_th005'\n                      ,'snd_wd_10000_len_th005'\n                     ,'snd_wd_10000_mean_th025','snd_wd_10000_sd_th025','snd_wd_10000_max_th025'\n                      ,'snd_wd_10000_len_th025'\n\n                     ,'sln_wd_2000_mean_th001','sln_wd_2000_sd_th001','sln_wd_2000_max_th001'\n                      ,'sln_wd_2000_len_th001'\n                     ,'sln_wd_2000_mean_th005','sln_wd_2000_sd_th005','sln_wd_2000_max_th005'\n                      ,'sln_wd_2000_len_th005'\n                     ,'sln_wd_10000_mean_th001','sln_wd_10000_sd_th001','sln_wd_10000_max_th001'\n                      ,'sln_wd_10000_len_th001'\n                     ,'sln_wd_10000_mean_th005','sln_wd_10000_sd_th005','sln_wd_10000_max_th005'\n                      ,'sln_wd_10000_len_th005'\n\n                     , 'time_to_attack_2000', 'time_to_relax_2000'\n                     , 'time_to_attack_10000', 'time_to_relax_10000'\n                     , 'zerocross_tot','zerocross_prop'\n                     ]\n    df_peaks.columns = ['f_sc','n10','n20','n50','n90','n99','lead_min'\n                        ,'high_freq_sum_50','high_freq_sum_90'\n                       ,'w10_mean','w10_sd','w10_max','w10_len'\n                       ,'w20_mean','w20_sd','w20_max','w20_len'\n                       ,'w50_mean','w50_sd','w50_max','w50_len'\n                       ,'w90_mean','w90_sd','w90_max','w90_len'\n                       ,'w99_mean','w99_sd','w99_max','w99_len'\n                       ,'coef_deg1','coef_deg2','coef_deg0'\n                       ,'S','S_n','S2','S2_n'\n                       ,'coef_deg1_w','coef_deg2_w','coef_deg0_w']    \n        \n    return df_peaks, df_sig, df_mfcc, df_rm, df_m, df_sc, df_fbank, df_ssc, le, labels_processed, filenames_processed\n","execution_count":18,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"fdebc7c86cda44777184363a00148a23d047b2e6"},"cell_type":"markdown","source":"As it takes quite a lot of time to load all files and create all features, I don't do it here. Change parameters to get training and scoring datasets."},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"12f9afe68405b8c7727437e3cc66b5944789bc8b"},"cell_type":"code","source":"get_data = 0\n\nif get_data == 1:\n    dataframe_list_training = get_training_dataset()\n    df_peaks, df_sig, df_mfcc, df_rm, df_m, df_sc, df_fbank, df_ssc, le, labels, files = dataframe_list_training\n    df_result = pd.concat([df_peaks, df_sig, df_mfcc], axis=1, ignore_index=True)\n    df_result.columns = list(df_peaks.columns) + list(df_sig.columns) + list(df_mfcc.columns)\n\n    Xall = np.array(df_result)\n    Xall = Xall.reshape(df_result.shape)\n    yall = np.array(labels)\n    Xall, yall = shuffle(Xall, yall, random_state=0)\n\n    df_result.to_csv('.../df_result_ens.csv')\n    np.save('.../Xalle_ens.npy', Xall)\n    np.save('.../yall_ens.npy', yall)\n    np.save('.../le_ens.npy', le)\n","execution_count":24,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1d923eef99898287186eae546c3df465a8049a7b","collapsed":true},"cell_type":"code","source":"get_data = 0\nif get_data == 1:\n    dataframe_list_training = get_training_dataset(False,'.../audio_test/')\n    df_peaks_scoring, df_sig_scoring, df_mfcc_scoring, df_rm_scoring, df_m_scoring, df_sc_scoring, _, _, _, _, files = dataframe_list_scoring\n    df_result_scoring = pd.concat([df_peaks_scoring, df_sig_scoring, df_mfcc_scoring], axis=1, ignore_index=True)\n    \n    X_scoring = np.array(df_result_scoring)\n    X_scoring = X_scoring.reshape(df_result_scoring.shape)\n    \n    np.save('.../X_scoring_ens.npy', X_scoring)\n    \n    score_filelist = [str(x) for x in os.listdir('.../audio_test/')]\n    np.save('.../score_filelist.npy',np.array(score_filelist))\n","execution_count":25,"outputs":[]},{"metadata":{"_uuid":"657ae22a0d6cbe7b2f62cfcd1aaacdea44a8d94f"},"cell_type":"markdown","source":"## Loading data"},{"metadata":{"trusted":true,"_uuid":"73c0025bfbe7a3e58e93711f7003e4643fc6c083","collapsed":true},"cell_type":"code","source":"Xall = np.load('../input/ensemble/Xall_ens.npy')\nyall = np.load('../input/ensemble/yall_ens.npy')\nX_scoring = np.load('../input/ensemble/X_scoring_ens.npy')\nle = np.load('../input/ensemble/le_ens.npy')\nscore_filelist = np.load('../input/score-filelist/score_filelist.npy')","execution_count":26,"outputs":[]},{"metadata":{"_uuid":"6d7df56bc7da8969f3a89f57af569186df364f65"},"cell_type":"markdown","source":"## EDA\nAs my first goal was to classify files using RF and xgb I created more classic dataset if compared with CNN or similar approaches, i.e. I had dataset with one row per observation. Interesting thing to do is to visualize the data using standard dimensionality reduction approaches like Spectral embedding, tSNE, MDS, etc. \n\nI scale data before applying all these techniques.\n"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"85c9f23db0366a274daa1aec51a60085b4b1fff2"},"cell_type":"code","source":"X_eda = Xall[0:2000,:]\ny_eda = yall[0:2000]\nX_eda = scale(X_eda)","execution_count":27,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ba9889ac041329b8f5bfd33b7603a839d50c583b","collapsed":true},"cell_type":"code","source":"mds = manifold.MDS(2, max_iter=100, n_init=1)\nY = mds.fit_transform(X_eda)\nplt.figure(1)\nplt.figure(figsize=(5,5))\nfor (j,cl) in enumerate([29,38,12,1,9,35]):\n    cname = le.tolist().inverse_transform([cl])[0]\n    plt.subplot(3,2,j+1)\n    plt.scatter(Y[y_eda!=cl][:,0], Y[y_eda!=cl][:,1], c='blue', alpha=0.75)\n    plt.scatter(Y[y_eda==cl][:,0], Y[y_eda==cl][:,1], c='red', alpha=0.75)\n    plt.title(cname+' MDS')\n    plt.axis('tight')\nplt.tight_layout() \nplt.show()\n   \n","execution_count":28,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fc8a0b28c1ccc90cae70756a65191811b4e79e9a","collapsed":true},"cell_type":"code","source":"for (i,prplx) in enumerate([5,20,100]):\n    tsne = manifold.TSNE(n_components=2, init='pca', random_state=0, perplexity = prplx)\n    Y = tsne.fit_transform(X_eda)\n    plt.figure(1)\n    plt.figure(figsize=(5,5))\n    for (j,cl) in enumerate([29,38,12,1,9,35]):\n        cname = le.tolist().inverse_transform([cl])[0]\n        plt.subplot(3,2,j+1)\n        plt.scatter(Y[y_eda!=cl][:,0], Y[y_eda!=cl][:,1], c='blue', alpha=0.75)\n        plt.scatter(Y[y_eda==cl][:,0], Y[y_eda==cl][:,1], c='red', alpha=0.75)\n        plt.title(cname+' t-SNE perplexity: '+str(prplx))\n        plt.axis('tight')\n    plt.tight_layout() \n    plt.show()\n   ","execution_count":29,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fe04be5879d76b20f46e4a80daa3ac11121fd244","collapsed":true},"cell_type":"code","source":"for (i,nneigh) in enumerate([20,50,100]):\n    se = manifold.SpectralEmbedding(n_components=2, n_neighbors=nneigh)\n    Y = se.fit_transform(X_eda)\n    plt.figure(1)\n    plt.figure(figsize=(10,10))\n    for (j,cl) in enumerate([29,38,12,1,9,35]):\n        cname = le.tolist().inverse_transform([cl])[0]\n        plt.subplot(3,2,j+1)\n        plt.scatter(Y[y_eda!=cl][:,0], Y[y_eda!=cl][:,1], c='blue', alpha=0.75)\n        plt.scatter(Y[y_eda==cl][:,0], Y[y_eda==cl][:,1], c='red', alpha=0.75)\n        plt.title(cname+' Spect Embed nniegh: '+str(nneigh))\n        plt.axis('tight')\n    plt.tight_layout() \n    plt.show()\n   ","execution_count":30,"outputs":[]},{"metadata":{"_uuid":"e97175d9f12aa8b6426d32da20c2c70555be98eb"},"cell_type":"markdown","source":"It seems that MDS does not give nice results in reducing dimensionality of data. \n\ntSNE and spectral embedding, on the other hand, give some indication that observations of each of the 'Applause', 'Oboe' and 'Tambourine' classes might form its own clusters even in reduced space. "},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"e37017f48d5cf31906dae1c6e9d3d9c65415e39b"},"cell_type":"markdown","source":"## Training classifiers\nOn my laptop it takes around 20 minutes to train random forest (with 1000 trees) and around one hour to train xgb (also with 1000 trees). I train on the whole dataset which gives better resultsif trained on stratified 10-folds with 10% holdout.\n\nHere I set both parameters equal to 10 to decrease computational time."},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"34edbde89181f47d15301b9d3821096cb7533a8e"},"cell_type":"code","source":"def get_proba(clf,X,y,Xtest):\n    clf.fit(X, y)\n    pred_clf_proba = clf.predict_proba(Xtest)\n    pred_clf_classes = [list(clf.classes_[np.argsort(x).tolist()[::-1]]) for x in pred_clf_proba] \n    return pred_clf_proba, pred_clf_classes\n\nclf = RandomForestClassifier(n_estimators=10, class_weight = 'balanced', random_state = 7)\npred_rf_proba, pred_rf_classes = get_proba(clf, Xall, yall, X_scoring)\n\nclf = xgb.XGBClassifier(n_estimators=10, learning_rate=0.05, max_depth=2)\npred_xgb_proba, pred_xgb_classes = get_proba(clf, Xall, yall, X_scoring)\n","execution_count":31,"outputs":[]},{"metadata":{"_uuid":"9ca38d0480f6d1cd459427c6b97f53c881a2f173"},"cell_type":"markdown","source":"## Ensembling predictions\nNext I take geometric mean of predictions produced by xgb and rf to get final ones."},{"metadata":{"trusted":true,"_uuid":"580b2543f7b3a311ff7debb7016d99eb9936d55c","collapsed":true},"cell_type":"code","source":"pred_ens_proba = (pred_rf_proba * pred_xgb_proba ) ** (1/2)\npred_ens_classes = [list(np.argsort(x).tolist()[::-1]) for x in pred_ens_proba] \n\ndf_output = pd.DataFrame(pred_ens_classes)\ndf_output = list(df_output.apply(lambda x: list(le.tolist().inverse_transform(x[0:3])), axis=1))\ndf_output = pd.DataFrame(df_output)\ndf_output['fname'] = score_filelist \ndf_output['label'] = df_output.apply(lambda x: str(x[0]) + ' ' + str(x[1]) +' '+ str(x[2]) , axis=1)\n\nnp.save('pred_ens_proba_simple_ens_newfeatures_libr.npy', pred_ens_proba)\ndf_output[['fname','label']].to_csv('pred_classes_rfxgb_ens.csv', index = False) # gives 0.844 on LB if trained on 1000 trees for both rf and xgb\n\n","execution_count":32,"outputs":[]},{"metadata":{"_uuid":"681ef5818638160a50d4207e5a55a815c1bca95a"},"cell_type":"markdown","source":"# Part 4: CNN\nAs I had zero experience in building CNNs, I checked available kernels and started with Zafar's solution [Beginner's Guide to Audio Data\n](https://www.kaggle.com/fizzbuzz/beginner-s-guide-to-audio-data), which was a great starting point.\n\nI did some meta parameter adjustment (e.g. adding dropouts and increasing number of filters helped to improve accuracy and mapk if compared with initial Zafar's solution).\n"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"23f424e130822126abcf8c64a9997dbd12c87311"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn.model_selection import StratifiedShuffleSplit","execution_count":33,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"996cb66932e83e49ca03231998c4b97de52e05f0","collapsed":true},"cell_type":"code","source":"from keras.models import Model, Sequential # basic class for specifying and training a neural network\nfrom keras.layers import Input, Convolution2D, MaxPooling2D, Dense, Dropout, Flatten, BatchNormalization, Activation, MaxPool2D\nfrom keras.utils import np_utils # utilities for one-hot encoding of ground truth values\nfrom keras.models import save_model, load_model\nfrom keras.callbacks import Callback\nfrom keras import losses, models, optimizers\nfrom keras import backend as K\nfrom keras.models import load_model\nimport h5py as h5py\nfrom keras.callbacks import ModelCheckpoint\n","execution_count":34,"outputs":[]},{"metadata":{"_uuid":"ca21d86eea94c444677a8b49b2c90f0adce50e13"},"cell_type":"markdown","source":"we can calculate mapk after each epoch"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"abe0e99ed87547ede4c3a59c44bb1e2cceb173d2"},"cell_type":"code","source":"class mapk_callback(Callback):\n    def __init__(self,training_data,validation_data):\n        \n        self.x_trn = training_data[0]\n        self.y_trn = training_data[1]\n        self.x_val = validation_data[0]\n        self.y_val = validation_data[1]        \n    \n    def on_train_begin(self, logs={}):\n        return\n    def on_train_end(self, logs={}):\n        return\n    def on_epoch_begin(self, epoch, logs={}):\n        return\n    def on_epoch_end(self, epoch, logs={}):        \n        y_pred_trn = self.model.predict(self.x_trn)\n        y_pred_top3_classes_trn = [np.argsort(x).tolist()[::-1][0:3] for x in y_pred_trn]\n        obs_y_trn = [[x] for x in self.y_trn.tolist()]\n        mapk_score_trn = mapk(obs_y_trn, y_pred_top3_classes_trn,k=3)\n\n        y_pred_val = self.model.predict(self.x_val)\n        y_pred_top3_classes_val = [np.argsort(x).tolist()[::-1][0:3] for x in y_pred_val]\n        obs_y_val = [[x] for x in self.y_val.tolist()]\n        mapk_score_val = mapk(obs_y_val, y_pred_top3_classes_val,k=3)\n        \n        print('\\rmapk: %s - mapk_val: %s' % (str(round(mapk_score_trn,4)),str(round(mapk_score_val,4))),end=100*' '+'\\n')\n        return\n    \n    def on_batch_begin(self, batch, logs={}):\n        return\n\n    def on_batch_end(self, batch, logs={}):\n        return   \n","execution_count":35,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"4ff691e6cd335e96889dd12c80f377732a407fe3"},"cell_type":"code","source":"load = 0\nif load == 1:\n    Xall = np.load('.../Xall_cnn_libr.npy')\n    yall = np.load('.../yall_cnn_libr.npy')\n    le = np.load('.../le_cnn_libr.npy')\n    X_scoring = np.load('.../X_scoring_cnn_libr.npy')\n\n","execution_count":36,"outputs":[]},{"metadata":{"_uuid":"b8b180ff2f21504b18fb67c6c4661d005e1d1004"},"cell_type":"markdown","source":"Next function defines and compiles model. "},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"8259b8442fa99e2197a04c4981a8ae0653df1f8c"},"cell_type":"code","source":"# (9)\n\ndef get_compiled_cnn_model():\n    \n    inp = Input(shape=(173, 40, 1))\n\n    lr = Convolution2D(32, (11,5), padding=\"same\")(inp)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(64, (7,5), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(64, (7,5), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(128, (5,3), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Flatten()(lr)\n    lr = Dense(128)(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    out = Dense(num_classes, activation='softmax')(lr)\n\n    model = Model(inputs=inp, outputs=out)\n    \n    model.compile(loss='categorical_crossentropy', # using the cross-entropy loss function\n              optimizer='adam', # using the Adam optimiser\n              metrics=['accuracy'],\n              )\n    \n    return model\n\n","execution_count":37,"outputs":[]},{"metadata":{"_uuid":"aad2da0d65e64ba433672b9092f8b368899ce493"},"cell_type":"markdown","source":"Next we train the model.\n\nTraining on all data (without validation) 5-10 times and aveaging results using geometric mean helped to improve performance. "},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"3981fb951fad3bff3355e58a1505e023c4892fe1"},"cell_type":"code","source":"\nnum_classes = np.unique(yall).shape[0]\npred_proba_cnn_strat = np.ones(shape=(X_scoring.shape[0],num_classes))\n\nmy_n_splits = 5\nstrat_split = StratifiedShuffleSplit(n_splits=my_n_splits, test_size=0.1, random_state=2)\ni = 0\nrun_training = 0\n\nif run_training == 1:\n    for ii, (train_index, val_index) in enumerate(strat_split.split(Xall, yall)):\n\n        # K.clear_session()\n        print()\n        print('-------------------------- Strata',ii,' --------------------------')\n        print()\n\n        X_train = Xall[train_index]\n        y_train = yall[train_index]\n        X_val = Xall[val_index]\n        y_val = yall[val_index]\n\n        mean = np.mean(X_train, axis=0)\n        std = np.std(X_train, axis=0)\n        X_train = (X_train - mean)/std\n        X_val = (X_val - mean)/std\n        X_scoring_iter = (X_scoring - mean)/std\n\n        y_train_1dim = y_train\n        y_train = np_utils.to_categorical(y_train, num_classes)\n\n        y_val_1dim = y_val\n        y_val = np_utils.to_categorical(y_val, num_classes)\n\n        cnn_model = get_compiled_cnn_model()\n\n        callbacks = [mapk_callback(training_data=(X_train,y_train_1dim),validation_data=(X_val, y_val_1dim))\n                     , ModelCheckpoint('.../weights_valid.model',\n                                        monitor='val_acc',\n                                        verbose=1,\n                                        save_best_only=True,\n                                        save_weights_only=True)\n        ]\n\n        history = cnn_model.fit(X_train, y_train\n                                , validation_data=(X_val, y_val)\n                                , callbacks = callbacks\n                                , verbose=1\n                                , batch_size = 64\n                                , epochs = 1)\n\n        pred_iter = cnn_model.predict(X_scoring_iter)\n        np.save('.../pred_cnn_proba_strat_%d.npy'%ii, pred_iter)\n        pred_proba_cnn_strat = pred_proba_cnn_strat * pred_iter\n\n    df_output = pd.DataFrame(pred_cnn_classes)\n    df_output = list(df_output.apply(lambda x: list(le.tolist().inverse_transform(x[0:3])), axis=1))\n    df_output = pd.DataFrame(df_output)\n    df_output['fname'] = score_filelist \n    df_output['label'] = df_output.apply(lambda x: str(x[0]) + ' ' + str(x[1]) +' '+ str(x[2]) , axis=1)\n\n    np.save('.../predictions_cnn_final.npy', pred_cnn_proba)\n    df_output[['fname','label']].to_csv('.../predictions_cnn_strat.csv', index = False)\n\n        ","execution_count":38,"outputs":[]},{"metadata":{"_uuid":"f0f8deb366bd9a7f17ed2b4a73e890d9831d9ed4"},"cell_type":"markdown","source":"# Part 5: CNN with two inputs\n\nAnother interesting thing which I tried to use (and it helped to increase accuracy and mapk) with to use xgb/rf dataset as another input into CNN with more or less same architecture as above.\n\nAs training process is the same as above, I show here only architecture.\n"},{"metadata":{"trusted":true,"_uuid":"5c811c7feefc4590f3823cefdb208a5da7c36fc7","collapsed":true},"cell_type":"code","source":"load_data = 0\n\nif load_data == 1:\n    Xall = np.load('.../Xall_cnn_libr.npy')\n    yall = np.load('.../yall_cnn_libr.npy')\n    le = np.load('.../le_cnn_libr.npy')\n    X_scoring = np.load('.../X_scoring_cnn_libr.npy')\n\n    Zall = np.load('.../Xall_ens.npy')\n    Z_scoring = np.load('.../X_ens.npy')\n","execution_count":39,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"75dc9a332af9f109c8a7156d5dab5529a0720c88"},"cell_type":"code","source":"def get_compiled_CNN_model():\n    inp1 = Input(shape=(173,40,1))\n\n    lr = Convolution2D(32, (10,4), padding=\"same\")(inp1)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(32, (10,4), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(32, (10,4), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    lr = Convolution2D(32, (10,4), padding=\"same\")(lr)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    lr = MaxPool2D()(lr)\n    lr = Dropout(0.1)(lr)\n\n    flat = Flatten()(lr)\n\n    inp2 = Input(shape = (Z_train.shape[1], ))\n    den2 = Dense(64)(inp2)\n    concatFeatures = Concatenate(axis = -1)([flat, den2])\n\n    lr = Dense(256)(concatFeatures)\n    lr = BatchNormalization()(lr)\n    lr = Activation(\"relu\")(lr)\n    out = Dense(num_classes, activation='softmax')(lr)\n\n    model = Model(inputs=[inp1,inp2], outputs=out)\n    \n    model.compile(loss='categorical_crossentropy', # using the cross-entropy loss function\n                  optimizer=optimizers.Adam(lr=0.001),\n                  metrics=['accuracy'],\n                  )    \n    \n    return model\n","execution_count":40,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"5096b86e0d4974633487be1f49a1166aab2da13e"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"7b8378289891efdb97103e660e7014e6b2da8718"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}