{"cells":[{"metadata":{},"cell_type":"markdown","source":" 3 questions:\n\n# What if treat it as sound?\n\n# What if we denoise the signal?\n\n# Can we speed up the computation (even more)\n\n\n## Aditional : automated stacking with vecstack"},{"metadata":{},"cell_type":"markdown","source":"First question:\n\nIts time series, why not treat it as sound and do feature engineering of audio data.\nVery nice paper if you want a comprehensive survey https://www.cambridge.org/core/services/aop-cambridge-core/content/view/S2048770314000122\n\n\n\nFeatures I used:\nstft,mfccs,chroma,mel,tonnetz. What are they?\n\n1. stft- Short-time Fourier transform. -Fourier-related transform used to determine the sinusoidal frequency and phase content of local sections of a signal as it changes over time\n\n2. MFCCS- feature in speech recognition https://en.wikipedia.org/wiki/Mel-frequency_cepstrum\n\n3. chroma https://en.wikipedia.org/wiki/Chroma_feature\n\n4. mel- https://librosa.github.io/librosa/generated/librosa.feature.melspectrogram.html\n\n5. tonnetz- https://sites.google.com/site/tonalintervalspace/\n"},{"metadata":{},"cell_type":"markdown","source":"Second question from https://www.kaggle.com/jackvial/dwt-signal-denoising:\n\nOK, so we could seperate signal from the noise and let the algos focus on the essentials without extracting the essentials trough new features\n\n\nBUT, the key is. Are we doing it correctly? In the sence that are we also removing the important signals from the data along the way?\n\n\nWell from local CV they are very close. I did not even modify (expect for sampling rate) the denoising functions. I would assume It could be improved. In the worst case extract the most useful engineered features from the DENOISED data set and add them to the new features"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Third question:\n    \nWhy do we even care about speeding it up more? \nDataset is huge, in order to iterate fast and try new features, techniques etc... we should speed up the code.\n\nFirst of all THANK YOU Ashi for the nice class in the discussions. That made me wonder what are some other techniques (parallelisation in python) that could be of help. I made a [compilation of that subject](https://www.kaggle.com/zikazika/parallelisation-in-python).\n\nWhere I investigated:\n\n1. multiprocessing library- Difference Pool and Process\n2. Numba decorators\n3. Joblib parallelisation (there is an option for threading)\n4. Some relevant topics and questions that came along\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T12:13:54.488961Z","start_time":"2019-04-19T12:13:46.577734Z"},"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.metrics import mean_absolute_error\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport xgboost as xgb\nimport time\nimport datetime\nfrom catboost import CatBoostRegressor\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.linear_model import LinearRegression\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\nfrom scipy import stats\nfrom sklearn.kernel_ridge import KernelRidge\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nfrom joblib import Parallel, delayed\nimport librosa\n\nimport pandas as pd\nfrom sklearn.datasets import load_iris\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom xgboost import XGBRegressor\nfrom catboost import CatBoostRegressor\nfrom vecstack import stacking\nfrom lightgbm import LGBMRegressor\n\n\nimport keras\nfrom keras.layers import Dense\nfrom keras.layers import BatchNormalization\nfrom keras.layers import Activation\nfrom keras.layers import Dropout\nfrom keras.models import Sequential\nfrom keras import optimizers\n\n\nfrom scipy import stats\nfrom sklearn.svm import NuSVR, SVR\nfrom scipy.stats import kurtosis\nfrom keras.wrappers.scikit_learn import KerasRegressor\nfrom numba import jit, int32\n\nimport pyarrow.parquet as pq\nimport gc\nimport pywt\nfrom statsmodels.robust import mad\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import butter\n\nfrom sklearn.preprocessing import MinMaxScaler","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T12:14:46.199919Z","start_time":"2019-04-19T12:14:46.191931Z"},"trusted":true},"cell_type":"code","source":"def maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef high_pass_filter(x, low_cutoff=1000, sample_rate=4000000):\n    \"\"\"\n    From @randxie https://github.com/randxie/Kaggle-VSB-Baseline/blob/master/src/utils/util_signal.py\n    Modified to work with scipy version 1.1.0 which does not have the fs parameter\n    \"\"\"\n    \n    # nyquist frequency is half the sample rate https://en.wikipedia.org/wiki/Nyquist_frequency\n    nyquist = 0.5 * sample_rate\n    norm_low_cutoff = low_cutoff / nyquist\n    \n    # Fault pattern usually exists in high frequency band. According to literature, the pattern is visible above 10^4 Hz.\n    # scipy version 1.2.0\n    #sos = butter(10, low_freq, btype='hp', fs=sample_fs, output='sos')\n    \n    # scipy version 1.1.0\n    sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n    filtered_sig = signal.sosfilt(sos, x)\n\n    return filtered_sig\n\ndef denoise_signal( x, wavelet='db4', level=1):\n    \"\"\"\n    1. Adapted from waveletSmooth function found here:\n    http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n    2. Threshold equation and using hard mode in threshold as mentioned\n    in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:\n    http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    \"\"\"\n    \n    # Decompose to get the wavelet coefficients\n    coeff = pywt.wavedec( x, wavelet, mode=\"per\" )\n    \n    # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation\n    sigma = (1/0.6745) * maddest( coeff[-level] )\n\n    # Calculte the univeral threshold\n    uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) )\n    coeff[1:] = ( pywt.threshold( i, value=uthresh, mode='hard' ) for i in coeff[1:] )\n    \n    # Reconstruct the signal using the thresholded coefficients\n    return pywt.waverec( coeff, wavelet, mode='per' )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T12:14:58.811382Z","start_time":"2019-04-19T12:14:58.759706Z"},"code_folding":[],"trusted":true},"cell_type":"code","source":"class FeatureGenerator(object):\n    def __init__(self, dtype, n_jobs=1, chunk_size=None):\n        self.chunk_size = chunk_size\n        self.dtype = dtype\n        self.filename = None\n        self.n_jobs = n_jobs\n        self.test_files = []\n        if self.dtype == 'train':\n            self.filename = '../input/train.csv'\n            self.total_data = int(629145481 / self.chunk_size)\n        else:\n            submission = pd.read_csv('../input/sample_submission.csv')\n            for seg_id in submission.seg_id.values:\n                self.test_files.append((seg_id, '../input/test/' + seg_id + '.csv'))\n            self.total_data = int(len(submission))\n\n    def read_chunks(self):\n        if self.dtype == 'train':\n            iter_df = pd.read_csv(self.filename, iterator=True, chunksize=self.chunk_size,\n                                  dtype={'acoustic_data': np.float64, 'time_to_failure': np.float64})\n            for counter, df in enumerate(iter_df):\n                x = df.acoustic_data.values\n                y = df.time_to_failure.values[-1]\n                seg_id = 'train_' + str(counter)\n                yield seg_id, x, y # its a generator\n        else:\n            for seg_id, f in self.test_files:\n                df = pd.read_csv(f, dtype={'acoustic_data': np.float64})\n                x = df.acoustic_data.values\n                yield seg_id, x, -999\n\n    def features(self, x, y, seg_id):\n        \n        x_hp = high_pass_filter(x, low_cutoff=10000, sample_rate=4000000)\n    \n\n        x = denoise_signal(x_hp, wavelet='haar', level=1)\n        feature_dict = dict()\n        feature_dict['target'] = y\n        feature_dict['seg_id'] = seg_id\n        #audio\n  \n        sample_rate=4000000\n        feature_dict['stft1']=np.abs(librosa.stft(x))\n        feature_dict['stft']= np.mean(np.abs(librosa.stft(x)))\n\n        feature_dict['mfccs']=np.mean(np.mean(librosa.feature.mfcc(y=x, sr=sample_rate, n_mfcc=40).T,axis=0))\n\n        feature_dict['chroma']=np.mean(np.mean((librosa.feature.chroma_stft(S=feature_dict['stft1'], sr=sample_rate).T)))\n\n        feature_dict['mel']=np.mean(np.mean(librosa.feature.melspectrogram(x, sr=sample_rate).T))\n\n        feature_dict['tonnetz']= np.mean(np.mean(librosa.feature.tonnetz(y=librosa.effects.harmonic(x),sr=sample_rate).T))\n\n        x=pd.Series(x)\n        \n        #some of mine some of andrews features\n\n        def calc_change_rate(x):\n            change = (np.diff(x) / x[:-1])\n            change = change[np.nonzero(change)[0]]\n            change = change[~np.isnan(change)]\n            change = change[change != -np.inf]\n            change = change[change != np.inf]\n            return np.mean(change)\n\n\n        def 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]\n\n        def classic_sta_lta(x, length_sta, length_lta):\n\n            sta = np.cumsum(x ** 2)\n\n            # Convert to float\n            sta = np.require(sta, dtype=np.float)\n\n            # Copy for LTA\n            lta = sta.copy()\n\n            # Compute the STA and the LTA\n            sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n            sta /= length_sta\n            lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n            lta /= length_lta\n\n            # Pad zeros\n            sta[:length_lta - 1] = 0\n\n            # Avoid division by zero by setting zero values to tiny float\n            dtiny = np.finfo(0.0).tiny\n            idx = lta < dtiny\n            lta[idx] = dtiny\n\n            return sta / lta\n        \n\n        feature_dict['mean'] = x.mean()\n        feature_dict['std'] = x.std()\n        feature_dict['max'] = x.max()\n        feature_dict['min'] = x.min()\n\n        feature_dict['mean_change_abs'] = np.mean(np.diff(x))\n        feature_dict['mean_change_rate'] = calc_change_rate(x)\n        feature_dict['abs_max'] = np.abs(x).max()\n        feature_dict['abs_min'] = np.abs(x).min()\n\n        feature_dict['std_first_50000'] = x[:50000].std()\n        feature_dict['std_last_50000'] = x[100000:].std()\n        feature_dict['std_first_10000'] = x[:10000].std()\n\n\n        feature_dict['avg_first_50000'] = x[:50000].mean()\n\n        feature_dict['avg_first_10000'] = x[:10000].mean()\n\n\n        feature_dict['min_first_50000'] = x[:50000].min()\n\n        feature_dict['min_first_10000'] = x[:10000].min()\n        feature_dict['min_last_10000'] = x[-10000:].min()\n\n        feature_dict['max_first_50000'] = x[:50000].max()\n\n        feature_dict['max_first_10000'] = x[:10000].max()\n        feature_dict['max_last_10000'] = x[-10000:].max()\n\n        feature_dict['max_to_min'] = x.max() / np.abs(x.min())\n        feature_dict['max_to_min_diff'] = x.max() - np.abs(x.min())\n        feature_dict['count_big'] = len(x[np.abs(x) > 500])\n        feature_dict['sum'] = x.sum()\n\n        feature_dict['mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n\n        feature_dict['mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n        feature_dict['mean_change_rate_last_10000'] = calc_change_rate(x[:-10000])\n\n        feature_dict['q95'] = np.quantile(x, 0.95)\n        feature_dict['q99'] = np.quantile(x, 0.99)\n        feature_dict['q05'] = np.quantile(x, 0.05)\n        feature_dict['q01'] = np.quantile(x, 0.01)\n\n        feature_dict['abs_q95'] = np.quantile(np.abs(x), 0.95)\n        feature_dict['abs_q99'] = np.quantile(np.abs(x), 0.99)\n        feature_dict['abs_q05'] = np.quantile(np.abs(x), 0.05)\n        feature_dict['abs_q01'] = np.quantile(np.abs(x), 0.01)\n\n        feature_dict['trend'] = add_trend_feature(x)\n        feature_dict['abs_trend'] = add_trend_feature(x, abs_values=True)\n        feature_dict['abs_mean'] = np.abs(x).mean()\n        feature_dict['abs_std'] = np.abs(x).std()\n\n        feature_dict['kurt'] = x.kurtosis()\n        feature_dict['skew'] = x.skew()\n        feature_dict['med'] = x.median()\n\n        feature_dict['Hilbert_mean'] = np.abs(hilbert(x)).mean()\n        feature_dict['Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n        feature_dict['classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n        feature_dict['classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n        feature_dict['classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n        feature_dict['classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n        feature_dict['classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n        feature_dict['classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n        feature_dict['classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n        feature_dict['classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n        feature_dict['Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n        ewma = pd.Series.ewm\n        feature_dict['exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n        feature_dict['exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n        feature_dict['exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n        no_of_std = 3\n        feature_dict['MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n        feature_dict['MA_700MA_BB_high_mean'] = (feature_dict['Moving_average_700_mean'] + no_of_std * feature_dict['Moving_average_700_mean'])\n        feature_dict['MA_700MA_BB_low_mean'] = (feature_dict['Moving_average_700_mean'] - no_of_std * feature_dict['Moving_average_700_mean'])\n        feature_dict['MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n        feature_dict['MA_400MA_BB_high_mean'] = (feature_dict['MA_400MA_std_mean'] + no_of_std * feature_dict['MA_400MA_std_mean'])\n        feature_dict['MA_400MA_BB_low_mean'] = (feature_dict['MA_400MA_std_mean'] - no_of_std * feature_dict['MA_400MA_std_mean'])\n        feature_dict['MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n\n\n        feature_dict['iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n        feature_dict['q999'] = np.quantile(x,0.999)\n        feature_dict['q001'] = np.quantile(x,0.001)\n        feature_dict['ave10'] = stats.trim_mean(x, 0.1)\n\n\n\n\n\n\n\n\n        \n        return feature_dict\n\n    def generate(self):\n        feature_list = []\n        res =Parallel(n_jobs=self.n_jobs,\n                       backend='threading')(delayed(self.features)(x, y, s)  for s, x, y in tqdm(self.read_chunks(), total=self.total_data)) # its a generator now, we can loop over it\n        for r in res:\n            feature_list.append(r)\n        return pd.DataFrame(feature_list)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets compare it now with a one without DWT denoising---https://www.kaggle.com/jackvial/dwt-signal-denoising"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T12:15:02.366956Z","start_time":"2019-04-19T12:15:02.310620Z"},"trusted":true},"cell_type":"code","source":"class FeatureGenerator1(object):\n    def __init__(self, dtype, n_jobs=1, chunk_size=None):\n        self.chunk_size = chunk_size\n        self.dtype = dtype\n        self.filename = None\n        self.n_jobs = n_jobs\n        self.test_files = []\n        if self.dtype == 'train':\n            self.filename = '../input/train.csv'\n            self.total_data = int(629145481 / self.chunk_size)\n        else:\n            submission = pd.read_csv('../input/sample_submission.csv')\n            for seg_id in submission.seg_id.values:\n                self.test_files.append((seg_id, '../input/test/' + seg_id + '.csv'))\n            self.total_data = int(len(submission))\n\n    def read_chunks(self):\n        if self.dtype == 'train':\n            iter_df = pd.read_csv(self.filename, iterator=True, chunksize=self.chunk_size,\n                                  dtype={'acoustic_data': np.float64, 'time_to_failure': np.float64})\n            for counter, df in enumerate(iter_df):\n                x = df.acoustic_data.values\n                y = df.time_to_failure.values[-1]\n                seg_id = 'train_' + str(counter)\n                yield seg_id, x, y # its a generator\n        else:\n            for seg_id, f in self.test_files:\n                df = pd.read_csv(f, dtype={'acoustic_data': np.float64})\n                x = df.acoustic_data.values\n                yield seg_id, x, -999\n\n    def features(self, x, y, seg_id):\n        \n        \n        feature_dict = dict()\n        feature_dict['target'] = y\n        feature_dict['seg_id'] = seg_id\n        #audio\n  \n        sample_rate=4000000\n        feature_dict['stft1']=np.abs(librosa.stft(x))\n        feature_dict['stft']= np.mean(np.abs(librosa.stft(x)))\n\n        feature_dict['mfccs']=np.mean(np.mean(librosa.feature.mfcc(y=x, sr=sample_rate, n_mfcc=40).T,axis=0))\n\n        feature_dict['chroma']=np.mean(np.mean((librosa.feature.chroma_stft(S=feature_dict['stft1'], sr=sample_rate).T)))\n\n        feature_dict['mel']=np.mean(np.mean(librosa.feature.melspectrogram(x, sr=sample_rate).T))\n\n        feature_dict['tonnetz']= np.mean(np.mean(librosa.feature.tonnetz(y=librosa.effects.harmonic(x),sr=sample_rate).T))\n\n        x=pd.Series(x)\n        #some of mine some of andrews features\n\n        def calc_change_rate(x):\n            change = (np.diff(x) / x[:-1])\n            change = change[np.nonzero(change)[0]]\n            change = change[~np.isnan(change)]\n            change = change[change != -np.inf]\n            change = change[change != np.inf]\n            return np.mean(change)\n\n\n        def 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]\n\n        def classic_sta_lta(x, length_sta, length_lta):\n\n            sta = np.cumsum(x ** 2)\n\n            # Convert to float\n            sta = np.require(sta, dtype=np.float)\n\n            # Copy for LTA\n            lta = sta.copy()\n\n            # Compute the STA and the LTA\n            sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n            sta /= length_sta\n            lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n            lta /= length_lta\n\n            # Pad zeros\n            sta[:length_lta - 1] = 0\n\n            # Avoid division by zero by setting zero values to tiny float\n            dtiny = np.finfo(0.0).tiny\n            idx = lta < dtiny\n            lta[idx] = dtiny\n\n            return sta / lta\n        \n\n        feature_dict['mean'] = x.mean()\n        feature_dict['std'] = x.std()\n        feature_dict['max'] = x.max()\n        feature_dict['min'] = x.min()\n\n        feature_dict['mean_change_abs'] = np.mean(np.diff(x))\n        feature_dict['mean_change_rate'] = calc_change_rate(x)\n        feature_dict['abs_max'] = np.abs(x).max()\n        feature_dict['abs_min'] = np.abs(x).min()\n\n        feature_dict['std_first_50000'] = x[:50000].std()\n        feature_dict['std_last_50000'] = x[100000:].std()\n        feature_dict['std_first_10000'] = x[:10000].std()\n\n\n        feature_dict['avg_first_50000'] = x[:50000].mean()\n\n        feature_dict['avg_first_10000'] = x[:10000].mean()\n\n\n        feature_dict['min_first_50000'] = x[:50000].min()\n\n        feature_dict['min_first_10000'] = x[:10000].min()\n        feature_dict['min_last_10000'] = x[-10000:].min()\n\n        feature_dict['max_first_50000'] = x[:50000].max()\n\n        feature_dict['max_first_10000'] = x[:10000].max()\n        feature_dict['max_last_10000'] = x[-10000:].max()\n\n        feature_dict['max_to_min'] = x.max() / np.abs(x.min())\n        feature_dict['max_to_min_diff'] = x.max() - np.abs(x.min())\n        feature_dict['count_big'] = len(x[np.abs(x) > 500])\n        feature_dict['sum'] = x.sum()\n\n        feature_dict['mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n\n        feature_dict['mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n        feature_dict['mean_change_rate_last_10000'] = calc_change_rate(x[:-10000])\n\n        feature_dict['q95'] = np.quantile(x, 0.95)\n        feature_dict['q99'] = np.quantile(x, 0.99)\n        feature_dict['q05'] = np.quantile(x, 0.05)\n        feature_dict['q01'] = np.quantile(x, 0.01)\n\n        feature_dict['abs_q95'] = np.quantile(np.abs(x), 0.95)\n        feature_dict['abs_q99'] = np.quantile(np.abs(x), 0.99)\n        feature_dict['abs_q05'] = np.quantile(np.abs(x), 0.05)\n        feature_dict['abs_q01'] = np.quantile(np.abs(x), 0.01)\n\n        feature_dict['trend'] = add_trend_feature(x)\n        feature_dict['abs_trend'] = add_trend_feature(x, abs_values=True)\n        feature_dict['abs_mean'] = np.abs(x).mean()\n        feature_dict['abs_std'] = np.abs(x).std()\n\n        feature_dict['kurt'] = x.kurtosis()\n        feature_dict['skew'] = x.skew()\n        feature_dict['med'] = x.median()\n\n        feature_dict['Hilbert_mean'] = np.abs(hilbert(x)).mean()\n        feature_dict['Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n        feature_dict['classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n        feature_dict['classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n        feature_dict['classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n        feature_dict['classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n        feature_dict['classic_sta_lta5_mean'] = classic_sta_lta(x, 50, 1000).mean()\n        feature_dict['classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n        feature_dict['classic_sta_lta7_mean'] = classic_sta_lta(x, 333, 666).mean()\n        feature_dict['classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n        feature_dict['Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n        ewma = pd.Series.ewm\n        feature_dict['exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n        feature_dict['exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n        feature_dict['exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n        no_of_std = 3\n        feature_dict['MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n        feature_dict['MA_700MA_BB_high_mean'] = (feature_dict['Moving_average_700_mean'] + no_of_std * feature_dict['Moving_average_700_mean'])\n        feature_dict['MA_700MA_BB_low_mean'] = (feature_dict['Moving_average_700_mean'] - no_of_std * feature_dict['Moving_average_700_mean'])\n        feature_dict['MA_400MA_std_mean'] = x.rolling(window=400).std().mean()\n        feature_dict['MA_400MA_BB_high_mean'] = (feature_dict['MA_400MA_std_mean'] + no_of_std * feature_dict['MA_400MA_std_mean'])\n        feature_dict['MA_400MA_BB_low_mean'] = (feature_dict['MA_400MA_std_mean'] - no_of_std * feature_dict['MA_400MA_std_mean'])\n        feature_dict['MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n\n\n        feature_dict['iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n        feature_dict['q999'] = np.quantile(x,0.999)\n        feature_dict['q001'] = np.quantile(x,0.001)\n        feature_dict['ave10'] = stats.trim_mean(x, 0.1)\n\n\n\n\n\n\n\n\n        \n        return feature_dict\n\n    def generate(self):\n        feature_list = []\n        res =Parallel(n_jobs=self.n_jobs,\n                       backend='threading')(delayed(self.features)(x, y, s)  for s, x, y in tqdm(self.read_chunks(), total=self.total_data)) # its a generator now, we can loop over it\n        for r in res:\n            feature_list.append(r)\n        return pd.DataFrame(feature_list)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T13:05:54.345061Z","start_time":"2019-04-19T12:15:06.358673Z"},"trusted":true},"cell_type":"code","source":"%time\n\ntraining_fg = FeatureGenerator(dtype='train', n_jobs=4, chunk_size=150000)\ntraining_data = training_fg.generate()\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train=training_data.target\ntraining_data.drop([\"stft1\",\"target\",\"seg_id\"],inplace=True,axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = MinMaxScaler()\ntraining_data = pd.DataFrame(scaler.fit_transform(training_data), columns=training_data.columns)\n","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T13:15:51.137804Z","start_time":"2019-04-19T13:15:51.129815Z"}},"cell_type":"markdown","source":"Dataset without DWT"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T14:04:23.460938Z","start_time":"2019-04-19T13:15:59.667984Z"},"trusted":true},"cell_type":"code","source":"%time\n\ntraining_fg_noDWT = FeatureGenerator1(dtype='train', n_jobs=4, chunk_size=150000)\ntraining_data_noDWT = training_fg_noDWT.generate()\ny_train_noDWT=training_data_noDWT.target\ntraining_data_noDWT.drop([\"stft1\",\"target\",\"seg_id\"],inplace=True,axis=1)\nscaler = MinMaxScaler()\ntraining_data_noDWT = pd.DataFrame(scaler.fit_transform(training_data_noDWT), columns=training_data_noDWT.columns)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Set data without DWT"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T14:58:02.977868Z","start_time":"2019-04-19T14:27:25.737893Z"},"trusted":true},"cell_type":"code","source":"%time\ntest_fg_noDWT = FeatureGenerator1(dtype='test', n_jobs=4, chunk_size=None)\ntest_data_noDWT = test_fg_noDWT.generate()\ntest_data_noDWT.drop([\"stft1\",\"target\",\"seg_id\"],inplace=True,axis=1)\nscaler = MinMaxScaler()\ntest_data_noDWT = pd.DataFrame(scaler.fit_transform(test_data_noDWT), columns=test_data_noDWT.columns)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"denoised test data:"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T15:55:38.547501Z","start_time":"2019-04-19T15:24:32.087523Z"},"trusted":true},"cell_type":"code","source":"%time\ntest_fg = FeatureGenerator(dtype='test', n_jobs=4, chunk_size=None)\ntest_data = test_fg.generate()\ntest_data.drop([\"stft1\",\"target\",\"seg_id\"],inplace=True,axis=1)\nscaler = MinMaxScaler()\ntest_data = pd.DataFrame(scaler.fit_transform(test_data), columns=test_data.columns)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:08:46.890775Z","start_time":"2019-04-19T16:08:46.598549Z"}},"cell_type":"markdown","source":"There is 1 missing value:"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:37:25.565893Z","start_time":"2019-04-19T16:37:25.550092Z"},"trusted":true},"cell_type":"code","source":"training_data.fillna(training_data.mean_change_rate_first_10000.mean(), inplace=True)\ntraining_data_noDWT.fillna(training_data_noDWT.mean_change_rate_first_10000.mean(), inplace=True)\n","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:37:26.542292Z","start_time":"2019-04-19T16:37:26.525752Z"},"trusted":true},"cell_type":"code","source":"test_data.fillna(test_data.mean_change_rate_first_10000.mean(), inplace=True)\ntest_data_noDWT.fillna(training_data_noDWT.mean_change_rate_first_10000.mean(), inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Vecstack, just playing around with it. No special hyperparameters or nothing"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:12:57.424495Z","start_time":"2019-04-19T16:12:57.405470Z"},"trusted":true},"cell_type":"code","source":"models = [\n    KNeighborsRegressor(n_neighbors=5,\n                        n_jobs=-1),\n        \n    RandomForestRegressor(random_state=0, n_jobs=-1, \n                           n_estimators=100, max_depth=3),\n        \n    XGBRegressor(random_state=0, n_jobs=-1, learning_rate=0.1, \n                  n_estimators=100, max_depth=3),\n    LinearRegression(),\n    CatBoostRegressor(random_state=0),\n    NuSVR(gamma='scale', nu=0.9, C=10.0, tol=0.01)\n    \n   \n]","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:13:02.828158Z","start_time":"2019-04-19T16:13:02.465259Z"},"trusted":true,"collapsed":true},"cell_type":"code","source":"training_data.std_last_50000[4194]=4\ntraining_data_noDWT.std_last_50000[4194]=4\n\n","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:27:53.112206Z","start_time":"2019-04-19T16:13:03.576080Z"},"trusted":true,"collapsed":true},"cell_type":"code","source":"S_train, S_test = stacking(models,                   \n                           training_data, y_train, test_data,   \n                           regression=True, \n     \n                           mode='oof_pred_bag', \n       \n                           needs_proba=False,\n         \n                           save_dir=None, \n            \n                           metric=mean_absolute_error, \n    \n                           n_folds=10, \n                 \n                           stratified=True,\n            \n                           shuffle=True,  \n            \n                           random_state=0,    \n         \n                           verbose=2)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:34:48.262000Z","start_time":"2019-04-19T16:34:48.255278Z"},"trusted":true},"cell_type":"code","source":"params = {'num_leaves': 128,\n          'min_data_in_leaf': 79,\n          'objective': 'huber',\n          'max_depth': -1,\n          'learning_rate': 0.01,\n          \"boosting\": \"gbdt\",\n          \"bagging_freq\": 5,\n          \"bagging_fraction\": 0.8126672064208567,\n          \"bagging_seed\": 11,\n          \"metric\": 'mae',\n          \"verbosity\": -1,\n          'reg_alpha': 0.1302650970728192,\n          'reg_lambda': 0.3603427518866501\n         }","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"CV with stacking, 3*5 CV\n\nNow I do realise that CV will be inflated because we used the same y_train data for the weak learners but one can use it as a benchmark for trying other things, not necessarily to correlate it to LB"},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:34:49.488066Z","start_time":"2019-04-19T16:34:49.476429Z"},"trusted":true},"cell_type":"code","source":"skf = KFold(n_splits=5, shuffle=True, random_state=123)\noof = pd.DataFrame(y_train)\noof['predict'] = 0\n\nval_mae = []","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:34:50.501349Z","start_time":"2019-04-19T16:34:50.494118Z"},"trusted":true},"cell_type":"code","source":"S_train=pd.DataFrame(S_train)\ny_train1=pd.Series(y_train)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:35:13.465809Z","start_time":"2019-04-19T16:34:52.597566Z"},"trusted":true,"collapsed":true},"cell_type":"code","source":"for fold, (trn_idx, val_idx) in enumerate(skf.split(S_train, y_train1)):\n    X_train, y_train = S_train.iloc[trn_idx], y_train1.iloc[trn_idx]\n    X_valid, y_valid = S_train.iloc[val_idx], y_train1.iloc[val_idx]\n    \n    N = 3\n    p_valid,yp = 0,0\n    for i in range(N):\n    \n        trn_data = lgb.Dataset(X_train, label=y_train)\n        val_data = lgb.Dataset(X_valid, label=y_valid)\n        evals_result = {}\n        lgb_clf = lgb.train(params,trn_data,1000,valid_sets = [trn_data, val_data],early_stopping_rounds=100,verbose_eval = 50,evals_result=evals_result)\n        p_valid += lgb_clf.predict(X_valid)\n\n    \n    oof['predict'][val_idx] = p_valid/N\n    mae = mean_absolute_error(y_valid, p_valid)\n    val_mae.append(mae)\n\n   \n    \n    ","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:35:17.305481Z","start_time":"2019-04-19T16:35:17.286009Z"},"trusted":true},"cell_type":"code","source":"mae1 = mean_absolute_error(oof['target'], oof['predict'])\nprint(\"local mae1 = {}\".format(mae1))\n      ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok lets try without denoising of data."},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:54:12.672567Z","start_time":"2019-04-19T16:41:05.381044Z"},"trusted":true,"collapsed":true},"cell_type":"code","source":"S_train_noDWT, S_test_noDWT = stacking(models,                   \n                           training_data_noDWT, y_train_noDWT, test_data_noDWT,   \n                           regression=True, \n     \n                           mode='oof_pred_bag', \n       \n                           needs_proba=False,\n         \n                           save_dir=None, \n            \n                           metric=mean_absolute_error, \n    \n                           n_folds=10, \n                 \n                           stratified=True,\n            \n                           shuffle=True,  \n            \n                           random_state=0,    \n         \n                           verbose=2)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:56:23.100186Z","start_time":"2019-04-19T16:56:23.080152Z"},"trusted":true},"cell_type":"code","source":"\noof_noDWT = pd.DataFrame(y_train_noDWT)\noof_noDWT['predict'] = 0\n\nval_mae_noDWT = []","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T16:57:00.951874Z","start_time":"2019-04-19T16:57:00.939041Z"},"trusted":true},"cell_type":"code","source":"S_train_noDWT=pd.DataFrame(S_train_noDWT)\ny_train_noDWT=pd.Series(y_train_noDWT)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T17:17:49.284463Z","start_time":"2019-04-19T17:17:27.163637Z"},"trusted":true,"collapsed":true},"cell_type":"code","source":"for fold, (trn_idx, val_idx) in enumerate(skf.split(S_train_noDWT, y_train_noDWT)):\n    X_train1, y_train1 = S_train_noDWT.iloc[trn_idx], y_train_noDWT.iloc[trn_idx]\n    X_valid1, y_valid1 = S_train_noDWT.iloc[val_idx], y_train_noDWT.iloc[val_idx]\n    \n    N = 3\n    p_valid,yp = 0,0\n    for i in range(N):\n    \n        trn_data1 = lgb.Dataset(X_train1, label=y_train1)\n        val_data1 = lgb.Dataset(X_valid1, label=y_valid1)\n        evals_result1 = {}\n        lgb_clf_noDWT = lgb.train(params,trn_data1,1000,valid_sets = [trn_data1, val_data1],early_stopping_rounds=100,verbose_eval = 50,evals_result=evals_result1)\n        p_valid += lgb_clf_noDWT.predict(X_valid1)\n\n    \n    oof_noDWT['predict'][val_idx] = p_valid/N\n    mae1 = mean_absolute_error(y_valid1, p_valid)\n    val_mae_noDWT.append(mae1)\n\n   \n    \n    ","execution_count":1,"outputs":[{"output_type":"error","ename":"NameError","evalue":"name 'skf' is not defined","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m<ipython-input-1-ac56ee2752c8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mfold\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtrn_idx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mval_idx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mskf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mS_train_noDWT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train_noDWT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m     \u001b[0mX_train\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mS_train_noDWT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrn_idx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train_noDWT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrn_idx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mX_valid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_valid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mS_train_noDWT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mval_idx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train_noDWT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mval_idx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mNameError\u001b[0m: name 'skf' is not defined"]}]},{"metadata":{"ExecuteTime":{"end_time":"2019-04-19T17:18:02.028702Z","start_time":"2019-04-19T17:18:02.013121Z"},"trusted":true},"cell_type":"code","source":"mae2 = mean_absolute_error(oof_noDWT['target'], oof_noDWT['predict'])\nprint(\"local mae1 = {}\".format(mae2))\n      ","execution_count":null,"outputs":[]}],"metadata":{"hide_input":false,"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"},"toc":{"base_numbering":1,"nav_menu":{},"number_sections":true,"sideBar":true,"skip_h1_title":false,"title_cell":"Table of Contents","title_sidebar":"Contents","toc_cell":false,"toc_position":{},"toc_section_display":true,"toc_window_display":false},"varInspector":{"cols":{"lenName":16,"lenType":16,"lenVar":40},"kernels_config":{"python":{"delete_cmd_postfix":"","delete_cmd_prefix":"del ","library":"var_list.py","varRefreshCmd":"print(var_dic_list())"},"r":{"delete_cmd_postfix":") ","delete_cmd_prefix":"rm(","library":"var_list.r","varRefreshCmd":"cat(var_dic_list()) "}},"types_to_exclude":["module","function","builtin_function_or_method","instance","_Feature"],"window_display":false}},"nbformat":4,"nbformat_minor":1}