{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Volcanic experiments"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom colorama import Fore, Style","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets check how many train files we have"},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_FOLDER = '../input/predict-volcanic-eruptions-ingv-oe/train/'\nTEST_FOLDER = '../input/predict-volcanic-eruptions-ingv-oe/test/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How does a train file look like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/train.csv')\ntrain_labels = np.expand_dims(np.array(train['time_to_eruption']), axis=-1)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's have a look on how a sensor data (from a csv file). W peek a random train file."},{"metadata":{"trusted":true},"cell_type":"code","source":"random_sensor_file = pd.read_csv(os.path.join(TRAIN_FOLDER, random.choice(os.listdir(TRAIN_FOLDER))))\nrandom_sensor_file.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axis = plt.subplots(2, 5, figsize=(25, 10))\n\nfor i, sensor_name in enumerate(random_sensor_file):\n    sensor_data = np.nan_to_num(np.array(random_sensor_file[sensor_name], dtype='float'))\n    axis[i % 2, i % 5].plot(range(len(sensor_data)), sensor_data, '.')\n    axis[i % 2, i % 5].grid(True)\n    axis[i % 2, i % 5].legend([sensor_name])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So as to encrease speed of convergence, we need to scale our data. Pandas provides a useful method, which returns a lot of features, which perfectly describe the data you have."},{"metadata":{},"cell_type":"markdown","source":"### Extracting features\nHere we need to build an extractor for our noisy dataset: 60000 valuse of signals take too much memory, which is not as useful as it could be for that much memory. So we can extract some mathematic features like mean, std, variance, max, min, skew, kurtosis etc. These will give a model much more knowladge for less memory."},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_features(signal, ts, sensor_id):\n    X = pd.DataFrame()\n    f = np.fft.fft(signal)\n    f_real = np.real(f)\n    X.loc[ts, f'{sensor_id}_sum']       = signal.sum()\n    X.loc[ts, f'{sensor_id}_mean']      = signal.mean()\n    X.loc[ts, f'{sensor_id}_std']       = signal.std()\n    X.loc[ts, f'{sensor_id}_var']       = signal.var() \n    X.loc[ts, f'{sensor_id}_max']       = signal.max()\n    X.loc[ts, f'{sensor_id}_min']       = signal.min()\n    X.loc[ts, f'{sensor_id}_skew']      = signal.skew()\n    X.loc[ts, f'{sensor_id}_mad']       = signal.mad()\n    X.loc[ts, f'{sensor_id}_kurtosis']  = signal.kurtosis()\n    X.loc[ts, f'{sensor_id}_quantile99']= np.quantile(signal, 0.99)\n    X.loc[ts, f'{sensor_id}_quantile95']= np.quantile(signal, 0.95)\n    X.loc[ts, f'{sensor_id}_quantile85']= np.quantile(signal, 0.85)\n    X.loc[ts, f'{sensor_id}_quantile75']= np.quantile(signal, 0.75)\n    X.loc[ts, f'{sensor_id}_quantile55']= np.quantile(signal, 0.55)\n    X.loc[ts, f'{sensor_id}_quantile45']= np.quantile(signal, 0.45) \n    X.loc[ts, f'{sensor_id}_quantile25']= np.quantile(signal, 0.25) \n    X.loc[ts, f'{sensor_id}_quantile15']= np.quantile(signal, 0.15) \n    X.loc[ts, f'{sensor_id}_quantile05']= np.quantile(signal, 0.05)\n    X.loc[ts, f'{sensor_id}_quantile01']= np.quantile(signal, 0.01)\n    X.loc[ts, f'{sensor_id}_fft_real_mean']= f_real.mean()\n    X.loc[ts, f'{sensor_id}_fft_real_std'] = f_real.std()\n    X.loc[ts, f'{sensor_id}_fft_real_max'] = f_real.max()\n    X.loc[ts, f'{sensor_id}_fft_real_min'] = f_real.min()\n\n    return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_subtrain_items(all_items, part, shuffle=False):\n    if shuffle:\n        random.shuffle(all_items)\n    split_idx = int(len(all_items) * part)\n    all_items = all_items[:split_idx]\n    return all_items","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Split train set\nWe need to have a smaller set due to huge size of our train dataset. Huge values are useful for training a model on different examples, but we are testing different perfomances, so we train (really test) on a smaller dataset (we're taking a subset)."},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_train_ids = get_subtrain_items(train.segment_id, 0.35, True)\nprint(f'sub train items {Fore.CYAN}{len(sub_train_ids)}{Style.RESET_ALL} of {Fore.BLUE}{len(train.segment_id)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_set = list()\nj=0\nfor seg in sub_train_ids:\n    signals = pd.read_csv(f'/kaggle/input/predict-volcanic-eruptions-ingv-oe/train/{seg}.csv')\n    train_row = []\n    if j%500 == 0:\n        print(j)\n    for i in range(0, 10):\n        sensor_id = f'sensor_{i+1}'\n        train_row.append(build_features(signals[sensor_id].fillna(0), seg, sensor_id))\n    train_row = pd.concat(train_row, axis=1)\n    train_set.append(train_row)\n    j+=1\ntrain_set = pd.concat(train_set)\ntrain_set = train_set.reset_index()\ntrain_set = train_set.rename(columns={'index': 'segment_id'})\ntrain_set = pd.merge(train_set, train, on='segment_id')\ntrain = train_set.drop(['segment_id', 'time_to_eruption'], axis=1)\ny = np.expand_dims(train_set['time_to_eruption'], axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"VULCAN_INPUT_SIZE = (23,10, 1)\nNUM_EXAMPLES = train.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = np.array(y)\ndl_train_dict = {'X': np.resize(np.array(train), (NUM_EXAMPLES,) + VULCAN_INPUT_SIZE), 'y': np.expand_dims(y, axis=-1)}\nml_train_dict = {'X': np.array(train), 'y': y}\nprint(f'trained data shape (already squared for deep learning)={Fore.CYAN}',dl_train_dict['X'].shape,f'{Style.RESET_ALL}')\nprint(f'trained data shape (for machine learning)={Fore.BLUE}',ml_train_dict['X'].shape,f'{Style.RESET_ALL}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need to remove train files, which have bug sensors, which didnt detect any data"},{"metadata":{},"cell_type":"markdown","source":"let's make a specific DataGenerator of volcanic data as our RAM isn't as big as train set, if ones wants to load the whole train dataset in memory. Our DataGenerator must contain all default methods for a classic DataGenerator"},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nimport sklearn\nimport xgboost","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Validation** is very important to check on unseen data your model. So I'd like to to check my model by the next sets"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_train_test_sets(full_data, split_part):\n    assert len(full_data['X']) == len(full_data['y'])\n    split_index = int(split_part * full_data['X'].shape[0])\n    #split full data by index into test part\n    test_data = full_data['X'][:split_index]\n    test_labels = full_data['y'][:split_index]\n    #now get train data\n    train_data = full_data['X'][split_index:]\n    train_labels = full_data['y'][split_index:]\n    return {'X': train_data, 'y':train_labels}, {'X': test_data, 'y': test_labels}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"val_part = 0.1\n#-----------DEEP LEARNING--------------#\nsplit_algo_train_dict = {'ML': {'X': None, 'y':None},'DL':{'X': None, 'y':None}}\nsplit_algo_valid_dict = {'ML': {'X': None, 'y':None},'DL':{'X': None, 'y':None}}\nfor name, algo_dict in [('DL', dl_train_dict), ('ML', ml_train_dict)]:\n    split_algo_train_dict[name], split_algo_valid_dict[name] = get_train_test_sets(algo_dict, val_part)\n    print(f'{name}: there are{Fore.BLUE}', len(split_algo_train_dict[name]['X']),\n          f'{Style.RESET_ALL}train examples and{Fore.GREEN}', len(split_algo_valid_dict[name]['X']),\n          f'{Style.RESET_ALL}valid examples')\n    print('of shape', split_algo_train_dict[name]['X'][0].shape)\n    print()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Models to experiment with:\n- small DNN (~10 000 hidden neuron values) with ReLU in hidden layers. **253925043208192.0** \n- medium (~200 000 hidden neuron values) with ReLU in hidden layers. **307476037632000.0**\n- medium DNN with convolution layers (~450 000 hidden neuron values) with mixed ReLU, LeakyReLU activations **200124135374848.0**\n- RandomTreeForest **None SCORE**\n- XGBoost **None SCORE**"},{"metadata":{},"cell_type":"markdown","source":"### Deep Learning experiments"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_small_dnn_model(name):\n    model = tf.keras.models.Sequential([\n        tf.keras.layers.Flatten(input_shape=VULCAN_INPUT_SIZE),\n        tf.keras.layers.Dense(64),\n        tf.keras.layers.Dense(128),\n        tf.keras.layers.Dense(1)\n    ], name=name)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_medium_dnn_model(name):\n    model = tf.keras.models.Sequential([\n        tf.keras.layers.Flatten(input_shape=VULCAN_INPUT_SIZE),\n        tf.keras.layers.Dense(64),\n        tf.keras.layers.Dense(128),\n        tf.keras.layers.Dense(128),\n        tf.keras.layers.Dense(256),\n        tf.keras.layers.Dense(256),\n        tf.keras.layers.Dense(1)\n    ], name=name)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_medium_cnn_model(name):\n    model = tf.keras.models.Sequential([\n        tf.keras.layers.Conv2D(64, (3, 3), padding='same', input_shape=VULCAN_INPUT_SIZE),\n        tf.keras.layers.MaxPooling2D(2, 2),\n        tf.keras.layers.Conv2D(128, (3, 3), padding='same', input_shape=VULCAN_INPUT_SIZE),\n        tf.keras.layers.AveragePooling2D(2, 2),\n        tf.keras.layers.Conv2D(128, (3, 3), padding='same', input_shape=VULCAN_INPUT_SIZE),\n        tf.keras.layers.MaxPooling2D(2, 2),\n        tf.keras.layers.Flatten(),\n        tf.keras.layers.Dense(256),        \n        tf.keras.layers.Dense(512),        \n        tf.keras.layers.Dense(1)\n    ], name=name)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"models_to_try = [get_small_dnn_model('small_dnn_model'),\n                 get_medium_dnn_model('medium_dnn_model'),\n                 get_medium_cnn_model('medium_cnn_model')]\nmodels_history = []\nfor cur_model in models_to_try:\n    cur_model.compile(optimizer=tf.keras.optimizers.Adam(), loss='mse')\n    cur_model.summary()\n    cur_model_history = cur_model.fit(x=split_algo_train_dict['DL']['X'], y=split_algo_train_dict['DL']['y'],\n                                      epochs=40, verbose=0,\n                                      validation_data=(split_algo_valid_dict['DL']['X'], split_algo_valid_dict['DL']['y']))\n    models_history.append((cur_model.name, cur_model_history))\n    print(f'final loss={Fore.RED}', \n          cur_model.evaluate(split_algo_valid_dict['DL']['X'], \n                             split_algo_valid_dict['DL']['y'], verbose=0),\n          f'{Style.RESET_ALL}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"colors = ['blue', 'red', 'orange']\nlegend = []\nplt.figure(figsize=(10, 10))\nfor i, model_hist in enumerate(models_history):\n    epochs = range(len(model_hist[1].history['loss']))\n    plt.plot(epochs, model_hist[1].history['loss'], color=colors[i])\n    legend.append((model_hist[0] + \"loss\"))\nplt.legend(legend)\nplt.grid()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, neural networks **dont** come to a nice solution at all. So it's proved, that using ML algorithms for this problem is a lot better."},{"metadata":{},"cell_type":"markdown","source":"### Machine learning experiments"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\ndef get_RFR_model(args):/\n    model = RandomForestRegressor(**args)\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_XGBoost_model(args):\n    model = xgboost.XGBRegressor(**args) \n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets train on default values of our ML models, and than choose what we can use for model accurate modeling"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_squared_error","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"good_args = {\n    'max_depth': 10, \n    'n_estimators': 349, \n    'learning_rate': 0.03463499307472963, \n    'gamma': 0.2584804713489563,\n    'random_state': 666\n}\nml_models_to_try = [('random forest regressor', get_RFR_model({'n_estimators': 200})), \n                    ('XGBoost regressor', get_XGBoost_model({'n_estimators': 200})),\n                   ('XGBoost regressor', get_XGBoost_model(good_args))]\nml_models_metrics = []\nfor model_name, ml_model in ml_models_to_try:\n    print('current model:', model_name)\n    ml_model.fit(split_algo_train_dict['ML']['X'], split_algo_train_dict['ML']['y'])\n    val_prediction = ml_model.predict(split_algo_valid_dict['ML']['X'])\n    mse = mean_squared_error(val_prediction, split_algo_valid_dict['ML']['y'])\n    ml_models_metrics.append((model_name, mse))\n    print(f'final loss={Fore.RED}', mse,f'{Style.RESET_ALL}')\n    print()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_estimators = range(10, 150, 10)\nhistory = []\nfor n_estims in test_estimators:\n    print('n estimators:', n_estims)\n    ml_model = get_RFR_model(n_estims)\n    ml_model.fit(split_algo_train_dict['ML']['X'], split_algo_train_dict['ML']['y'])\n    val_prediction = ml_model.predict(split_algo_valid_dict['ML']['X'])\n    mse = mean_squared_error(val_prediction, split_algo_valid_dict['ML']['y'])\n    ml_models_metrics.append((model_name, mse))\n    print(f'final loss={Fore.RED}', mse,f'{Style.RESET_ALL}')\n    history.append(mse)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"optimized_xgb = XGBR\n'max_depth': 10, \n    'n_estimators': 349, \n    'learning_rate': 0.03463499307472963, \n    'gamma': 0.2584804713489563,\n    'random_state': 666","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(test_estimators, history)\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def predict_test_data(model, test_folder, test_csv, example_shape, x_rescale, y_rescale, verbose=1):\n    print(test_folder, test_csv)\n    test_indexes = pd.read_csv(os.path.join('../input/predict-volcanic-eruptions-ingv-oe/',test_csv))['segment_id']\n    if verbose:\n        print(f'there are{Fore.GREEN}', len(train_file_names), f'{Style.RESET_ALL}test files')    \n    predictions = {}\n    num_test_files = len(test_indexes)\n    for i, ID in enumerate(test_indexes):\n        # Store sample\n        X = np.empty((*example_shape, 1))\n        X = get_test_data(test_folder, ID, example_shape) / x_rescale\n        predictions[ID] = model.predict(X) * y_rescale\n        if verbose:\n            if i % 20 == 0:\n                print(i,'/', len(test_indexes))\n    return predictions","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Predict test data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.read_csv('../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv').head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = predict_test_data(model=current_model, test_folder=TEST_FOLDER, \n                                test_csv='sample_submission.csv',\n                                example_shape=VULCAN_INPUT_SHAPE,\n                                x_rescale=sensor_abs_max_value,\n                                y_rescale=label_abs_max_value,\n                                verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predict_frame = pd.DataFrame({'segment_id': list(predictions.keys()), 'time_to_eruption': [int(val) for val in predictions.values()]})\npredict_frame.to_csv(path_or_buf='../working/submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predict_frame.head()","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}