{"cells":[{"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 pyarrow.parquet as pq\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy.signal import *\nimport gc\nfrom sklearn.feature_selection import f_classif\nimport lightgbm as lgbm\nfrom sklearn.model_selection import RandomizedSearchCV\nfrom scipy.stats import expon, uniform, norm\nfrom scipy.stats import randint, poisson\nfrom sklearn.metrics import confusion_matrix, make_scorer\n\nsns.set(style=\"darkgrid\", context=\"notebook\")\nrand_seed = 135\nnp.random.seed(rand_seed)\nxsize = 12.0\nysize = 8.0\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true,"_uuid":"774abb22b584cf423055462dbdfe548267d42737"},"cell_type":"code","source":"def reduce_mem_usage(df, verbose=True):\n    numerics = [\"int16\", \"int32\", \"int64\", \"float16\", \"float32\", \"float64\"]\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)    \n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print(\"Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)\".format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%%time\n\ntrain_meta_df = pd.read_csv(\"../input/metadata_train.csv\")\ntrain_df = pq.read_pandas(\"../input/train.parquet\").to_pandas()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"51116a22a9dd7f5ab74c154d5a873d8ea8a30cde"},"cell_type":"code","source":"%%time\n\ntrain_meta_df = reduce_mem_usage(train_meta_df)\ntrain_df = reduce_mem_usage(train_df)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1c6c57a67b65bfb45cd19f266e1bfeb71dfca713"},"cell_type":"code","source":"train_meta_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e2086636b37d115157833d424d356a69dd8c5899"},"cell_type":"code","source":"train_meta_df.head(6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"77751a49ac04b3d27c6e750b5bbbae2985a721ba"},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1182f8ca6f1403877a9630811872b1a5e2cab347"},"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2)\nfig.set_size_inches(xsize, 2.0*ysize)\n\nsns.countplot(x=\"phase\", data=train_meta_df, ax=axes[0])\n\nsns.countplot(x=\"target\", data=train_meta_df, ax=axes[1])\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2d33ecc99e93e5fd992eaac57c39e3d191837d76"},"cell_type":"code","source":"fig, ax = plt.subplots()\nfig.set_size_inches(xsize, ysize)\n\nsns.countplot(x=\"phase\", hue=\"target\", data=train_meta_df, ax=ax)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ab49ac8bebe33ccd425dc9f9d3259b09eccdf881"},"cell_type":"markdown","source":"So the phase counts are all equal, so this will not be a useful variable on its own for detecting a fault. Furthermore, it's interesting to not that the target much more likely to be 0, or the line has no fault, by default. This might might make models difficult to calibrate later on, but that's a later issue.\n\nNow let's take a look at some of these signals."},{"metadata":{"trusted":true,"_uuid":"77a870cd372f0c15bd76fb642b383f641cb71967"},"cell_type":"code","source":"fig, axes = plt.subplots(nrows=3, ncols=2)\nfig.set_size_inches(2.0*xsize, 2.0*ysize)\naxes = axes.flatten()\n\naxes[0].plot(train_df[\"0\"].values, marker=\"o\", linestyle=\"none\")\naxes[0].set_title(\"Signal ID: 0\")\n\naxes[1].plot(train_df[\"2\"].values, marker=\"o\", linestyle=\"none\")\naxes[1].set_title(\"Signal ID: 1\")\n\naxes[2].plot(train_df[\"3\"].values, marker=\"o\", linestyle=\"none\")\naxes[2].set_title(\"Signal ID: 2\")\n\naxes[3].plot(train_df[\"4\"].values, marker=\"o\", linestyle=\"none\")\naxes[3].set_title(\"Signal ID: 3\")\n\naxes[4].plot(train_df[\"5\"].values, marker=\"o\", linestyle=\"none\")\naxes[4].set_title(\"Signal ID: 4\")\n\naxes[5].plot(train_df[\"6\"].values, marker=\"o\", linestyle=\"none\")\naxes[5].set_title(\"Signal ID: 5\")\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d6492a6cef612c7c42d4dafd1fb46e83c3ac07d7"},"cell_type":"markdown","source":"Note signals 0, 1, and 2 are not faulty and signals 3, 4, and 5 are faulty. They're messy, noisy, and not obviously periodic, oh boy. However, there are quite a few signal processing techniques that can be used anyways. Speaking of which, it's time for some feature engineering. Starting with some basic aggregations."},{"metadata":{"trusted":true,"_uuid":"914e11872c89d4fa9e8535135107c7392cb7f925"},"cell_type":"code","source":"%%time\n\ntrain_meta_df[\"signal_mean\"] = train_df.agg(np.mean).values\ntrain_meta_df[\"signal_sum\"] = train_df.agg(np.sum).values\ntrain_meta_df[\"signal_std\"] = train_df.agg(np.std).values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a5e1393ea1185743d381b7b83e18e113f5adf970"},"cell_type":"code","source":"train_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e881d67dfb5018c76045f47d6f75393914b66021"},"cell_type":"markdown","source":"Now to look into some power spectrums since this is a signal processing challenge after all."},{"metadata":{"trusted":true,"_uuid":"1a2e73fc4704530930475200a5dc236c68ba326d"},"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=2)\nfig.set_size_inches(2.0*xsize, 2.0*ysize)\naxes = axes.flatten()\n\nf, Pxx = welch(train_df[\"0\"].values)\naxes[0].plot(f, Pxx, marker=\"o\", linestyle=\"none\")\naxes[0].set_title(\"Signal ID: 0\")\naxes[0].axhline(y=2.5, color=\"k\", linestyle=\"--\")\n\nf, Pxx = welch(train_df[\"1\"].values)\naxes[1].plot(f, Pxx, marker=\"o\", linestyle=\"none\")\naxes[1].set_title(\"Signal ID: 1\")\naxes[1].axhline(y=2.5, color=\"k\", linestyle=\"--\")\n\nf, Pxx = welch(train_df[\"2\"].values)\naxes[2].plot(f, Pxx, marker=\"o\", linestyle=\"none\")\naxes[2].set_title(\"Signal ID: 2\")\naxes[2].axhline(y=2.5, color=\"k\", linestyle=\"--\")\n\nf, Pxx = welch(train_df[\"3\"].values)\naxes[3].plot(f, Pxx, marker=\"o\", linestyle=\"none\")\naxes[3].set_title(\"Signal ID: 3\")\naxes[3].axhline(y=2.5, color=\"k\", linestyle=\"--\")\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b1805fb8ccc75cfba84d49cd72e04130803679b9"},"cell_type":"code","source":"%%time\n\ndef welch_max_power_and_frequency(signal):\n    f, Pxx = welch(signal)\n    ix = np.argmax(Pxx)\n    strong_count = np.sum(Pxx>2.5)\n    avg_amp = np.mean(Pxx)\n    sum_amp = np.sum(Pxx)\n    std_amp = np.std(Pxx)\n    median_amp = np.median(Pxx)\n    return [Pxx[ix], f[ix], strong_count, avg_amp, sum_amp, std_amp, median_amp]\n\npower_spectrum_summary = train_df.apply(welch_max_power_and_frequency, result_type=\"expand\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e1f81943c683e8f2a28845c456f6e26a9843e8f8"},"cell_type":"code","source":"power_spectrum_summary = power_spectrum_summary.T.rename(columns={0:\"max_amp\", 1:\"max_freq\", 2:\"strong_amp_count\", 3:\"avg_amp\", \n                                                                  4:\"sum_amp\", 5:\"std_amp\", 6:\"median_amp\"})\npower_spectrum_summary.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c850353413a07efc0c9b44554891784bf1c45769"},"cell_type":"code","source":"power_spectrum_summary.index = power_spectrum_summary.index.astype(int)\ntrain_meta_df = train_meta_df.merge(power_spectrum_summary, left_on=\"signal_id\", right_index=True)\ntrain_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"11b072a913435ae33a90146a1d05a4563a7fe688"},"cell_type":"code","source":"X_cols = [\"phase\"] + train_meta_df.columns[4:].tolist()\nX_cols","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0d9bcb358edaee3b24cb946839d3531395969fcc"},"cell_type":"code","source":"Fvals, pvals = f_classif(train_meta_df[X_cols], train_meta_df[\"target\"])\n\nprint(\"F-value | P-value | Feature Name\")\nprint(\"--------------------------------\")\n\nfor i, col in enumerate(X_cols):\n    print(\"%.4f\"%Fvals[i]+\" | \"+\"%.4f\"%pvals[i]+\" | \"+col)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"46b017b4380d0034ccc23e04ac5d727d5c724bfa"},"cell_type":"markdown","source":"So as expected phase is a useless feature on its own, but interestingly std_amp, median_amp, signal_std, max_amp may not be extremely useful variables because we cannot reject the null with a significance of 0.01 for these. However the features signal_mean, signal_sum, max_freq, strong_amp_count, avg_amp, and sum_amp all look like very useful features, even on their own."},{"metadata":{"trusted":true,"_uuid":"656914dcff816d875c37b686978b19a6b320aae9"},"cell_type":"code","source":"def mcc(y_true, y_pred, labels=None, sample_weight=None):\n    tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=labels, sample_weight=sample_weight).ravel()\n    mcc = (tp*tn - fp*fn)/np.sqrt((tp + fp)*(tp + fn)*(tn + fp)*(tn + fn))\n    return mcc\n\nmcc_scorer = make_scorer(mcc)\n\nlgbm_classifier = lgbm.LGBMClassifier(boosting_type='gbdt', max_depth=-1, subsample_for_bin=200000, objective=\"binary\", \n                                      class_weight=None, min_split_gain=0.0, min_child_weight=0.001, subsample=1.0, \n                                      subsample_freq=0, random_state=rand_seed, n_jobs=1, silent=True, importance_type='split')\n\nparam_distributions = {\n    \"num_leaves\": randint(16, 48),\n    \"learning_rate\": expon(),\n    \"reg_alpha\": expon(),\n    \"reg_lambda\": expon(),\n    \"colsample_bytree\": uniform(0.25, 1.0),\n    \"min_child_samples\": randint(10, 30),\n    \"n_estimators\": randint(50, 250)\n}\n\nclf = RandomizedSearchCV(lgbm_classifier, param_distributions, n_iter=100, scoring=mcc_scorer, fit_params=None, n_jobs=1, iid=True, \n                         refit=True, cv=5, verbose=1, random_state=rand_seed, error_score=-1.0, return_train_score=True)\nclf.fit(train_meta_df[X_cols], train_meta_df[\"target\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3e42e72888ed75373e2c6e6e73b56dc3d239fa0b"},"cell_type":"code","source":"print(clf.best_score_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22c6ec1dedb3bae2f3819d8924fc6d6467b0f3ec"},"cell_type":"code","source":"clf.best_estimator_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"48eed98d5f1af65aba78237d3e247551f9de893b"},"cell_type":"code","source":"fig, ax = plt.subplots()\nfig.set_size_inches(xsize, ysize)\n\nlgbm.plot_importance(clf.best_estimator_, ax=ax)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8791108fa6673562643d4913f7a05a213da7a2f7"},"cell_type":"markdown","source":"These results are interesting. The features signal_std, signal_mean, and avg_amp seem to be the most important. This makes sense intuitively because a faulty line will have more noise in its signal than a non faulty line, so for a faulty line we would expect a abnormally large signal_std, a signal_mean that is outside of the normal range due to outliers, and abnormally large avg_amp that results from lower frequencies becoming more present due to the noise from a faulty line. The next set of important features, max_amp, strong_amp_count, std_amp, and median_amp while not as important still support the current hypothesis of what the lgbm model is capturing. Finally signal_sum, sum_amp, max_freq are not important features because sum and median are robust to large outliers hence why they are not important and as determined earlier phase is not an important feature at all. \n\n\nSo the moral of this brief EDA is that we should look for features that quantify the abnormal \"noise of a signal,\" and we expect faulty lines to have large amounts of noise and not faulty lines to have low amounts of noise. Thanks for reading this kernel and good luck in detecting faulty power lines!\n\n\n*Correction*: Versions 1 and 2 of this kernel incorrectly calculated median_amp by instead calculating the std of power spectrum amplitudes (effectively creating two std_amp features). Fixing this goof doesn't change my hypothesis about what kinds of features will do well in this competition, but it does even out the distribution of feature importance a bit. If you see any other issues with the kernel let me know or if you have any questions or discussion bits let me know."},{"metadata":{"trusted":true,"_uuid":"488a4590b0005c921c2d20b31a20150218918c39"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}