{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Attempt at Genetic Programming with GPLearn"},{"metadata":{},"cell_type":"markdown","source":"#### Vettejeep"},{"metadata":{},"cell_type":"markdown","source":"### Introduction"},{"metadata":{},"cell_type":"markdown","source":"This is an attempt at genetic programming with gplearn in Python.  It is a learning exercise, not a serious attempt at the leader board.  The leader board score was 1.561, so there is a long way to go before I have enough knowledge of genetic programming to consider my efforts with gplearn being competitive with boosted trees.  Still, 1.561 seems to be a reasonable start and clearly indicates that the alorithm is learning.  Hopefully this can serve as a community learning exercise for interested Kagglers.  Maybe we can improve it together.  All sugestion for improvement is requested and appreciated.  Please feel free to fork and re-post with your improvements. "},{"metadata":{},"cell_type":"markdown","source":"Switched this over to [\"Andrews Features\"](https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples).  Still dicey as whether it will train in 9 hours, so the training code is commented out.  I do not have a leader board score with this yet."},{"metadata":{},"cell_type":"markdown","source":"### Notes"},{"metadata":{},"cell_type":"markdown","source":"The code was originally run with a feature set generated by code from my Masters Degree project.  Please see my [Masters Final Project – Model; LB~1.392](https://www.kaggle.com/vettejeep/masters-final-project-model-lb-1-392) for that feature generation code.  This notebook is not a part of my degree project, it is just an interest area."},{"metadata":{},"cell_type":"markdown","source":"The gplearn library is at https://gplearn.readthedocs.io/en/stable/.  It is a simple pip install & appears to be on Kaggle already."},{"metadata":{},"cell_type":"markdown","source":"### Impressions of gplearn"},{"metadata":{},"cell_type":"markdown","source":"Trains slowly but seems to consume few resources.  Locally I trained three instances at once without a problem."},{"metadata":{},"cell_type":"markdown","source":"Lacks built-in 'tanh' function, but it is easy to make our own.  Tanh has helped toy regression models that I have tried out."},{"metadata":{},"cell_type":"markdown","source":"The gplearn implementation of early stopping seems useless, so the call to fit() is implemented in a loop with warm start and a simple early stop based upon a failure to improve metric."},{"metadata":{},"cell_type":"markdown","source":"Have no clue how to match the quality of results obtained by Scirpus in genetic programming!  Tutorials for gplearn appear to be sadly lacking on the web.  It is not so clear on how to best change hyperparameters."},{"metadata":{},"cell_type":"markdown","source":"### Coding Attempt"},{"metadata":{},"cell_type":"markdown","source":"Get gplearn, install if needed.  Kaggle appears to have it installed."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":199},"colab_type":"code","executionInfo":{"elapsed":4052,"status":"ok","timestamp":1556718983617,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"vB-S33R1Vaa4","outputId":"48ea5d38-5689-4510-d497-3a69528cb6e1","trusted":true},"cell_type":"code","source":"# !pip install gplearn\nimport gplearn\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Imports."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":35},"colab_type":"code","executionInfo":{"elapsed":927,"status":"ok","timestamp":1556718987647,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"HU3nR6bfVk42","outputId":"2eb0a47c-5c64-4407-cd10-b267a19e601e","trusted":true},"cell_type":"code","source":"import os\nimport shutil\nimport copy\nimport numpy as np\nimport pandas as pd\nfrom sklearn.metrics import mean_absolute_error\nfrom gplearn.functions import make_function\nfrom gplearn.genetic import SymbolicRegressor\nfrom sklearn.model_selection import KFold\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir('../..'))\nprint(os.listdir('../input'))\nprint(os.listdir('../input/gppyfiles'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load data files."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":35},"colab_type":"code","executionInfo":{"elapsed":14010,"status":"ok","timestamp":1556719053894,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"cvI7cvq_YNU6","outputId":"820c7938-7fa3-4c3f-8fd5-3e66e93f31a3","trusted":true},"cell_type":"code","source":"DATA_DIR = '../input/gplearn-data'\nsubmission = pd.read_csv(os.path.join('../input/LANL-Earthquake-Prediction', 'sample_submission.csv'), index_col='seg_id')\nscaled_train_X = pd.read_csv(os.path.join(DATA_DIR, 'scaled_train_X_AF0.csv'))\nscaled_test_X = pd.read_csv(os.path.join(DATA_DIR, 'scaled_test_X_AF0.csv'))\ntrain_y = pd.read_csv(os.path.join(DATA_DIR, 'train_y_AF0.csv'))\npredictions = np.zeros(len(scaled_test_X))\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Optional: Reduce feature set by eliminating apparently useless features using Pearsons."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":92},"colab_type":"code","executionInfo":{"elapsed":947,"status":"ok","timestamp":1556719101049,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"McKU-WcQjIAP","outputId":"d71f3fb3-04b6-4126-d7d1-6016604f96be","trusted":true},"cell_type":"code","source":"from scipy.stats import pearsonr\ny = train_y['time_to_failure'].values\npcol = []\npcor = []\npval = []\nfor col in scaled_train_X.columns:\n    pcol.append(col)\n    pcor.append(abs(pearsonr(scaled_train_X[col], y)[0]))\n    pval.append(abs(pearsonr(scaled_train_X[col], y)[1]))\n\ndf = pd.DataFrame(data={'col': pcol, 'cor': pcor, 'pval': pval}, index=range(len(pcol)))\ndf.sort_values(by=['cor', 'pval'], inplace=True)\ndf.dropna(inplace=True)\ndf = df.loc[df['pval'] <= 0.05]\n\ndrop_cols = []\n\nfor col in scaled_train_X.columns:\n    if col not in df['col'].tolist():\n        drop_cols.append(col)\n\nscaled_train_X.drop(labels=drop_cols, axis=1, inplace=True)\nscaled_test_X.drop(labels=drop_cols, axis=1, inplace=True)\nprint(scaled_train_X.shape, scaled_test_X.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Make a 'tanh' function for gplearn.  I suspect it needs to be in a top level module to pickle properly in the notebook, so create as a separate file (gplearn_tanh.py) and import. Here is code for the import file. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# in an IDE this works as a lambda - not so in colab - have not tried a lambda on kaggle\nfrom gplearn.functions import make_function\n\ndef th(x):\n    return np.tanh(x)\n\ngptanh = make_function(th, 'tanh', 1)\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import our tanh function and create a custom function set for our gplearn session that includes it."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":316},"colab_type":"code","executionInfo":{"elapsed":260,"status":"error","timestamp":1556719171272,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"KqFzIPecmVb0","outputId":"838a3f8d-898d-41aa-bc42-9fc4bc2bfdc0","trusted":true},"cell_type":"code","source":"# import os\n# import sys\n# sys.path.append('../input/gppyfiles')\n# from gppyfiles.gplearn_tanh import gptanh\nfunction_set = ['add', 'sub', 'mul', 'div', 'inv', 'abs', 'neg', 'max', 'min', gptanh]  # 'sqrt', 'log', \nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Set up the folds for cross validation."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":244},"colab_type":"code","executionInfo":{"elapsed":413,"status":"error","timestamp":1556718969982,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"P7iHBQ0CbeVn","outputId":"43301b21-7c08-42fb-87f0-b7d4251daf59","trusted":true},"cell_type":"code","source":"predictions = np.zeros(len(scaled_test_X))\n\nn_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\n\nfold_importance_df = pd.DataFrame()\nfold_importance_df[\"Feature\"] = scaled_train_X.columns\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Optional sample weighting.  I tried to up-weight the samples above 10 seconds because they are uncommon, in the hope the model would fit them better."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":126},"colab_type":"code","executionInfo":{"elapsed":278,"status":"ok","timestamp":1556685210866,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"n8qo8w3Ok-x_","outputId":"0f23ec2f-fae1-422d-d126-f70d1486b2ca","trusted":true},"cell_type":"code","source":"sample_wts = np.sqrt(np.array([x - 10.0 if x > 10.0 else 0 for x in y]) + 1.0)\nprint(y[0:16])\nprint(sample_wts[-8:])\nprint('ok')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Genetic program model, main code loop."},{"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":72635},"colab_type":"code","executionInfo":{"elapsed":3911518,"status":"ok","timestamp":1556674131729,"user":{"displayName":"Kevin Maher","photoUrl":"","userId":"00562981053856066430"},"user_tz":360},"id":"Tg3X7MQ0aIk2","outputId":"ed8bd188-a766-4fe2-81b2-3ff25c59bc48","trusted":true},"cell_type":"code","source":"# still needs cleanup\n# GENS = 500\n# MAE_THRESH = 2.5\n# MAX_NO_IMPROVE = 50\n# np.random.seed(666)\n# maes = []\n# gens = []\n\n# for fold_, (trn_idx, val_idx) in enumerate(folds.split(scaled_train_X, train_y.values)):\n#   print('working fold %d' % fold_)\n#   X_tr, X_val = scaled_train_X.iloc[trn_idx], scaled_train_X.iloc[val_idx]\n#   y_tr, y_val = train_y['time_to_failure'].values[trn_idx].ravel(), train_y['time_to_failure'].values[val_idx].ravel()\n#   sample_wts_tr = sample_wts[trn_idx]\n#   np.random.seed(5591 + fold_)\n#   best = 1e10\n#   count = 1\n#   imp_count = 0\n#   best_mdl = None\n#   best_iter = 0\n  \n#   gp = SymbolicRegressor(population_size=2000,\n#                        generations=count,\n#                        tournament_size=50,  # consider 20, was 50\n#                        parsimony_coefficient=0.0001,  # oops: 0.0001?\n#                        const_range=(-16, 16),  # consider +/-20, was 100\n#                        function_set=function_set,\n#                        # stopping_criteria=1.0,\n#                        # p_hoist_mutation=0.05,\n#                        # max_samples=.875,  # was in\n#                        # p_crossover=0.7,\n#                        # p_subtree_mutation=0.1,\n#                        # p_point_mutation=0.1,\n#                        init_depth=(6, 16),\n#                        warm_start=True,\n#                        metric='mean absolute error', verbose=1, random_state=42, n_jobs=-1, low_memory=True)\n\n#   for run in range(GENS):\n#       mdl = gp.fit(X_tr, y_tr, sample_weight=sample_wts_tr)\n#       pred = gp.predict(X_val)\n#       mae = np.sqrt(mean_absolute_error(y_val, pred))\n\n#       if mae < best and imp_count < MAX_NO_IMPROVE:\n#           best = mae\n#           count += 1\n#           gp.set_params(generations=count, warm_start=True)\n#           imp_count = 0\n#           best_iter = run\n#           if mae < MAE_THRESH:\n#               best_mdl = copy.deepcopy(mdl)\n#       elif imp_count < MAX_NO_IMPROVE:\n#           count += 1\n#           gp.set_params(generations=count, warm_start=True)\n#           imp_count += 1\n#       else:\n#           break\n\n#       print('GP MAE: %.4f, Run: %d, Best Run: %d, Fold: %d' % (mae, run, best_iter, fold_))\n\n#   maes.append(best)\n#   gens.append(run)\n      \n#   print('Finish - GP MAE: %.4f, Run: %d, Best Run: %d' % (mae, run, best_iter))\n          \n#   preds = best_mdl.predict(scaled_test_X)\n#   print(preds[0:12])\n#   predictions += preds / folds.n_splits\n\n# try:\n#     print(maes)\n#     print(np.mean(maes))\n#     print(gens)\n# except:\n#     print('oops')\n# submission.time_to_failure = predictions\n# submission.to_csv('submission_gplearn_AF0_1.csv')\n# print(submission.head(12))\nprint('ok')","execution_count":null,"outputs":[]}],"metadata":{"colab":{"collapsed_sections":[],"name":"eq_gplearn.ipynb","provenance":[],"version":"0.3.2"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.6"}},"nbformat":4,"nbformat_minor":1}