{"cells":[{"metadata":{},"cell_type":"markdown","source":"## General Infomation\n\nCorrect prediction and warning on earthquake indirectly saves lifes from damage of buildings.\n\nIn this competition, we train on large dataset having several laboratory earthquakes continuously,\nand predict `time_to_failure` in each last sample of segment to next laboratory earthquake.\nEach segment we usually have about 150,000 continous samples.\n\n## Summary on EDA\n\n* Sampling rate is about 3.853 MHz\n\n* `time_to_failure` within a lab earthquake is in stairs-like descending order\n\n* About 0.31sec after each huge fluctuations, next laboratory earthquake comes.\n\n* After Fast Fourier Transformation, the main difference between huge fluctuation chunks and others is their distribution of real numbers."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport gc # garbage collection\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# change plot style and default figure size\n# plt.style.use('seaborn')\nplt.rc(\"font\", size=13)\nplt.rc(\"figure\", figsize=(14.4, 8.1), dpi=72)\n# plt.rc(\"savefig\", dpi=72)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(\"../input/train.csv\", dtype={'acoustic_data': 'int16', 'time_to_failure': 'float64'})\ngc.collect()\n# print(len(train_df))  # >> length of full_dataframe: 629,145,480","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot slice in dataframe\ndef df_slice_plot(df, start=0, stop=None, step=1, figsize=None):\n\n    if not stop:\n        start, stop = 0, start\n    if start == 0 and start == stop:\n        stop = len(df)\n\n    ps = \"\"\n    if step > 1:\n        ps = \"(down sampled)\"\n\n    df_slice = df.iloc[start:stop:step]\n\n    df_slice.acoustic_data.plot(kind='line', label=\"Signals\", figsize=figsize, legend=True)\n    df_slice.time_to_failure.plot(kind='line', label=\"TimeToFailure\", figsize=figsize, legend=True, secondary_y=True)\n    plt.title(\"Signals and time_to_failure within a slice from {:,d} to {:,d} {}\".format(start, stop, ps))\n    plt.show()\n\n    gc.collect()\n\n\ndef before_next_lab_earthquake(df, start=0, stop=None):\n    \n    if not stop:\n        start, stop = 0, start\n    \n    df_slice = df.iloc[start:stop].copy().reset_index()\n        \n    sig_top_at = df_slice.acoustic_data.idxmax()\n    sig_bot_at = df_slice.acoustic_data.idxmin()\n    t_2_fail_min_at = df_slice.time_to_failure.idxmin()\n    \n    top_2_closest = df_slice.time_to_failure.iloc[sig_top_at] - df_slice.time_to_failure.min()\n    bot_2_closest = df_slice.time_to_failure.iloc[sig_bot_at] - df_slice.time_to_failure.min()\n    \n    top_2_closest_x = t_2_fail_min_at - sig_top_at\n    bot_2_closest_x = t_2_fail_min_at - sig_bot_at\n    \n    sample_rate = top_2_closest_x / top_2_closest\n    \n    print((\n        \"In slice from {:,d} to {:,d}:\\n\"\n        \"Time from top value of signal to the bottom value of time_to_failure: {:.8f} sec.\\n\"\n        \"Time from bottom value of signal to the bottom value of time_to_failure: {:.8f} sec.\").format(\n            start, stop, top_2_closest, bot_2_closest))\n    print((\n        \"The location of top val is {:,d} samples to the bottom value of time_to_failure.\\n\"\n        \"The location of bot val is {:,d} samples to the bottom value of time_to_failure.\").format(\n            top_2_closest_x, bot_2_closest_x))\n    print(\"The sample rate is about: {} Hz.\\n\".format(sample_rate))\n\n\ndef plot_transformed(signals):\n\n    zc = np.fft.fft(signals)\n    freq = np.fft.fftfreq(signals.shape[-1])\n\n    fig, axes = plt.subplots(2, 1, figsize=(14.4, 16.2))\n\n    axes[0].plot(zc.real, zc.imag)\n    axes[0].set_aspect(\"equal\")\n    axes[0].set_title(\"Fourier transformed signals\")\n\n    axes[1].plot(freq, zc.real, label=\"Real\", alpha=0.6)\n    axes[1].plot(freq, zc.imag, label=\"Image\", alpha=0.5)\n    axes[1].legend()\n    axes[1].set_title(\"Fourier transformed signals in freq\")\n\n    plt.show()\n    \n    print(\"Mean in real: {:f}, in image: {:.16f}\".format(zc.real.mean(), zc.imag.mean()))\n    print(\"Std in real: {:f}, in image: {:.16f}\".format(zc.real.std(), zc.imag.std()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore time_to_failure"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.time_to_failure.iloc[:50_000].plot(kind='line', title=\"time_to_failure within a lab earthquake (stairs-like)\")\nplt.show()\n\ntrain_df.time_to_failure.iloc[:4_001].plot(kind='line', title=\"time_to_failure within a stair\")\nplt.show()\n\ngc.collect()\n\ntrain_df.time_to_failure.iloc[:50_000].diff().plot(kind='line', title=\"Diffs of time_to_failure within a lab earthquake\")\nplt.show()\n\ntrain_df.time_to_failure.iloc[5_656_570:5_656_580].plot(\n    kind='bar', logy=True,\n    title=\"Logarithmic(base 10) time_to_failure between 2 laboratory earthquakes\")\nplt.show()\n\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore the acoustic_data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Down sampling is risky to loss peak in huge fluctuations about 0.31sec before next laboratory earthquake (compared with fig1 and fig4)\ndf_slice_plot(train_df, 10_000_000)  # fig1\ndf_slice_plot(train_df, 10_000_000, 30_000_000)  # fig2\ndf_slice_plot(train_df, 30_000_000, 60_000_000)  # fig3\ndf_slice_plot(train_df, step=20, figsize=(24, 10))  # fig4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_signal = train_df.acoustic_data.iloc[4_435_000:4_445_000]\nslice_signal.plot(kind='line', title=\"Signals within huge fluctuations.\")\nplt.show()\n\nplot_transformed(slice_signal)\n\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Short before next lab earthquake"},{"metadata":{"trusted":true},"cell_type":"code","source":"before_next_lab_earthquake(train_df, 4_000_000, 6_000_000)\nbefore_next_lab_earthquake(train_df, 47_500_000, 51_000_000)\n\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore first 5 test segments"},{"metadata":{"trusted":true},"cell_type":"code","source":"seg_lst = list(os.listdir(\"../input/test/\"))\nfor f in seg_lst[:5]:\n    seg_df = pd.read_csv(f\"../input/test/{f}\")\n    seg_df.plot()\n    plt.title(f\"{f}\")\n    plt.show()\n    plot_transformed(seg_df.acoustic_data)\n    del seg_df\n    \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore test segments having huge fluctuations"},{"metadata":{"trusted":true},"cell_type":"code","source":"c = 0\nfor f in seg_lst:\n    seg_df = pd.read_csv(f\"../input/test/{f}\")\n    if seg_df.max().values > 2000:\n        print(f\"Huge fluctuations in {f}\")\n        seg_df.plot()\n        plt.title(f\"{f}\")\n        plt.show()\n        plot_transformed(seg_df.acoustic_data)\n        c += 1\n    del seg_df\n    if c >= 5:\n        break\n        \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Is it time to build models?\n\n* It's such a huge amount of samples that machine cannot learn quickly\n\n* So we aggregate on features to make derived dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"# drop the dataframe\ndel train_df\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering\n\n* Using threshold which does not gain feature importance does no helps.\n\n* `min`, `max` have extremely high correlation coefficient with `std`."},{"metadata":{"trusted":true},"cell_type":"code","source":"derived_trn = pd.DataFrame()\ntrain_rd = pd.read_csv(\"../input/train.csv\", iterator=True, chunksize=150000)\n\n\n# Derive from `acoustic_data`\ndef get_features(chunk):\n\n    curr_df = pd.DataFrame()\n\n    curr_df[\"mean\"] = [chunk.acoustic_data.mean()]\n    curr_df[\"std\"] = [chunk.acoustic_data.std()]\n    curr_df[\"kurtosis\"] = [chunk.acoustic_data.kurtosis()]\n    curr_df[\"skew\"] = [chunk.acoustic_data.skew()]\n    curr_df[\"quantile_05\"], curr_df[\"quantile_25\"], curr_df[\"quantile_75\"], curr_df[\"quantile_95\"] = \\\n        [[x] for x in chunk.acoustic_data.quantile([0.05, 0.25, 0.75, 0.95])]\n\n    # Attempt to apply window function on `acoustic_data`\n    win_width = 100\n    slepian_width = 51\n    windowed = chunk.acoustic_data.rolling(win_width, win_type='slepian').mean(width=slepian_width).dropna()\n    curr_df[f\"window_{win_width}_mean\"] = [windowed.mean()]\n    curr_df[f\"window_[win_width]_std\"] = [windowed.std()]\n\n    # Fast Fourier Transformed arrays\n    fft = np.fft.fft(chunk.acoustic_data)\n    curr_df[\"fft_real_mean\"] = [fft.real.mean()]\n    curr_df[\"fft_image_mean\"] = [fft.imag.mean()]\n    curr_df[\"fft_real_std\"] = [fft.real.std()]\n    curr_df[\"fft_image_std\"] = [fft.imag.std()]\n\n    # Aggregate on beginning number of samples.\n    # If a chunk (the last chunk) is not more than the num of samples, use all samples in the chunk.\n    sig_first = chunk.acoustic_data.iloc[:20000]\n    curr_df[\"mean_first\"] = [sig_first.mean()]\n    curr_df[\"std_first\"] = [sig_first.std()]\n    curr_df[\"kurtosis_first\"] = [sig_first.kurtosis()]\n    curr_df[\"skew_first\"] = [sig_first.skew()]\n    curr_df[\"quantile_05_first\"], curr_df[\"quantile_25_first\"], \\\n    curr_df[\"quantile_75_first\"], curr_df[\"quantile_95_first\"] = [[x] for x in sig_first.quantile([0.05, 0.25, 0.75, 0.95])]\n    \n    windowed_first = sig_first.rolling(win_width, win_type='slepian').mean(width=slepian_width).dropna()\n    curr_df[f\"window_{win_width}_mean_first\"] = [windowed_first.mean()]\n    curr_df[f\"window_{win_width}_std_first\"] = [windowed_first.std()]\n\n    fft_first = np.fft.fft(sig_first)\n    curr_df[\"fft_real_mean_first\"] = [fft_first.real.mean()]\n    curr_df[\"fft_image_mean_first\"] = [fft_first.imag.mean()]\n    curr_df[\"fft_real_std_first\"] = [fft_first.real.std()]\n    curr_df[\"fft_image_std_first\"] = [fft_first.imag.std()]\n\n    # Aggregate on last number of samples.\n    sig_last = chunk.acoustic_data.iloc[-20000:]\n    curr_df[\"mean_last\"] = [sig_last.mean()]\n    curr_df[\"std_last\"] = [sig_last.std()]\n    curr_df[\"kurtosis_last\"] = [sig_last.kurtosis()]\n    curr_df[\"skew_last\"] = [sig_last.skew()]\n    curr_df[\"quantile_05_last\"], curr_df[\"quantile_25_last\"], \\\n    curr_df[\"quantile_75_last\"], curr_df[\"quantile_95_last\"] = [[x] for x in sig_last.quantile([0.05, 0.25, 0.75, 0.95])]\n    \n    windowed_last = sig_last.rolling(win_width, win_type='slepian').mean(width=slepian_width).dropna()\n    curr_df[f\"window_{win_width}_mean_last\"] = [windowed_last.mean()]\n    curr_df[f\"window_{win_width}_std_last\"] = [windowed_last.std()]\n\n    fft_last = np.fft.fft(sig_last)\n    curr_df[\"fft_real_mean_last\"] = [fft_last.real.mean()]\n    curr_df[\"fft_image_mean_last\"] = [fft_last.imag.mean()]\n    curr_df[\"fft_real_std_last\"] = [fft_last.real.std()]\n    curr_df[\"fft_image_std_last\"] = [fft_last.imag.std()]\n\n    return curr_df\n\n\nfor chunk in train_rd:\n    curr_df = get_features(chunk)\n    curr_df[\"target\"] = chunk.time_to_failure.iloc[-1]\n    derived_trn = pd.concat([derived_trn, curr_df], ignore_index=True)\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"derived_trn.to_csv(\"derived_train.csv\", index=False)\ndisplay(derived_trn.head())\ndisplay(derived_trn.loc[(derived_trn.target <= 0.33) & (derived_trn.target >= 0.3)].head())\nfeatures = [col for col in derived_trn.columns if col != \"target\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_seg_list = list(os.listdir(\"../input/test/\"))\nderived_test = pd.DataFrame()\n\nfor csv_file in test_seg_list:\n    curr_test_df = pd.read_csv(f\"../input/test/{csv_file}\", dtype={\"acoustic_data\": \"int16\"})\n\n    curr_df = get_features(curr_test_df)\n    curr_df[\"seg_id\"] = [csv_file.split(\".\")[0]]\n\n    derived_test = pd.concat([derived_test, curr_df], ignore_index=True)\n\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"derived_test.to_csv(\"derived_test.csv\", index=False)\ndisplay(derived_test.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Exploration on Derived datasets -- Is it worthy to train? Why?"},{"metadata":{"trusted":true},"cell_type":"code","source":"display(derived_trn.describe())\ndisplay(derived_test.describe())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nbins=None\nthreshold=0.35\n    \nfor col in features:\n    fig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\n    sns.distplot(derived_trn[col].loc[derived_trn[\"target\"]>=threshold], bins=nbins, label=r\"`target`>={}\".format(threshold), ax=axes[0])\n    sns.distplot(derived_trn[col].loc[derived_trn[\"target\"]<threshold], bins=nbins, label=r\"`target`<{}\".format(threshold), ax=axes[0])\n    axes[0].set_title(\"derived_trn\")\n    axes[0].legend()\n\n    sns.distplot(derived_test[col], bins=nbins, ax=axes[1])\n    axes[1].set_title(\"derived_test\")\n    plt.tight_layout()\n    plt.show()\n\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(2, 1, figsize=(22, 40))\nsns.heatmap(data=derived_trn.corr().abs().round(2), ax=axes[0], annot=False)\naxes[0].set_title(\"Correlation coefficients in derived_trn\")\nsns.heatmap(data=derived_test.corr().abs().round(2), ax=axes[1], annot=False)\naxes[1].set_title(\"Correlation coefficients in derived_test\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model building and training -- Don't overfit"},{"metadata":{"trusted":true},"cell_type":"code","source":"from lightgbm import LGBMRegressor\nfrom sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold\nfrom sklearn.metrics import mean_absolute_error\n\nparams = dict(\n    objective='regression_l1',  # loss func, candidates: huber, regression_l1, regression_l2, fair\n    n_estimators=5000, learning_rate=0.01,\n    num_leaves=31, max_depth=-1,\n#     min_split_gain=10,\n#     min_child_weight=1,\n#     min_child_samples=90,\n#     subsample=1., subsample_freq=0,\n#     colsample_bytree=1.,\n    reg_alpha=0.3, reg_lambda=0.5,\n#     ramdom_state=42,\n    metric='mae',\n)\n\n\n# samples augmentation\ndef augment(df, threshold=None, repeat_times=1):\n    \n    df = df.copy()\n    \n    for _ in range(repeat_times):\n        df_copy = (df.loc[df[\"target\"] <= threshold].copy()\n                   if threshold is not None\n                   else df.copy())\n        df = pd.concat([df, df_copy], ignore_index=True)\n    \n    df = df.sample(frac=1, random_state=42)  # shuffle samples\n    return df","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"param_grid = dict(\n    n_estimators=[1000], learning_rate=[0.01],\n#     min_split_gain=[0.001, 10.],\n)\n\n\ndef search_params(estimator, param_grid, params):\n    \n    ori_n_estimators = params[\"n_estimators\"]\n    \n    gscv = GridSearchCV(estimator, cv=5, param_grid=param_grid, scoring='neg_mean_absolute_error', verbose=3)\n    best = gscv.fit(derived_trn[features], derived_trn[\"target\"])\n    \n    params.update(gscv.best_params_)\n    params[\"n_estimators\"] = ori_n_estimators  # keep n_estimators not changed\n    print(gscv.best_params_)\n    \n    return params\n\n\n# params = search_params(LGBMRegressor(**params), param_grid, params)  # No output if this line is commented","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"n_splits = 5\nn_repeats = 25\n\nrskf = RepeatedStratifiedKFold(n_repeats=n_repeats, n_splits=n_splits, random_state=42)\n\npredictions = pd.DataFrame()\nfeature_importances = pd.DataFrame()\noof = np.zeros((len(derived_trn), n_repeats))\n\nfor nf, (trn_idx, val_idx) in enumerate(rskf.split(derived_trn, (derived_trn['target'] <= 1.2))):\n    \n#     augmented_df = augment(derived_trn.iloc[trn_idx], 1, 2)\n#     trn_features, trn_targets = augmented_df[features], augmented_df[\"target\"]\n    trn_features, trn_targets = derived_trn[features].iloc[trn_idx], derived_trn[\"target\"].iloc[trn_idx]\n    val_features, val_targets = derived_trn[features].iloc[val_idx], derived_trn[\"target\"].iloc[val_idx]\n    \n    print(\"\\nNf: {}, train on {} samples, val on {} samples\".format(nf, len(trn_targets), len(val_targets)))\n    \n    model = LGBMRegressor(**params)\n    gc.collect()\n\n    fit_params = dict(\n        eval_set=[(trn_features, trn_targets), (val_features, val_targets)],\n        verbose=1000,\n        early_stopping_rounds=1500,\n    )\n\n    model.fit(trn_features, trn_targets, **fit_params)\n    \n    curr_importance = pd.DataFrame({\"feature\": features, \"importance\": model.feature_importances_})\n    feature_importances = pd.concat([feature_importances, curr_importance])\n    \n    oof[val_idx, nf // n_splits] = model.predict(val_features)\n    \n    curr_test_pred = model.predict(derived_test[features])\n    predictions[f\"pred_{nf}\"] = curr_test_pred","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## MAE in trainset and Feature importances"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"OOF's mae: {}\".format(mean_absolute_error(derived_trn.target, oof.mean(axis=1))))\nfeature_importances.groupby(\"feature\").mean().sort_values(by=\"importance\", ascending=True).plot(kind=\"barh\", figsize=(14.4, 10.8))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_y = predictions.mean(axis=1).values\nsubmission = pd.DataFrame({\"seg_id\": derived_test[\"seg_id\"], \"time_to_failure\": submission_y})\nsubmission.to_csv(\"submission.csv\", index=False)\ndisplay(submission.head())\ndisplay(submission.tail())\nprint(\"Minimum time_to_failure in submission:\", submission[\"time_to_failure\"].min(), \"sec.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Discussion\n\n* How to exactly apply window function in signal processing?(I doubt that I did not exactly apply)\n\n* Would smaller chunksize in trainset help to improve model?"}],"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}