{"cells":[{"metadata":{},"cell_type":"markdown","source":"# One Feature, No ML, Gold Medal Range"},{"metadata":{},"cell_type":"markdown","source":"#### This kernel demonstrates the ability to achieve a high score (private leaderboard: 2.33037; 9th place) in the LANL Earthquake Prediction competition with only _one_ well-crafted feature and information from the test set data leak. This kernel was inspired by my [Three Keys to This Competition](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/94355) discussion topic."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nfrom sklearn.model_selection import KFold\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import mean_absolute_error\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm_notebook\nimport pandas as pd\nfrom numba import njit\nimport gc","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Define a few functions"},{"metadata":{},"cell_type":"markdown","source":"`numba`tized version of information entropy. This is faster than the version in `scipy`."},{"metadata":{"trusted":true},"cell_type":"code","source":"@njit\ndef entropy_fast(vec, bins):\n    h = np.histogram(vec, bins=bins)[0] + 1E-15\n    h = h / np.sum(h)\n    return -np.sum(h * np.log(h))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Replaces peaks above the absolute value of a threshold with the average of nearby peaks."},{"metadata":{"trusted":true},"cell_type":"code","source":"def replace_peaks(X, win=1000, threshold=500):\n    np.random.seed(0)\n    \n    for idx in range(X.shape[0]):\n        max_idx = np.argmax(np.abs(X[idx]))\n        \n        while np.abs(X[idx,max_idx]) > threshold:\n            temp_win = win + np.random.randint(-500, 501)\n            X[idx,max_idx] = 0.5 * (X[idx, np.maximum(0, max_idx-temp_win)] + X[idx, np.minimum(X[idx].shape[0] - 1, max_idx+temp_win)])\n            max_idx = np.argmax(np.abs(X[idx]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndf = pd.read_csv(os.path.join('..','input','train.csv'), dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_samples = 150000\nn_segments = int(np.floor(df.shape[0] / n_samples))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = np.zeros((n_segments, n_samples))\ny = np.zeros(n_segments)\n\nfor seg_id in range(n_segments):\n    seg = df.iloc[seg_id*n_samples:seg_id*n_samples+n_samples]\n    X[seg_id] = seg['acoustic_data'].values\n    y[seg_id] = seg['time_to_failure'].values[-1]\n\ndel df\ngc.collect();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Replace peaks"},{"metadata":{"trusted":true},"cell_type":"code","source":"replace_peaks(X)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create information entropy feature"},{"metadata":{},"cell_type":"markdown","source":"Split the segment into _n_ parts. Calculate the entropy on each part and take the mean to reduce noise."},{"metadata":{"trusted":true},"cell_type":"code","source":"ent = np.zeros(X.shape[0])\n\nn=2000\nent_temp = np.zeros(n)\ncv = KFold(n, shuffle=False)\n\nfor idx in tqdm_notebook(range(X.shape[0])):\n    for idx2, (train_idx, test_idx) in enumerate(cv.split(X[idx])):\n        ent_temp[idx2] = entropy_fast(X[idx,test_idx], 300)\n    \n    ent[idx] = np.mean(ent_temp)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot the entropy."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,5))\nplt.plot(ent)\nplt.xlabel('time')\nplt.ylabel('entropy');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Do some clipping and linearly transform the feature to best match the TTF. The optimal values for clipping were determined through some iterative hand tuning (not shown here) to arrive at an MAE of 2.112. Some people may say that linear regression is machine learning. However, I would claim that with one variable, it's a simple univariate linear transformation. "},{"metadata":{"trusted":true},"cell_type":"code","source":"feature = np.clip(ent, a_min=2.37, a_max=2.66)\n\nmodel = LinearRegression(n_jobs=-1)\nmodel.fit(feature.reshape(-1,1), y.reshape(-1,1))\nfeature = model.predict(feature.reshape(-1,1))\n\nprint('MAE: ', mean_absolute_error(feature, y))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"View the feature and TTF together."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,5))\nplt.plot(feature,'r')\nplt.plot(y,'k')\nplt.ylabel('TTF')\nplt.xlabel('time')\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Load testing data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(os.path.join('..','input','sample_submission.csv'), index_col='seg_id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_test_segments = submission.shape[0]\n\nX_test = np.zeros((n_test_segments, n_samples))\n\nfor seg_idx, seg_id in enumerate(tqdm_notebook(submission.index)):\n    seg = pd.read_csv(os.path.join('..','input','test', seg_id + '.csv'))\n    X_test[seg_idx] = seg['acoustic_data'].values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Replace peaks and compute entropy feature for testing data."},{"metadata":{"trusted":true},"cell_type":"code","source":"replace_peaks(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ent = np.zeros(X_test.shape[0])\n\nn=2000\nent_temp = np.zeros(n)\ncv = KFold(n, shuffle=False)\n\nfor idx in tqdm_notebook(range(X_test.shape[0])):\n    for idx2, (train_idx, test_idx) in enumerate(cv.split(X_test[idx])):\n        ent_temp[idx2] = entropy_fast(X_test[idx,test_idx], 300)\n    \n    ent[idx] = np.mean(ent_temp)\n\ntest_feature = np.clip(ent, a_min=2.37, a_max=2.66)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Take advantage of test data set leak"},{"metadata":{},"cell_type":"markdown","source":"The peaks of the TTF in the test set were determined from the figures in [this discussion post](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/94086)."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_set_TTF_peaks = [11, 11, 11, 8, 11, 16, 9, 11, 16]\ntest_set_TTF_mean_peak = np.mean(test_set_TTF_peaks)\nscaling_factor = test_set_TTF_mean_peak / np.max(feature)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Multiply predictions by the ratio of the mean peak TTF values in the test set to the peak TTF value in the predictions from the one entropy feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = scaling_factor * model.predict(test_feature.reshape(-1,1)).flatten()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.time_to_failure = predictions\nsubmission.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.5.2"}},"nbformat":4,"nbformat_minor":1}