{"cells":[{"metadata":{},"cell_type":"markdown","source":" **Volcanic eruption Introduction**\n\n\n**Detecting volcanic eruptions before they happen is an important problem that has historically proven to be a very difficult. This competition provides you with readings from several seismic sensors around a volcano and challenges you to estimate how long it will be until the next eruption. The data represent a classic signal processing setup that has resisted traditional methods.**\n\nVolcanologists attempt to forecast volcanic eruptions, but this has proven to be nearly as difficult as predicting an earthquake. Many pieces of evidence can mean that a volcano is about to erupt, but the time and magnitude of the eruption are difficult to pin down. This evidence includes the history of previous volcanic activity, earthquakes, slope deformation, and gas emissions.\n\n# Earthquakes due to activity of volcano\n\n**Moving magma shakes the ground, so the number and size of earthquakes increases before an eruption. A volcano that is about to erupt may produce a sequence of earthquakes. Scientists use seismographs that record the length and strength of each earthquake to try to determine if an eruption is imminent. Magma and gas can push the volcanoβs slope upward. Most ground deformation is subtle and can only be detected by tiltmeters, which are instruments that measure the angle of the slope of a volcano. But ground swelling may sometimes create huge changes in the shape of a volcano. Mount St. Helens grew a bulge on its north side before its 1980 eruption. Ground swelling may also increase rockfalls and landslides.**\n\nThese earthquakes can be treated as vibrations and our sensor data has recordings of them \n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os # file handling\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 # Data Visuaization \nfrom tqdm import tqdm # Progress Bar ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_data_items = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\nprint(train_data_items.head(),'\\n\\n')\nprint(train_data_items.describe(), \"\\n\\n\")\nprint(train_data_items.info(), \"\\n\\n\")\nprint(train_data_items.dtypes, '\\n\\n')\nprint(train_data_items.count(), '\\n\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Showing the number of Segment_ids for Training**"},{"metadata":{"trusted":true},"cell_type":"code","source":"listing_training_sensors = os.listdir('../input/predict-volcanic-eruptions-ingv-oe/train')\nprint(pd.Series(listing_training_sensors).shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#### Checking if data of some training segments is missing ####\n\n\n# function to remove \".csv\" from a string\n\ndef removing_CSV(arr):\n return arr[0:len(arr)-4]\n\n##########################\n\nlisting_training_sensors = list(map(removing_CSV, listing_training_sensors))\n\nif listing_training_sensors.sort() == list(train_data_items.segment_id).sort():\n print(\"Yes we have data for all the training segment_ids\")\n \ntrain_data_items.sort_values(by = 'time_to_eruption', inplace =True) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_items.nunique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Creating a hexbin plot to check the randomness in Time to eruption"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_items.plot.hexbin(x= \"segment_id\",y =\"time_to_eruption\",gridsize=25, figsize = (10,10))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Hexbin Plot is clearly showing that the time of erruption is quite random and flat distributed so now let's dig into the data of sensors \n\n**Some information about Sensors:**\n1. generally sensors are placed in an array in different locations near the volcano\n2. these sensors are seismic sensors and it might be possible that they would have recorded noise due to external factors \n3. Last point I would like to mention is that a higher magnitude of sensor reading might account for less time remaining for net erruption\n \nOur furthur findings will be based on these points \n\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Now Looking into Data of Sensors"},{"metadata":{"trusted":true},"cell_type":"code","source":"# taking a sample file to see how it looks like \n\nSegment_1 = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/1000015382.csv\")\nSegment_1.plot(subplots = True, layout = (5,2), figsize = (20,15))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Now Let's try to Visualize the difference between values of sensors in segment ids which are having the maximum and the minimum value of time_of_eruption"},{"metadata":{"trusted":true},"cell_type":"code","source":"# segment_id_min_time = Segment id with minimum time to erruption in train.csv\n# segment_id_max_time = Segment id with maximum time to erruption in train.csv\n\n\n\nsegment_id_min_time = train_data_items[train_data_items[\"time_to_eruption\"] == min(train_data_items[\"time_to_eruption\"])]['segment_id']\nsegment_id_max_time = train_data_items[train_data_items[\"time_to_eruption\"] == max(train_data_items[\"time_to_eruption\"])]['segment_id']\nprint(f\"Segment_id with maximum time ({ max( train_data_items.time_to_eruption )}) to erruption is {int(segment_id_max_time)}\")\nprint(f\"Segment_id with maximum time ({min(train_data_items.time_to_eruption)}) to erruption is {int(segment_id_min_time)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" \nI am assuming that big differeces between the values of sensor data with segment_Id of minimum and maximum time to erruption will be seen by me πππ\n lets check this hypothesis\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sensor_max_time = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/\"+ str(int(segment_id_max_time)) +\".csv\")\nsensor_max_time.plot(subplots = True, layout = (5,2), figsize = (20,15), title = 'Plotting Sensor Data of Segment_id with maximum time to eruption')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sensor_min_time = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/\"+ str(int(segment_id_min_time)) +\".csv\")\n\nsensor_min_time.plot(subplots = True, layout = (5,2), figsize = (20,15),title = 'Plotting Sensor Data of Segment_id with maximum time to eruption')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.rcParams.update({'font.size': 10})\nsensor_min_time.plot(subplots = True, layout = (1,10), figsize = (60,5),title = 'Minimum Time Segment')\nsensor_max_time.plot(subplots = True, layout = (1,10), figsize = (60,5),title = \"Maximum Time Segment\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It can be clearly seen that these segment IDs have some missing sensors so we can not propagate this way ,\n**But the important thing is that data has some peaks, it may be noise also so let's try to clip and clean the sensors data by clipping through 99 percentile of the recorded readings**"},{"metadata":{},"cell_type":"markdown","source":"# Finding the outliers and Noise in Data Of Sensors\n\nFirst lets consider only the sensor_1 data of the segment_id with minimum time_to_erruption\nto Get the things done in proper way ,\nNOTE : we need to calculate separate percentiles for the positive and negative readings of sensors\n\n **let me show you what i mean** \n\n\nI have separated the postive and negative values for the percentile analysis and their plots are shown below\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"positive_values = sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1']\nnegative_values = sensor_min_time[sensor_min_time['sensor_1'] <= 0]['sensor_1']\nspliter_df = pd.DataFrame({'positive_values' : positive_values,\n 'negative_values' : negative_values})\nspliter_df[\"Actual_sensor_data\"] = sensor_min_time['sensor_1']\n\nspliter_df.Actual_sensor_data.plot(figsize = (20,5), ylabel = \"Actual values\",fontsize = 20,title = \"Actual values of Sensor 1\",color = 'y')\nplt.show()\nspliter_df[['negative_values','positive_values']].plot(figsize = (20,5), ylabel = \"Negative values\",fontsize = 20,title = \"Negative values of Sensor 1\",colormap = 'winter')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# finding different percentile values for our Data in sensor 1\npercent_dict = dict()\nfor i in range(1,101):\n# print(f\"{i} percentile value :\",end = ' ')\n# print(\"{0:.2f}\".format(np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],i)))\n percent_dict[i] = np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],i)\n\nplt.figure(figsize = (30,10), facecolor = 'w',edgecolor = 'g') \nplt.grid(True)\nplt.bar(list(percent_dict.keys()),list(percent_dict.values()), tick_label = list(range(1,101)), color = 'y')\nplt.show()\n\n\n# counting the number of values under each percentile\ncheck_df = pd.DataFrame()\nvalues = []\nfor i in range(1,101):\n values.append(len(sensor_min_time[sensor_min_time['sensor_1'] > percent_dict[i]]['sensor_1']))\n \ncheck_df = pd.DataFrame({'Percentile': list(range(1,101)),\n 'No of values ahead': values})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets see how many values are greater than 99 percentile\nprint(\"99 percentile : \", np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],99))\nprint(\"Total postive values are : \", sensor_min_time[sensor_min_time['sensor_1'] > 0][\"sensor_1\"].count())\nprint(\"the number of values greater than 99 percentile are :\", sensor_min_time[sensor_min_time['sensor_1'] > np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],99)]['sensor_1'].count())\nprint(\"ratio :\", (sensor_min_time[sensor_min_time['sensor_1'] > 0][\"sensor_1\"].count()/sensor_min_time[sensor_min_time['sensor_1'] > np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],99)]['sensor_1'].count())**-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 300 values can spike up the central tendency of the data a lot !!!!\n\n\nWe have assumed that these immediate spikes are due to some frequent explosions inside the active volcano and the prediction of time_of_erruption should totally depend on the continuous vibrations produced by the active volcano the standard deviation and mean values of the readings from each sensor must be\nstable and for a volcano having low value of time_of_erruption the mean and standard deviation values must be higher since high amplitude waves will be generated from a active volcano\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nprint(\"Displaying differences between successive percentile values: \\n\".upper() ) \nfor i in range(80,101):\n print(f\"{i} percentile - {i-1} percentile :\",end = ' ')\n diff = np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],i) - np.percentile(sensor_min_time[sensor_min_time['sensor_1'] > 0]['sensor_1'],i-1)\n print(\"{0:.2f}\".format(diff))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Lets write logic to clip these spike values from the sensor data "},{"metadata":{"trusted":true},"cell_type":"code","source":"def clipped(arr, after = 99): \n upper_bound = np.percentile(arr[arr > 0],after)\n lower_bound = -np.percentile(-arr[arr <= 0], after)\n return np.clip(arr, a_min = lower_bound, a_max = upper_bound)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\nLet's visualize the outputs of this clipper function on some Samples\n\n\n\nYou can clearly see the difference in data after clipping "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"arr = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/1002275321.csv\")\n\nfor i in list(arr.columns):\n if not arr[i].isnull().values.all():\n arr[i+'clipped'] = clipped(arr[i])\n else:\n arr[i+'clipped'] = arr[i]\n \ncolumn_list = ['sensor_1', 'sensor_1clipped',\n 'sensor_2', 'sensor_2clipped', 'sensor_3', 'sensor_3clipped',\n 'sensor_4', 'sensor_4clipped', 'sensor_5', 'sensor_5clipped',\n 'sensor_6', 'sensor_6clipped', 'sensor_7', 'sensor_7clipped',\n 'sensor_8', 'sensor_8clipped', 'sensor_9', 'sensor_9clipped',\n 'sensor_10', 'sensor_10clipped',]\n\narr = arr[column_list]\n\narr.plot(subplots = True, layout = (10,2),figsize =(20,60), sharex = True,sharey =True ,colormap = 'winter')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Making the Dataset Developing Function "},{"metadata":{"trusted":true},"cell_type":"code","source":"def fft_analysis(arr):\n inspection = arr\n transformed = np.fft.fft(inspection,len(inspection))\n psd = transformed*np.conjugate(transformed)/len(transformed)\n frequencies = np.arange(len(transformed))\n border = np.ones(len(frequencies)) * np.percentile(psd**2,99)\n \n index = psd < border\n rectified_fft = transformed*index\n inversed = np.fft.ifft(rectified_fft)\n return inversed.real,rectified_fft\n\n\n\ndef develop_data(arr,segment_id, appending = pd.DataFrame()):\n appending['segment_id'] = [int(segment_id)]\n arr.fillna(0)\n for i in list(arr.columns):\n if (not arr[i].isnull().values.all()) and (arr[i].any()):\n arr[i] = clipped(arr[i])\n rect_1,fft_1 = fft_analysis(arr[i])\n \n fft_1_real = pd.Series(fft_1.real)\n fft_1_imag = pd.Series(fft_1.imag)\n rect_1 = pd.Series(rect_1)\n \n \n appending['mean_' + i] = [float(arr[i].sum())]\n appending['median_' + i] = [float(arr[i].median())]\n appending['skew_' + i] = [float(arr[i].skew())]\n appending['mad_' + i] = [float(arr[i].mad())]\n appending['min_' + i] = [float(arr[i].min())]\n appending['max_' + i] = [float(arr[i].max())]\n appending['kurtosis_' + i]= [float(arr[i].kurtosis())]\n \n appending['rectified_mean_'+i] = [float(rect_1.sum())]\n appending['rectified_median_'+i] = [float(rect_1.median())]\n appending['rectified_skew_'+i] = [float(rect_1.median())]\n appending['rectified_mad_'+i] = [float(rect_1.mad())]\n appending['rectified_min_'+i] = [float(rect_1.min())]\n appending['rectified_max_'+i] = [float(rect_1.max())]\n appending['rectified_kurtosis_'+i] = [float(rect_1.kurtosis())]\n \n appending['fft_1_real_mean_'+i] = [float(fft_1_real.sum())]\n appending['fft_1_real_median_'+i] = [float(fft_1_real.median())]\n appending['fft_1_real_skew_'+i] = [float(fft_1_real.median())]\n appending['fft_1_real_mad_'+i] = [float(fft_1_real.mad())]\n appending['fft_1_real_min_'+i] = [float(fft_1_real.min())]\n appending['fft_1_real_max_'+i] = [float(fft_1_real.max())]\n appending['fft_1_real_kurtosis_'+i] = [float(fft_1_real.kurtosis())]\n\n\n appending['fft_1_imag_mean_'+i] = [float(fft_1_imag.sum())]\n appending['fft_1_imag_median_'+i] = [float(fft_1_imag.median())]\n appending['fft_1_imag_skew_'+i] = [float(fft_1_imag.median())]\n appending['fft_1_imag_mad_'+i] = [float(fft_1_imag.mad())]\n appending['fft_1_imag_min_'+i] = [float(fft_1_imag.min())]\n appending['fft_1_imag_max_'+i] = [float(fft_1_imag.max())]\n appending['fft_1_imag_kurtosis_'+i] = [float(fft_1_imag.kurtosis())]\n \n\n return appending \n\ntraining_data_frame = pd.DataFrame() \nfor i in tqdm(list(train_data_items.segment_id) , desc = \"Prepairing Training Data\" , total = len(list(train_data_items.segment_id))):\n segment_id = i\n raw = pd.read_csv(f\"../input/predict-volcanic-eruptions-ingv-oe/train/{i}.csv\")\n training_data_frame = training_data_frame.append(develop_data(raw,segment_id), ignore_index=True)\n \ntesting_data_frame = pd.DataFrame()\nfor i in tqdm(os.listdir('../input/predict-volcanic-eruptions-ingv-oe/test'), desc = \"Preparing Testing Data\", total = len(os.listdir('../input/predict-volcanic-eruptions-ingv-oe/test'))):\n segment_id = int(i[0:len(i)-4])\n raw = pd.read_csv(f\"../input/predict-volcanic-eruptions-ingv-oe/test/{i}\")\n testing_data_frame = testing_data_frame.append(develop_data(raw,segment_id), ignore_index=True)\n\ntraining_data_frame.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data_frame.to_csv(\"./training_data_frame.csv\")\ntesting_data_frame.to_csv(\"./testing_data_frame.csv\")\ntraining_data_frame.head(5)\nprint(testing_data_frame.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_data_frame = training_data_frame.fillna(0)\ntraining_data_frame = pd.merge(training_data_frame ,train_data_items , on = 'segment_id')\ntraining_data_frame","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = training_data_frame.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"testing_data_frame = testing_data_frame.fillna(0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preparing the lightGBM with MAE(Mean Absolute Error) and training the model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgbm\nfrom sklearn.model_selection import KFold\nimport gc\n\nn_fold = 7\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=101)\nfeatures = list(data.drop([\"segment_id\", \"time_to_eruption\"], axis=1).columns)\ntarget_name = [\"time_to_eruption\"]\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\": 101,\n \"nthread\": -1,\n \"device\": \"cpu\",\n}\n\nsub_preds = np.zeros(testing_data_frame.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 print(f\"Fold {n_fold}:\")\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=200, early_stopping_rounds=50\n )\n\n feature_importance.iloc[n_fold, :] = model.feature_importances_\n sub_preds += model.predict(testing_data_frame[features], num_iteration=model.best_iteration_) / folds.n_splits","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Creating Submission file"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame()\nsubmission['segment_id'] = testing_data_frame[\"segment_id\"]\nsubmission['time_to_eruption'] = sub_preds\nsubmission.to_csv('submission.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport base64\nfrom IPython.display import HTML\nimport pandas as pd\n\n# function that takes in a dataframe and creates a text link to \n# download it (will only work for files < 2MB or so)\ndef create_download_link(df, title = \"Download CSV file\", filename = \"data.csv\"): \n csv = df.to_csv()\n b64 = base64.b64encode(csv.encode())\n payload = b64.decode()\n html = '{title}'\n html = html.format(payload=payload,title=title,filename=filename)\n return HTML(html)\n\n# create a random sample dataframe\ndf = pd.DataFrame(np.random.randn(50, 4), columns=list('ABCD'))\n\n# create a link to download the dataframe\ncreate_download_link(df)\n\n# β β β Yay, download link! β β β \ncreate_download_link(submission, filename = \"submissions.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}