{"cells":[{"metadata":{},"cell_type":"markdown","source":"# ****Volcanic Explorer****"},{"metadata":{},"cell_type":"markdown","source":"### This notebook is in its preliminary phase of data analysis. So, I will be performing some basic signal visualiztion and EDA."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import pandas  as pd\nimport numpy   as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nfrom tqdm import tqdm\nfrom pathlib import Path\n\nplt.rcParams['figure.figsize'] = (15, 10)\n\nrandom_state = 10\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So, let's begin with with some data loading."},{"metadata":{"trusted":true},"cell_type":"code","source":"train    = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\nsample   = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')\n\n#mean = train['time_to_eruption'].mean()\n\n#sample.iloc[:,1:] = mean\n#sample.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The 'train' metadeta comprises 4431 stations, each stations containing 10 sensors for collecting seismic events. Small scale earthquakes in the vicinity of volcanic area can help to guess if there is an imminent volcanic activities. "},{"metadata":{},"cell_type":"markdown","source":"### Now, lets take a look at an actual reading from one of the station."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_csvs = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train/1000015382.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train_csvs.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_csvs.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Looks like one of the sensor from this station might have missing values. "},{"metadata":{"trusted":true},"cell_type":"code","source":"train_csvs.isna().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Lets take a look at the seismograms.\nA seismogram is a graphical display of the seismic measurement by a seismograph or geophones. It capture the motion of the ground whenever a seismic wave travels through it. Usually, it is a measure of acceleration. The actual measurement consists of time in one axis and amplitude in the other axis. The ***large spikes*** in the seismogram display in the below image tells certain activity at a given point of time. Larger the amplitude, stronger the event is. The weaker or low amplitude measurements may contain a lot of noises such as ground rolls, winds, or even electric interference depending upon what kind of sensor it is used."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_csvs.plot()\nplt.show","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### From the wiggle display, it seems that the the not all the sensor have same amplitude characteristics. I would expect if all of these sensors to have similar amplitude range. It is possible that these sensors might be far apart and therefore different characteristics. Lets look at them side by side."},{"metadata":{"trusted":true},"cell_type":"code","source":"curves = train_csvs.columns\nnum_curves = len(train_csvs.columns)\n\nf, ax = plt.subplots(nrows=1, ncols = num_curves)\n\nfor ic, col in enumerate(curves):\n    if np.all(np.isnan(train_csvs[col])):\n        curve = np.empty(train_csvs[col].values.shape)\n        curve[:] = np.nan\n    else:\n        curve = train_csvs[col]\n        \n    ax[ic].plot(curve, curve.index)\n    ax[ic].set_xlabel(col)\n    ax[ic].invert_yaxis()\n    #ax[ic].set_xlim(1,60000)\n    #ax[ic].set_yticklabels([])   \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### It looks like these sensors are not located in the same geographical vicinity. As you can see some sensors dont have the spikes at the same time. If a sensor is farther from the source/epicenter of an Earthquake, the longer it takes the seismic waves to travel and reach the sensor. The image below (example from notebook by [Jasper Dramsch](https://www.kaggle.com/jesperdramsch/introduction-to-volcanology-seismograms-and-lgbm)) gives you an idea on how sensors are laid out in the field:\nDigital elevation model from Etna from ([Bonaccorso 2011]\n![Bonaccorso 2011 DEM of Etna](https://i.imgur.com/2b99LHc.jpg)\n(https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2010GC003480))\n\n\n\n\n### Looking at the seismogram display earlier, sensors 2,3 & 10 are quite similar, sensor 1, 4 & 9 are quate similar. Other sensors, especially 6 and 7 seems to pick up quite a lot of noise as compared to other sensors which responds to the stronger event nicely. \n\nTo be continued..."},{"metadata":{},"cell_type":"markdown","source":"# ****Prepare train and test set for simple learning"},{"metadata":{"trusted":true},"cell_type":"code","source":"def agg_stats(df, idx):\n    df = df.agg(['sum','min', 'mean', 'std', 'median', 'skew', 'kurtosis'])\n    df_flat = df.stack()\n    df_flat.index = df_flat.index.map('{0[1]}_{0[0]}'.format)\n    df_out = df_flat.to_frame().T\n    df_out[\"segment_id\"] = int(idx)\n    return df_out","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"summary_stats = pd.DataFrame()\nfor csv in tqdm(Path(\"../input/predict-volcanic-eruptions-ingv-oe/train/\").glob(\"**/*.csv\"), total=4501):\n    df = pd.read_csv(csv)\n    summary_stats = summary_stats.append(agg_stats(df, csv.stem))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data = pd.DataFrame()\nfor csv in tqdm(Path(\"../input/predict-volcanic-eruptions-ingv-oe/test/\").glob(\"**/*.csv\"), total=4501):\n    df = pd.read_csv(csv)\n    test_data = test_data.append(agg_stats(df, csv.stem))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features = list(summary_stats.drop([\"segment_id\"], axis=1).columns)\ntarget_name = [\"time_to_eruption\"]\nsummary_stats = summary_stats.merge(train, on=\"segment_id\")\nsummary_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"summary_stats.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Training with LGBM\nimport lightgbm as lgbm\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import roc_auc_score\nimport gc\n%%timeit\n\nn_fold = 2\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=random_state)\n\ndata = summary_stats\n\nparams = {\n    \"n_estimators\": 100,\n    \"boosting_type\": \"gbdt\",\n    \"metric\": \"mae\",\n    \"num_leaves\": 66,\n    \"learning_rate\": 0.1,\n    \"feature_fraction\": 0.9,\n    \"bagging_fraction\": 0.8,\n    \"agging_freq\": 3,\n    \"max_bins\": 2048,\n    \"verbose\": 0,\n    \"random_state\": random_state,\n    \"nthread\": -1,\n    #\"device\": \"gpu\",\n    }\n\noof_preds = np.zeros(data.shape[0])\nsub_preds = np.zeros(test_data.shape[0])\nfeature_importance = pd.DataFrame(index=list(range(n_fold)), columns=features)\n\nfor n_fold, (trn_idx, val_idx) in enumerate(folds.split(data)):\n    X_train, y_train = data[features].iloc[trn_idx], data[target_name].iloc[trn_idx]\n    X_val, y_val = data[features].iloc[val_idx], data[target_name].iloc[val_idx]\n    \n    model = lgbm.LGBMRegressor(**params)\n    \n    model.fit(X_train, y_train,  \n            eval_set= [(X_train, y_train), (X_val, y_val)], \n            eval_metric=\"mae\", verbose=0, early_stopping_rounds=150)\n    \n    feature_importance.iloc[n_fold, :] = model.feature_importances_\n    \n    oof_preds[val_idx] = model.predict(X_val, num_iteration = model.best_iteration_)\n    sub_preds += model.predict(test_data[features], num_iteration = model.best_iteration_) / folds.n_splits\n    \n    \n    print('Fold %2d AUC: %.6f' % (n_fold+1, roc_auc_score(val_y, oof_preds[val_idx])))\n    del clf, trn_x, trn_y, val_x, val_y\n    gc.collect()\n    \nprint('Full AUC score %.6f' % roc_auc_score(y,oof_preds))\n       \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best = feature_importance.mean().sort_values(ascending=False)\nbest_idx = best[best > 5].index\n\nplt.figure(figsize=(14,26))\nsns.boxplot(data=feature_importance[best_idx], orient=\"h\")\nplt.title(\"Features Importance per Fold\")\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submit Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame() \nsubmission['segment_id'] = test_data[\"segment_id\"] \nsubmission['time_to_eruption'] = sub_preds \nsubmission.to_csv('submission.csv', header = True, index = False)","execution_count":null,"outputs":[]}],"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"}},"nbformat":4,"nbformat_minor":4}