{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# LANL: Review\n***\nIn May/June 2019, I completed my first Kaggle competition - the Earthquake Prediction challenge hosted by Los Alamos National Laboratory (LANL). It was a great chance to take some rather freshly-acquired Python and ML knowledge out for a spin. Admitedly, I did it in a bit of a rush as I joined the competition late. Also, I didn't really take the time to look back and document the lessons learnt (and make sense of all the mess I created as a newbie). \n\nSo 8 months later, here I am... (I am also yet to complete another competition and thought this would be a good way to get back into it!)\n\n## What I did ... in a nutshell\n\nFor me, this competition was mostly about feature engineering - I knew almost nothing about seismology and acoustic data processing, but at least I had a kitchen sink full of statistical functions to throw at the data. I submitted a very simple ML model - an Elastic Net with minimal hyperparameter tuning. I did play around with a neural network, but couldn't really improve performance all that much and the deadline was looming...\n\nI did get in the top 10% and nabbed a Bronze medal, jumping up 2875 spots all the way up to 443 when the leaderboard was finalised with all the test data. Maybe the fact I joined late worked in my favour as I would have probably ended up overfitting my model if I had more time.\n\nAnyway, here's a review of what I did (quite a bit more readable than my actual submision here: https://www.kaggle.com/slashie/lanl-final), plus a little bonus XGBoost exploration at the end, because I've always been itching to know if XGBoost could've really *boosted* my score :)\n\n(I also placed some of the more detailed code in this utility script here https://www.kaggle.com/slashie/lanl-udf to keep this notebook relatively tidy)"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Bread and butter\nimport os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nplt.style.use('seaborn-darkgrid')\n%matplotlib inline\n\n# ML modelling\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.linear_model import LinearRegression, ElasticNet\nfrom sklearn.metrics import mean_absolute_error\nimport xgboost as xgb\n\n# Utility module: my user-defined functions | see: https://www.kaggle.com/slashie/lanl-udf\n#import lanl_udf\nfrom lanl_udf import *\n#help(lanl_udf)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Exploring the data\n***\nFirst, let's see what's living in the data directory (I've noticed the Kaggle directory setup has changed slightly since last time...):"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = '/kaggle/input/LANL-Earthquake-Prediction'\nprint(os.listdir(data_dir)) # let's see what's in the data directory!","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Inspecting training data\n\nThe training dataset in `'train.csv'` is very large and should be extracted in batches/chunks.\n\n*I still remember trying to import the whole training dataset in one go and running out of memory! Luckily `pandas.read_csv()` has pretty good documentation and so I eventually figured out to limit the number of rows extracted! #noobmoment*"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_obs = 150000 # number of obs to extract in one go (corresponds with how test data is set up - see test data section below)\nn_skip = 4 * (10 ** 5) # number of rows to skip (i.e. to be able to look at different data sections, not just one!)\ntrain_path = os.path.join(data_dir,'train.csv')\nsample = pd.read_csv(train_path, \n                     nrows=n_obs, header=None, skiprows=n_skip) # header set to None, else values will be set as column names when skipping\nsample.columns = ['acoustic_data','time_to_failure']\nsample.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At first glance of the dataframe head, we see some variation in `acoustic_data`, but none in `time_to_failure`. \n\nWhen we plot it out below, we can see that `acoustic_data` is quite volatile, whereas the `time_to_failure` decreases in a consistent step-wise linear manner. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=2, facecolor='white', figsize=(14,7))\nsample.acoustic_data.plot(linewidth=0.5, ax=ax[0])\nax[0].set_ylabel('Acoustic Data',fontsize=12)\nsample.time_to_failure.plot(linewidth=1.5, ax=ax[1])\nax[1].set_ylabel('Time to Failure (seconds)',fontsize=12)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Inspecting test files\nBy looking at files in the `test` folder, it becomes clear that training data should be split into chunks of 150000. This way the training data is chunked into files that correspond in size to the test data files. For each of these test files, we eventually have to generate a *single* prediction of `time_to_failure` using the patterns in the `acoustic_data`.\n\nGiven we have to generate only one prediction for each segment in the `test` folder, that means we don't need to retain all the `time_to_failure` values from each training data segment as a training y-variable. We just need the last `time_to_failure` value.\n\n*When I was doing this the first time, it took me a while to realise this, but I eventually understood after looking through some of the public notebooks for tips (they are a great learning resource for beginners)... #noobmoment*"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_path = os.path.join(data_dir,'test')\ntest_files = os.listdir(test_path)\ntest_file_num = 28\ntest_sample = pd.read_csv(os.path.join(data_dir,'test',test_files[test_file_num]))\nprint(\"There are %d rows in each test file\" %(test_sample.shape[0]))\ndisplay(test_sample.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training Data Setup\n***\n### Feature Generating Functions\n\nMy general approach to feature engineering was to summarise the `acoustic_data` values in a given 150000 segment using statistical features such as *autocorrelation*, *volatility*, *concentration* etc. I also attempted to read up on frequency measures from the field signal processing and found an approach to compute the sine-based frequency of a series using `numpy.fft()`.\n\nMy user-defined functions are in the `lanl_udf.py` module, which is imported at the beggining. Using my functions, as well as some `scipy` summary statistics functions, I generated the following 14 features in each 150000-long segment of `acoustic_data`:\n1. *Mean*\n2. *Standard deviation*\n3. *Autocorrelation (first order)*\n4. *Log of skewness squared*\n5. *Log of kurtosis*\n6. *Autocorrelation (first order) of first differences*\n7. *Mean of absolute deviations*\n8. *Geometric mean of absolute deviations*\n9. *Harmonic mean of absolute deviations*\n10. *Fraction of sum of absolute deviations in top 500 observations*\n11. *Fraction of sum of absolute deviations in top 25000 observations*\n12. *Fraction of absolute deviations above 750*\n13. *Fraction of observations equal to the mode*\n14. *Wave frequency measure*\n\nWith most of these features I tried to capitalise on the fact that the autocorrelation of `acoustic_data` tended to fall and volatility (as well as extreme values) would spike as `time_to_failure` approached 0 (as will be shown shortly). I experimented with a lot more than this, but many other features were often almost perfectly correlated with one of the 14 features above and thus redundant. \n\n### Training Data Extraction\nTo increase the number of observations in my training data, I extracted overlapping (rather than mutually exclusive) training segments. To implement this, I started iterating through the training dataset at different points, using the `skiprows` argument of `pandas.read_csv()`. Therefore, I iterated through the training data multiple times (5 to be exact), each time starting at a different point, so that the segments are different. This is implemented below by the `gen_training_data()` function from the `lanl_udf.py` module. "},{"metadata":{"trusted":true},"cell_type":"code","source":"X_names = ['mean','stdev','AC(1)','log(skew^2)','log(kurt)','AC(1)_diff',\n           'mean(abs_dev)','gmean(abs_dev)','hmean(abs_dev)', 'frac_top500', \n           'frac_top25000', 'frac_dev>750', 'frac_eq_mode', 'wave_freq']\ny_name = 'time_to_failure'\ntry:\n    try:\n        df_train = pd.read_csv('df_train.csv')\n    except:\n        df_train = pd.read_csv('/kaggle/input/lanl-review/df_train.csv')\n    X_train = df_train[X_names].values\n    y_train = df_train[y_name].values\nexcept:\n    X_train, y_train = gen_training_data(train_path)\n    df_train = pd.DataFrame(X_train)\n    df_train.columns = X_names\n    df_train.loc[:,y_name] = y_train\ndf_train.to_csv('df_train.csv', index=False)\ndf_train.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visual Analysis\n***\n### Correlation of each feature with the target\n\nAs I trialled various features, I looked at each of their individual correlations with the target in the search of something that would provide some predictive power. "},{"metadata":{"trusted":true},"cell_type":"code","source":"corr_Xy = [np.corrcoef(y_train,X_train[:,i])[0,1] for i in range(X_train.shape[1])]\ndf = pd.DataFrame({'feature':X_names, 'corr_w_target':corr_Xy}).set_index('feature')\nfig = plt.figure(facecolor='white', figsize=(14,7))\ndf['corr_w_target'].plot(kind='bar', fontsize=12)\nplt.xlabel(None)\nplt.xticks(rotation=45)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Feature Correologram\n\nI was also worried about multi-collinearity across my features, so I regularly checked the feature correologram as I added new features. In other words, I wanted each new feature I experiemented with to bring *new information* to the table i.e. it would not be too strongly correlated with other features (but still correlated somewhat with the y-variable). \n\nTo be honest, looking at some of the correlations in the figure below, I'm not sure if I succeeded. What I should have done is to extract as many features as possible and apply Principal Component Analysis to extract the key components of variation. Instead, I just went ahead and hoped some form of regularization would take care of redundant features later on..."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_df = df_train[X_names]\nfig = plt.figure(facecolor='white', figsize=(12,10))\nhm = sns.heatmap(X_df.corr(), cmap='viridis')\nhm.tick_params(labelsize=12)\nplt.xticks(rotation=45)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## The Model\n***\n### Evaluating model with train-test splits\n\nInstead of setting up a cross-validation K-fold type environment, I went for a less elegant approach to test the model's predictive stability. I basically ran it multiple times with different `random_state` arguments in `train_test_split()`. I also did it manually back then (not even using the for-loop approach below), picking whichever number came to my head (not the best, I know). The things I was looking out for:\n* The R squared\n* The Mean Absolute Error - which was the competition's evaluation metric\n* Parameter Estimates - must be my econometrics background, but I really cared about those! \n\nAnyway, I manually tried different linear models with different hyperparameters, dabbled unsuccessfully with a neural network, but in the end, my *super-scientfic* methods led me to the model below!\n\nThe `ElasticNet` was quite good because it muted the large coefficient estimates on some of the parameters that I got when I ran a simple `LinearRegression`. This likely happened because of multi-collinearity i.e. very high correlation among some of my features, as can be seen in the correlogram above. Principal Components Analysis (PCA) could have helped, but I went ahead with the `ElasticNet` as a quicker solution given the time pressure. You can see the `0` coeffcient values where the `ElasticNet` did its job and thereby improved predictive power by making the model more generalizable. "},{"metadata":{"trusted":true},"cell_type":"code","source":"for state in [7, 42, 88, 101]:\n    print(\"Random state is %d\"%state,\"\\n\",\"-\"*30)\n    X_fit, X_eval, y_fit, y_eval = train_test_split(X_train, y_train, test_size = 0.3, random_state=state)\n    steps = [('scaler', StandardScaler()),\n            ('reg', ElasticNet(alpha=0.01))]\n    pipeline = Pipeline(steps)\n    pipeline.fit(X_fit, y_fit)\n    y_pred = pipeline.predict(X_eval)\n    print(\"R^2: {}\".format(pipeline.score(X_eval, y_eval)))\n    MAE = mean_absolute_error(y_eval,y_pred)\n    print(\"Mean Absolute Error: {}\".format(MAE))\n    print(pipeline.steps[1][1].coef_,'\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Training the full model\n\nThis was the best I could do, so I proceeded to fit the model on all the training data, taking a peek at the fitted parameters to make sure they didn't change in an unexpectedly crazy way. They did not."},{"metadata":{"trusted":true},"cell_type":"code","source":"pipeline.fit(X_train, y_train)\nprint(pipeline.steps[1][1].coef_)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Making and Submitting Predictions\n***\n### Generating test data\n\nHere, the `gen_test_data()` function from my utility module `lanl_udf` is applied to extract the same set of features as above from the `acoustic_data` in each of the test segments. This works very similarly to the `gen_train_data()` function, except that, being *unseen* test data, we  don't have any `time_to_failure` y-data to extract. Instead, we need to keep track of the name of each test segment file using `seg_id`, so that we have an identifier for the Kaggle submission file."},{"metadata":{"trusted":true},"cell_type":"code","source":"seg_id, X_test = gen_test_data(test_path)\nprint(\"Generated features for %d test segments\"%len(seg_id))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Predictions\n\nNow we can use the model to generate predictions of `time_to_failure` for each test segment and construct the Kaggle `submission.csv` file. I also adjusted any negative predictions to equal 0, as negative `time_to_failure` does not make sense."},{"metadata":{"trusted":true},"cell_type":"code","source":"elnet_pred = pipeline.predict(X_test)\nelnet_pred[elnet_pred<0] = 0\nelnet_submit_df = pd.DataFrame({'seg_id': seg_id, 'time_to_failure': elnet_pred})\n#elnet_submit_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Trying out XGBoost\n\nWith this competition, I felt like more than 90% of my time went into the feature engineering part and I didn't get to play around with other ML models - particularly `XGBoost`, which is what all the cool kids use. I always wondered whether using `XGBoost` could have gotten me a much better score, holding my training data and feature set fixed. Well... let's find out!\n\n***\n### Quick parameter tune-up\n\nI'm not going to spend a lot of time fine-tuning here since the comeptition is long done, but I still want to use something that's reasonable. The parameters I will *roughly* fine-tune are: `max_depth`, `eta` (learning rate) and `num_boost_round`, using the MAE score from the test sample across 4 cross-validation folds. "},{"metadata":{"trusted":true},"cell_type":"code","source":"xgb_train = xgb.DMatrix(data=X_train, label=y_train)\nall_results = {'MAE-test':[], 'MAE-train':[], 'max_depth':[], 'eta':[], 'num_boost_round':[]}\nfor max_depth in [10, 12, 14]:\n    for eta in [0.2, 0.25, 0.3]:\n        for num_boost in [9, 11, 13]:\n            params = {\"objective\":\"reg:squarederror\", \"max_depth\":max_depth, \"eta\":eta}\n            cv_results = xgb.cv(dtrain=xgb_train, params=params, nfold=4, num_boost_round=num_boost, metrics=\"mae\", seed=42)\n            all_results['max_depth'].append(max_depth)\n            all_results['eta'].append(eta)\n            all_results['num_boost_round'].append(num_boost)\n            all_results['MAE-test'].append(cv_results['test-mae-mean'].values[-1])\n            all_results['MAE-train'].append(cv_results['train-mae-mean'].values[-1])\nall_results_df = pd.DataFrame(all_results).sort_values(by='MAE-test')\nall_results_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Hyperparameter choice and prediction\n\nI am not going to use the top 2 models from above, because I am worried they are over-fitted given the significantly lower training MAE. To me, the third and fourth models seem most resonable and they are quite similar. I am actually going to go with the fourth model, because it uses fewer boosting rounds and thus potentially less prone to overfitting.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {\"objective\":\"reg:squarederror\", \"max_depth\":10, \"eta\":0.25}\nxg_reg = xgb.train(params=params, dtrain=xgb_train, num_boost_round=11)\nxgb_test = xgb.DMatrix(data=X_test, label=np.zeros([X_test.shape[0],]))\nxgb_pred = xg_reg.predict(xgb_test)\nxgb_submit_df = pd.DataFrame({'seg_id': seg_id, 'time_to_failure': xgb_pred})\nxgb_submit_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### ElasticNet vs XGBoost\n\nOut of curiosity, I wanted to see whether the predictions of the `XGBoost` model are that different from the `ElasticNet` model. At this point, I have not submitted the new XGBoost predictions yet, but, given the similarity of predictions, my initial instinct is that `XGBoost` would not have gotten me any prize money :) A lot more work would've had to be done! "},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Correlation between XGBoost and ElasticNet predictions is %.3f\"%np.corrcoef(xgb_pred, elnet_pred)[0,1])\nfig = plt.figure(facecolor='white', figsize=(12,9))\nplt.scatter(xgb_pred, elnet_pred)\nplt.xlabel('XGBoost', fontsize=12)\nplt.ylabel('ElasticNet', fontsize=12)\nplt.show()","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":1}