{"cells":[{"metadata":{},"cell_type":"markdown","source":"# What is Volcano Forecasting?\nThis is an introduction to understanding the data and the problem.\n\n**Volcanoes are awesome!** both in the California way and the \"instilling awe\" way.\n\n![Mount Etna CC-BY Kuhnmi Flickr](https://i.imgur.com/TZUE3ht.jpg)\nMount Etna on a calm day. (CC-BY Kuhnmi [Flickr](https://www.flickr.com/photos/31176607@N05/25720452925))\n\nVolcano monitoring is important for both the inhabitants on and next to volcanoes, but also globally, as seen with the [Eyjafjallajökull eruption](https://en.wikipedia.org/wiki/2010_eruptions_of_Eyjafjallaj%C3%B6kull) disrupting air travel ten years ago. Geophysics is the field that largely works with active volcanoes and their activity, measuring earthquakes, tiltmeters, changes in gravimetry, etc. Specifically, seismologists record the rumble of the Earth when the magma forces its way upward. Geochemists collect data on degassing on Volcanos. And geologists are looking how the rocks form from lava. \n\nIn the spirit of the competition, we should not find additional metadata, however, it's important to understand the kind of data we are dealing with. A Nature article ([Hall 2018](https://www.nature.com/articles/d41586-018-07420-y#ref-CR1)) describes the world's first automatic volcano forecast system on Mount Etna. The data for this is mostly of acoustic nature, specifically infrasound ([Ripepe et al. 2018](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018JB015561)). Seeing that we have 10 sensors with timeseries, it's pretty safe to make the assumption that we are dealing with a seismological problem. These tend to be the most reliable, as compared to e.g. gas analyzers. Gas analyzers have to be right on top of a vent (which tend to rebuild in the lifetime of volcanoes), whereas the seismological stations just have to listen and not run out of battery.\n\nSo let's have a look at the actual data!\n\n# Data\nFiles\ntrain.csv Metadata for the train files.\n\n`segment_id`: ID code for the data segment. Matches the name of the associated data file.\n\n`time_to_eruption`: The target value, the time until the next eruption.\n\n`[train|test]/*.csv`: the data files. Each file contains ten minutes of logs from ten different sensors arrayed around a volcano. The readings have been normalized within each segment, in part to ensure that the readings fall within the range of int16 values. If you are using the Pandas library you may find that you still need to load the data as float32 due to the presence of some nulls."},{"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\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt # plots\nimport seaborn as sns\n\nfrom pathlib import Path\nfrom tqdm import tqdm\n\nrandom_state = 42\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train.csv\")\ntrain.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sequence = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/1000015382.csv\", dtype=\"Int16\")\nsequence.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sequence.tail()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Each sequence is 10 minutes long with 600001 samples. The data is `int16` but contains nan's, luckily with Pandas 1.0 they introduced the [nullable integer datatype](https://pandas.pydata.org/pandas-docs/stable/user_guide/integer_na.html), just make sure to actually call it capitalized `\"Int16\"` so Pandas knows."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"sequence.fillna(0).plot(subplots=True, figsize=(25, 10))\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These stations clearly have similar data, but are shifted in time. Why?\n\nLook at this digital elevation model from Etna from ([Bonaccorso 2011](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2010GC003480))\n\n![Bonaccorso 2011 DEM of Etna](https://i.imgur.com/2b99LHc.jpg)\n\nThe different stations are located all around the crater. That means the further away a station is from the magma or hypocentre of an Earthquake or other accoustic signals like streets, or [humans](https://www.weforum.org/agenda/2020/07/seismic-anthropogenic-noise-lockdown-covid19/), the longer the acoustic wave has to travel, arriving later at the station. So we can assume the with regards to that one signel event above, station 9 and 10 are almost identical in how far they are from it. Sensor 8 is a bit further.\n\nAdditionally, you can see the different noise levels of the data. Sensor 5 is fantastically quiet, only responding to big events, while Sensor 6 has some extremely periodic noise on it. It'd be worth to investigate if this is always the case.\n\nOverall, the sensors will have different lags, depending where the volcanic activity is happening. So when we build a model that uses the time series, this is something to keep in mind.\n\n# Play With the Data\nTime series are hard. Can we get away with initially playing with some aggregate statistics?"},{"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":{},"cell_type":"markdown","source":"# Train a LightGBM Regressor\n\nUse [Cross Validation](https://scikit-learn.org/stable/modules/cross_validation.html), because if you don't [shakeup](https://www.kaggle.com/jtrotman/meta-kaggle-competition-shake-up) will not be your friend."},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgbm\nfrom sklearn.model_selection import KFold\nimport gc\n\n\nn_fold = 7\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=random_state)\n\ndata = summary_stats\n\nparams = {\n    \"n_estimators\": 2000,\n    \"boosting_type\": \"gbdt\",\n    \"metric\": \"mae\",\n    \"num_leaves\": 66,\n    \"learning_rate\": 0.005,\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\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    trn_x, trn_y = data[features].iloc[trn_idx], data[target_name].iloc[trn_idx]\n    val_x, val_y = data[features].iloc[val_idx], data[target_name].iloc[val_idx]\n    \n    model = lgbm.LGBMRegressor(**params)\n    \n    model.fit(trn_x, trn_y, \n            eval_set= [(trn_x, trn_y), (val_x, val_y)], \n            eval_metric=\"mae\", verbose=0, early_stopping_rounds=150\n           )\n\n    feature_importance.iloc[n_fold, :] = model.feature_importances_\n    \n    sub_preds += model.predict(test_data[features], num_iteration=model.best_iteration_) / folds.n_splits\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":"So that gives us a nice idea which data actually matters for the next iteration.\n\n# Submit Prediction\nLet's build a csv for submission."},{"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":{},"cell_type":"markdown","source":"# How Far Can One Go?\nNeed Inspiration?\n\nOf course someone tried Transformers on Earthquake time series, to detect the event and different phases. I love multi tasking like that!\n![Nature Paper on Earthquake Transformers](https://i.imgur.com/KvRcjLh.png)\nSeismogram tagging on full sequences with an Earthquake Transformer in Nature ([Mousavi 2020](https://www.nature.com/articles/s41467-020-17591-w))."}],"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}