{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":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 gc\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nfrom sklearn.linear_model import SGDRegressor\nfrom sklearn.linear_model import Lasso\nfrom sklearn.linear_model import HuberRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nimport xgboost as xgb\n\nimport os\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n    #for filename in filenames:\n        #print(os.path.join(dirname, filename))\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n        \nfrom sklearn.decomposition import PCA\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nfrom scipy.ndimage import maximum_filter1d\nfrom scipy.ndimage import minimum_filter1d\n\nfrom datetime import datetime\n\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.preprocessing import StandardScaler\nStSc = StandardScaler()\nMMS = MinMaxScaler()\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":{},"cell_type":"markdown","source":"# Note: I have shared the TRAIN data feed notebook so that the settings / code used are visible. \n\nhttps://www.kaggle.com/davidedwards1/volcano-train-fts-gen-v1/notebook\n\n# This takes a while to run. If you have any ideas to make it faster and to get better features, feel free to share..."},{"metadata":{},"cell_type":"markdown","source":"# This is a rough first attempt based on summary data for each train dataframe\n# #So far, lower cv = lower LB, when the CV has reduced significantly. BUT the cv remains consistently lower than Lb\n# \n# Comments\n# Created a ton of features (rolling stats, summaries of rolling stats)\n# So far, more features generally works better\n# # I'm loading these from another workbook/data feed as they take a while to calculate\n# Ive not spent tons of time fine tuning. just some basic directional testing\n# Huber/SGD etc regression dont seem to work\n# Tree models seem to work better at this stage\n# I've chosen XGB as it has GPU option and can handle missing values\n# Running multiple seeds\n# \n# Version updates\n# Replaced zeros with np.nan (missing)\n# Removed scaling (not needed for xgb)\n# Testing dropping features with a lot more NAs in test data set to see if this reduces CV absolute error compared to LB absolute error\n# #I have added a correlation feed - this quite simply uses the corr() function on each set of data to provide a table (10 x 10) of correlations between each sensor for the whole 60001 rows. this did not help.\n# I tried a stronger cutoff for features in test missing values. this did not help.\n\n#added feature importance analysis at the end"},{"metadata":{"trusted":true},"cell_type":"code","source":"V_PATH = '/kaggle/input/predict-volcanic-eruptions-ingv-oe/'\nTRAIN_PATH = V_PATH + 'train/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SENSOR_COLS = ['sensor_1', 'sensor_2', 'sensor_3', 'sensor_4', 'sensor_5', 'sensor_6',\n       'sensor_7', 'sensor_8', 'sensor_9', 'sensor_10']\n\n\nSENSOR_RMEANS = [x+'_rmin' for x in SENSOR_COLS] \nSENSOR_RSTDS = [x+'_rstd' for x in SENSOR_COLS] \nSENSOR_RMINS = [x+'_rmin' for x in SENSOR_COLS] \nSENSOR_RMAXES = [x+'_rmax' for x in SENSOR_COLS]\nSENSOR_RGRADMEAN = [x+'_grad_rmean' for x in SENSOR_COLS]\nSENSOR_RGRADSTD = [x+'_grad_rstd' for x in SENSOR_COLS]\n\nSENSOR_RSTATS = [SENSOR_RMEANS, SENSOR_RSTDS, SENSOR_RMINS, SENSOR_RMAXES,\n               SENSOR_RGRADMEAN,  SENSOR_RGRADSTD]\n\nROLL_DESCR = ['rmin', 'rstd', 'rmin', 'rmax', 'grad_rmean','grad_rstd']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(V_PATH+'train.csv')\nprint(train.shape)\nprint(train.columns)\ntrain.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"total_rows_estimate = 60001 * len(train) / 1000000\nprint('estimate of total TRAIN rows (millions)',total_rows_estimate)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv(V_PATH+'sample_submission.csv')\nprint(sample_submission.shape)\nprint(sample_submission.columns)\n\ntotal_rows_estimate = 60001 * len(sample_submission) / 1000000\nprint('estimate of total TEST rows (millions)',total_rows_estimate)\n\nsample_submission.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#examine the distribution of time until eruption\n\nsns.kdeplot(train['time_to_eruption'] / 1000000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train['time_to_eruption'].min(), train['time_to_eruption'].max(), train['time_to_eruption'].mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sz = train['time_to_eruption'].size-1\ntrain['PCNT_TIME'] = train['time_to_eruption'].rank(method='max').apply(lambda x: 1.0*(x-1)/sz)\ntrain.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_rolling(df, cols, window=50):\n    for col in cols:\n        df[col+'_grad'] = np.gradient(df[col])\n        df[col+'_grad'] = df[col+'_grad'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_grad_abs'] = np.gradient(np.abs(df[col]))\n        df[col+'_grad'] = df[col+'_grad'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_rmin'] = minimum_filter1d(df[col].values, size=window)\n        df[col+'_rmax'] = maximum_filter1d(df[col].values, size=window)\n        \n        df[col+'_rmin'] = df[col+'_rmin'].fillna(method='bfill').fillna(method='ffill')\n        df[col+'_rmax'] = df[col+'_rmax'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_rmean'] = df[col].rolling(window=window, center=True).mean().fillna(method='bfill').fillna(method='ffill')\n        df[col+'_rstd'] = df[col].rolling(window=window, center=True).std().fillna(method='bfill').fillna(method='ffill')\n        \n        #add also for gradients\n        df[col+'_grad_rmin'] = minimum_filter1d(df[col+'_grad_abs'].values, size=window)\n        df[col+'_grad_rmax'] = maximum_filter1d(df[col+'_grad_abs'].values, size=window)\n        \n        df[col+'_grad_rmin'] = df[col+'_grad_rmin'].fillna(method='bfill').fillna(method='ffill')\n        df[col+'_grad_rmax'] = df[col+'_grad_rmax'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_grad_rmean'] = df[col+'_grad_abs'].rolling(window=window, center=True).mean().fillna(method='bfill').fillna(method='ffill')\n        df[col+'_grad_rstd'] = df[col+'_grad_abs'].rolling(window=window, center=True).std().fillna(method='bfill').fillna(method='ffill')\n        \n    return df\n\ndef get_stats(df, sensor_cols=SENSOR_COLS, rolling_cols=SENSOR_RSTATS):\n    #we create the min max etc of original sensor columns\n    df['max'] = df[sensor_cols].max(axis=1)\n    df['min'] = df[sensor_cols].min(axis=1)\n    df['std'] = df[sensor_cols].std(axis=1)\n    \n    #and with absolute values\n    df['max_abs'] = np.abs(df[sensor_cols]).max(axis=1)\n    df['min_abs'] = np.abs(df[sensor_cols]).min(axis=1)\n    df['std_abs'] = np.abs(df[sensor_cols]).std(axis=1)\n    \n    #we take mins and maxes of groups of rolling columns\n    for count,rc in enumerate(rolling_cols): #this takes a SINGLE mean, max across each GROUP of rolling\n        #columns - e.g. the max of all rolling mins\n        df[ROLL_DESCR[count]+'_max'] = df[rolling_cols[count]].max(axis=1)\n        df[ROLL_DESCR[count]+'_min'] = df[rolling_cols[count]].min(axis=1)\n        df[ROLL_DESCR[count]+'_std'] = df[rolling_cols[count]].std(axis=1)\n        df[ROLL_DESCR[count]+'_mean'] = df[rolling_cols[count]].mean(axis=1)    \n   \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#lets drop the rolling mean - does not seem that useful\n\n\ndef get_all_stats(df, cols, rolling_cols, window=50):\n    \n    df = get_rolling(df, cols, window=window)\n    df = get_stats(df, sensor_cols=cols, rolling_cols=rolling_cols)\n    df = df.groupby(['segment'])[[x for x in df.columns if x != 'segment']].agg(['mean',\n                                                                                'max','min','std'])\n    df.columns=[a+b for a,b in df.columns]\n    return df\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loaded_dfs = pd.read_csv('/kaggle/input/volcano-train-fts/volcano_train_fts.csv',index_col=0)\nprint(loaded_dfs.shape)\nloaded_dfs.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_dfs = pd.read_csv('/kaggle/input/volcano-test-features/volcano_test_fts.csv',index_col=0)\nprint(test_dfs.shape)\ntest_dfs.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NON_FTS = ['time_to_eruption', 'segment']\nLABEL = 'time_to_eruptionmean'\n\nREGRESSION_FTS = [x for x in loaded_dfs.columns if 'time_to_eruption' not in x]\nREGRESSION_FTS = [x for x in REGRESSION_FTS if 'segment' not in x]\nprint('Number of features,', len(REGRESSION_FTS))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#it looks from the EDA like some sensors may be quite reasonably correlated\n#lets try to create some features by examining differences between stats of sensor 2 and sensor 4\n#these sensors (based on some limited sample data) looked much better correlated close to eruptions\n\ns1 = 'sensor_2'\ns2 = 'sensor_4'\n\ns1_feats = [x for x in REGRESSION_FTS if s1 in x]\ns2_feats = [x for x in REGRESSION_FTS if s2 in x]\n\nprint(s1_feats[0:10])\nprint(s2_feats[0:10])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for sd1, sd2 in zip(s1_feats, s2_feats):\n    loaded_dfs[sd1+'_delta_'+s2] = loaded_dfs[sd1] - loaded_dfs[sd2]\n    test_dfs[sd1+'_delta_'+s2] = test_dfs[sd1] - test_dfs[sd2]\n    \n    REGRESSION_FTS+=[sd1+'_delta_'+s2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loaded_dfs['time_to_eruptionmean'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"LABEL = 'time_to_eruptionmean'\n\nsns.kdeplot(loaded_dfs[LABEL])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loaded_dfs = loaded_dfs.fillna(value=0)\ntest_dfs = test_dfs.fillna(value=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"regression_importance = pd.Series(index=REGRESSION_FTS, data=0.0)\nfrom sklearn.feature_selection import mutual_info_regression\nfrom sklearn.feature_selection import f_regression\n#for RF in REGRESSION_FTS:\nregression_importance[:] = f_regression(loaded_dfs[REGRESSION_FTS], loaded_dfs[LABEL])[0]\n\nCUTOFF = regression_importance.quantile(0.25)\nprint('Number of features over cutoff', sum(regression_importance>CUTOFF))\nSEL_FTS = regression_importance.index[regression_importance>CUTOFF]\n\nsns.kdeplot(regression_importance)\nregression_importance.sort_values(ascending=False).head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FILL_ZEROS=True\nif FILL_ZEROS==True:\n    loaded_dfs[SEL_FTS] = loaded_dfs[SEL_FTS].replace({0: np.nan})\n    test_dfs[SEL_FTS] = test_dfs[SEL_FTS].replace({0: np.nan})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"na_analysis = pd.DataFrame(index=SEL_FTS,\n                          data=0.0, columns=['Train', 'Test'])\n\nna_analysis['Train'] = loaded_dfs[SEL_FTS].isna().sum().values / len(train)\nna_analysis['Test'] = test_dfs[SEL_FTS].isna().sum().values / len(test_dfs)\nna_analysis['Delta'] = na_analysis['Test'] - na_analysis['Train']\n\nfig,axes=plt.subplots(figsize=(10,4))\nsns.kdeplot(na_analysis['Train'], color='Green')\nsns.kdeplot(na_analysis['Test'], color='Red')\nsns.kdeplot(na_analysis['Delta'], color='Blue')\naxes.set_title('Distribution of Zeros/NAs')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"na_analysis['Delta'].sort_values()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"drop_fts = [x for x in na_analysis[na_analysis['Delta']>0.2].index]\nprint(len(drop_fts))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"drop_fts","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SEL_FTS = [x for x in SEL_FTS if x not in drop_fts]\nprint(len(SEL_FTS))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loaded_dfs['label_strat'] = np.round(loaded_dfs[LABEL] * 20, 0)\nloaded_dfs['label_strat'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import StratifiedKFold\nNFOLDS=10\nskf5 = StratifiedKFold(n_splits=NFOLDS)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loaded_dfs.columns[~loaded_dfs.columns.isin(test_dfs.columns)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#run xgb with multiple seeds and gpu support\n\nbaseline_error = mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                  np.full((len(loaded_dfs),), loaded_dfs['time_to_eruptionmean'].mean()))\n\nprint('baseline error', baseline_error)\n\npredictions = np.zeros((len(loaded_dfs),))\ntest_predictions = np.zeros((len(test_dfs),))\n\nft_imps=pd.Series(index=SEL_FTS,\n                  data=0.0)\n\nrslist=range(20)\nrs_errors=[]\n\nfor count1, RS in enumerate(rslist):\n    xgbr = xgb.XGBRegressor(random_state=RS,\n                           tree_method='gpu_hist' ,\n                            colsample_bytree=0.5,\n                            reg_alpha=0.1,\n                            missing =np.nan,\n                            subsample=0.75\n                       )\n    models = [xgbr]\n    for count,mod in enumerate(models):\n        #print(mod)\n\n        for trn_idx, val_idx in skf5.split(loaded_dfs[SEL_FTS], loaded_dfs['label_strat']):\n            print('run fold')\n            mod.fit(loaded_dfs.loc[trn_idx, SEL_FTS].values, \n                      loaded_dfs.loc[trn_idx,'time_to_eruptionmean'].values)\n\n            predictions[val_idx] +=mod.predict(loaded_dfs.loc[val_idx, SEL_FTS].values)\n\n            print('Fold val Error',mean_absolute_error(loaded_dfs['time_to_eruptionmean'][val_idx],\n                          predictions[val_idx]/((count+1) * (count1+1))))\n\n            test_predictions += mod.predict(test_dfs[SEL_FTS].values)\n            \n            ft_imps+=xgbr.feature_importances_\n\n        print('Error end of model run',mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                          predictions/((count+1) * (count1+1))))\n        \n    rs_errors+=[mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                          predictions/((count+1) * (count1+1)))]\n\npredictions = predictions/(len(models)* len(rslist))\ntest_predictions = test_predictions / (len(models)*NFOLDS * len(rslist))\n\npredictions = np.where(predictions<0, 0, predictions)\ntest_predictions = np.where(test_predictions<0, 0, test_predictions)\n\nprint(mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                      predictions))\n\nprint('Scaled CV error',mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                      predictions) * train['time_to_eruption'].max())\n\n\nfig,axes=plt.subplots(nrows=1,ncols=2,figsize=(18,6))\naxes[0].scatter(x=loaded_dfs['time_to_eruptionmean'],\n           y=predictions, color='Red')\nsns.lineplot(x=range(len(rs_errors)),\n           y=np.array(rs_errors), ax=axes[1])\n\naxes[0].set_title('CV predictions vs actual time to eruption')\naxes[1].set_title('CV error vs random seed cycle')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Scaled CV error',mean_absolute_error(loaded_dfs['time_to_eruptionmean'],\n                      predictions) * train['time_to_eruption'].max())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.kdeplot(predictions, color='Green')\nsns.kdeplot(test_predictions, color='Red')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# let's check feature importance\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.kdeplot(ft_imps, color='Green')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(figsize=(8,20))\nft_imps=ft_imps.sort_values(ascending=False)\naxes.barh(y=ft_imps.index[0:20], width=ft_imps[0:20])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(test_predictions), len(sample_submission))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission['time_to_eruption'] = test_predictions * train['time_to_eruption'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.kdeplot(train['time_to_eruption'], color='Green')\nsns.kdeplot(sample_submission['time_to_eruption'], color='Red')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(datetime.now())","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}