{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Import statements"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#!pip install category_encoders\n#!pip install LightGBM\n\nimport os\nimport gc\nimport math\nimport numpy as np\nimport pandas as pd\nimport category_encoders as ce\nimport lightgbm as lgbm\nimport matplotlib.pyplot as plt\n\nfrom copy import copy\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error as mae","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting to display more columns than default\n\npd.set_option(\"display.max_rows\", 1000)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Loading"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Data Loading\ntrain = pd.read_csv(\"../input/train.csv\")\ntest = pd.read_csv(\"../input/test.csv\")\nsample_sub = pd.read_csv(\"../input/sample_submission.csv\")\nstructures = pd.read_csv(\"../input/structures.csv\")\n\n# Original size\nprint(f\"train.shape: {train.shape}\")\nprint(f\"test.shape: {test.shape}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Making traning, target and testing dataframes"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = train.drop(columns=['scalar_coupling_constant']).copy()\ny_train = train['scalar_coupling_constant'].copy()\nX_test = test.copy()\n\nX_train = X_train.drop(columns = ['id'])\nX_test = X_test.drop(columns = ['id'])\n\n# This is done because some steps in data pre-processing can chnage the order of rows in X_train\n# and order of rows in target may not be chnaged becuase there is not processing done on it. So to maintain which row in X_train corresponds to \n# target the index is reset.\nX_train = X_train.reset_index()\nX_test = X_test.reset_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Baseline Model"},{"metadata":{},"cell_type":"markdown","source":"First defining some functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# function to convert object columns to category\n\ndef convert_object_to_categories(X_train, X_test):\n    for col in X_train.columns:\n        if X_train[col].dtype == 'O':\n            X_train[col] = X_train[col].astype('category')\n            X_test[col] = X_test[col].astype('category')\n    return X_train, X_test\n\nX_train, X_test = convert_object_to_categories(X_train, X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The score function\n\ndef calc_score(X_train, y_train, y_val):\n    X_train_new = X_train.copy()\n    y_train_new = y_train.copy()\n    y_val_new = y_val.copy()\n    y_val_new = pd.Series(y_val_new)\n    X_train_new = X_train_new.reset_index(drop=True)\n    y_train_new = y_train_new.reset_index(drop=True)\n    X_train_new = X_train_new.merge(pd.DataFrame(y_train_new, columns=['scalar_coupling_constant']), left_index=True, right_index=True)\n    X_train_new = X_train_new.merge(pd.DataFrame(y_val_new, columns=['y_val']), left_index=True, right_index=True)\n    X_train_new['error'] = (X_train_new['scalar_coupling_constant'] - X_train_new['y_val']).abs()\n    X_train_new['count'] = 1\n    score_df = X_train_new.groupby(by = ['type']).agg({'count': 'count', 'error': 'sum'})\n    score_df['error'] = (score_df['error']/score_df['count']).apply(np.log, dtype=float)\n    score = (1/score_df.shape[0])*(score_df['error'].sum())\n    return score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Cross Validation\n\ndef cross_val(X, y):\n    kf = KFold(n_splits=5, shuffle=True, random_state=10)\n    fold = 0\n    for train_index, val_index in kf.split(X):\n        fold +=1\n        lgbm_model = lgbm.LGBMRegressor(random_state=10, n_estimators=1000)\n        lgbm_model.fit(X.loc[train_index,:], y[train_index])\n        y_val = lgbm_model.predict(X.loc[val_index,:])\n        print(f\"fold{fold} score: {calc_score(X.loc[val_index,:],y[val_index],y_val)}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Baseline cross validation score\n# cross_val(X_train, y_train)\n\n# score: 1.13 approx","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Baseline model\n# lgbm_model = lgbm.LGBMRegressor()\n# lgbm_model.fit(X_train, y_train)\n# y_predict = lgbm_model.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Joining training/test and structures\n\nX_train = X_train.merge(structures, left_on = ['molecule_name','atom_index_0'], right_on = ['molecule_name', 'atom_index'], sort = True)\nX_train = X_train.rename(columns={'atom_index': 'atom_index_0_0', 'x':'atom_index_0_x', 'y':'atom_index_0_y', 'z':'atom_index_0_z', 'atom': 'atom_0'})\nX_test = X_test.merge(structures, left_on = ['molecule_name','atom_index_0'], right_on = ['molecule_name', 'atom_index'], sort = True)\nX_test = X_test.rename(columns={'atom_index': 'atom_index_0_0', 'x':'atom_index_0_x', 'y':'atom_index_0_y', 'z':'atom_index_0_z', 'atom': 'atom_0'})\n\n\nX_train = X_train.merge(structures, left_on = ['molecule_name','atom_index_1'], right_on = ['molecule_name', 'atom_index'], sort = True)\nX_train = X_train.rename(columns={'atom_index': 'atom_index_1_1', 'x':'atom_index_1_x', 'y':'atom_index_1_y', 'z':'atom_index_1_z', 'atom': 'atom_1'})\nX_test = X_test.merge(structures, left_on = ['molecule_name','atom_index_1'], right_on = ['molecule_name', 'atom_index'], sort = True)\nX_test = X_test.rename(columns={'atom_index': 'atom_index_1_1', 'x':'atom_index_1_x', 'y':'atom_index_1_y', 'z':'atom_index_1_z', 'atom': 'atom_1'})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Dropping redundant columns\n\nX_train = X_train.drop(columns=['atom_index_0_0','atom_index_1_1'])\nX_test = X_test.drop(columns=['atom_index_0_0','atom_index_1_1'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distance between atoms\n\nX_train['distance'] = ((X_train['atom_index_0_x'] - X_train['atom_index_1_x'])**2 + (X_train['atom_index_0_y'] - X_train['atom_index_1_y'])**2 + (X_train['atom_index_0_z'] - X_train['atom_index_1_z'])**2)**0.5\nX_test['distance'] = ((X_test['atom_index_0_x'] - X_test['atom_index_1_x'])**2 + (X_test['atom_index_0_y'] - X_test['atom_index_1_y'])**2 + (X_test['atom_index_0_z'] - X_test['atom_index_1_z'])**2)**0.5  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Making column based in join type (eg. 2JHC is 2J type of join)\n\nX_train['join_type'] = X_train['type'].str.slice(0,2)\nX_test['join_type'] = X_test['type'].str.slice(0,2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(X_train['atom_0'].unique())\nprint(X_test['atom_0'].unique())\n\n# So basically atom 0 in all molecules is always hydrogen. Hence will not hep in predicting target so we can drop it\n\nX_train = X_train.drop(columns=['atom_0'])\nX_test = X_test.drop(columns=['atom_0'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Number of atoms in molecule\n\nX_train['num_atoms']=X_train.groupby(['molecule_name'])['atom_index_0'].transform('max') + 1\nX_test['num_atoms']=X_test.groupby(['molecule_name'])['atom_index_0'].transform('max') + 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Lets check correlation between all the variables\n\ndf = X_train.set_index(keys='index', drop=False).merge(pd.DataFrame(y_train, columns=['scalar_coupling_constant']), left_index=True, right_index=True)\ndf.corr()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note: There is a high negative correlation between distance and scalar_coupling_constant"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test = convert_object_to_categories(X_train, X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting index column as index so that cross validation function can pick matching rows from X_train_new and y_train. This is the reason \n# index was reset earlier. Reseting and setting index for test is not important would have made no chnage but I did because code looks more \n# consistent and also doing it is no harm\n\nX_train_new = X_train.set_index(keys='index')\nX_test_new = X_test.set_index(keys='index')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# cross_val(X_train_new, y_train)\n# score is 0.7 approx","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" So new variables especially distance is very important they are decreasing error"},{"metadata":{"trusted":true},"cell_type":"code","source":"# So that dataframe looks fine.\n\nX_train_new = X_train_new.sort_index(axis=0)\nX_test_new = X_test_new.sort_index(axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Making num_of_bonds in between atoms\n# This is just an intuition 1J is H to direct carbon(ie one bond), 2J- 'HH' hydrogens attached to a common carbon atom so (2 bonds in between) or HC means hydrogen and carbvon but different carbon which is two bonds away   \n# 3J 'HH' mean hydrgen attached to diferent carbons (which are 3 bonds away)\n\nX_train_new['num_bonds'] = X_train_new['join_type'].str.slice(0,1)\nX_test_new['num_bonds'] = X_test_new['join_type'].str.slice(0,1)\n\n# 'num_bonds' should be integer\nX_train_new['num_bonds'] = X_train_new['num_bonds'].astype('int')\nX_test_new['num_bonds'] = X_test_new['num_bonds'].astype('int')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Making angle between atoms, may be an important feature engineering. Angle will be in radians\n\ndef angle_between_vectors(df):\n    dot_products = (df['atom_index_0_x']*df['atom_index_1_x'] + df['atom_index_0_y']*df['atom_index_1_y'] + df['atom_index_0_z']*df['atom_index_1_z'])\n    magnitudes_product = (df['atom_index_0_x']**2+df['atom_index_0_y']**2+df['atom_index_0_z']**2)**0.5*(df['atom_index_1_x']**2+df['atom_index_1_y']**2+df['atom_index_1_z']**2)**0.5\n    df['angle'] = np.arccos(dot_products/magnitudes_product)\n    return df\n\nX_train_new = angle_between_vectors(X_train_new)\nX_test_new = angle_between_vectors(X_test_new)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# cross_val(X_train_new, y_train)\n\n# score: 0.66 approx","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Training error\n\n# lgbm_model = lgbm.LGBMRegressor(random_state=10, n_estimators=1000)\n# lgbm_model.fit(X_train_new, y_train)\n# y_predict = lgbm_model.predict(X_train_new)\n# print(f\"training score: {calc_score(X_train_new, y_train, y_predict)}\")\n\n# training score: 0.6372931683989951","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Modelling. I tried to see if increasing n_estomators improves score. It does so I am taking n_estomator=1000 \n# More estimators also takes a lot of time to compute.\n\nlgbm_model = lgbm.LGBMRegressor(random_state=10, n_estimators=1000)\nlgbm_model.fit(X_train_new, y_train)\ny_predict = lgbm_model.predict(X_test_new)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot importance"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(12,9))\nlgbm.plot_importance(lgbm_model, ax)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see num_bonds is not very important feature. There can be two reasons either the intuition of making num_bonds is wrong or simply num_bonds is not enough information i.e. for for two rows with num_bonds=2, scalar_coupling constant can be very different since it can also depend upon whether it is between two hydrogen atoms or one hydrogen and one carbon. Also worth noting that the feature 'type' and 'join_type' both are more important than num_bonds but 'type' is much more important than both of them."},{"metadata":{},"cell_type":"markdown","source":"Making Submission file."},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_sub['scalar_coupling_constant'] = list(y_predict)\nsample_sub.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### References"},{"metadata":{},"cell_type":"markdown","source":"I watched this video to get an intuition for making more features\n1. https://www.youtube.com/watch?v=CUI9bWH1i1Y "}],"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}