{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\nCan you measure the magnetic interactions between a pair of atoms?\n\nIn this competition, you will develop an algorithm that can predict the magnetic interaction between two atoms in a molecule (i.e., the `scalar_coupling_constant`).\n https://www.kaggle.com/c/champs-scalar-coupling/overview/description. \n \n## General information\n* based on : https://www.kaggle.com/robikscube/exploring-molecular-properties-data\n* Forked/based on: https://www.kaggle.com/artgor/molecular-properties-eda-and-models\n* https://www.kaggle.com/robikscube/exploring-molecular-properties-data\n* https://www.kaggle.com/tarunpaparaju/champs-competition-chemistry-background-and-eda\n\n* Feature engineering brute force: https://www.kaggle.com/artgor/brute-force-feature-engineering\n* https://www.kaggle.com/adrianoavelar/bond-calculation-lb-0-82  , Chemical Bond Calculation \n\n* https://www.kaggle.com/buchan/a-neural-network-approach  (Includes covalent bond calc - https://www.kaggle.com/adrianoavelar/bond-calculation-lb-0-82 - Important feature!)\n* https://www.kaggle.com/borisdee/predicting-mulliken-charges-with-acsf-descriptors  - Using external libraries + calc the features given for the train\n\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport matplotlib.pylab as plt\nimport seaborn as sns\nfrom sklearn import metrics\nfrom sklearn.model_selection import GroupKFold\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n# plt.style.use('ggplot')\ncolor_pal = [x['color'] for x in plt.rcParams['axes.prop_cycle']]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# The Data\n\nIn this competition, you will be predicting the scalar_coupling_constant between atom pairs in molecules, given the two atom types (e.g., C and H), the coupling type (e.g., 2JHC), and any features you are able to create from the molecule structure (xyz) files.\n\nFor this competition, you will not be predicting all the atom pairs in each molecule rather, you will only need to predict the pairs that are explicitly listed in the train and test files. For example, some molecules contain Fluorine (F), but you will not be predicting the scalar coupling constant for any pair that includes F.\n\nThe training and test splits are by molecule, so that no molecule in the training data is found in the test data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# # Show how the files appear in the input folder\n# !ls -GFlash --color ../input","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## train.csv and test.csv\nThe training set, where the first column `molecule_name` is the name of the molecule where the coupling constant originates (the corresponding XYZ file is located at ./structures/.xyz), the second `atom_index_0` and third column `atom_index_1` is the atom indices of the atom-pair creating the coupling and the fourth column `scalar_coupling_constant` is the scalar coupling constant that we want to be able to predict.\n\n#### Important: \n* Since this is a DEMO we will only read in a sample of the train data\n    * for realistic purposes, better would be to read in all the data, then"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('../input/train.csv').sample(frac=0.04)\nprint(\"Train data shape\",train_df.shape)\ntest_df = pd.read_csv('../input/test.csv')\nprint(\"Test shape\",test_df.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see each observation provides the:\n- molecule_name\n- atom_index_0 - index of the first atom pair\n- atom_index_1 - index of the second atom pair\n- scalar_coupling_constant (target)"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The training set is larger than the test set.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('The training set has shape {}'.format(train_df.shape))\nprint('The test set has shape {}'.format(test_df.shape))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The target distribution is pretty interesting! Spikes near zero, and -20. There is also a good bit around 80. We will return to these files later."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distribution of the target\ntrain_df['scalar_coupling_constant'].plot(kind='hist', figsize=(20, 5), bins=700, title='Distribution of the target scalar coupling constant')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Number of of atoms in molecule\nfig, ax = plt.subplots(1, 2)\ntrain_df.groupby('molecule_name').count().sort_values('id')['id'].plot(kind='hist',\n                                                                       bins=25,\n                                                                       color=color_pal[6],\n                                                                      figsize=(20, 5),\n                                                                      title='# of Atoms in Molecule (Train Set)',\n                                                                      ax=ax[0])\ntest_df.groupby('molecule_name').count().sort_values('id')['id'].plot(kind='hist',\n                                                                       bins=25,\n                                                                       color=color_pal[2],\n                                                                      figsize=(20, 5),\n                                                                      title='# of Atoms in Molecule (Test Set)',\n                                                                     ax=ax[1])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Get list of molecule types\n* We'll see the target distribution is very different for the different types: We may want to make a model for each type seperately! \n* We also have relatively few types, so this is easily tractionable"},{"metadata":{"trusted":true},"cell_type":"code","source":"typelist = train_df['type'].unique()\nprint(typelist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.groupby(\"type\")['scalar_coupling_constant'].plot(kind='hist', figsize=(18, 5), bins=500, title='Distribution of the target coupling constant, given type');\n\n# plt.figure(figsize=(26, 24))\n# for i, col in enumerate(typelist):\n#     plt.subplot(4,2, i + 1)\n#     sns.distplot(train_df[train_df['type']==col]['scalar_coupling_constant'],color ='indigo')\n#     plt.title(col)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## structures.zip annd structures csv files.\nfolder containing molecular structure (xyz) files, where: \n- the first line is the number of atoms in the molecule,\n- followed by a blank line\n- and then a line for every atom, where the first column contains the atomic element (H for hydrogen, C for carbon etc.) and the remaining columns contain the X, Y and Z cartesian coordinates (a standard format for chemists and molecular visualization programs)\n\n\n...lets have a look at the first example from the training set!"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head(1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Look at the xyz file for this example\n- The first number `5` is the number of atoms in this molecule.\n- A blank line....\n- Each following line is the element and their cartesian coordinates. So in this examples we have one Carbon atom and four Hydrogen atoms!"},{"metadata":{"trusted":true},"cell_type":"code","source":"! cat ../input/structures/dsgdb9nsd_000001.xyz","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Structures.csv\nThis file contains the same information as the individual xyz structure files, but in a single file.\n\nThis csv is a lot more useable than the `xyz` files (for most usages, but not all!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')\nstructures.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 3D Plot!\nfrom mpl_toolkits.mplot3d import Axes3D\nfig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\nexample = structures.loc[structures['molecule_name'] == 'dsgdb9nsd_000001']\nax.scatter(xs=example['x'], ys=example['y'], zs=example['z'], s=100)\nplt.suptitle('dsgdb9nsd_000001')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Additional Data\n*NOTE: additional data is provided for the molecules in Train only!*\n\n* If we want to use it for our model, we'll need to extract it externally and/or predict it!\n\n* Good 101 on the pre-provided values: https://www.kaggle.com/tarunpaparaju/champs-competition-chemistry-background-and-eda\n\n"},{"metadata":{},"cell_type":"markdown","source":"## dipole_moments.csv\n- contains the molecular electric dipole moments. These are three dimensional vectors that indicate the charge distribution in the molecule. The first column (molecule_name) are the names of the molecule, the second to fourth column are the X, Y and Z components respectively of the dipole moment."},{"metadata":{"trusted":true},"cell_type":"code","source":"dm = pd.read_csv('../input/dipole_moments.csv')\ndm.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## magnetic_shielding_tensors.csv\n- contains the magnetic shielding tensors for all atoms in the molecules. The first column (molecule_name) contains the molecule name, the second column (atom_index) contains the index of the atom in the molecule, the third to eleventh columns contain the XX, YX, ZX, XY, YY, ZY, XZ, YZ and ZZ elements of the tensor/matrix respectively."},{"metadata":{"trusted":true},"cell_type":"code","source":"mst = pd.read_csv('../input/magnetic_shielding_tensors.csv')\nmst.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## mulliken_charges.csv\n- contains the mulliken charges for all atoms in the molecules. The first column (molecule_name) contains the name of the molecule, the second column (atom_index) contains the index of the atom in the molecule, the third column (mulliken_charge) contains the mulliken charge of the atom."},{"metadata":{"trusted":true},"cell_type":"code","source":"mul = pd.read_csv('../input/mulliken_charges.csv')\nmul.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot the distribution of mulliken_charges\nmul['mulliken_charge'].plot(kind='hist', figsize=(15, 5), bins=500, title='Distribution of Mulliken Charges')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## potential_energy.csv\n- contains the potential energy of the molecules. The first column (molecule_name) contains the name of the molecule, the second column (potential_energy) contains the potential energy of the molecule."},{"metadata":{"trusted":true},"cell_type":"code","source":"pote = pd.read_csv('../input/potential_energy.csv')\npote.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot the distribution of potential_energy\npote['potential_energy'].plot(kind='hist',\n                              figsize=(15, 5),\n                              bins=500,\n                              title='Distribution of Potential Energy',\n                              color='b')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## scalar_coupling_contributions.csv\n- The scalar coupling constants in train.csv (or corresponding files) are a sum of four terms. scalar_coupling_contributions.csv contain all these terms.\n    - The first column (molecule_name) are the **name of the molecule**,\n    - the second **(atom_index_0)** and\n    - third column **(atom_index_1)** are the atom indices of the atom-pair,\n    - the fourth column indicates the **type of coupling**,\n    - the fifth column (fc) is the **Fermi Contact contribution**,\n    - the sixth column (sd) is the **Spin-dipolar contribution**,\n    - the seventh column (pso) is the **Paramagnetic spin-orbit contribution** and\n    - the eighth column (dso) is the **Diamagnetic spin-orbit contribution**."},{"metadata":{"trusted":true},"cell_type":"code","source":"scc = pd.read_csv('../input/scalar_coupling_contributions.csv')\nscc.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scc.groupby('type').count()['molecule_name'].sort_values().plot(kind='barh',\n                                                                color='grey',\n                                                               figsize=(15, 5),\n                                                               title='Count of Coupling Type in Train Set')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2, 2, figsize=(20, 10))\nscc['fc'].plot(kind='hist', ax=ax.flat[0], bins=500, title='Fermi Contact contribution', color=color_pal[0])\nscc['sd'].plot(kind='hist', ax=ax.flat[1], bins=500, title='Spin-dipolar contribution', color=color_pal[1])\nscc['pso'].plot(kind='hist', ax=ax.flat[2], bins=500, title='Paramagnetic spin-orbit contribution', color=color_pal[2])\nscc['dso'].plot(kind='hist', ax=ax.flat[3], bins=500, title='Diamagnetic spin-orbit contribution', color=color_pal[3])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Relationship between Target and Features\n** Keep in mind these features are provided for the training data ONLY**"},{"metadata":{"trusted":true},"cell_type":"code","source":"scc = scc.merge(train_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These plots are beautiful. It's a shame we don't have this data for the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Downsample to speed up plot time.\nsns.pairplot(data=scc.sample(5000), hue='type', vars=['fc','sd','pso','dso','scalar_coupling_constant'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Target vs. Atom Count"},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_count_dict = structures.groupby('molecule_name').count()['atom_index'].to_dict()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['atom_count'] = train_df['molecule_name'].map(atom_count_dict)\ntest_df['atom_count'] = test_df['molecule_name'].map(atom_count_dict)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When we look at the target `scalar_coupling_constant` in relation to the `atom_count` - there visually appears to be a relationship. We notice the gap in coupling constant values, between ~25 and ~75. It is rare to see a value within this range. Could this be a good case for a classification problem between the two clusters?"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.sample(600).plot(x='atom_count',\n                           y='scalar_coupling_constant',\n                           kind='scatter',\n                           color=color_pal[0],\n                           figsize=(20, 5),\n                           alpha=0.5)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Super Simple Baseline Model [1.239 Public LB]\nThe second simplest thing we can do as a model is predict that the target is the **average** value that we observe for that **type** in the training set!"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.groupby('type')['scalar_coupling_constant'].mean().plot(kind='barh',\n                                                                 figsize=(15, 5),\n                                                                title='Average Scalar Coupling Constant by Type')\nplt.show()\ntype_mean_dict = train_df.groupby('type')['scalar_coupling_constant'].mean().to_dict()\ntest_df['scalar_coupling_constant'] = test_df['type'].map(type_mean_dict)\ntest_df[['id','scalar_coupling_constant']].to_csv('super_simple_submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluation Metric\n\nSubmissions are evaluated on the Log of the Mean Absolute Error, calculated for each scalar coupling type, and then averaged across types, so that a 1% decrease in MAE for one type provides the same improvement in score as a 1% decrease for another type.\n\n![Eval Metric](https://i.imgur.com/AK6z3Dn.png)\n\nWhere:\n\n- `T` is the number of scalar coupling types\n- `nt` is the number of observations of type t\n- `yi` is the actual scalar coupling constant for the observation\n- `yi^` is the predicted scalar coupling constant for the observation\n\nFor this metric, the MAE for any group has a floor of 1e-9, so that the minimum (best) possible score for perfect predictions is approximately -20.7232."},{"metadata":{},"cell_type":"markdown","source":"Evaluation metric is important to understand as it determines how your model will be scored. Ideally we will set the loss function of our machine learning algorithm to use this metric so we can minimize the specific type of error.\n\nCheck out this kernel by `@abhishek` with code for the evaluation metric: https://www.kaggle.com/abhishek/competition-metric\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def metric(df, preds):\n    df[\"prediction\"] = preds\n    maes = []\n    for t in df.type.unique():\n        y_true = df[df.type==t].scalar_coupling_constant.values\n        y_pred = df[df.type==t].prediction.values\n        mae = np.log(metrics.mean_absolute_error(y_true, y_pred))\n        maes.append(mae)\n    return np.mean(maes)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Distance Feature Creation\nThis feature was found from `@inversion` 's kernel here: https://www.kaggle.com/inversion/atomic-distance-benchmark/output\nThe code was then made faster by `@seriousran` here: https://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark"},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n    df = pd.merge(df, structures, how = 'left',\n                  left_on  = ['molecule_name', f'atom_index_{atom_idx}'],\n                  right_on = ['molecule_name',  'atom_index'])\n    \n    df = df.drop('atom_index', axis=1)\n    df = df.rename(columns={'atom': f'atom_{atom_idx}',\n                            'x': f'x_{atom_idx}',\n                            'y': f'y_{atom_idx}',\n                            'z': f'z_{atom_idx}'})\n    return df\n\ntrain_df = map_atom_info(train_df, 0)\ntrain_df = map_atom_info(train_df, 1)\n\ntest_df = map_atom_info(test_df, 0)\ntest_df = map_atom_info(test_df, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark\ntrain_p_0 = train_df[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train_df[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test_df[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test_df[['x_1', 'y_1', 'z_1']].values\n\ntrain_df['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest_df['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# old, base from : @artgor's  kernel - https://www.kaggle.com/artgor/molecular-properties-eda-and-models\n# train_df['dist_to_type_mean'] = train_df['dist'] / train_df.groupby('type')['dist'].transform('mean')\n# test_df['dist_to_type_mean'] = test_df['dist'] / test_df.groupby('type')['dist'].transform('mean')\n\n### I made my own Versions: \n# train_df['dist_to_type_mean'] = train_df['dist'] / train_df.groupby('type')['dist'].transform('mean')\n# test_df['dist_to_type_mean'] = test_df['dist'] / test_df.groupby('type')['dist'].transform('mean')\n\n# train_df['dist_to_type_0_mean'] = train_df['dist'] / train_df.groupby('atom_0')['dist'].transform('mean')\n# test_df['dist_to_type_0_mean'] = test_df['dist'] / test_df.groupby('atom_0')['dist'].transform('mean')\n\n# train_df['dist_to_type_1_mean'] = train_df['dist'] / train_df.groupby('atom_1')['dist'].transform('mean')\n# test_df['dist_to_type_1_mean'] = test_df['dist'] / test_df.groupby('atom_1')['dist'].transform('mean')\n\n# train_df['molecule_type_dist_mean'] = train_df.groupby([ 'type'])['dist'].transform('mean')\n# test_df['molecule_type_dist_mean'] = test_df.groupby(['type'])['dist'].transform('mean')\n\n\ntrain_df['dist_to_type_0_mean'] = train_df['dist'] / train_df.groupby(['type','atom_0'])['dist'].transform('mean')\ntest_df['dist_to_type_0_mean'] = test_df['dist'] / test_df.groupby(['type','atom_0'])['dist'].transform('mean')\n\ntrain_df['dist_to_type_1_mean'] = train_df['dist'] / train_df.groupby(['type','atom_1'])['dist'].transform('mean')\ntest_df['dist_to_type_1_mean'] = test_df['dist'] / test_df.groupby(['type','atom_1'])['dist'].transform('mean')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# make categorical variables\natom_map = {'H': 0,\n            'C': 1,\n            'N': 2}\ntrain_df['atom_0_cat'] = train_df['atom_0'].map(atom_map).astype('int')\ntrain_df['atom_1_cat'] = train_df['atom_1'].map(atom_map).astype('int')\ntest_df['atom_0_cat'] = test_df['atom_0'].map(atom_map).astype('int')\ntest_df['atom_1_cat'] = test_df['atom_1'].map(atom_map).astype('int')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# One Hot Encode the Type\ntrain_df = pd.concat([train_df, pd.get_dummies(train_df['type'])], axis=1)\ntest_df = pd.concat([test_df, pd.get_dummies(test_df['type'])], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"color_index = 0\naxes_index = 0\nfig, axes = plt.subplots(8, 1, figsize=(20, 20), sharex=True)\nfor mtype, d in train_df.groupby('type'):\n    d['dist'].plot(kind='hist',\n                  bins=1000,\n                  title='Distribution of Distance Feature for {}'.format(mtype),\n                  color=color_pal[color_index],\n                  ax=axes[axes_index])\n    if color_index == 6:\n        color_index = 0\n    else:\n        color_index += 1\n    axes_index += 1\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Baseline Models\n- using **atom_count** and **type** as categorical features\n- Run this in a notebook to see the interactive plot of training and test error metrics."},{"metadata":{},"cell_type":"markdown","source":"# LightGBM - CV"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Configurables\nFEATURES = ['atom_index_0', 'atom_index_1',\n            'atom_0_cat',\n            'x_0', 'y_0', 'z_0',\n            'atom_1_cat', \n            'x_1', 'y_1', 'z_1', 'dist', \n#             'dist_to_type_mean',\n            'atom_count',\n            '1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN'\n            ,'dist_to_type_0_mean',\n       'dist_to_type_1_mean'\n           ]\n\n# # instead of whitelist, blacklist: # broken in feat importance part\n# DROP_FEATS = ['id', 'molecule_name', 'type','scalar_coupling_constant', 'atom_0','atom_1',]\n\nTARGET = 'scalar_coupling_constant'\nCAT_FEATS = ['atom_0','atom_1']\n## ORIG: \n# N_ESTIMATORS = 2000\n# VERBOSE = 500\n# EARLY_STOPPING_ROUNDS = 200\n# RANDOM_STATE = 529\n\n# faster:\nN_ESTIMATORS = 100\nVERBOSE = 50\nEARLY_STOPPING_ROUNDS = 5\nRANDOM_STATE = 529\n\n# ## whitelist feats: \nX = train_df[FEATURES]\nX_test = test_df[FEATURES]\n\n# ## exclude cols:\n# X = train_df.drop(DROP_FEATS,axis=1)\n# X_test = test_df.drop(DROP_FEATS,axis=1)\n\ny = train_df[TARGET]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error\nimport lightgbm as lgb\n\nlgb_params = {\n#     'num_leaves': 128, #  orig\n    'num_leaves': 64,\n#               'min_child_samples': 64, # orig\n              'min_child_samples': 32,\n              'objective': 'regression',\n#               'max_depth': 6, # ORIG\n                            'max_depth': 5,\n#               'learning_rate': 0.1, # orig\n              'learning_rate': 0.2,\n              \"boosting_type\": \"gbdt\",\n              \"subsample_freq\": 1,\n              \"subsample\": 0.9,\n              \"bagging_seed\": 11,\n              \"metric\": 'mae',\n              \"verbosity\": -1,\n              'reg_alpha': 0.1,\n              'reg_lambda': 0.4,\n              'colsample_bytree': 0.9\n         }\n\nRUN_LGB = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if RUN_LGB:\n    n_fold = 3\n    folds = KFold(n_splits=n_fold, shuffle=True, random_state=RANDOM_STATE)\n\n    # Setup arrays for storing results\n    oof = np.zeros(len(X))\n    prediction = np.zeros(len(X_test))\n    scores = []\n    feature_importance = pd.DataFrame()\n\n    # Train the model\n    for fold_n, (train_idx, valid_idx) in enumerate(folds.split(X)):\n        X_train, X_valid = X.iloc[train_idx], X.iloc[valid_idx]\n        y_train, y_valid = y.iloc[train_idx], y.iloc[valid_idx]\n        model = lgb.LGBMRegressor(**lgb_params, n_estimators = N_ESTIMATORS, n_jobs = -1)\n        model.fit(X_train, y_train,\n                  eval_set=[(X_train, y_train), (X_valid, y_valid)],\n                  eval_metric='mae',\n                  verbose=VERBOSE,\n                  early_stopping_rounds=EARLY_STOPPING_ROUNDS)\n\n        y_pred_valid = model.predict(X_valid)\n        y_pred = model.predict(X_test, num_iteration=model.best_iteration_)\n\n        # feature importance\n        fold_importance = pd.DataFrame()\n        fold_importance[\"feature\"] = FEATURES\n        fold_importance[\"importance\"] = model.feature_importances_\n        fold_importance[\"fold\"] = fold_n + 1\n        feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n        prediction /= folds.n_splits\n        scores.append(mean_absolute_error(y_valid, y_pred_valid))\n        print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n        oof[valid_idx] = y_pred_valid.reshape(-1,)\n        scores.append(mean_absolute_error(y_valid, y_pred_valid))\n        prediction += y_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importance.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save LGB Results, OOF, and Feature Importance\nIt's always a good idea to save your OOF, predictions and feature importances. You never know when they will come in handy in the future.\n\nWe'll save the Number of folds and CV score in the filename."},{"metadata":{"trusted":true},"cell_type":"code","source":"if RUN_LGB:\n    # Save Prediction and name appropriately\n    submission_csv_name = 'submission_lgb_{}folds_{}CV.csv'.format(n_fold, np.mean(scores))\n    oof_csv_name = 'oof_lgb_{}folds_{}CV.csv'.format(n_fold, np.mean(scores))\n    fi_csv_name = 'fi_lgb_{}folds_{}CV.csv'.format(n_fold, np.mean(scores))\n\n    print('Saving LGB Submission as:')\n    print(submission_csv_name)\n    ss = pd.read_csv('../input/sample_submission.csv')\n    ss['scalar_coupling_constant'] = prediction\n    ss.to_csv(submission_csv_name, index=False)\n    ss.head()\n    \n    # OOF\n    oof_df = train_df[['id','molecule_name','scalar_coupling_constant']].copy()\n    oof_df['oof_pred'] = oof\n    oof_df.to_csv(oof_csv_name, index=False)\n    \n    # Feature Importance\n    feature_importance.to_csv(fi_csv_name, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n        by=\"importance\", ascending=False)[:40].plot(kind=\"bar\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if RUN_LGB:\n    # Plot feature importance as done in https://www.kaggle.com/artgor/artgor-utils\n    feature_importance[\"importance\"] /= folds.n_splits\n    cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n        by=\"importance\", ascending=False)[:50].index\n\n    best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n    plt.figure(figsize=(15, 20));\n    ax = sns.barplot(x=\"importance\",\n                y=\"feature\",\n                hue='fold',\n                data=best_features.sort_values(by=\"importance\", ascending=False));\n    plt.title('LGB Features (avg over folds)');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Catboost\n* Similar to LGBM and XGBoost. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# from catboost import Pool, cv\n\n# RUN_CATBOOST_CV = False\n\n# if RUN_CATBOOST_CV:\n#     labels = train_df['scalar_coupling_constant'].values\n#     cat_features = ['type','atom_count','atom_0','atom_1']\n#     cv_data = train_df[['type','atom_count','atom_0','atom_1',\n#                         'x_0','y_0','z_0','x_1','y_1','z_1','dist']]\n#     cv_dataset = Pool(data=cv_data,\n#                       label=labels,\n#                       cat_features=cat_features)\n\n# ##     ITERATIONS = 100000 # ORIG\n#     ITERATIONS = 1234\n#     params = {\"iterations\": ITERATIONS,\n#               \"learning_rate\" : 0.02,\n#               \"depth\": 7,\n#               \"loss_function\": \"MAE\",\n#               \"verbose\": False,\n#               \"task_type\" : \"GPU\"}\n\n#     scores = cv(cv_dataset,\n#                 params,\n#                 fold_count=5, \n#                 plot=\"True\")\n    \n#     scores['iterations'] = scores['iterations'].astype('int')\n#     scores.set_index('iterations')[['test-MAE-mean','train-MAE-mean']].plot(figsize=(15, 5), title='CV (MAE) Score by iteration (5 Folds)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# from catboost import CatBoostRegressor, Pool\n\n# # ITERATIONS = 200000 # Default\n# ITERATIONS = 1234 # faster\n\n# FEATURES = [#'atom_index_0',\n#             'atom_index_1',\n#             'atom_0',\n#             'x_0', 'y_0', 'z_0',\n#             'atom_1', \n#             'x_1', 'y_1', 'z_1',\n#             'dist', 'dist_to_type_mean',\n#             'atom_count',\n#             'type']\n# TARGET = 'scalar_coupling_constant'\n# CAT_FEATS = ['atom_0','atom_1','type']\n\n# train_dataset = Pool(data=train_df[FEATURES],\n#                   label=train_df['scalar_coupling_constant'].values,\n#                   cat_features=CAT_FEATS)\n\n# cb_model = CatBoostRegressor(iterations=ITERATIONS,\n#                              learning_rate=0.2,\n#                              depth=7,\n#                              eval_metric='MAE',\n#                              random_seed = 529,\n#                              task_type=\"GPU\")\n\n# # Fit the model\n# cb_model.fit(train_dataset, verbose=500)\n\n# # Predict\n# test_data = test_df[FEATURES]\n\n# test_dataset = Pool(data=test_data,\n#                     cat_features=CAT_FEATS)\n\n# ss = pd.read_csv('../input/sample_submission.csv')\n# ss['scalar_coupling_constant'] = cb_model.predict(test_dataset)\n# ss.to_csv('basline_catboost_submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}