{"cells":[{"metadata":{},"cell_type":"markdown","source":"Importing libraries."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":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\nfrom itertools import product\n\nfrom tsfresh.feature_extraction import feature_calculators\nfrom joblib import Parallel, delayed","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preparing functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create a training file with simple derived features\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]\n\ndef 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\ndef calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\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)","execution_count":2,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Main class for feature generation\n\nFrom abhishek kernel"},{"metadata":{"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                del df\n                yield seg_id, x, y\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[-self.chunk_size:]\n                del df\n                yield seg_id, x, -999\n    \n    def get_features(self, x, y, seg_id):\n        \"\"\"\n        Gets three groups of features: from original data and from reald and imaginary parts of FFT.\n        \"\"\"\n        \n        x = pd.Series(x)\n    \n        zc = np.fft.fft(x)\n        realFFT = pd.Series(np.real(zc))\n        imagFFT = pd.Series(np.imag(zc))\n        \n        main_dict = self.features(x, y, seg_id)\n        r_dict = self.features(realFFT, y, seg_id)\n        i_dict = self.features(imagFFT, y, seg_id)\n        \n        for k, v in r_dict.items():\n            if k not in ['target', 'seg_id']:\n                main_dict[f'fftr_{k}'] = v\n                \n        for k, v in i_dict.items():\n            if k not in ['target', 'seg_id']:\n                main_dict[f'ffti_{k}'] = v\n        \n        return main_dict\n        \n    \n    def features(self, x, y, seg_id):\n        feature_dict = dict()\n        feature_dict['target'] = y\n        feature_dict['seg_id'] = seg_id\n\n        # create features here\n\n        # lists with parameters to iterate over them\n        percentiles = [1, 5, 10, 20, 25, 30, 40, 50, 60, 70, 75, 80, 90, 95, 99]\n        hann_windows = [50, 150, 1500, 15000]\n        spans = [300, 3000, 30000, 50000]\n        windows = [10, 50, 100, 500, 1000, 10000]\n        borders = list(range(-4000, 4001, 1000))\n        peaks = [10, 20, 50, 100]\n        coefs = [1, 5, 10, 50, 100]\n        lags = [10, 100, 1000, 10000]\n        autocorr_lags = [5, 10, 50, 100, 500, 1000, 5000, 10000]\n\n        # basic stats\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        # basic stats on absolute values\n        feature_dict['mean_change_abs'] = np.mean(np.diff(x))\n        feature_dict['abs_max'] = np.abs(x).max()\n        feature_dict['abs_mean'] = np.abs(x).mean()\n        feature_dict['abs_std'] = np.abs(x).std()\n\n        # geometric and harminic means\n        feature_dict['hmean'] = stats.hmean(np.abs(x[np.nonzero(x)[0]]))\n        feature_dict['gmean'] = stats.gmean(np.abs(x[np.nonzero(x)[0]])) \n\n        # k-statistic and moments\n        for i in range(1, 5):\n            feature_dict[f'kstat_{i}'] = stats.kstat(x, i)\n            feature_dict[f'moment_{i}'] = stats.moment(x, i)\n\n        for i in [1, 2]:\n            feature_dict[f'kstatvar_{i}'] = stats.kstatvar(x, i)\n\n        # aggregations on various slices of data\n        for agg_type, slice_length, direction in product(['std', 'min', 'max', 'mean'], [1000, 10000, 50000], ['first', 'last']):\n            if direction == 'first':\n                feature_dict[f'{agg_type}_{direction}_{slice_length}'] = x[:slice_length].agg(agg_type)\n            elif direction == 'last':\n                feature_dict[f'{agg_type}_{direction}_{slice_length}'] = x[-slice_length:].agg(agg_type)\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'] = calc_change_rate(x)\n        # calc_change_rate on slices of data\n        for slice_length, direction in product([1000, 10000, 50000], ['first', 'last']):\n            if direction == 'first':\n                feature_dict[f'mean_change_rate_{direction}_{slice_length}'] = calc_change_rate(x[:slice_length])\n            elif direction == 'last':\n                feature_dict[f'mean_change_rate_{direction}_{slice_length}'] = calc_change_rate(x[-slice_length:])\n\n        # percentiles on original and absolute values\n        for p in percentiles:\n            feature_dict[f'percentile_{p}'] = np.percentile(x, p)\n            feature_dict[f'abs_percentile_{p}'] = np.percentile(np.abs(x), p)\n\n        feature_dict['trend'] = add_trend_feature(x)\n        feature_dict['abs_trend'] = add_trend_feature(x, abs_values=True)\n\n        feature_dict['mad'] = x.mad()\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\n        for hw in hann_windows:\n            feature_dict[f'Hann_window_mean_{hw}'] = (convolve(x, hann(hw), mode='same') / sum(hann(hw))).mean()\n\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\n        # exponential rolling statistics\n        ewma = pd.Series.ewm\n        for s in spans:\n            feature_dict[f'exp_Moving_average_{s}_mean'] = (ewma(x, span=s).mean(skipna=True)).mean(skipna=True)\n            feature_dict[f'exp_Moving_average_{s}_std'] = (ewma(x, span=s).mean(skipna=True)).std(skipna=True)\n            feature_dict[f'exp_Moving_std_{s}_mean'] = (ewma(x, span=s).std(skipna=True)).mean(skipna=True)\n            feature_dict[f'exp_Moving_std_{s}_std'] = (ewma(x, span=s).std(skipna=True)).std(skipna=True)\n\n        feature_dict['iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n        feature_dict['iqr1'] = np.subtract(*np.percentile(x, [95, 5]))\n        feature_dict['ave10'] = stats.trim_mean(x, 0.1)\n        \n        for slice_length, threshold in product([50000, 100000, 150000],\n                                                     [5, 10, 20, 50, 100]):\n            feature_dict[f'count_big_{slice_length}_threshold_{threshold}'] = (np.abs(x[-slice_length:]) > threshold).sum()\n            feature_dict[f'count_big_{slice_length}_less_threshold_{threshold}'] = (np.abs(x[-slice_length:]) < threshold).sum()\n\n        # tfresh features take too long to calculate, so I comment them for now\n\n#         feature_dict['abs_energy'] = feature_calculators.abs_energy(x)\n#         feature_dict['abs_sum_of_changes'] = feature_calculators.absolute_sum_of_changes(x)\n#         feature_dict['count_above_mean'] = feature_calculators.count_above_mean(x)\n#         feature_dict['count_below_mean'] = feature_calculators.count_below_mean(x)\n#         feature_dict['mean_abs_change'] = feature_calculators.mean_abs_change(x)\n#         feature_dict['mean_change'] = feature_calculators.mean_change(x)\n#         feature_dict['var_larger_than_std_dev'] = feature_calculators.variance_larger_than_standard_deviation(x)\n        feature_dict['range_minf_m4000'] = feature_calculators.range_count(x, -np.inf, -4000)\n        feature_dict['range_p4000_pinf'] = feature_calculators.range_count(x, 4000, np.inf)\n\n        for i, j in zip(borders, borders[1:]):\n            feature_dict[f'range_{i}_{j}'] = feature_calculators.range_count(x, i, j)\n\n#         feature_dict['ratio_unique_values'] = feature_calculators.ratio_value_number_to_time_series_length(x)\n#         feature_dict['first_loc_min'] = feature_calculators.first_location_of_minimum(x)\n#         feature_dict['first_loc_max'] = feature_calculators.first_location_of_maximum(x)\n#         feature_dict['last_loc_min'] = feature_calculators.last_location_of_minimum(x)\n#         feature_dict['last_loc_max'] = feature_calculators.last_location_of_maximum(x)\n\n#         for lag in lags:\n#             feature_dict[f'time_rev_asym_stat_{lag}'] = feature_calculators.time_reversal_asymmetry_statistic(x, lag)\n        for autocorr_lag in autocorr_lags:\n            feature_dict[f'autocorrelation_{autocorr_lag}'] = feature_calculators.autocorrelation(x, autocorr_lag)\n            feature_dict[f'c3_{autocorr_lag}'] = feature_calculators.c3(x, autocorr_lag)\n\n#         for coeff, attr in product([1, 2, 3, 4, 5], ['real', 'imag', 'angle']):\n#             feature_dict[f'fft_{coeff}_{attr}'] = list(feature_calculators.fft_coefficient(x, [{'coeff': coeff, 'attr': attr}]))[0][1]\n\n#         feature_dict['long_strk_above_mean'] = feature_calculators.longest_strike_above_mean(x)\n#         feature_dict['long_strk_below_mean'] = feature_calculators.longest_strike_below_mean(x)\n#         feature_dict['cid_ce_0'] = feature_calculators.cid_ce(x, 0)\n#         feature_dict['cid_ce_1'] = feature_calculators.cid_ce(x, 1)\n\n        for p in percentiles:\n            feature_dict[f'binned_entropy_{p}'] = feature_calculators.binned_entropy(x, p)\n\n        feature_dict['num_crossing_0'] = feature_calculators.number_crossing_m(x, 0)\n\n        for peak in peaks:\n            feature_dict[f'num_peaks_{peak}'] = feature_calculators.number_peaks(x, peak)\n\n        for c in coefs:\n            feature_dict[f'spkt_welch_density_{c}'] = list(feature_calculators.spkt_welch_density(x, [{'coeff': c}]))[0][1]\n            feature_dict[f'time_rev_asym_stat_{c}'] = feature_calculators.time_reversal_asymmetry_statistic(x, c)  \n\n        # statistics on rolling windows of various sizes\n        for w in windows:\n            x_roll_std = x.rolling(w).std().dropna().values\n            x_roll_mean = x.rolling(w).mean().dropna().values\n\n            feature_dict[f'ave_roll_std_{w}'] = x_roll_std.mean()\n            feature_dict[f'std_roll_std_{w}'] = x_roll_std.std()\n            feature_dict[f'max_roll_std_{w}'] = x_roll_std.max()\n            feature_dict[f'min_roll_std_{w}'] = x_roll_std.min()\n\n            for p in percentiles:\n                feature_dict[f'percentile_roll_std_{p}_window_{w}'] = np.percentile(x_roll_std, p)\n\n            feature_dict[f'av_change_abs_roll_std_{w}'] = np.mean(np.diff(x_roll_std))\n            feature_dict[f'av_change_rate_roll_std_{w}'] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n            feature_dict[f'abs_max_roll_std_{w}'] = np.abs(x_roll_std).max()\n\n            feature_dict[f'ave_roll_mean_{w}'] = x_roll_mean.mean()\n            feature_dict[f'std_roll_mean_{w}'] = x_roll_mean.std()\n            feature_dict[f'max_roll_mean_{w}'] = x_roll_mean.max()\n            feature_dict[f'min_roll_mean_{w}'] = x_roll_mean.min()\n\n            for p in percentiles:\n                feature_dict[f'percentile_roll_mean_{p}_window_{w}'] = np.percentile(x_roll_mean, p)\n\n            feature_dict[f'av_change_abs_roll_mean_{w}'] = np.mean(np.diff(x_roll_mean))\n            feature_dict[f'av_change_rate_roll_mean_{w}'] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n            feature_dict[f'abs_max_roll_mean_{w}'] = np.abs(x_roll_mean).max()       \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.get_features)(x, y, s)\n                                            for s, x, y in tqdm_notebook(self.read_chunks(), total=self.total_data))\n        for r in res:\n            feature_list.append(r)\n        return pd.DataFrame(feature_list)","execution_count":3,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Generating features\n\nI calculate statistics on original values and on real and imaginary parts of FFT."},{"metadata":{"trusted":true},"cell_type":"code","source":"training_fg = FeatureGenerator(dtype='train', n_jobs=20, chunk_size=150000)\ntraining_data = training_fg.generate()\n\ntest_fg = FeatureGenerator(dtype='test', n_jobs=20, chunk_size=150000)\ntest_data = test_fg.generate()\n\nX = training_data.drop(['target', 'seg_id'], axis=1)\nX_test = test_data.drop(['target', 'seg_id'], axis=1)\ntest_segs = test_data.seg_id\ny = training_data.target","execution_count":null,"outputs":[{"output_type":"display_data","data":{"text/plain":"HBox(children=(IntProgress(value=0, max=4194), HTML(value='')))","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"463dd2a49ffd4b788b60fdec01b87198"}},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"Now we have much more features instead of 138 in my previous kernels.\nLets's see which of them are the most correlated with the target!"},{"metadata":{"trusted":true},"cell_type":"code","source":"X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.abs(X.corrwith(y)).sort_values(ascending=False).head(12)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Fixing missing values"},{"metadata":{"trusted":true},"cell_type":"code","source":"means_dict = {}\nfor col in X.columns:\n    if X[col].isnull().any():\n        print(col)\n        mean_value = X.loc[X[col] != -np.inf, col].mean()\n        X.loc[X[col] == -np.inf, col] = mean_value\n        X[col] = X[col].fillna(mean_value)\n        means_dict[col] = mean_value","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for col in X_test.columns:\n    if X_test[col].isnull().any():\n        X_test.loc[X_test[col] == -np.inf, col] = means_dict[col]\n        X_test[col] = X_test[col].fillna(means_dict[col])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X.to_csv('train_features.csv', index=False)\nX_test.to_csv('test_features.csv', index=False)\npd.DataFrame(y).to_csv('y.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}