{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Libraries used"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd","execution_count":null,"outputs":[]},{"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 in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"## Reading Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\nsub = pd.read_csv('../input/sample_submission.csv')\nstructures = pd.read_csv('../input/structures.csv')\nscalar_coupling_contributions = pd.read_csv('../input/scalar_coupling_contributions.csv')\n\nprint('Train dataset shape is -> rows: {} cols:{}'.format(train.shape[0],train.shape[1]))\nprint('Test dataset shape is  -> rows: {} cols:{}'.format(test.shape[0],test.shape[1]))\nprint('Sub dataset shape is  -> rows: {} cols:{}'.format(sub.shape[0],sub.shape[1]))\nprint('Structures dataset shape is  -> rows: {} cols:{}'.format(structures.shape[0],structures.shape[1]))\nprint('Scalar_coupling_contributions dataset shape is  -> rows: {} cols:{}'.format(scalar_coupling_contributions.shape[0],\n                                                                                   scalar_coupling_contributions.shape[1]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data merge"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.merge(train, scalar_coupling_contributions, how = 'left',\n                  left_on  = ['molecule_name', 'atom_index_0', 'atom_index_1', 'type'],\n                  right_on = ['molecule_name', 'atom_index_0', 'atom_index_1', 'type'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Exploration"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('The training set has shape {}'.format(train.shape))\nprint('The test set has shape {}'.format(test.shape))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distribution of the target\ntrain['scalar_coupling_constant'].plot(kind='hist', figsize=(20, 5), bins=1000, 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.groupby('molecule_name').count().sort_values('id')['id'].plot(kind='hist',\n                                                                       bins=25,\n                                                                      figsize=(20, 5),\n                                                                      title='# of Atoms in Molecule (Train Set)',\n                                                                      ax=ax[0])\ntest.groupby('molecule_name').count().sort_values('id')['id'].plot(kind='hist',\n                                                                       bins=25,\n                                                                      figsize=(20, 5),\n                                                                      title='# of Atoms in Molecule (Test Set)',\n                                                                     ax=ax[1])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Additional Data\nmst = pd.read_csv('../input/magnetic_shielding_tensors.csv')\nmul = pd.read_csv('../input/mulliken_charges.csv')\npote = pd.read_csv('../input/potential_energy.csv')\nscc = pd.read_csv('../input/scalar_coupling_contributions.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mst.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mul.head(3)","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":{"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":{"trusted":true},"cell_type":"code","source":"scalar_coupling_contributions.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='red',\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=\"red\")\nscc['sd'].plot(kind='hist', ax=ax.flat[1], bins=500, title='Spin-dipolar contribution', color=\"blue\")\nscc['pso'].plot(kind='hist', ax=ax.flat[2], bins=500, title='Paramagnetic spin-orbit contribution', color=\"green\")\nscc['dso'].plot(kind='hist', ax=ax.flat[3], bins=500, title='Diamagnetic spin-orbit contribution', color=\"orange\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import seaborn as sns\nsns.pairplot(data=train.sample(5000), hue='type', vars=['fc','sd','pso','dso','scalar_coupling_constant'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Target vs atom count\natom_count_dict = structures.groupby('molecule_name').count()['atom_index'].to_dict()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['atom_count'] = train['molecule_name'].map(atom_count_dict)\ntest['atom_count'] = test['molecule_name'].map(atom_count_dict)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Creation"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Map the atom structure data into train and test files\n\ndef 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 = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 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[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)","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['atom_0_cat'] = train['atom_0'].map(atom_map).astype('int')\ntrain['atom_1_cat'] = train['atom_1'].map(atom_map).astype('int')\ntest['atom_0_cat'] = test['atom_0'].map(atom_map).astype('int')\ntest['atom_1_cat'] = test['atom_1'].map(atom_map).astype('int')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# One Hot Encode the Type\ntrain = pd.concat([train, pd.get_dummies(train['type'])], axis=1)\ntest = pd.concat([test, pd.get_dummies(test['type'])], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['dist_to_type_mean'] = train['dist'] / train.groupby('type')['dist'].transform('mean')\ntest['dist_to_type_mean'] = test['dist'] / test.groupby('type')['dist'].transform('mean')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Baseline models"},{"metadata":{},"cell_type":"markdown","source":"# LIGHTGBM -5 Fold Cross validation"},{"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', 'dist_to_type_mean',\n            'atom_count',\n            '1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN'\n           ]\nTARGET = 'scalar_coupling_constant'\nCAT_FEATS = ['atom_0','atom_1']\nN_ESTIMATORS = 2000\nVERBOSE = 500\nEARLY_STOPPING_ROUNDS = 200\nRANDOM_STATE = 529\n\nX = train[FEATURES]\nX_test = test[FEATURES]\ny = train[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 = {'num_leaves': 128,\n              'min_child_samples': 64,\n              'objective': 'regression',\n              'max_depth': 6,\n              'learning_rate': 0.1,\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': 1.0\n         }\nRUN_LGB = True\nif RUN_LGB:\n    n_fold = 5\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        # 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":"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":{"trusted":true},"cell_type":"code","source":"prediction_final = model.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction_final","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction_final.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub[\"scalar_coupling_constant\"] = prediction_final\nsub.to_csv('submission.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}