{"cells":[{"metadata":{"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\nfrom itertools import product\n\nfrom tsfresh.feature_extraction import feature_calculators\nfrom joblib import Parallel, delayed","execution_count":null,"outputs":[]},{"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":null,"outputs":[]},{"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 = np.sign(df.acoustic_data.values)*np.log1p(np.abs(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                x = np.sign(x)*np.log1p(np.abs(x))\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 in ['classic_sta_lta2_mean',\n                     'exp_Moving_std_3000_mean',\n                     'classic_sta_lta1_mean']:\n                main_dict[f'fftr_{k}'] = v\n                \n        for k, v in i_dict.items():\n            if k in ['classic_sta_lta2_mean']:\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        for p in [50]:\n            feature_dict[f'abs_percentile_{p}'] = np.percentile(np.abs(x), p)\n        for autocorr_lag in [5]:\n            feature_dict[f'autocorrelation_{autocorr_lag}'] = feature_calculators.autocorrelation(x, autocorr_lag)\n        for p in [95,99]:\n            feature_dict[f'binned_entropy_{p}'] = feature_calculators.binned_entropy(x, p)\n        \n        # calc_change_rate on slices of data\n        for slice_length, direction in product([50000], ['first']):\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        for peak in [10]:\n            feature_dict[f'num_peaks_{peak}'] = feature_calculators.number_peaks(x, peak)\n        \n        for p in [20]:\n            feature_dict[f'percentile_{p}'] = np.percentile(x, p)\n        \n        x_roll_std = x.rolling(1000).std().dropna().values\n        feature_dict[f'percentile_roll_std_30_window_1000'] = np.percentile(x_roll_std, 30)   \n        \n        x_roll_std = x.rolling(500).std().dropna().values\n        feature_dict[f'percentile_roll_std_75_window_500'] = np.percentile(x_roll_std, 75)   \n        feature_dict[f'percentile_roll_std_80_window_500'] = np.percentile(x_roll_std, 80)    \n        \n        ewma = pd.Series.ewm\n        feature_dict[f'exp_Moving_std_3000_mean'] = (ewma(x, span=3000).std(skipna=True)).mean(skipna=True)\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['skew'] = x.skew()    \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":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_fg = FeatureGenerator(dtype='train', n_jobs=4, chunk_size=150000)\ntraining_data = training_fg.generate()\ntraining_data.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"testing_fg = FeatureGenerator(dtype='test', n_jobs=4, chunk_size=150000)\ntest_data = testing_fg.generate()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = 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":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nss = StandardScaler()\nfeats = X_test.columns\nss.fit(pd.concat([X[feats],X_test[feats]]))\nX[feats] = ss.transform(X[feats])\nX_test[feats] = ss.transform(X_test[feats])\nX.insert(0,'seg_id',training_data.seg_id.values)\nX['time_to_failure'] = y.values\nX_test.insert(0,'seg_id',test_segs.values)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def GPI(data):\n    return (5.683668 +\n            1.0*np.tanh(((((data[\"abs_percentile_50\"]) - (((data[\"percentile_roll_std_75_window_500\"]) * 2.0)))) - (((((((((((data[\"percentile_roll_std_75_window_500\"]) + (data[\"num_peaks_10\"]))) * 2.0)) + (data[\"abs_percentile_50\"]))) * 2.0)) * 2.0)))) +\n            1.0*np.tanh(((-1.0) - (((((data[\"percentile_roll_std_30_window_1000\"]) + ((((data[\"binned_entropy_95\"]) + (data[\"num_peaks_10\"]))/2.0)))) + (data[\"percentile_roll_std_75_window_500\"]))))) +\n            1.0*np.tanh(((((((((data[\"abs_percentile_50\"]) + (data[\"binned_entropy_99\"]))) * ((((data[\"percentile_roll_std_75_window_500\"]) + (data[\"num_peaks_10\"]))/2.0)))) * ((-1.0*((((data[\"percentile_roll_std_75_window_500\"]) + (data[\"percentile_roll_std_30_window_1000\"])))))))) * 2.0)) +\n            1.0*np.tanh((((-1.0*((((((data[\"fftr_exp_Moving_std_3000_mean\"]) * (((1.0) + (((data[\"percentile_roll_std_30_window_1000\"]) * 2.0)))))) / 2.0))))) * ((((data[\"num_peaks_10\"]) + (data[\"abs_percentile_50\"]))/2.0)))) +\n            1.0*np.tanh(((((data[\"exp_Moving_std_3000_mean\"]) - (data[\"percentile_roll_std_75_window_500\"]))) + (((data[\"num_peaks_10\"]) * (((data[\"fftr_classic_sta_lta1_mean\"]) - (data[\"fftr_exp_Moving_std_3000_mean\"]))))))))\n\ndef GPII(data):\n    return (5.683668 +\n            1.0*np.tanh((((((((((-1.0*((((data[\"num_peaks_10\"]) + (data[\"percentile_roll_std_75_window_500\"])))))) * 2.0)) - (((data[\"num_peaks_10\"]) + (data[\"abs_percentile_50\"]))))) * 2.0)) * 2.0)) +\n            1.0*np.tanh((((((((data[\"exp_Moving_std_3000_mean\"]) + (-2.0))) + (((((((data[\"percentile_roll_std_80_window_500\"]) + (data[\"autocorrelation_5\"]))) / 2.0)) / 2.0)))/2.0)) - (((((data[\"percentile_roll_std_75_window_500\"]) * 2.0)) * 2.0)))) +\n            1.0*np.tanh((((-1.0*((data[\"fftr_exp_Moving_std_3000_mean\"])))) * (((((data[\"num_peaks_10\"]) + (data[\"abs_percentile_50\"]))) * (((((data[\"percentile_roll_std_30_window_1000\"]) * 2.0)) + (((1.0) + (data[\"percentile_roll_std_30_window_1000\"]))))))))) +\n            1.0*np.tanh(((data[\"exp_Moving_std_3000_mean\"]) - (((data[\"percentile_roll_std_75_window_500\"]) + ((((((((data[\"abs_percentile_50\"]) + (data[\"percentile_roll_std_75_window_500\"]))/2.0)) * ((((data[\"abs_percentile_50\"]) + (data[\"percentile_roll_std_75_window_500\"]))/2.0)))) * (data[\"fftr_exp_Moving_std_3000_mean\"]))))))) +\n            1.0*np.tanh((((((((((data[\"percentile_roll_std_80_window_500\"]) - (data[\"fftr_exp_Moving_std_3000_mean\"]))) + (((data[\"fftr_classic_sta_lta1_mean\"]) - (data[\"fftr_classic_sta_lta2_mean\"]))))/2.0)) + (((data[\"fftr_classic_sta_lta1_mean\"]) - (data[\"fftr_exp_Moving_std_3000_mean\"]))))) * (data[\"num_peaks_10\"]))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(feats)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"mean_absolute_error(X.time_to_failure,GPI(X))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_absolute_error(X.time_to_failure,GPII(X))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test['time_to_failure'] = (GPI(X_test)+GPII(X_test))/2.\nX_test[['seg_id','time_to_failure']].to_csv('submission.csv',index=False)","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.7.1"}},"nbformat":4,"nbformat_minor":1}