{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Imports and loading data"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom scipy import stats, signal, fft\nfrom sklearn.linear_model import Lasso\nfrom sklearn.svm import LinearSVR\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.model_selection import RepeatedKFold\nfrom sklearn.compose import TransformedTargetRegressor\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom multiprocessing import Pool\nfrom astropy.stats import biweight_location, biweight_scale\nimport random\nimport os","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Locking seeds"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(1)\nrandom.seed(1)\nos.environ['PYTHONHASHSEED'] = '0'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading the data"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/train.csv')\ntest = pd.read_csv('/kaggle/input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature Extraction"},{"metadata":{},"cell_type":"markdown","source":"Functions to build our training and test sets"},{"metadata":{"trusted":true},"cell_type":"code","source":"def df_parallelize_run(func, t_split, cores=4):\n    \"\"\"\n    inspired by https://www.kaggle.com/kyakovlev/m5-three-shades-of-dark-darker-magic\n    \"\"\"\n    num_cores = np.min([cores, len(t_split)])\n    pool = Pool(processes=num_cores)\n    df = pd.concat(pool.starmap(func, t_split), axis=0)\n    pool.close()\n    pool.join()\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_basic_features(df, direct='train'):\n    \n    result = pd.DataFrame(dtype=np.float32)\n    result['segment_id'] = df['segment_id']\n    result['time_to_eruption'] = df['time_to_eruption']\n    result.set_index('segment_id', inplace=True)\n    \n    for ids in df['segment_id']:\n        \n        f = pd.read_csv(f'/kaggle/input/predict-volcanic-eruptions-ingv-oe/{direct}/{ids}.csv')\n        f['sensor_all'] = f.sum(1)                          # Adding another read that will be the sum of all sensors\n                                                            # Or you can use other statistics as well (min, max, etc)\n        for sensor in f.columns:\n            \n            s_ = f[sensor].ffill().fillna(0).values         # filling nans\n            raw_s = s_                                      # raw signal\n            filtered_s = signal.medfilt(s_, 5)              # median filtered signal with window=5\n            cos_s = fft.dct(s_)                             # direct cosine transformed signal\n\n            for num, data in enumerate([raw_s, cos_s, filtered_s]):\n                \n                mean = np.mean(data)                        # mean\n                std = np.std(data, ddof=0)                  # std\n                min_ = np.min(data)                         # minimum\n                max_ = np.max(data)                         # maximum\n                mr = (min_ + max_)*0.5                      # midrange\n                ptp = max_ - min_                           # range\n                q1, q3 = np.quantile(data, q=[0.25, 0.75])  # 1st & 3rd quartile\n                mh = (q3 + q1)*0.5                          # midhinge\n                iqr = q3 - q1                               # IQR\n                mad_1 = np.median(np.abs(data - mean))      # Median deviation from the mean\n                mad_2 = np.mean(np.abs(data - mean))        # Mean deviation from the mean\n                bwl = biweight_location(data)               # Biweight location\n                bws = biweight_scale(data)                  # Biweight scale\n                \n                result.loc[ids, f'{sensor}_mean_{num}'] = mean\n                result.loc[ids, f'{sensor}_std_{num}'] = std\n                result.loc[ids, f'{sensor}_midrange_{num}'] = mr\n                result.loc[ids, f'{sensor}_midhinge_{num}'] = mh\n                result.loc[ids, f'{sensor}_min_{num}'] = min_\n                result.loc[ids, f'{sensor}_max_{num}'] = max_\n                result.loc[ids, f'{sensor}_ptp_{num}'] = ptp\n                result.loc[ids, f'{sensor}_q1_{num}'] = q1\n                result.loc[ids, f'{sensor}_q3_{num}'] = q3\n                result.loc[ids, f'{sensor}_iqr_{num}'] = iqr\n                result.loc[ids, f'{sensor}_med_abs_dev_mean_{num}'] = mad_1\n                result.loc[ids, f'{sensor}_mean_abs_dev_mean_{num}'] = mad_2\n                result.loc[ids, f'{sensor}_biweight_location_{num}'] = bwl\n                result.loc[ids, f'{sensor}_biweight_scale_{num}'] = bws\n            \n    return result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What else could be added here?\n\n1. Add different modification of signal: abs signal, FFT, DST, MFCC/Cepstrum, different filters with different windows, STA/LTA (1s and 60s default, with 100Hz sampling frequency it will be 100 and 6000 samples), 1st/2nd derivative and their abs, 1st/2nd percentage change and their abs, cumulative and window statistics;\n2. Add other statistics to compute. I removed a lot just because I do not want to impute NaN's that will emerge, but these features also could give additional information. Examples will be kurtosis, skewness, CV, informational entropy. You can also add other TS features to the mix (from *tsfresh* or *tslearn* libs);\n3. Watch out for computation times, number of features and collinearity! Use L1/L2 regulariazations (and feature subsampling if you are planing to use DT-based models). Blindly adding additional features and signals will lead to bloated dataset which will eat your RAM alive. Do not throw everything in the mixer."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_CORES = 4\nchunk_size = train.shape[0] // N_CORES\nsplits_train = [train[:chunk_size], train[chunk_size:2*chunk_size], train[2*chunk_size:3*chunk_size], train[3*chunk_size:]]\nsplits_test = [test[:chunk_size], test[chunk_size:2*chunk_size], test[2*chunk_size:3*chunk_size], test[3*chunk_size:]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\ntrain_df = df_parallelize_run(get_basic_features, \n                              [(chunk, 'train') for chunk in splits_train],\n                              N_CORES)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\ntest_df = df_parallelize_run(get_basic_features, \n                             [(chunk, 'test') for chunk in splits_test],\n                             N_CORES)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It will take a good amount of time, so do not overbear the CPU with all the features you want to compute. Or compute them on your local machine (if it has more cores) and then upload *train* and *test* dataframes to Kaggle for reuse."},{"metadata":{},"cell_type":"markdown","source":"Some sanity checks"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Features containing nans in train are {train_df.columns[(train_df.replace((np.inf, -np.inf), np.nan).isna().sum() > 0)].tolist()}')\nprint(f'Features containing nans in test are {test_df.columns[(test_df.replace((np.inf, -np.inf), np.nan).isna().sum() > 0)].tolist()}')\nprint(f'Low variation features in train are {train_df.columns[train_df.var() <= 1].tolist()}')\nprint(f'Low variation features in test are {test_df.columns[test_df.var() <= 1].tolist()}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train_df.drop('time_to_eruption', axis=1)\nX_test = test_df.drop('time_to_eruption', axis=1)\ny = train_df['time_to_eruption']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will transform our target to:\n1. Reduce the range of possible values;\n2. Get an \"easier\" to predict distribution than uniform."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20,6));\naxes[0].set_title('Raw target');\ny.hist(bins='fd', density=True, ax=axes[0]);\naxes[1].set_title('Target with sqrt transform');\nnp.sqrt(y).hist(bins='fd', density=True, ax=axes[1]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will also check correlation for different targets"},{"metadata":{"trusted":true},"cell_type":"code","source":"X.corrwith(y).abs().sort_values(ascending=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X.corrwith(np.sqrt(y)).abs().sort_values(ascending=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since we will use linear model correlation is an OK way to look at usefulness of our features. For rule-based models (like Decision Trees and all the variations) is better to use \n> method='spearman'\n\nor\n\n> method='kendall'\n\noptinos in *corrwith* method. They are rank-based correlation and thus a bit better suited for rule-based models."},{"metadata":{},"cell_type":"markdown","source":"And a check for distribution shift between train and test features with 2 sample two-sided Kolmogorov-Smirnov test."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfor col in X.columns:\n    p = stats.ks_2samp(X[col], X_test[col], mode='exact', alternative='two-sided')[1]\n    if p < 0.05:\n        print(col)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Modeling"},{"metadata":{},"cell_type":"markdown","source":"First of all we will use Lasso to select features.\n\nWe will need *StandardScaler* to scale our data, *Lasso* to get coefs, *TargetTransformer* to transform our targets and cross-validation to tune hyperparameters."},{"metadata":{},"cell_type":"markdown","source":"We will use repeated K fold validation scheme with 3 splits and 4 repetitions. Number of splits were chosen based on factorization of number of rows in train: 4431 = 3 ∙ 7 ∙ 211. Number of reprtitions were chosen to get number of folds divisible by 4 (number of cores) to get (hopefully) the most out of parallel computations."},{"metadata":{"trusted":true},"cell_type":"code","source":"rkf = RepeatedKFold(n_splits=3, n_repeats=4, random_state=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I am already using good hyperparameters for my notebook, but with different features it could vary."},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler_lasso = StandardScaler()\nlasso = Lasso(alpha=2.5, fit_intercept=True, normalize=False, max_iter=20_000, tol=1e-3, positive=False, random_state=1, selection='random')\ntt_lasso = TransformedTargetRegressor(regressor=lasso, func=np.sqrt, inverse_func=np.square)\npipe_lasso = Pipeline([('scaler', scaler_lasso), ('reg', tt_lasso)])\ncv_lasso = cross_val_score(pipe_lasso, \n                           X, \n                           y, \n                           scoring='neg_mean_absolute_error', \n                           cv=rkf, \n                           n_jobs=4, \n                           verbose=True, \n                           error_score='raise')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Mean MAE is {(cv_lasso*-1).mean()}')\nprint(f'STD MAE is {(cv_lasso*-1).std()}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Fitting our pipeline and getting coefficients for feature selection."},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe_lasso.fit(X, y)\nlasso_coefs = pipe_lasso.named_steps['reg'].regressor_.coef_\nimportant_cols = X.columns[lasso_coefs != 0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now our main model. It will be *LinearSVR* with target transformation."},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler_svr = StandardScaler()\nlin_svr = LinearSVR(epsilon=0, tol=1e-2, C=333, loss='epsilon_insensitive', fit_intercept=True, intercept_scaling=1, dual=True, verbose=0, random_state=1, max_iter=1000)\ntt_svr = TransformedTargetRegressor(regressor=lin_svr, func=np.sqrt, inverse_func=np.square)\npipe_svr = Pipeline([('scaler', scaler_svr), ('reg', tt_svr)])\ncv_svr = cross_val_score(pipe_svr, \n                         X[important_cols], \n                         y, \n                         scoring='neg_mean_absolute_error', \n                         cv=rkf, \n                         n_jobs=4, \n                         verbose=True, \n                         error_score='raise')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Mean MAE is {(cv_svr*-1).mean()}')\nprint(f'STD MAE is {(cv_svr*-1).std()}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pipe_svr.fit(X[important_cols], y)\npredictions = pipe_svr.predict(X_test[important_cols])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What else could be done here?\n\n1. You can use boosting. I personally do not recommend it. GBDT can only interpolate, it could not extrapolate. So, if test set has target values that are out of range of train target values, you won't be able to predict them properly. Linear models or NNs can. You can use blend of GBDT and LR, this will be the best combination. Another option will be XGBoost with boosting='gblinear';\n2. You can use different scalers and different target transform functions (log, 1/x, cbrt);\n3. You can predict hours, minutes, seconds and mseconds, then blend results together. Or you can split target into hours-minutes-seconds-ms, predict them and then sum results."},{"metadata":{},"cell_type":"markdown","source":"# Submission"},{"metadata":{},"cell_type":"markdown","source":"We won't doing any postprocessing here, but the options would be floor/ceil/round, clip to get rid of negative or/and too big values."},{"metadata":{"trusted":true},"cell_type":"code","source":"test['time_to_eruption'] = predictions\ntest","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.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}