{"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_minor":4,"nbformat":4,"cells":[{"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\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\nplt.style.use('ggplot')\ncolor_pal = [x['color'] for x in plt.rcParams['axes.prop_cycle']]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show how the files appear in the input folder\n!ls -GFlash --color ../input","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv('../input/train.csv')\ntest_df = pd.read_csv('../input/test.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('The training set has shape {}'.format(train_df.shape))\nprint('The test set has shape {}'.format(test_df.shape))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Distribution of the target\ntrain_df['scalar_coupling_constant'].plot(kind='hist', figsize=(20, 5), bins=1000, title='Distribution of the target scalar coupling constant')\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ss = pd.read_csv('../input/sample_submission.csv')\nss.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.head(1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! cat ../input/structures/dsgdb9nsd_000001.xyz","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')\nstructures.head()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dm = pd.read_csv('../input/dipole_moments.csv')\ndm.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mst = 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')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mst.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mul.head()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pote.head()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scc.head()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scc = scc.merge(train_df)","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"atom_count_dict = structures.groupby('molecule_name').count()['atom_index'].to_dict()","metadata":{},"execution_count":null,"outputs":[]},{"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)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.sample(1000).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()","metadata":{},"execution_count":null,"outputs":[]},{"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()\n# THIS IS A MODEL!!! This is a model??\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)","metadata":{},"execution_count":null,"outputs":[]},{"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)","metadata":{},"execution_count":null,"outputs":[]},{"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)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_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)","metadata":{},"execution_count":null,"outputs":[]},{"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')","metadata":{},"execution_count":null,"outputs":[]},{"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)","metadata":{},"execution_count":null,"outputs":[]},{"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()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df['dist_to_type_mean'] = train_df['dist'] / train_df.groupby('type')['dist'].transform('mean')\ntest_df['dist_to_type_mean'] = test_df['dist'] / test_df.groupby('type')['dist'].transform('mean')","metadata":{},"execution_count":null,"outputs":[]},{"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_df[FEATURES]\nX_test = test_df[FEATURES]\ny = train_df[TARGET]","metadata":{},"execution_count":null,"outputs":[]},{"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         }\n\nRUN_LGB = False\n\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    # 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","metadata":{},"execution_count":null,"outputs":[]},{"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    # 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    # Feature Importance\n    feature_importance.to_csv(fi_csv_name, index=False)","metadata":{},"execution_count":null,"outputs":[]},{"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)');","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from catboost import Pool, cv\n\nRUN_CATBOOST_CV = False\n\nif 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\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)')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from catboost import CatBoostRegressor, Pool\n\nITERATIONS = 200000\n\nFEATURES = [#'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']\nTARGET = 'scalar_coupling_constant'\nCAT_FEATS = ['atom_0','atom_1','type']\n\ntrain_dataset = Pool(data=train_df[FEATURES],\n                  label=train_df['scalar_coupling_constant'].values,\n                  cat_features=CAT_FEATS)\n\ncb_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\ncb_model.fit(train_dataset, verbose=1000)\n\n# Predict\ntest_data = test_df[FEATURES]\n\ntest_dataset = Pool(data=test_data,\n                    cat_features=CAT_FEATS)\n\nss = pd.read_csv('../input/sample_submission.csv')\nss['scalar_coupling_constant'] = cb_model.predict(test_dataset)\nss.to_csv('basline_catboost_submission.csv', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}