{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":11000,"databundleVersionId":875412,"sourceType":"competition"}],"dockerImageVersionId":30527,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<h1><center><font size=\"6\">LANL Earthquake EDA and Prediction</font></center></h1>\n\n<h2><center><font size=\"4\">Dataset used: LANL Earthquake Prediction</font></center></h2>\n\n<br>\n\n# <a id='0'>Content</a>\n\n- <a href='#1'>Introduction</a>  \n- <a href='#2'>Prepare the data analysis</a>  \n- <a href='#3'>Data exploration</a>   \n- <a href='#4'>Feature engineering</a>\n- <a href='#5'>Model</a>\n- <a href='#6'>Submission</a>  \n- <a href='#7'>References</a>","metadata":{"_uuid":"3b7efa15014e822ca40ee55d2e7e712a6ca63563"}},{"cell_type":"markdown","source":"# <a id='1'>Introduction</a>  \n\nIn Earth science, predicting earthquakes is equaly important and challenging. The current research related to forecasting earthquakes focuses on three aspects:\n- when such an event will occur\n- where this will happen (position, depth)\n- how large it will be (in terms of energy)\n\nThe current results of research are still far from allowing scientists to predict with accuracy even only one of these dimmensions.\n\n## Simulated earthquake experiment\n\nThe data are from an experiment conducted on rock in a double direct shear geometry subjected to bi-axial loading, a classic laboratory earthquake model. Two fault gouge layers are sheared simultaneously while subjected to a constant normal load and a prescribed shear velocity. The laboratory faults fail in repetitive cycles of stick and slip that is meant to mimic the cycle of loading and failure on tectonic faults. While the experiment is considerably simpler than a fault in Earth, it shares many physical characteristics. Los Alamos' initial work showed that the prediction of laboratory earthquakes from continuous seismic data is possible in the case of quasi-periodic laboratory seismic cycles.   \n\n## Competition \n\nIn this competition, the objective is to predict when an earthquake will take place. More precisely, competitors are asked to predict the time remaining before laboratory earthquake occur from real-time seismic data.\n\n## Data  \nThe training data is a single, continous segment of experiment data. The test data consist of a folder containing many small segments. The data within each test file is continous, but the test file do not represent a continous segment of the experiment. For each data segment in test data (identified by `seg_id`), the competitors should predict a single `time_to_failure` - this correspond to the time between the last row of the segment and the next laboratory earthquake.\n\n### Data files:\n* train.csv - single, continous training segment of experimental data.\n* test - a folder with many segments of test data.\n* sample_submission.csv - a sample submission ile in the correct format.\n\n## Kernel\n\nThis solution uses  Andrew's Data Munging Notebook plus some inspiration from Scirpus's [Notebook](https://www.kaggle.com/scirpus/andrews-script-plus-a-genetic-program-model/).\n\n","metadata":{"_uuid":"0bf2c2c00334bdc9f8f3d78d10d958bc521b4c7a"}},{"cell_type":"markdown","source":"# <a id='1'>Prepare the data analysis</a>\n\n## Load packages\n\nHere we define the packages for data manipulation, feature engineering and model training.","metadata":{"_uuid":"3b198f4c2141e0e94e0853cea5798c0e83343716"}},{"cell_type":"code","source":"import gc\nimport os\nimport time\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport xgboost as xgb\nimport lightgbm as lgb\nfrom scipy import stats\nfrom scipy.signal import hann\nfrom tqdm import tqdm_notebook\nimport matplotlib.pyplot as plt\nfrom scipy.signal import hilbert\nfrom scipy.signal import convolve\nfrom sklearn.svm import NuSVR, SVR\nfrom catboost import CatBoostRegressor\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nwarnings.filterwarnings(\"ignore\")","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2023-07-31T17:29:49.266944Z","iopub.execute_input":"2023-07-31T17:29:49.267572Z","iopub.status.idle":"2023-07-31T17:29:52.687726Z","shell.execute_reply.started":"2023-07-31T17:29:49.267538Z","shell.execute_reply":"2023-07-31T17:29:52.686820Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load the data\n\nLet's see first what files we have in input directory.","metadata":{"_uuid":"f7dc473a43cbd799a184ca71a27391bb6faa9984"}},{"cell_type":"code","source":"IS_LOCAL = False\nif(IS_LOCAL):\n    PATH=\"../input/LANL/\"\nelse:\n    PATH=\"../input/LANL-Earthquake-Prediction\"\nos.listdir(PATH)","metadata":{"_kg_hide-input":true,"_uuid":"45e5498b89764854036baa88eda6ead21d41f382","execution":{"iopub.status.busy":"2023-07-31T17:29:52.691473Z","iopub.execute_input":"2023-07-31T17:29:52.691811Z","iopub.status.idle":"2023-07-31T17:29:52.703869Z","shell.execute_reply.started":"2023-07-31T17:29:52.691782Z","shell.execute_reply":"2023-07-31T17:29:52.702623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have two files in the **input** directory and another directory, with the **test** data.  \n\nLet's see how many files are in **test** folder.","metadata":{"_uuid":"95ef0a76216010678fb760844307b24719eb4bec"}},{"cell_type":"code","source":"print(\"There are {} files in test folder\".format(len(os.listdir(os.path.join(PATH, 'test' )))))","metadata":{"_kg_hide-input":true,"_uuid":"91bda65accbd9172a37799acd2e9198a3477060e","execution":{"iopub.status.busy":"2023-07-31T17:29:52.705435Z","iopub.execute_input":"2023-07-31T17:29:52.705937Z","iopub.status.idle":"2023-07-31T17:29:53.132451Z","shell.execute_reply.started":"2023-07-31T17:29:52.705898Z","shell.execute_reply":"2023-07-31T17:29:53.131375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n\nLet's load the train file.\nTo reduce the size of data, we are loading the dataset in csv format with reduced size for each column.","metadata":{"_uuid":"68ced52082fc77c768db9ff14fed7cbe43e097d1"}},{"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(os.path.join(PATH,'train.csv'), dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","metadata":{"_kg_hide-input":true,"_uuid":"bb8916c762c26e98e31a721e4aae673aa12e6c5f","execution":{"iopub.status.busy":"2023-07-31T17:29:53.135093Z","iopub.execute_input":"2023-07-31T17:29:53.137304Z","iopub.status.idle":"2023-07-31T17:33:37.074283Z","shell.execute_reply.started":"2023-07-31T17:29:53.137253Z","shell.execute_reply":"2023-07-31T17:33:37.072431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check the data imported.","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a"}},{"cell_type":"code","source":"print(\"Train: rows:{} cols:{}\".format(train_df.shape[0], train_df.shape[1]))","metadata":{"_kg_hide-input":true,"_uuid":"40c15befb5782fe1d671ae4ad0ca592aff316073","execution":{"iopub.status.busy":"2023-07-31T17:33:37.076654Z","iopub.execute_input":"2023-07-31T17:33:37.077818Z","iopub.status.idle":"2023-07-31T17:33:37.085185Z","shell.execute_reply.started":"2023-07-31T17:33:37.077771Z","shell.execute_reply":"2023-07-31T17:33:37.084052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.options.display.precision = 15\ntrain_df.head(10)","metadata":{"_kg_hide-input":true,"_uuid":"0ee349f0d7f482a0cc4827a3b7515a915e4dd37f","execution":{"iopub.status.busy":"2023-07-31T17:33:37.086822Z","iopub.execute_input":"2023-07-31T17:33:37.087320Z","iopub.status.idle":"2023-07-31T17:33:37.118577Z","shell.execute_reply.started":"2023-07-31T17:33:37.087281Z","shell.execute_reply":"2023-07-31T17:33:37.117321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <a id='3'>Data exploration</a>  \n\nThe dimmension of the data is quite large, in excess of 600 millions rows of data.  \nThe two columns in the train dataset have the following meaning:   \n*  accoustic_data: is the accoustic signal measured in the laboratory experiment;  \n* time to failure: this gives the time until a failure will occurs.\n\nLet's plot 1% of the data. For this we will sample every 100 points of data.  ","metadata":{"_uuid":"877ce7e2c62dc0c7b9a86c081754a7f3b19bb364"}},{"cell_type":"code","source":"def plot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% sampled data\"):\n    \"\"\"\n    Plot accoustic data and time to failure\n    Args:\n        train_ad_sample_df: train accoustic data sample\n        train_ttf_sample_df: train time to failure data sample\n        title: title of the plot\n    Returns:\n        None\n    \"\"\"\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(title)\n    plt.plot(idx, train_ad_sample_df, color='r')\n    ax1.set_ylabel('acoustic data', color='r')\n    plt.legend(['acoustic data'], loc=(0.01, 0.95))\n    ax2 = ax1.twinx()\n    plt.plot(idx, train_ttf_sample_df, color='b')\n    ax2.set_ylabel('time to failure', color='b')\n    plt.legend(['time to failure'], loc=(0.01, 0.9))\n    plt.grid(True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-07-31T17:33:37.119765Z","iopub.execute_input":"2023-07-31T17:33:37.120080Z","iopub.status.idle":"2023-07-31T17:33:37.128238Z","shell.execute_reply.started":"2023-07-31T17:33:37.120038Z","shell.execute_reply":"2023-07-31T17:33:37.127097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[::100]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[::100]\nidx = train_df.index[::100]\nplot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df)\ndel idx\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"_kg_hide-input":true,"_uuid":"206af56767581d8e98777da3b0670e11d4e444bd","scrolled":true,"execution":{"iopub.status.busy":"2023-07-31T17:33:37.129549Z","iopub.execute_input":"2023-07-31T17:33:37.129930Z","iopub.status.idle":"2023-07-31T17:33:40.095665Z","shell.execute_reply.started":"2023-07-31T17:33:37.129899Z","shell.execute_reply":"2023-07-31T17:33:40.094756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*Discussion*\n\nThe plot shows only 1% of the full data. \nThe acoustic data shows complex oscilations with variable amplitude. Just before each failure there is an increase in the amplitude of the acoustic data. We see that large amplitudes are also obtained at different moments in time (for example about the mid-time between two succesive failures).  \n\nLet's plot as well the first 1% of the data.","metadata":{"_uuid":"36c363a5f2e56e54c8ffddd802aa86d91c61d97f"}},{"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[:6291455]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[:6291455]\nidx = train_df.index[:6291455]\nplot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: first 1% of data\")\ndel idx\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"_kg_hide-input":true,"_uuid":"845e5a241e545d73cfb1dbaea2cfa464faf8148a","execution":{"iopub.status.busy":"2023-07-31T17:33:40.096979Z","iopub.execute_input":"2023-07-31T17:33:40.097509Z","iopub.status.idle":"2023-07-31T17:33:42.649701Z","shell.execute_reply.started":"2023-07-31T17:33:40.097478Z","shell.execute_reply":"2023-07-31T17:33:42.647898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"On this zoomed-in-time plot we can see that actually the large oscilation before the failure is not quite in the last moment. There are also trains of intense oscilations preceeding the large one and also some oscilations with smaller peaks after the large one. Then, after some minor oscilations, the failure occurs.","metadata":{"_uuid":"78948b35358020156422a76c8d60dc132f17417b"}},{"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[6291455:2*6291455]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[6291455:2*6291455]\nidx = train_df.index[6291455:2*6291455]\nplot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 2nd 1% of data\")\ndel idx\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-07-31T17:33:42.653847Z","iopub.execute_input":"2023-07-31T17:33:42.654948Z","iopub.status.idle":"2023-07-31T17:33:45.448344Z","shell.execute_reply.started":"2023-07-31T17:33:42.654904Z","shell.execute_reply":"2023-07-31T17:33:45.447247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can observe that for this time interval, there was not a complete failure. We can see there are multiple small peaks during this time interval.","metadata":{}},{"cell_type":"markdown","source":"Let's select a larger time interval, to capture 2 failures. We will also sample every 10 values.","metadata":{}},{"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[:10*6291455][::10]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[:10*6291455][::10]\nidx = train_df.index[:10*6291455][::10]\nplot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: first 10% of data\")\ndel idx\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-07-31T17:33:45.449978Z","iopub.execute_input":"2023-07-31T17:33:45.450767Z","iopub.status.idle":"2023-07-31T17:33:48.066889Z","shell.execute_reply.started":"2023-07-31T17:33:45.450727Z","shell.execute_reply":"2023-07-31T17:33:48.065852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can observe the large peaks before the failure. Due to the 1/10 sampling rate, we didn't catched the maximum value of peaks. We also can see that there are multiple acoustic data peaks before the actual large peak that preceeds the failure.\n\nLet's also look more close to the last part of the train data.","metadata":{}},{"cell_type":"code","source":"train_ad_sample_df = train_df['acoustic_data'].values[-5*6291455:]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[-5*6291455:]\nidx = train_df.index[-5*6291455:]\nplot_acc_ttf_data(idx, train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: last 5% of data\")\ndel idx\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-07-31T17:33:48.068264Z","iopub.execute_input":"2023-07-31T17:33:48.068650Z","iopub.status.idle":"2023-07-31T17:33:58.725992Z","shell.execute_reply.started":"2023-07-31T17:33:48.068617Z","shell.execute_reply":"2023-07-31T17:33:58.724872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot also the acoustic signals from few test segments.","metadata":{}},{"cell_type":"code","source":"def plot_acc_data(test_sample_df, segment_name):\n    \"\"\"\n    Plot accoustic data for a train segment\n    Args:\n        test_sample_df: test accoustic data sample\n        segment_name: title of the plot\n    Returns:\n        None\n    \"\"\"\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(f\"Test segment: {segment_name}\")\n    plt.plot(test_sample_df, color='r')\n    ax1.set_ylabel('acoustic data', color='r')\n    plt.legend([f\"acoustic data: {segment_name}\"], loc=(0.01, 0.95))\n    plt.grid(True)","metadata":{"execution":{"iopub.status.busy":"2023-07-31T17:39:38.734529Z","iopub.execute_input":"2023-07-31T17:39:38.734942Z","iopub.status.idle":"2023-07-31T17:39:38.741458Z","shell.execute_reply.started":"2023-07-31T17:39:38.734911Z","shell.execute_reply":"2023-07-31T17:39:38.740277Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_sample_name = \"seg_00030f\"\ntest_sample_path = \"/kaggle/input/LANL-Earthquake-Prediction/test/\" + test_sample_name + \".csv\"\ntest_sample_df = pd.read_csv(test_sample_path)\nplot_acc_data(test_sample_df, test_sample_name)\ndel test_sample_df","metadata":{"execution":{"iopub.status.busy":"2023-07-31T17:40:08.377981Z","iopub.execute_input":"2023-07-31T17:40:08.378400Z","iopub.status.idle":"2023-07-31T17:40:09.028088Z","shell.execute_reply.started":"2023-07-31T17:40:08.378366Z","shell.execute_reply":"2023-07-31T17:40:09.026910Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_sample_name = \"seg_0012b5\"\ntest_sample_path = \"/kaggle/input/LANL-Earthquake-Prediction/test/\" + test_sample_name + \".csv\"\ntest_sample_df = pd.read_csv(test_sample_path)\nplot_acc_data(test_sample_df, test_sample_name)\ndel test_sample_df","metadata":{"execution":{"iopub.status.busy":"2023-07-31T17:40:33.490557Z","iopub.execute_input":"2023-07-31T17:40:33.490948Z","iopub.status.idle":"2023-07-31T17:40:34.055691Z","shell.execute_reply.started":"2023-07-31T17:40:33.490917Z","shell.execute_reply":"2023-07-31T17:40:34.054661Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_sample_name = \"seg_0042cc\"\ntest_sample_path = \"/kaggle/input/LANL-Earthquake-Prediction/test/\" + test_sample_name + \".csv\"\ntest_sample_df = pd.read_csv(test_sample_path)\nplot_acc_data(test_sample_df, test_sample_name)\ndel test_sample_df","metadata":{"execution":{"iopub.status.busy":"2023-07-31T17:41:00.071644Z","iopub.execute_input":"2023-07-31T17:41:00.072290Z","iopub.status.idle":"2023-07-31T17:41:00.622977Z","shell.execute_reply.started":"2023-07-31T17:41:00.072258Z","shell.execute_reply":"2023-07-31T17:41:00.621810Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_sample_name = \"seg_02042f\"\ntest_sample_path = \"/kaggle/input/LANL-Earthquake-Prediction/test/\" + test_sample_name + \".csv\"\ntest_sample_df = pd.read_csv(test_sample_path)\nplot_acc_data(test_sample_df, test_sample_name)\ndel test_sample_df","metadata":{"execution":{"iopub.status.busy":"2023-07-31T17:41:40.365756Z","iopub.execute_input":"2023-07-31T17:41:40.366168Z","iopub.status.idle":"2023-07-31T17:41:40.953039Z","shell.execute_reply.started":"2023-07-31T17:41:40.366140Z","shell.execute_reply":"2023-07-31T17:41:40.952080Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <a id='4'>Features engineering</a>  \n\nThe test segments are 150,000 data points each.   \nWe split the train data in segments of the same dimmension with the test sets.\n\nWe will create additional aggregation features, calculated on the segments. \n","metadata":{"_uuid":"8dee473309c39d74020b37f2e4863094191c121f"}},{"cell_type":"code","source":"rows = 150000\nsegments = int(np.floor(train_df.shape[0] / rows))\nprint(\"Number of segments: \", segments)","metadata":{"_kg_hide-input":true,"_uuid":"abcf886f07660f12d170a2182675e15efffefde0","execution":{"iopub.status.busy":"2023-07-30T12:40:58.540018Z","iopub.execute_input":"2023-07-30T12:40:58.540458Z","iopub.status.idle":"2023-07-30T12:40:58.547764Z","shell.execute_reply.started":"2023-07-30T12:40:58.540422Z","shell.execute_reply":"2023-07-30T12:40:58.546524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's define some computation helper functions.\n\nWe define a function to detect/extract trend feature and one to calculate classic STA/LTA factor.\n\nSTa/LTA is the ratio between the amplitude (calculated as average of absolute values) of the signal on a short time window of length STA and on a long time window of length LTA.","metadata":{"_uuid":"d107cd036c9a931d1e85795f11e5853811a44246"}},{"cell_type":"code","source":"def add_trend_feature(arr, abs_values=False):\n    \"\"\"\n    Calculate trend features\n    Uses a linear regression algorithm to extract the trend\n    from the list of values in the array (arr)\n    Args:\n        arr: array of values\n        abs_values: flag if to use abs values, default is False\n    \n    Returns:\n        trend feature\n    \"\"\"\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    Calculate classic STA/LTA\n    STA/LTA represents the ratio between amplitude of the\n    signal on a short time window of length LTA and on a\n    long time window LTA\n    Args:\n        length_sta: length of short time average window\n        length_lta: length of long time average window\n    Returns:\n        STA/LTA\n    \"\"\"\n    sta = np.cumsum(x ** 2)\n    # Convert to float\n    sta = np.require(sta, dtype=np.float)\n    # Copy for LTA\n    lta = sta.copy()\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    # Pad zeros\n    sta[:length_lta - 1] = 0\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    return sta / lta","metadata":{"_kg_hide-input":true,"_uuid":"8b46e3d6267988a865b69efc1f92ba90a20ca7b7","execution":{"iopub.status.busy":"2023-07-30T12:50:38.011809Z","iopub.execute_input":"2023-07-30T12:50:38.012200Z","iopub.status.idle":"2023-07-30T12:50:38.021224Z","shell.execute_reply.started":"2023-07-30T12:50:38.012170Z","shell.execute_reply":"2023-07-30T12:50:38.020091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Process train file\n\nNow let's calculate the aggregated functions for train set.  We start by calculating:\n- mean\n- standard deviation\n- maximum\n- minimum\n- sum\n- absolute sum\n","metadata":{"_uuid":"0cec83d75c3498355fb232c33ef43dcd2017f24e"}},{"cell_type":"code","source":"train_X = pd.DataFrame(index=range(segments), dtype=np.float64)\ntrain_y = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\ntotal_mean = train_df['acoustic_data'].mean()\ntotal_std = train_df['acoustic_data'].std()\ntotal_max = train_df['acoustic_data'].max()\ntotal_min = train_df['acoustic_data'].min()\ntotal_sum = train_df['acoustic_data'].sum()\ntotal_abs_sum = np.abs(train_df['acoustic_data']).sum()","metadata":{"_kg_hide-input":true,"_uuid":"c65536bf1bf32037175518c3cc8fc7439e217e9d","execution":{"iopub.status.busy":"2023-07-30T13:27:18.750987Z","iopub.execute_input":"2023-07-30T13:27:18.751384Z","iopub.status.idle":"2023-07-30T13:27:26.419943Z","shell.execute_reply.started":"2023-07-30T13:27:18.751355Z","shell.execute_reply":"2023-07-30T13:27:26.417743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we create signal-processing specific features per each segment, as following:\n\n* Fourier fast transform (FFT) applied to entire segment\n* Mean, Standard deviation, max and min applied to the entire segment\n* Extract real and imaginary part of FFT transform\n* Calculate mean, standard deviation, min and max for both real and imaginary part of FFT transform\n* Calculate mean, standard deviation, min and max for both real and imaginary part of FFT transform for a 5K and 15K data points from the end of segment\n* Mean change for entire segment\n* Absolute change for entire segment\n* Max and min of absolute values (per entire segment)\n* Standard deviation for first and last 50K and 10K values per segment\n* Average value for first and last 50K and 10K values per segment\n* Minimum values for first and last 50K and 10K values per segment\n* Maximum values for first and last 50K and 10K values per segment\n* Ratio of maximum to minimum values for entire segment\n* Difference between max and minimum values for entire segment\n* Mean change rate (excluding nulls) for first and last 10K and 50K data points per segment\n* 01, 05, 95 and 99 quantile values for entire segment\n* absolute values for 01, 05, 95 and 99 quantile values\n* trend values (calculated with `add_trend_values` with absolute flag off)\n* absolute trend values (calculated with `add_trend_values` with absolute flag on)\n* Median absolute deviation (mad), Kurtosis, Skew and Median values\n* Hilbert mean\n* Hann window mean\n* Classical STA/LTA mean for 500-10K, 5K-100K, 3333-6666, and 10K-25K STA/LTA windows\n* Moving averages means for 700, 1.5K, 3K, 6K windows (and excluding nans)\n* Exponentially weighted moving average with spans of 300, 3K and 6K\n* Average standard deviation moving average over 700 and 400 windows\n* Moving averages means for 700 windown plus or minus 2 x average standard deviation moving average over same size window\n* IQR, 001 and 999 quantiles\n* For windows in 10, 100 and 1000, calculate the moving average and moving standard deviation. With these values then calculate  min, max, mean, standard deviation, average absolute and relative change, 01, 05, 95, and 99 quantiles, and absolute max rol.","metadata":{}},{"cell_type":"code","source":"def nonzero(x):\n    \"\"\"\n    Utility function to simplify call of numpy `nonzero` function\n    \"\"\"\n    return np.nonzero(np.atleast_1d(x))\n\ndef create_features(seg_id, seg, X):\n    \"\"\"\n    Create features\n    Args:\n        seg_id: the id of current data segment to process\n        seg: the current selected segment data\n        X: transformed train data\n    Returns:\n        None\n    \"\"\"\n    xc = pd.Series(seg['acoustic_data'].values)\n    zc = np.fft.fft(xc)\n    \n    X.loc[seg_id, 'mean'] = xc.mean()\n    X.loc[seg_id, 'std'] = xc.std()\n    X.loc[seg_id, 'max'] = xc.max()\n    X.loc[seg_id, 'min'] = xc.min()\n    \n    #FFT transform values\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    X.loc[seg_id, 'Rmean'] = realFFT.mean()\n    X.loc[seg_id, 'Rstd'] = realFFT.std()\n    X.loc[seg_id, 'Rmax'] = realFFT.max()\n    X.loc[seg_id, 'Rmin'] = realFFT.min()\n    X.loc[seg_id, 'Imean'] = imagFFT.mean()\n    X.loc[seg_id, 'Istd'] = imagFFT.std()\n    X.loc[seg_id, 'Imax'] = imagFFT.max()\n    X.loc[seg_id, 'Imin'] = imagFFT.min()\n    X.loc[seg_id, 'Rmean_last_5000'] = realFFT[-5000:].mean()\n    X.loc[seg_id, 'Rstd__last_5000'] = realFFT[-5000:].std()\n    X.loc[seg_id, 'Rmax_last_5000'] = realFFT[-5000:].max()\n    X.loc[seg_id, 'Rmin_last_5000'] = realFFT[-5000:].min()\n    X.loc[seg_id, 'Rmean_last_15000'] = realFFT[-15000:].mean()\n    X.loc[seg_id, 'Rstd_last_15000'] = realFFT[-15000:].std()\n    X.loc[seg_id, 'Rmax_last_15000'] = realFFT[-15000:].max()\n    X.loc[seg_id, 'Rmin_last_15000'] = realFFT[-15000:].min()\n    \n    X.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(xc))\n    X.loc[seg_id, 'mean_change_rate'] = np.mean(nonzero(np.diff(xc) / xc[:-1])[0])\n\n    X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n    X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\n    \n    X.loc[seg_id, 'std_first_50000'] = xc[:50000].std()\n    X.loc[seg_id, 'std_last_50000'] = xc[-50000:].std()\n    X.loc[seg_id, 'std_first_10000'] = xc[:10000].std()\n    X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()\n    \n    X.loc[seg_id, 'avg_first_50000'] = xc[:50000].mean()\n    X.loc[seg_id, 'avg_last_50000'] = xc[-50000:].mean()\n    X.loc[seg_id, 'avg_first_10000'] = xc[:10000].mean()\n    X.loc[seg_id, 'avg_last_10000'] = xc[-10000:].mean()\n    \n    X.loc[seg_id, 'min_first_50000'] = xc[:50000].min()\n    X.loc[seg_id, 'min_last_50000'] = xc[-50000:].min()\n    X.loc[seg_id, 'min_first_10000'] = xc[:10000].min()\n    X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()\n    \n    X.loc[seg_id, 'max_first_50000'] = xc[:50000].max()\n    X.loc[seg_id, 'max_last_50000'] = xc[-50000:].max()\n    X.loc[seg_id, 'max_first_10000'] = xc[:10000].max()\n    X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max()\n    \n    X.loc[seg_id, 'max_to_min'] = xc.max() / np.abs(xc.min())\n    X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())\n    X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > 500])\n    X.loc[seg_id, 'sum'] = xc.sum()\n    \n    X.loc[seg_id, 'mean_change_rate_first_50000'] = np.mean(nonzero((np.diff(xc[:50000]) / xc[:50000][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_last_50000'] = np.mean(nonzero((np.diff(xc[-50000:]) / xc[-50000:][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_first_10000'] = np.mean(nonzero((np.diff(xc[:10000]) / xc[:10000][:-1]))[0])\n    X.loc[seg_id, 'mean_change_rate_last_10000'] = np.mean(nonzero((np.diff(xc[-10000:]) / xc[-10000:][:-1]))[0])\n    \n    X.loc[seg_id, 'q95'] = np.quantile(xc, 0.95)\n    X.loc[seg_id, 'q99'] = np.quantile(xc, 0.99)\n    X.loc[seg_id, 'q05'] = np.quantile(xc, 0.05)\n    X.loc[seg_id, 'q01'] = np.quantile(xc, 0.01)\n    \n    X.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(xc), 0.95)\n    X.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(xc), 0.99)\n    X.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(xc), 0.05)\n    X.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(xc), 0.01)\n    \n    X.loc[seg_id, 'trend'] = add_trend_feature(xc)\n    X.loc[seg_id, 'abs_trend'] = add_trend_feature(xc, abs_values=True)\n    X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()\n    X.loc[seg_id, 'abs_std'] = np.abs(xc).std()\n    \n    X.loc[seg_id, 'mad'] = xc.mad()\n    X.loc[seg_id, 'kurt'] = xc.kurtosis()\n    X.loc[seg_id, 'skew'] = xc.skew()\n    X.loc[seg_id, 'med'] = xc.median()\n    \n    X.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(xc)).mean()\n    X.loc[seg_id, 'Hann_window_mean'] = (convolve(xc, hann(150), mode='same') / sum(hann(150))).mean()\n    X.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(xc, 500, 10000).mean()\n    X.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(xc, 5000, 100000).mean()\n    X.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(xc, 3333, 6666).mean()\n    X.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(xc, 10000, 25000).mean()\n    X.loc[seg_id, 'Moving_average_700_mean'] = xc.rolling(window=700).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_1500_mean'] = xc.rolling(window=1500).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_3000_mean'] = xc.rolling(window=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_6000_mean'] = xc.rolling(window=6000).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    X.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(xc, span=300).mean()).mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(xc, span=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(xc, span=6000).mean().mean(skipna=True)\n    no_of_std = 2\n    X.loc[seg_id, 'MA_700MA_std_mean'] = xc.rolling(window=700).std().mean()\n    X.loc[seg_id,'MA_700MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id,'MA_700MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_400MA_std_mean'] = xc.rolling(window=400).std().mean()\n    X.loc[seg_id,'MA_400MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id,'MA_400MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_1000MA_std_mean'] = xc.rolling(window=1000).std().mean()\n    \n    X.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(xc, [75, 25]))\n    X.loc[seg_id, 'q999'] = np.quantile(xc,0.999)\n    X.loc[seg_id, 'q001'] = np.quantile(xc,0.001)\n    X.loc[seg_id, 'ave10'] = stats.trim_mean(xc, 0.1)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = xc.rolling(windows).std().dropna().values\n        x_roll_mean = xc.rolling(windows).mean().dropna().values\n        \n        X.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean(nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n        X.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n        X.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"_kg_hide-input":true,"_uuid":"0745eb336b1313d3c46c957d7b18c19c98076743","execution":{"iopub.status.busy":"2023-07-30T13:49:06.852771Z","iopub.execute_input":"2023-07-30T13:49:06.853200Z","iopub.status.idle":"2023-07-30T13:49:06.899910Z","shell.execute_reply.started":"2023-07-30T13:49:06.853164Z","shell.execute_reply":"2023-07-30T13:49:06.898534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For each data segment of dimmension equal with the dimmension of the test data, we perform the following processing:\n- select the data segment\n- call the `create_features` function to process the data from current segment i.e. create the features\n- set the `time_to_failure` value which is the value for time to feature for the last value in the segment","metadata":{}},{"cell_type":"code","source":"# iterate over all segments\nfor seg_id in tqdm_notebook(range(segments)):\n    seg = train_df.iloc[seg_id*rows:seg_id*rows+rows]\n    create_features(seg_id, seg, train_X)\n    train_y.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1]","metadata":{"_kg_hide-input":true,"_uuid":"5a446450bfd25296d6b98dde3e750e11d96250a4","execution":{"iopub.status.busy":"2023-07-30T13:49:12.174180Z","iopub.execute_input":"2023-07-30T13:49:12.174632Z","iopub.status.idle":"2023-07-30T14:08:04.092741Z","shell.execute_reply.started":"2023-07-30T13:49:12.174595Z","shell.execute_reply":"2023-07-30T14:08:04.091454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check the result. We plot the shape and the head of train_X.","metadata":{"_uuid":"c2e064c25b00c08acd8f8973f565a0b29d09704f"}},{"cell_type":"code","source":"train_X.shape","metadata":{"_kg_hide-input":true,"_uuid":"ff98b8f65290cfa14a5fed507a5567a589109fc9","execution":{"iopub.status.busy":"2023-07-30T14:08:11.787436Z","iopub.execute_input":"2023-07-30T14:08:11.787879Z","iopub.status.idle":"2023-07-30T14:08:11.797561Z","shell.execute_reply.started":"2023-07-30T14:08:11.787845Z","shell.execute_reply":"2023-07-30T14:08:11.796261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_X.head(10)","metadata":{"_kg_hide-input":true,"_uuid":"dc6dd182bf1819c4dcfe823416efca51ff15d5bf","execution":{"iopub.status.busy":"2023-07-30T14:08:15.055072Z","iopub.execute_input":"2023-07-30T14:08:15.055491Z","iopub.status.idle":"2023-07-30T14:08:15.104463Z","shell.execute_reply.started":"2023-07-30T14:08:15.055454Z","shell.execute_reply":"2023-07-30T14:08:15.103192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We scale the data.","metadata":{"_uuid":"58196ccf11c59b8bca96279f6bdc7e66eaad6fbe"}},{"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(train_X)\nscaled_train_X = pd.DataFrame(scaler.transform(train_X), columns=train_X.columns)","metadata":{"_kg_hide-input":true,"_uuid":"5e9f4c996ac9842fd3243e4f9d3a6e77b8b65cf7","execution":{"iopub.status.busy":"2023-07-30T14:08:19.259024Z","iopub.execute_input":"2023-07-30T14:08:19.259445Z","iopub.status.idle":"2023-07-30T14:08:19.285509Z","shell.execute_reply.started":"2023-07-30T14:08:19.259403Z","shell.execute_reply":"2023-07-30T14:08:19.284467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check the obtained dataframe.","metadata":{"_uuid":"885370309aca6e12d22a811d3ccc126205de468d"}},{"cell_type":"code","source":"scaled_train_X.head(10)","metadata":{"_kg_hide-input":true,"_uuid":"1212dd6508332a32bfb4e337f3076ee9680584d4","execution":{"iopub.status.busy":"2023-07-30T14:08:31.817501Z","iopub.execute_input":"2023-07-30T14:08:31.817951Z","iopub.status.idle":"2023-07-30T14:08:31.846788Z","shell.execute_reply.started":"2023-07-30T14:08:31.817918Z","shell.execute_reply":"2023-07-30T14:08:31.845944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Process test data\n\nWe apply the same processing done for the training data to the test data.\n\nWe read the submission file and prepare the test file.","metadata":{"_uuid":"500232cf8ff46fd1d0a9a43088c815d9624a222a"}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\ntest_X = pd.DataFrame(columns=train_X.columns, dtype=np.float64, index=submission.index)","metadata":{"_kg_hide-input":true,"_uuid":"38e67f8e457d7c7ce4900b454d3d4a1ab5c42d04","execution":{"iopub.status.busy":"2023-07-30T14:09:16.746145Z","iopub.execute_input":"2023-07-30T14:09:16.746661Z","iopub.status.idle":"2023-07-30T14:09:16.787632Z","shell.execute_reply.started":"2023-07-30T14:09:16.746622Z","shell.execute_reply":"2023-07-30T14:09:16.786468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check the shape of the submission and test_X datasets.","metadata":{"_uuid":"b9b80332f6e9af2ee82f8e0feca4e9d1d18925ee"}},{"cell_type":"code","source":"submission.shape, test_X.shape","metadata":{"_kg_hide-input":true,"_uuid":"1f06fee3edb09ac05890fa109808487de395d1de","execution":{"iopub.status.busy":"2023-07-30T14:09:20.506269Z","iopub.execute_input":"2023-07-30T14:09:20.506695Z","iopub.status.idle":"2023-07-30T14:09:20.513624Z","shell.execute_reply.started":"2023-07-30T14:09:20.506662Z","shell.execute_reply":"2023-07-30T14:09:20.512477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for seg_id in tqdm_notebook(test_X.index):\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    create_features(seg_id, seg, test_X)","metadata":{"_kg_hide-input":true,"_uuid":"e3b54cd654438c2af63060b03ddb767eacbd70c6","execution":{"iopub.status.busy":"2023-07-30T14:09:28.974196Z","iopub.execute_input":"2023-07-30T14:09:28.975006Z","iopub.status.idle":"2023-07-30T14:22:33.356618Z","shell.execute_reply.started":"2023-07-30T14:09:28.974969Z","shell.execute_reply":"2023-07-30T14:22:33.355414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We scale also the test data.","metadata":{"_uuid":"5630b8d9234af698ed83af94980538e9ed8eaa59"}},{"cell_type":"code","source":"scaled_test_X = pd.DataFrame(scaler.transform(test_X), columns=test_X.columns)","metadata":{"_kg_hide-input":true,"_uuid":"a18c26409ced2a71c3b7b4926d05455afd5fa26a","execution":{"iopub.status.busy":"2023-07-30T14:22:33.359068Z","iopub.execute_input":"2023-07-30T14:22:33.360175Z","iopub.status.idle":"2023-07-30T14:22:33.371302Z","shell.execute_reply.started":"2023-07-30T14:22:33.360130Z","shell.execute_reply":"2023-07-30T14:22:33.370297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaled_test_X.shape","metadata":{"_kg_hide-input":true,"_uuid":"9c5701e9cc1c9076bde818af54fa67ec4ad3eb4b","execution":{"iopub.status.busy":"2023-07-30T14:22:33.373224Z","iopub.execute_input":"2023-07-30T14:22:33.374072Z","iopub.status.idle":"2023-07-30T14:22:33.384101Z","shell.execute_reply.started":"2023-07-30T14:22:33.374029Z","shell.execute_reply":"2023-07-30T14:22:33.382675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaled_test_X.tail(10)","metadata":{"_kg_hide-input":true,"_uuid":"7bf4d31bea2f1fd82fbfdc1a4a84f36834daf5c3","execution":{"iopub.status.busy":"2023-07-30T14:22:33.386933Z","iopub.execute_input":"2023-07-30T14:22:33.387640Z","iopub.status.idle":"2023-07-30T14:22:33.421228Z","shell.execute_reply.started":"2023-07-30T14:22:33.387596Z","shell.execute_reply":"2023-07-30T14:22:33.420118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# <a id='5'>Model</a>  \n\nLet's prepare the model.\n\n## Run model\n\nWe define the folds for cross-validation.","metadata":{"_uuid":"2bb7ba0d1ec723be641adf5f925c24237c9198c0"}},{"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\ntrain_columns = scaled_train_X.columns.values","metadata":{"_kg_hide-input":true,"_uuid":"99d43c789dde8cac37270dc16430cf98870788d3","execution":{"iopub.status.busy":"2023-07-30T14:22:33.423159Z","iopub.execute_input":"2023-07-30T14:22:33.423563Z","iopub.status.idle":"2023-07-30T14:22:33.434669Z","shell.execute_reply.started":"2023-07-30T14:22:33.423531Z","shell.execute_reply":"2023-07-30T14:22:33.433607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We define the model parameters.","metadata":{"_uuid":"d506ac30835fbfd5ed324bab3a3370e01739f907"}},{"cell_type":"code","source":"params = {'num_leaves': 51,\n         'min_data_in_leaf': 10, \n         'objective':'regression',\n         'max_depth': -1,\n         'learning_rate': 0.001,\n         \"boosting\": \"gbdt\",\n         \"feature_fraction\": 0.91,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.91,\n         \"bagging_seed\": 42,\n         \"metric\": 'mae',\n         \"lambda_l1\": 0.1,\n         \"verbosity\": -1,\n         \"nthread\": -1,\n         \"random_state\": 42}","metadata":{"_kg_hide-input":true,"_uuid":"b3fc19ff1bf1bb0897e29747421514f211ea0edc","execution":{"iopub.status.busy":"2023-07-30T14:22:33.436262Z","iopub.execute_input":"2023-07-30T14:22:33.436630Z","iopub.status.idle":"2023-07-30T14:22:33.455733Z","shell.execute_reply.started":"2023-07-30T14:22:33.436555Z","shell.execute_reply":"2023-07-30T14:22:33.454523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We run the model. During training for each fold, we validate using the validation set and also we predict using the current model for test set. The final result will be the average over the all folds for the predictions done at each fold training.","metadata":{"_uuid":"5f5e691c0c52f7adfb13064794c7bf2600f07f39"}},{"cell_type":"code","source":"oof = np.zeros(len(scaled_train_X))\npredictions = np.zeros(len(scaled_test_X))\nfeature_importance_df = pd.DataFrame()\n#run model\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_train_X,train_y.values)):\n    strLog = \"fold {}\".format(fold_)\n    print(strLog)\n    \n    X_tr, X_val = scaled_train_X.iloc[trn_idx], scaled_train_X.iloc[val_idx]\n    y_tr, y_val = train_y.iloc[trn_idx], train_y.iloc[val_idx]\n\n    model = lgb.LGBMRegressor(**params, n_estimators = 20000, n_jobs = -1)\n    model.fit(X_tr, \n              y_tr, \n              eval_set=[(X_tr, y_tr), (X_val, y_val)], \n              eval_metric='mae',\n              verbose=1000, \n              early_stopping_rounds=500)\n    oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)\n    #feature importance\n    fold_importance_df = pd.DataFrame()\n    fold_importance_df[\"Feature\"] = train_columns\n    fold_importance_df[\"importance\"] = model.feature_importances_[:len(train_columns)]\n    fold_importance_df[\"fold\"] = fold_ + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n    #predictions\n    predictions += model.predict(scaled_test_X, num_iteration=model.best_iteration_) / folds.n_splits","metadata":{"_kg_hide-input":true,"_uuid":"d0b278f0a14dbd3dee610e2e27a0cc6d55844678","execution":{"iopub.status.busy":"2023-07-30T14:29:18.150146Z","iopub.execute_input":"2023-07-30T14:29:18.150560Z","iopub.status.idle":"2023-07-30T14:39:53.635028Z","shell.execute_reply.started":"2023-07-30T14:29:18.150529Z","shell.execute_reply":"2023-07-30T14:39:53.634104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Features importance\n\nLet's print features importance.","metadata":{"_uuid":"cc8966810c47ae092a300dbb3a5ec582c9740b8f"}},{"cell_type":"code","source":"cols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:200].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,26))\nsns.barplot(x=\"importance\", y=\"Feature\", data=best_features.sort_values(by=\"importance\",ascending=False))\nplt.title('LightGBM Features (averaged over folds)')\nplt.tight_layout()\nplt.savefig('lgbm_importances_200.png')","metadata":{"_kg_hide-input":true,"_uuid":"0f64306227bc30ffa306b5821a63f04ab6a874c0","execution":{"iopub.status.busy":"2023-07-30T14:39:53.637286Z","iopub.execute_input":"2023-07-30T14:39:53.638199Z","iopub.status.idle":"2023-07-30T14:40:01.115495Z","shell.execute_reply.started":"2023-07-30T14:39:53.638161Z","shell.execute_reply":"2023-07-30T14:40:01.114411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:20].index)\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,10))\nsns.barplot(x=\"importance\", y=\"Feature\", data=best_features.sort_values(by=\"importance\",ascending=False))\nplt.title('Top 20 LightGBM Features (averaged over folds)')\nplt.tight_layout()\nplt.savefig('lgbm_importances_20.png')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <a id='6'>Submission</a>  \n\nWe set the predicted time to failure in the submission file.","metadata":{"_uuid":"0fce1845bdc118dea09c917d2322c840426e34e1"}},{"cell_type":"code","source":"submission.time_to_failure = predictions\nsubmission.to_csv('submission.csv',index=True)","metadata":{"_kg_hide-input":true,"_uuid":"641d4feeb57e87e24c3f4bc29863c35169f9302f","execution":{"iopub.status.busy":"2023-07-30T14:40:38.980109Z","iopub.execute_input":"2023-07-30T14:40:38.980493Z","iopub.status.idle":"2023-07-30T14:40:38.994121Z","shell.execute_reply.started":"2023-07-30T14:40:38.980465Z","shell.execute_reply":"2023-07-30T14:40:38.993098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <a id='7'>References</a>  \n\n[1] Fast Fourier Transform, https://en.wikipedia.org/wiki/Fast_Fourier_transform   \n[2] Shifting aperture, in Neural network for inverse mapping in eddy current testing, https://www.researchgate.net/publication/3839126_Neural_network_for_inverse_mapping_in_eddy_current_testing   \n[3] Andrews Script plus a Genetic Program Model, https://www.kaggle.com/scirpus/andrews-script-plus-a-genetic-program-model/","metadata":{"_uuid":"31de440320aa87a962deb1d166abba4ac82e4c04"}}]}