{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\n\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\n\nfrom matplotlib import pyplot as plt\n\nfrom sklearn.base import TransformerMixin\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.model_selection import train_test_split, KFold\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.inspection import permutation_importance\nfrom sklearn.feature_selection import chi2, SelectKBest\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.metrics import mean_absolute_error\n\nimport optuna\n\nfrom xgboost import XGBRegressor","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Задача состоит в том, чтобы по поведению временного ряда в 10-минутном интервале предсказать величину `time_to_eruption`.\nПри этом не требуется декомпозировать ряд или прогнозировать его изменение.\nВ таком случае важно вычислять такие характеристики, как наибольший выброс, коэффициенты эксцесса, ассиметрии и так далее."},{"metadata":{"trusted":true},"cell_type":"code","source":"data_folder = Path('../input/predict-volcanic-eruptions-ingv-oe/')\n\ndf = pd.read_csv(data_folder / 'train.csv')\n\ndf.head().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Рассмотрим, как выглядит временной ряд в 10-минутном промежутке.\nВыберем наблюдения для промежутков с наименьшим и наибольшим доступными показателями `time_to_eruption`."},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df.sort_values('time_to_eruption')\n\ntime_min = df.head(1).iloc[0]['segment_id']\ntime_max = df.tail(1).iloc[0]['segment_id']\n\ndf_min = pd.read_csv(data_folder / 'train' / f'{time_min}.csv')\ndf_max = pd.read_csv(data_folder / 'train' / f'{time_max}.csv')\n\ndf_min['time_to_eruption'] = df.head(1).iloc[0]['time_to_eruption']\ndf_max['time_to_eruption'] = df.tail(1).iloc[0]['time_to_eruption']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_min.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_max.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Заметим, что у ряда с минимальным `time_to_eruption` не работают 3 датчика, а у ряда с максимальным значением всего 2.\nМожем использовать число отсутствующих датчиков как признак при предсказании."},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(10):\n    sensor = f'sensor_{i + 1}'\n    if df_min[sensor].isnull().all() or df_max[sensor].isnull().all():\n        del df_min[sensor]\n        del df_max[sensor]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_sensors_data(ld, rd):\n    figure, axs = plt.subplots(5, 2, figsize=(16, 20))\n\n    for (i, c) in zip(range(5), df_min.columns):\n        axs[i, 0].plot(ld[c])\n        axs[i, 0].set_title(f'Minimal time to eruption, {c}')\n    \n        axs[i, 1].plot(rd[c])\n        axs[i, 1].set_title(f'Maximal time to eruption, {c}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sensors_data(df_min, df_max)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Видно, что при малом `time_to_eruption` наблюдается всплеск около `28000` на всех датчиках, особенно массивный на последних двух.\nНа `sensor_7` также есть заметный всплеск около `48000`.\n\nУменьшим масштаб и посмотрим, как ведет себя ряд на малых промежутках времени."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sensors_data(df_min[:100], df_max[:100])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Нетрудно заметить, что частоты отличаются.\nДля извлечения этой и других характеристик используем библиотеку `tsfresh`. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# !pip install tsfresh\n\nfrom tsfresh import extract_features\nfrom tsfresh.feature_extraction import MinimalFCParameters","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tsfresh_parameters = MinimalFCParameters()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Так как данных много, то вручную настраиваем параметры, которые хотим извлечь."},{"metadata":{"trusted":true},"cell_type":"code","source":"del tsfresh_parameters['length']\n\ntsfresh_parameters['skewness'] = None\ntsfresh_parameters['kurtosis'] = None\ntsfresh_parameters['last_location_of_maximum'] = None\ntsfresh_parameters['first_location_of_maximum'] = None\ntsfresh_parameters['last_location_of_minimum'] = None\ntsfresh_parameters['first_location_of_minimum'] = None\ntsfresh_parameters['first_location_of_minimum'] = None\ntsfresh_parameters['benford_correlation'] = None\ntsfresh_parameters['percentage_of_reoccurring_values_to_all_values'] = None\ntsfresh_parameters['percentage_of_reoccurring_datapoints_to_all_datapoints'] = None\n\ntsfresh_parameters['number_peaks'] =  [\n    {'n': 1}, \n    {'n': 3}, \n    {'n': 5}, \n    {'n': 10}, \n    {'n': 50}\n]\ntsfresh_parameters['binned_entropy']  = [\n    {'max_bins': 10}\n]\ntsfresh_parameters['fft_aggregated']  = [\n    {'aggtype': 'centroid'},\n    {'aggtype': 'variance'},\n    {'aggtype': 'skew'},\n    {'aggtype': 'kurtosis'}\n]\ntsfresh_parameters['autocorrelation'] = [\n    {'lag': 0},\n    {'lag': 1},\n    {'lag': 2},\n    {'lag': 3},\n    {'lag': 4},\n    {'lag': 5},\n    {'lag': 6},\n    {'lag': 7},\n    {'lag': 8},\n    {'lag': 9}\n]\ntsfresh_parameters['agg_autocorrelation'] = [\n    {'f_agg': 'mean', 'maxlag': 40},\n    {'f_agg': 'median', 'maxlag': 40},\n    {'f_agg': 'var', 'maxlag': 40}\n]\ntsfresh_parameters['friedrich_coefficients'] = [\n    {'coeff': 0, 'm': 3, 'r': 30},\n    {'coeff': 1, 'm': 3, 'r': 30},\n    {'coeff': 2, 'm': 3, 'r': 30},\n    {'coeff': 3, 'm': 3, 'r': 30}\n]\ntsfresh_parameters['count_above'] = [{'t': 0}]\ntsfresh_parameters['count_below'] = [{'t': 0}]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Читаем данные всех рядов."},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_timeseries(data, parameters, is_train=True):\n    df_result = None\n    \n    if is_train:\n        segments = df.iterrows()\n    else:\n        segments = enumerate(os.listdir(data_folder / 'test'))\n        \n    for idx, row in segments:\n        if is_train:\n            segment, time_to_eruption = row\n            segment_timeseries_path = data_folder / 'train/{segment}.csv'\n        else:\n            segment = row\n            segment_timeseries_path = data_folder / 'test/{segment}'\n\n        segment_timeseries = pd.read_csv(segment_timeseries_path)\n        segment_timeseries = segment_timeseries.fillna(0).reset_index()\n        segment_timeseries['id'] = idx\n        \n        # Извлекаем признаки из временного ряда\n        extracted_features = extract_features(segment_timeseries, \\\n                                              column_id='id', \\\n                                              column_sort='index', \\\n                                              disable_progressbar=True,\n                                              default_fc_parameters=parameters)\n        extracted_features['segment'] = segment\n        \n        if is_train:\n            extracted_features['time_to_eruption'] = time_to_eruption\n\n        if df_result is None:\n            df_result = extracted_features\n        else:\n            df_result = pd.concat([df_result, extracted_features], \\\n                                  axis=0, \\\n                                  ignore_index=True, \\\n                                  sort=True)\n            \n        print(f'Processed segment #{idx}')\n        \n    return df_result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Сохраняем полученные датафреймы, чтобы не пересчитывать вновь. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def save(parameters):\n    ts_train = pd.read_csv(data_folder / 'train.csv')\n    df_train = preprocess_timeseries(ts_train, parameters)\n    df_train.to_csv('train.csv', index=False)\n    \n    df_test = preprocess_timeseries(None, parameters, is_train=False)\n    df_test.to_csv('test.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Мы загрузили посчитанные значения в отдельный `kaggle`-датасет. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/volcanoes-ts-processed-with-tsfresh/train.csv')\ndf_train = df_train.dropna(axis='columns')\n\nfeatures = [c for c in df_train.columns if c not in ['time_to_eruption', 'segment']]\n\nX = df_train[features]\ny = df_train['time_to_eruption']\n\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Используем случайный лес, чтобы отсеять признаки с низкой информативностью."},{"metadata":{"trusted":true},"cell_type":"code","source":"class LowImportanceSelector(TransformerMixin):\n    def __init__(self, threshold, n_estimators=100):\n        self.features = None\n        self.threshold = threshold\n        self.n_estimators = n_estimators\n\n    def fit(self, X, y):\n        estimator = RandomForestRegressor(n_estimators=self.n_estimators)\n        estimator.fit(X, y)\n    \n        importances = pd.DataFrame({\n            'feature': X.columns,\n            'importance': estimator.feature_importances_\n        })\n        importances = importances[importances['importance'] > self.threshold]\n        \n        self.features = importances['feature']\n        \n        return self\n\n    def transform(self, X):\n        return X[self.features]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Удаляем признаки с высоким коэффициентом корреляции."},{"metadata":{"trusted":true},"cell_type":"code","source":"class CorrelationSelector(TransformerMixin):\n    def __init__(self, threshold):\n        self.columns = None\n        self.threshold = threshold\n   \n    def fit(self, X, y=None):\n        X = X.copy()\n        self.columns = set()\n        C = X.corr()\n        for i in range(len(C.columns)):\n            for j in range(i):\n                if (C.iloc[i, j] >= self.threshold) and (C.columns[j] not in self.columns):\n                    c = C.columns[i]\n                    self.columns.add(c)\n                    if c in X.columns:\n                        del X[c]\n        \n        return self\n\n    def transform(self, X, y=None):\n        X.drop(columns=list(self.columns)).shape\n        \n        return X.drop(columns=list(self.columns))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def objective(trial, data=X, target=y):\n    parameters = {\n        'tree_method': 'gpu_hist',\n        'lambda': trial.suggest_loguniform('lambda', 1e-3, 10.0),\n        'alpha': trial.suggest_loguniform('alpha', 1e-3, 10.0),\n        'colsample_bytree': trial.suggest_categorical('colsample_bytree', [0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,1.0]),\n        'subsample': trial.suggest_categorical('subsample', [0.4, 0.5, 0.6, 0.7, 0.8, 1.0]),\n        'learning_rate': trial.suggest_categorical('learning_rate', [0.008, 0.009, 0.01, 0.012, 0.014, 0.016, 0.018, 0.02]),\n        'n_estimators': 1000,\n        'max_depth': trial.suggest_categorical('max_depth', [5, 7, 9, 11, 13, 15, 17, 20]),\n        'random_state': trial.suggest_categorical('random_state', [24, 48, 2020]),\n        'min_child_weight': trial.suggest_int('min_child_weight', 1, 300),\n    }\n    \n    X_train, X_test, y_train, y_test = train_test_split(data, target, test_size=0.2, random_state=1337)\n    \n    model = XGBRegressor(**parameters)\n    model.fit(X_train, y_train, eval_set=[(X_test, y_test)], early_stopping_rounds=100, verbose=False)\n    \n    return mean_absolute_error(y_test, model.predict(X_test))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\nstudy = optuna.create_study(direction='minimize')\nstudy.optimize(objective, n_trials=100)\n\nprint(f'Number of finished trials: {len(study.trials)}')\nprint(f'Best trial: {study.best_trial.params}')\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"parameters = {\n    'lambda': 0.0020555245431348778, \n    'alpha': 0.11298627316540845, \n    'colsample_bytree': 0.6, \n    'subsample': 1.0, \n    'learning_rate': 0.01, \n    'max_depth': 20, \n    'random_state': 48, \n    'min_child_weight': 18\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe = Pipeline([\n    ('correlation', CorrelationSelector(threshold=0.85)),\n    # ('importance', LowImportanceSelector(threshold=1e-4)),\n    ('scaler', MinMaxScaler()),\n    ('xgboost', XGBRegressor(objective='reg:squarederror', n_estimators=1000, **parameters))\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1337)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe.score(X_test, y_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe.fit(X, y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_predict = pd.read_csv('../input/volcanoes-ts-processed-with-tsfresh/predict.csv')\ndf_predict = df_predict.dropna(axis='columns')\ndf_predict['segment'] = df_predict['segment'].apply(lambda s: s.split('.')[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"target = pipe.predict(df_predict.drop(columns='segment'))\ntarget = pd.DataFrame({\n    'segment_id': df_predict['segment'], 'time_to_eruption': target\n})\ntarget.to_csv('submission.csv', 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}