{"cells":[{"metadata":{},"cell_type":"markdown","source":"# LANL Earthquake Prediction"},{"metadata":{},"cell_type":"markdown","source":"**Note:** I had to chop the features engineered for the train dataset to just 100 segments so the kernel could finish running in under 9 hours.  So, the private and public score for this kernel when we engineer features for the complete 2000 segments is 2.54742 and 1.90667.  Which would place first on the private leaderboard."},{"metadata":{},"cell_type":"markdown","source":"\n### 1: Data Exploration/Manipulation & Feature Engineering"},{"metadata":{},"cell_type":"markdown","source":"Load data analysis modules we may require for the EDA:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\n#Load data analysis/plotting modules\nfrom tqdm import tqdm_notebook as tqdm\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport lightgbm as lgb\nimport pandas as pd\nimport seaborn as sns\nimport numpy as np\nimport matplotlib\nimport scipy\nimport numpy\nimport os\n\n\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading train data and the names of the files in the test data folder:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# only read first 300 million rows to conserve memory.\ndf_train = pd.read_csv(\"../input/train.csv\", nrows=150e6)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ntest_files = os.listdir(\"../input/test/\")\nprint(test_files[0:5])\n# load sample_submission.csv\ndf_samplesub = pd.read_csv(\"../input/sample_submission.csv\", index_col='seg_id')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are approximately 4194 'segment-sized' partitions in the train data.  This could be a good way of segmenting the training data to generate features. (as per the benchmark)."},{"metadata":{},"cell_type":"markdown","source":"Changing name of columns to simpler abbreviations and increase the number of decimals displayed by pandas."},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.options.display.precision = 15\n\ndf_train.rename({'acoustic_data': 'ad', 'time_to_failure': 'ttf'}, axis=1, inplace=True)\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can plot 1% of the train data to get a big picture of the acoustic data and time to failure behaviour.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# from \"LANL Earthquake EDA and Prediction\" by Gabriel Preda\n\ndf_train_ad_1per_all = df_train['ad'].values[::100]\ndf_train_ttf_1per_all = df_train['ttf'].values[::100]\n\nfig, ax1 = plt.subplots(figsize=(12, 8))\nplt.title(\"Acoustic Data and Time to Failure\")\nplt.plot(df_train_ad_1per_all, color='b')\nax1.set_ylabel('acoustic data', color='b')\nplt.legend(['acoustic data'], loc=(0.01, 0.95))\nax2 = ax1.twinx()\nplt.plot(df_train_ttf_1per_all, color='r')\nax2.set_ylabel('time to failure', color='r')\nplt.legend(['time to failure'], loc=(0.01, 0.9))\nplt.grid(True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The time to failure reaches 0 sixteen times in the train dataset, this corresponds to 16 lab earthquakes in total.  \n\nI've read a number of interesting kaggle kernels that study the properties of the signal and produce an assortment of statistical and rolling features for the signal.  In this kernel, I'd like to focus specifically on frequency domain analysis for EDA and feature engineering.  \n\nLet's take a \"segment-sized\" sample of the signal and study the frequency domain behaviour of the sample."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train_sample = df_train[0:150000].ad\ndf_train_sample.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The power spectrum for the sample signal represents the degree to which the frequencies, that the signal is composed of, contribute to the overall power of the signal.\n\nThere are a number of estimators for the power spectrum:  PSD (from DFT), spectrogram (Blackman-Tukey method), and the Welch method.  \n\nLet's compare the power spectra for each method."},{"metadata":{"trusted":true},"cell_type":"code","source":"#get the power spectrum for the sample signal\nimport scipy.signal\nfs = 4e6\ndef DFT(sig_in):\n\n    fs = 4e6 #4e6 is the signal sample rate\n    ps = np.abs(np.fft.fft(sig_in))**2\n    freqs = np.fft.fftfreq(sig_in.size, 1/fs)\n    return freqs, ps\n\nfreqs, ps = DFT(df_train_sample)\nfreqs_pdg, ps_pdg = scipy.signal.periodogram(df_train_sample.T, fs, scaling = 'spectrum', window = 'triang')\nfreqs_welch, ps_welch = scipy.signal.welch(df_train_sample.T, fs, scaling = 'spectrum')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DFT\nfig = plt.figure(figsize=(30, 20))\nax = fig.add_subplot(3, 1, 1)\nax.margins(x=0.003)\nplt.plot(freqs,ps)\nplt.xlim(left=0)\nplt.ylim(0,5e8)\nplt.title('Power Spectrum (DFT)', fontsize=24, loc='center')\nplt.xlabel('Frequency (Hz)', fontsize = 18)\nplt.ylabel('Power (signal strength)', fontsize = 18)\n\n# Spectrogram \n#fig = plt.figure(figsize=(30, 20))\nax = fig.add_subplot(3, 1, 2)\nax.margins(x=0.003)\nplt.plot(freqs_pdg, ps_pdg)\nplt.title('Power Spectrum (Periodogram)', fontsize=24, loc='center')\nplt.xlabel('Frequency (Hz)', fontsize = 18)\nplt.ylabel('Power (signal strength)', fontsize = 18)\n\n# Welch method\n#fig = plt.figure(figsize=(30, 20))\nax = fig.add_subplot(3, 1, 3)\nax.margins(x=0.003)\nplt.plot(freqs_welch, ps_welch)\nplt.title('Power Spectrum (Welch method)', fontsize=24, loc='center')\nplt.xlabel('Frequency (Hz)', fontsize = 18)\nplt.ylabel('Power (signal strength)', fontsize = 18)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For the sample signal, the power spectrum (PS) from the DFT shows a significant peak at ~70 kHz with a lot of the most powerful frequencies centered around it along with a number of scattered local peaks in the higher frequency range.\n\nThe DFT method for the power spectrum is susceptible to large variance in the PS coefficients, so there may be peaks in the plot that are due to noise rather than actual features of the signal.\n\nThe PS for the spectrogram/Welch method reduce the variance at the cost of detail/resolution in the PS. \n\nThe periodogram method seems to be a suitable compromise between the reduced variance vs reduced detail, but the DFT method has some interesting minor peaks in the high frequency range that might yield useful features.  We'll use the DFT method going forward for EDA and feature generation.\n\n"},{"metadata":{},"cell_type":"markdown","source":"Let's plot the power spectrum for a number of sampled signals from the training data and detect peaks using the PeakUtils library to see how they compare."},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install PeakUtils","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import peakutils\n\ni = 29\ndf_train_sampled = df_train[i*150000:(i*150000)+150000].ad\n\nfreqs, ps = DFT(df_train_sampled)\npeaks_ind1 = peakutils.indexes(ps, thres=0.0001, min_dist=200 )\nplt.margins(x=0.003)\nplt.plot(freqs, ps)\nfor p in peaks_ind1:\n    plt.scatter(freqs[p], ps[p], marker='s', color='red', label='v1')\nplt.xlim(left=0)\nplt.ylim([0,0.4e12])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"\ndef plot_ps():\n    sns.set_style('whitegrid')\n    plt.figure()\n    fig, ax = plt.subplots(3, 4, figsize=(15, 15))\n    fig.subplots_adjust(hspace=0.4, wspace=0.3)\n    \n    for i in tqdm(range(12)):\n        df_train_sampled = df_train[i*150000:(i*150000)+150000].ad\n        freqs, ps = DFT(df_train_sampled)\n        peaks_ind = peakutils.indexes(ps, thres=1e-4, min_dist=200)\n        plt.subplot(3, 4, i+1)\n        plt.margins(x=0.003)\n        plt.plot(freqs, ps)\n        for p in peaks_ind:\n            plt.scatter(freqs[p], ps[p], marker='s', color='red', label='v1')\n        plt.xlim(0,500000)\n        plt.ylim(0,1e9)\n    plt.show()\n    \nplot_ps()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The plots show a wide range of peaks at a number of different frequencies as the signal evolves in time.  \n\nLet's also plot the spectrogram for a sample signal.  The spectrogram is comprised of a series of periodograms that are evaluated at discrete periods of time to form a 3-dimensional plot of frequency vs. power with time along the x-axis."},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.signal import spectrogram\n\ndef make_spectrogram(sig_in):\n    nperseg = 512 # default 256\n    noverlap = nperseg // 4 # default: nperseg // 8\n    fs = 4e6  # raw signal sample rate is 4MHz\n    window = 'boxcar'\n    scaling = 'spectrum' # {'density', 'spectrum'}\n    detrend = 'linear' # {'linear', 'constant', False}\n    f, t, Sxx = spectrogram(sig_in.T, nperseg=nperseg, noverlap=noverlap,\n                                   fs=fs, window=window,\n                                   scaling=scaling, detrend=detrend)\n    return f, t, Sxx\n\n\nf, t, Sxx_out = make_spectrogram(df_train_sample)\n\nprint('Sxx_out:', Sxx_out.shape)\nprint('f:', f.shape)\nprint('t:', t.shape)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(30, 20))\nax = fig.add_subplot(4, 1, 1)\nax.margins(x=0.003)\nplt.plot(df_train_sample)\nplt.title('sample signal:', fontsize=18, loc='left')\n\n\nax = fig.add_subplot(4, 1, 2)\ncmap = plt.get_cmap('magma')\nspec = plt.pcolormesh(t, f, Sxx_out, cmap=cmap, norm = matplotlib.colors.Normalize(0,1))\nplt.title('normalized log spectrogram:',\n          fontsize=18, loc='left')\nax.set_ylim([0,5e5])\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's plot a slice of the spectrogram and detect peaks."},{"metadata":{"trusted":true},"cell_type":"code","source":"signal = Sxx_out[0:257,10]\npeaks_ind = peakutils.indexes(signal, thres=0.05, min_dist=1)\nplt.plot(f[0:257], signal)\nfor p in peaks_ind:\n    plt.scatter(f[p], signal[p], marker='s', color='red', label='v1')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Engineering Features"},{"metadata":{},"cell_type":"markdown","source":"Now, we can engineer features using the information we've gathered in the frequency domain.  \nLet's first create some computation functions for the frequency analysis."},{"metadata":{"trusted":true},"cell_type":"code","source":"  \n# getting errors during feature engineering where peaks couldn't be found for certain segments, \n# so we'll set the threshold really low to avoid that possibility.\n\ndef peak_finder(freqs, ps ,thres=0.00001, min_dist=200):\n    cb = np.array([d[0] for d in ps])\n    peaks_ind = peakutils.indexes(cb, thres=thres, min_dist=min_dist)\n    peaks_data = pd.DataFrame(index=range(peaks_ind.size), dtype=np.float64, columns=['peak_freq','peak_power'])\n    for p in range(peaks_ind.size):\n            peaks_data['peak_freq'].loc[p] = freqs[peaks_ind[p]]\n            peaks_data['peak_power'].loc[p] = ps[peaks_ind[p]][0]\n    peaks_data.replace([\"NaN\", 'NaT'], np.nan, inplace = True)\n    peaks_data = peaks_data.dropna()\n    return peaks_data\n\n\ndef specgram_peak_slice(Sxx_out, f, t, thres=0.00001, min_dist=1):\n    signal = Sxx_out[0][0:257,t]\n    peaks_ind = peakutils.indexes(signal, thres=thres, min_dist=min_dist)\n    peaks_data = pd.DataFrame(index=range(peaks_ind.size), dtype=np.float64, columns=['peak_freq','peak_power'])\n    for p in range(peaks_ind.size):\n        if f[peaks_ind[p]] > 0:\n            peaks_data['peak_freq'].loc[p] = f[peaks_ind[p]]\n            peaks_data['peak_power'].loc[p] = Sxx_out[0][peaks_ind[p]][t]\n            \n    peaks_data.replace([\"NaN\", 'NaT'], np.nan, inplace = True)\n    peaks_data = peaks_data.dropna()\n    return peaks_data\n\ndef specgram_temporal_features(power, f):\n    tmaxfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tminfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmeanfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tstdfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmaxpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tminpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmeanpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tstdpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmadfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tkurtfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tskewfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmedfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmadpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tkurtpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tskewpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmedpow = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tmaxpowfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    tminpowfreq = pd.DataFrame(index=range(power.shape[1]), dtype=np.float64, columns=['val'])\n    \n    for seg in range(power.shape[1]):\n        peak_slice = specgram_peak_slice(power, f, seg)\n        \n        tmaxfreq['val'].loc[seg] = peak_slice['peak_freq'].max()\n        tminfreq['val'].loc[seg] = peak_slice['peak_freq'].min()\n        tmeanfreq['val'].loc[seg] = peak_slice['peak_freq'].mean()\n        tstdfreq['val'].loc[seg] = peak_slice['peak_freq'].std()\n\n        tmaxpow['val'].loc[seg] = peak_slice['peak_power'].max()\n        tminpow['val'].loc[seg] = peak_slice['peak_power'].min()\n        tmeanpow['val'].loc[seg] = peak_slice['peak_power'].mean()\n        tstdpow['val'].loc[seg] = peak_slice['peak_power'].std()\n\n        tmadfreq['val'].loc[seg] = peak_slice['peak_freq'].mad()\n        tkurtfreq['val'].loc[seg] = peak_slice['peak_freq'].kurtosis()\n        tskewfreq['val'].loc[seg] = peak_slice['peak_freq'].skew()\n        tmedfreq['val'].loc[seg] = peak_slice['peak_freq'].median()\n\n        tmadpow['val'].loc[seg] = peak_slice['peak_power'].mad()\n        tkurtpow['val'].loc[seg] = peak_slice['peak_power'].kurtosis()\n        tskewpow['val'].loc[seg] = peak_slice['peak_power'].skew()\n        tmedpow['val'].loc[seg] = peak_slice['peak_power'].median()\n\n        tmaxpowfreq['val'].loc[seg] = peak_slice['peak_freq'][peak_slice['peak_power'].idxmax()]\n        tminpowfreq['val'].loc[seg] = peak_slice['peak_freq'][peak_slice['peak_power'].idxmin()]\n\n    return  tmaxfreq, tminfreq, tmeanfreq, tstdfreq, tmaxpow, tminpow, tmeanpow, tstdpow, tmadfreq, tkurtfreq, tskewfreq, tmedfreq, tmadpow, tkurtpow, tskewpow, tmedpow, tmaxpowfreq, tminpowfreq\n        \ndef add_trend_feature(arr, abs_values=False):\n    idx = np.array(range(len(arr)))\n    if abs_values:\n        arr = np.abs(arr)\n    lr = LinearRegression()\n    lr.fit(idx.reshape(-1, 1), arr)\n    return lr.coef_[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"segments = 500  #could only load 1000 segments worth of training data due to RAM restrictions.\nX_train = pd.DataFrame(index=range(segments), dtype=np.float64)\ny_train = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['ttf'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def feature_generation(seg_id, seg, X):\n        \n    xc = pd.DataFrame(seg['ad'])\n    freqs, ps = DFT(xc)\n    avg_peaks = peak_finder(freqs, ps) #sorted at low to high frequency\n    f, t, power = make_spectrogram(xc)\n    \n    \n#     average FFT peak features\n    \n    X.loc[seg_id, 'aFFT_max'] = avg_peaks['peak_freq'].max()\n    X.loc[seg_id, 'aFFT_min'] = avg_peaks['peak_freq'].min()\n    X.loc[seg_id, 'aFFT_mean'] = avg_peaks['peak_freq'].mean()\n    X.loc[seg_id, 'aFFT_std'] = avg_peaks['peak_freq'].std()\n    \n    \n    X.loc[seg_id, 'aFFT_max_pow'] = avg_peaks['peak_power'].max()\n    X.loc[seg_id, 'aFFT_min_pow'] = avg_peaks['peak_power'].min()\n    X.loc[seg_id, 'aFFT_mean_pow'] = avg_peaks['peak_power'].mean()\n    X.loc[seg_id, 'aFFT_std_pow'] = avg_peaks['peak_power'].std()\n    \n    X.loc[seg_id, 'aFFT_mad_freq'] = avg_peaks['peak_freq'].mad()\n    X.loc[seg_id, 'aFFT_kurt_freq'] = avg_peaks['peak_freq'].kurtosis()\n    X.loc[seg_id, 'aFFT_skew_freq'] = avg_peaks['peak_freq'].skew()\n    X.loc[seg_id, 'aFFT_med_freq'] = avg_peaks['peak_freq'].median()\n    \n    X.loc[seg_id, 'aFFT_mad_pow'] = avg_peaks['peak_power'].mad()\n    X.loc[seg_id, 'aFFT_kurt_pow'] = avg_peaks['peak_power'].kurtosis()\n    X.loc[seg_id, 'aFFT_skew_pow'] = avg_peaks['peak_power'].skew()\n    X.loc[seg_id, 'aFFT_med_pow'] = avg_peaks['peak_power'].median()\n    \n    \n    X.loc[seg_id, 'aFFT_max_pow_freq'] = avg_peaks['peak_freq'][avg_peaks['peak_power'].idxmax()]\n    X.loc[seg_id, 'aFFT_min_pow_freq'] = avg_peaks['peak_freq'][avg_peaks['peak_power'].idxmin()]\n    \n    X.loc[seg_id, 'aFFT_mean_pow_q99'] = np.quantile(avg_peaks['peak_power'], 0.99)\n    X.loc[seg_id, 'aFFT_mean_pow_q95'] = np.quantile(avg_peaks['peak_power'], 0.95)\n    X.loc[seg_id, 'aFFT_mean_pow_q05'] = np.quantile(avg_peaks['peak_power'], 0.05)\n    X.loc[seg_id, 'aFFT_mean_pow_q01'] = np.quantile(avg_peaks['peak_power'], 0.01)\n    \n    X.loc[seg_id, 'aFFT_mean_freq_q99'] = np.quantile(avg_peaks['peak_freq'], 0.99)\n    X.loc[seg_id, 'aFFT_mean_freq_q95'] = np.quantile(avg_peaks['peak_freq'], 0.95)\n    X.loc[seg_id, 'aFFT_mean_freq_q05'] = np.quantile(avg_peaks['peak_freq'], 0.05)\n    X.loc[seg_id, 'aFFT_mean_freq_q01'] = np.quantile(avg_peaks['peak_freq'], 0.01)\n    \n#     Spectrogram FFT features (focusing on temporal characteristics and trend features)\n\n    \n    tmaxfreq, tminfreq, tmeanfreq, tstdfreq, tmaxpow, tminpow, tmeanpow, tstdpow, tmadfreq, tkurtfreq, tskewfreq, tmedfreq, tmadpow, tkurtpow, tskewpow, tmedpow, tmaxpowfreq, tminpowfreq  = specgram_temporal_features(power, f)\n    \n    X.loc[seg_id, 'tmaxfreq_trend'] = add_trend_feature(tmaxfreq['val'])\n    X.loc[seg_id, 'tminfreq_trend'] = add_trend_feature(tminfreq['val'])\n    X.loc[seg_id, 'tmeanfreq_trend'] = add_trend_feature(tmeanfreq['val'])\n    X.loc[seg_id, 'tstdfreq_trend'] = add_trend_feature(tstdfreq['val'])\n    X.loc[seg_id, 'tmaxpow_trend'] = add_trend_feature(tmaxpow['val'])\n    X.loc[seg_id, 'tminpow_trend'] = add_trend_feature(tminpow['val'])\n    X.loc[seg_id, 'tmeanpow_trend'] = add_trend_feature(tmeanpow['val'])\n    X.loc[seg_id, 'tstdpow_trend'] = add_trend_feature(tstdpow['val'])\n    X.loc[seg_id, 'tmadfreq_trend'] = add_trend_feature(tmadfreq['val'])\n    X.loc[seg_id, 'tkurtfreq_trend'] = add_trend_feature(tkurtfreq['val'])\n    X.loc[seg_id, 'tskewfreq_trend'] = add_trend_feature(tskewfreq['val'])\n    X.loc[seg_id, 'tmedfreq_trend'] = add_trend_feature(tmedfreq['val'])\n    X.loc[seg_id, 'tmadpow_trend'] = add_trend_feature(tmadpow['val'])\n    X.loc[seg_id, 'tkurtpow_trend'] = add_trend_feature(tkurtpow['val'])\n    X.loc[seg_id, 'tskewpow_trend'] = add_trend_feature(tskewpow['val'])\n    X.loc[seg_id, 'tmedpow_trend'] = add_trend_feature(tmedpow['val'])\n    X.loc[seg_id, 'tmaxpowfreq_trend'] = add_trend_feature(tmaxpowfreq['val'])\n    X.loc[seg_id, 'tmaxpowfreq_trend'] = add_trend_feature(tminpowfreq['val'])\n    \n    X.loc[seg_id, 'tmaxfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmaxfreq['val'])))\n    X.loc[seg_id, 'tminfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tminfreq['val'])))\n    X.loc[seg_id, 'tmeanfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmeanfreq['val'])))\n    X.loc[seg_id, 'tstdfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tstdfreq['val'])))\n    X.loc[seg_id, 'tmaxpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmaxpow['val'])))\n    X.loc[seg_id, 'tminpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tminpow['val'])))\n    X.loc[seg_id, 'tmeanpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmeanpow['val'])))\n    X.loc[seg_id, 'tstdpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tstdpow['val'])))\n    X.loc[seg_id, 'tmadfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmadfreq['val'])))\n    X.loc[seg_id, 'tkurtfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tkurtfreq['val'])))\n    X.loc[seg_id, 'tskewfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tskewfreq['val'])))\n    X.loc[seg_id, 'tmedfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmedfreq['val'])))\n    X.loc[seg_id, 'tmadpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmadpow['val'])))\n    X.loc[seg_id, 'tkurtpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tkurtpow['val'])))\n    X.loc[seg_id, 'tskewpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tskewpow['val'])))\n    X.loc[seg_id, 'tmedpow_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmedpow['val'])))\n    X.loc[seg_id, 'tmaxpowfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tmaxpowfreq['val'])))\n    X.loc[seg_id, 'tmaxpowfreq_grad_skew'] = scipy.stats.skew(np.diff(np.diff(tminpowfreq['val'])))\n    \n    X.loc[seg_id, 'tmaxfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmaxfreq['val'])))\n    X.loc[seg_id, 'tminfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tminfreq['val'])))\n    X.loc[seg_id, 'tmeanfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmeanfreq['val'])))\n    X.loc[seg_id, 'tstdfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tstdfreq['val'])))\n    X.loc[seg_id, 'tmaxpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmaxpow['val'])))\n    X.loc[seg_id, 'tminpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tminpow['val'])))\n    X.loc[seg_id, 'tmeanpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmeanpow['val'])))\n    X.loc[seg_id, 'tstdpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tstdpow['val'])))\n    X.loc[seg_id, 'tmadfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmadfreq['val'])))\n    X.loc[seg_id, 'tkurtfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tkurtfreq['val'])))\n    X.loc[seg_id, 'tskewfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tskewfreq['val'])))\n    X.loc[seg_id, 'tmedfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmedfreq['val'])))\n    X.loc[seg_id, 'tmadpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmadpow['val'])))\n    X.loc[seg_id, 'tkurtpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tkurtpow['val'])))\n    X.loc[seg_id, 'tskewpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tskewpow['val'])))\n    X.loc[seg_id, 'tmedpow_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmedpow['val'])))\n    X.loc[seg_id, 'tmaxpowfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tmaxpowfreq['val'])))\n    X.loc[seg_id, 'tmaxpowfreq_grad_kurt'] = scipy.stats.kurtosis(np.diff(np.diff(tminpowfreq['val'])))\n    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The original "},{"metadata":{"trusted":true},"cell_type":"code","source":"# iterate over 100 segments\n\nrows = 150000\nfor seg_id in tqdm(range(100)):\n    seg = df_train.iloc[seg_id*rows:seg_id*rows+rows]\n    feature_generation(seg_id, seg, X_train)\n    y_train.loc[seg_id, 'ttf'] = seg['ttf'].values[-1]\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Scaling the data may help so we don't overfit for large powered frequencies."},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_train)\nscaled_X_train = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_X_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we generate features and scale the test data."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = pd.DataFrame(columns=X_train.columns, dtype=np.float64, index=df_samplesub.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_PATH = \"../input/test/\"\nfor seg_id in tqdm(X_test.index):\n    \n    seg = pd.read_csv(test_PATH + seg_id + '.csv')\n    seg.rename({'acoustic_data': 'ad'}, axis=1, inplace=True)\n    feature_generation(seg_id, seg, X_test)\n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_X_test = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_X_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model Train/Prediction and Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nn_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\ntrain_columns = scaled_X_train.columns.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'num_leaves': 51,\n         'min_data_in_leaf': 10, \n         'objective':'regression',\n         'max_depth': -1,\n         'learning_rate': 0.001,\n         \"boosting\": \"gbdt\",\n         \"feature_fraction\": 0.91,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.91,\n         \"bagging_seed\": 42,\n         \"metric\": 'mae',\n         \"lambda_l1\": 0.1,\n         \"verbosity\": -1,\n         \"nthread\": -1,\n         \"random_state\": 42}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"oof = np.zeros(len(scaled_X_train))\npredictions = np.zeros(len(scaled_X_test))\nfeature_importance_df = pd.DataFrame()\n#run model\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_X_train,y_train.values)):\n    strLog = \"fold {}\".format(fold_)\n    print(strLog)\n    \n    X_tr, X_val = scaled_X_train.iloc[trn_idx], scaled_X_train.iloc[val_idx]\n    y_tr, y_val = y_train.iloc[trn_idx], y_train.iloc[val_idx]\n\n    model = lgb.LGBMRegressor(**params, n_estimators = 20000, n_jobs = -1)\n    model.fit(X_tr, \n              y_tr, \n              eval_set=[(X_tr, y_tr), (X_val, y_val)], \n              eval_metric='mae',\n              verbose=1000, \n              early_stopping_rounds=500)\n    oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)\n    #feature importance\n    fold_importance_df = pd.DataFrame()\n    fold_importance_df[\"Feature\"] = train_columns\n    fold_importance_df[\"importance\"] = model.feature_importances_[:len(train_columns)]\n    fold_importance_df[\"fold\"] = fold_ + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n    #predictions\n    predictions += model.predict(scaled_X_test, num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:200].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,26))\nsns.barplot(x=\"importance\", y=\"Feature\", data=best_features.sort_values(by=\"importance\",ascending=False))\nplt.title('LightGBM Features (averaged over folds)')\nplt.tight_layout()\nplt.savefig('lgbm_importances.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Submission "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_samplesub.time_to_failure = predictions\ndf_samplesub.to_csv('submission.csv',index=True)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.8"}},"nbformat":4,"nbformat_minor":1}