{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this kernel, I tried to create a more accurate algorithm for detecting major and minor failures, which is based on cluster analysis and analysis of PCA components."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport gc\nimport time\nimport numpy as np\nimport pandas as pd\nfrom scipy import signal\nfrom scipy import stats\nfrom joblib import Parallel, delayed\nfrom tqdm import tqdm_notebook as tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tsfresh.feature_extraction import feature_calculators\nfrom scipy.signal import hilbert\nimport pywt \nfrom sklearn.cluster import DBSCAN\nfrom statsmodels.robust import mad\n%matplotlib inline","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## <center> Data Preparing"},{"metadata":{"trusted":true},"cell_type":"code","source":"def maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    \n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\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    # Calculate 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')\n","execution_count":2,"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 = 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        main_dict = self.features(x, y, seg_id)\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        x = pd.Series(denoise_signal(x, wavelet='db1', level=1))\n        #x = x - np.mean(x)\n    \n        zc = np.fft.fft(x)\n        zc = zc[:37500]\n\n        # FFT transform values\n        realFFT = np.real(zc)\n        imagFFT = np.imag(zc)\n\n        freq_bands = [x for x in range(0, 37500, 7500)]\n        magFFT = np.sqrt(realFFT ** 2 + imagFFT ** 2)\n        phzFFT = np.arctan(imagFFT / realFFT)\n        phzFFT[phzFFT == -np.inf] = -np.pi / 2.0\n        phzFFT[phzFFT == np.inf] = np.pi / 2.0\n        phzFFT = np.nan_to_num(phzFFT)\n\n        for freq in freq_bands:\n            if freq == 0:\n                continue\n            feature_dict['FFT_Mag_01q%d' % freq] = np.quantile(magFFT[freq: freq + 7500], 0.01)\n            feature_dict['FFT_Mag_10q%d' % freq] = np.quantile(magFFT[freq: freq + 7500], 0.1)\n            feature_dict['FFT_Mag_90q%d' % freq] = np.quantile(magFFT[freq: freq + 7500], 0.9)\n            feature_dict['FFT_Mag_99q%d' % freq] = np.quantile(magFFT[freq: freq + 7500], 0.99)\n            feature_dict['FFT_Mag_mean%d' % freq] = np.mean(magFFT[freq: freq + 7500])\n            feature_dict['FFT_Mag_std%d' % freq] = np.std(magFFT[freq: freq + 7500])\n            feature_dict['FFT_Mag_max%d' % freq] = np.max(magFFT[freq: freq + 7500])\n            \n        for p in [10]:\n            feature_dict[f'num_peaks_{p}'] = feature_calculators.number_peaks(x, 10)\n            \n        feature_dict['cid_ce'] = feature_calculators.cid_ce(x, normalize=True)\n            \n        for w in [5]:\n            feature_dict[f'autocorrelation_{w}'] = feature_calculators.autocorrelation(x, w)\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(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":{"trusted":true},"cell_type":"code","source":"training_fg = FeatureGenerator(dtype='train', n_jobs=-1, chunk_size=150000)\ntraining_data = training_fg.generate()\n\ntest_fg = FeatureGenerator(dtype='test', n_jobs=-1, 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":"3cf359507b6446e988d7fbb005e86df4"}},"metadata":{}},{"output_type":"stream","text":"/opt/conda/lib/python3.6/site-packages/pywt/_multilevel.py:148: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n  a = a[[slice(s) for s in d.shape]]\n","name":"stderr"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"del training_fg, training_data, test_fg, test_data; gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_dict = {}\nstd_dict = {}\nfor col in X.columns:\n    mean_value = X.loc[X[col] != -np.inf, col].mean()\n    std_value = X.loc[X[col] != -np.inf, col].std()\n    if X[col].isnull().any():\n        print(col)\n        X.loc[X[col] == -np.inf, col] = mean_value\n        X[col] = X[col].fillna(mean_value)\n    mean_dict[col] = mean_value\n    std_dict[col] = std_value\n    X[col] = (X[col] - mean_value)/std_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] = mean_dict[col]\n        X_test[col] = X_test[col].fillna(mean_dict[col])\n    X_test[col] = (X_test[col] - mean_dict[col])/std_dict[col]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## <center> Step by step"},{"metadata":{},"cell_type":"markdown","source":"Algorithm:\n* feature preparing\n* pca decomposition\n* split based on evristic major failures\n* split based on DBSCAN algorithm minor failures"},{"metadata":{},"cell_type":"markdown","source":"#### PCA Decomposition and Major failures splitting"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.decomposition import PCA\npca = PCA(n_components=15)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def major_failures(x):\n    return x > 45","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x_pca = pca.fit_transform(X.values)\nx_te_pca = pca.transform(X_test.values)\n\nmajor_failure = major_failures(x_pca[:, 0])\nmajor_failure_te = major_failures(x_te_pca[:, 0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 12))\n\nax1.plot(x_pca[major_failure, 0], x_pca[major_failure, 1], '.r', alpha=0.75, label='major failure')\nax1.plot(x_pca[~major_failure, 0], x_pca[~major_failure, 1], '.', color='g', alpha=0.75, label='other time domain')\nax1.axvline(45, color='b')\nax1.set_title(f'Train PCA (Major Failures {sum(major_failure)})')\nax1.legend()\n\nax2.plot(x_te_pca[major_failure_te, 0], x_te_pca[major_failure_te, 1], '.r', alpha=0.75, label='major failure')\nax2.plot(x_te_pca[~major_failure_te, 0], x_te_pca[~major_failure_te, 1], '.', color='g', alpha=0.75, label='other time domain')\nax2.axvline(45, color='b')\nax2.set_title(f'Test PCA (Major Failures {sum(major_failure_te)})')\nax2.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Minor failures splitting "},{"metadata":{"trusted":true},"cell_type":"code","source":"other_time_domain = ~major_failure\nother_te_time_domain = ~major_failure_te\n\nx_other = x_pca[other_time_domain]\nx_te_other = x_te_pca[other_te_time_domain]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dbscan = DBSCAN(eps=1.5)\nclust = dbscan.fit_predict(x_other[:, :3])\nclust_te = dbscan.fit_predict(x_te_other[:, :3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 12))\n\nax1.plot(x_other[clust==-1, 0], x_other[clust==-1, 1], '.', color='r', alpha=0.75, label='minor failure')\nax1.plot(x_other[clust!=-1, 0], x_other[clust!=-1, 1], '.', color='g', alpha=0.5, label='other time domain')\nax1.set_title(f'Train PCA (Minor Failures {sum(clust==-1)})')\nax1.legend()\n\nax2.plot(x_te_other[clust_te==-1, 0], x_te_other[clust_te==-1, 1], '.', color='r', alpha=0.75, label='minor failure')\nax2.plot(x_te_other[clust_te!=-1, 0], x_te_other[clust_te!=-1, 1], '.', color='g', alpha=0.5, label='other time domain')\nax2.set_title(f'Test PCA (Minor Failures {sum(clust_te==-1)})')\nax2.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Composed all below in one object"},{"metadata":{"trusted":true},"cell_type":"code","source":"class UFailureStateDetector:\n    def __init__(self, major_thr=45, minor_eps=1,n_components=5):\n        self.major_thr = major_thr\n        self.minor_eps = minor_eps\n        self.n_components = n_components \n        self.pca = PCA(n_components=0.99)\n        \n    def fit_predict(self, x):      \n        self.pca.fit(x)\n        return self.predict(x)\n    \n    def predict(self, x):\n        x = self.pca.transform(x)\n        \n        indexs = np.array([i for i in range(len(x))])\n        major = self.major_failures(x, self.major_thr)\n        major_indx = indexs[major]\n        \n        _x = x[~major]\n        _indexs = indexs[~major]\n        minor = self.minor_failures(_x, self.minor_eps, self.n_components)\n        minor_indx = _indexs[minor]\n        \n        other_indx = np.array([i for i in indexs if i not in major_indx and i not in minor_indx])\n        return other_indx, minor_indx, major_indx\n        \n    @staticmethod\n    def major_failures(x, thr):\n        return x[:, 0] > thr\n    \n    @staticmethod\n    def minor_failures(x, eps, n_components):\n        dbscan = DBSCAN(eps=eps)\n        clust = dbscan.fit_predict(x[:, :n_components])\n        minor = clust == -1\n        return minor","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ufsd = UFailureStateDetector(minor_eps=1.5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"other, minor, major = ufsd.fit_predict(X.values)\nother_te, minor_te, major_te = ufsd.predict(X_test.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 12))\n\nax1.plot(x_pca[major, 0], x_pca[major, 1], '.r', alpha=0.75, label='major failure')\nax1.plot(x_pca[minor, 0], x_pca[minor, 1], '.b', alpha=0.5, label='minor failure')\nax1.plot(x_pca[other, 0], x_pca[other, 1], '.', color='g', alpha=0.35, label='other time domain')\nax1.axvline(45, color='b')\nax1.set_title(f'Train PCA (Major Failures - {len(major)}, Minor Failures - {len(minor)})')\nax1.legend()\n\nax2.plot(x_te_pca[major_te, 0], x_te_pca[major_te, 1], '.r', alpha=0.75, label='major failure')\nax2.plot(x_te_pca[minor_te, 0], x_te_pca[minor_te, 1], '.b', alpha=0.5, label='minor failure')\nax2.plot(x_te_pca[other_te, 0], x_te_pca[other_te, 1], '.', color='g', alpha=0.35, label='other time domain')\nax2.axvline(45, color='b')\nax2.set_title(f'Test PCA (Major Failures - {len(major_te)}, Minor Failures - {len(minor_te)})')\nax2.legend()\nfig.savefig('pca_representation.png');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(24, 8))\nplt.plot(X.FFT_Mag_10q7500)\n\nfor m in major:\n    plt.axvspan(m-5, m+5, color='r', alpha=0.5)\n    \nfor m in minor:\n    plt.axvspan(m-1, m+1, color='g', alpha=0.5)\n\nfig.savefig('time_feature_representattion.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Minor in Train:', *minor, '\\n')\nprint('Major in Train:', *major, '\\n')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Minor in Test:', *test_segs.values[minor_te].tolist(), '\\n')\nprint('Major in Train:', *test_segs.values[major_te].tolist(), '\\n')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## <center> Conclusion\nThe current approach is more accurate and more powerful than the [previous one](https://www.kaggle.com/miklgr500/fast-failure-detector). It allows to detect not only major, but also minor failures. This algorithm was also able to detect two moments of failure, not detected in my previous core."},{"metadata":{},"cell_type":"markdown","source":"### <center> Reference\n*  https://www.kaggle.com/vettejeep/masters-final-project-model-lb-1-392\n*  https://www.kaggle.com/tarunpaparaju/lanl-earthquake-prediction-signal-denoising\n*  https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/80250#latest-532497\n*  https://www.kaggle.com/artgor/even-more-features"}],"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}