{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\n---\n\nThe goal of this competition is to find a **coupling constant** through intermolecular bonding. Coupling Constant is constant indicating the strength of the physical interaction (here between atoms), called a complete combination when the combined constant is 1.\n\nWe receive the coupling constant information from the *'train data'* by the coupling of two atoms. And from the *'Structures data'* we get the x, y, z axis information for each atom.  \n\n#### Training Strategy\nFor model training, I did not train the entire data, but I did the train by each **'type'** of molecular.\nAnd I used the **LightGBM** for model training. "},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\nimport lightgbm as lgb\n\n# path\npath_dir = '../input/champs-scalar-coupling/'\nfile_list = os.listdir(path_dir)\nfile_list","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. Load Train/Test Data\n**Columns**\n- molecule_name\n- atom_index_0 / atom_index_1\n- type\n- Coupling Constant : A constant indicating the strength of the physical interaction (here between atoms), called a complete combination when the combined constant is 1."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(path_dir+'train.csv')\ntest_df = pd.read_csv(path_dir+'test.csv')   # target = 'scalar_coupling_constant'\n\nprint('Length of train set: {}'.format(len(train_df)))\nprint('Length of test set: {}'.format(len(test_df)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Unique molecule of train set: {}'.format(len(train_df['molecule_name'].unique())))\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Unique molecule of test set: {}'.format(len(test_df['molecule_name'].unique())))\ntest_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. EDA"},{"metadata":{},"cell_type":"markdown","source":"#### 2.1 Distribution of Target ('scalar_coupling_constant')\n- Min Value : -36.2186\n- Max Value : 204.88\n- Most are between -20 and +20\n- Small distribution exists between 80 and 100"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distribution of target\nprint('Min Value of Target : {}'.format(train_df['scalar_coupling_constant'].min()))\nprint('Max Value of Target : {}'.format(train_df['scalar_coupling_constant'].max()))\n\nplt.figure(figsize=(11, 5))\nsns.distplot(train_df['scalar_coupling_constant'])\nplt.title('Distribution of scalar_coupling_constant')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 2.2 Distribution of 'scalar_coupling_constant' by type\n- '1JHC' type is distributed in a relatively high scalar coupling range(+66.6 ~ +204.8) \n- '2JHH' type is distributed in a relatively low scalar coupling range(-35.1 ~ +11.8"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distribution of 'scalar_coupling_constant' by type\nplt.figure(figsize=(14, 13))\nfor i, t in enumerate(train_df['type'].unique()):\n    plt.subplot(4,2, i+1)\n    sns.distplot(train_df[train_df['type'] == t]['scalar_coupling_constant'])\n    plt.title('Distribution of coupling constant by type '+ t)\n    plt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 2.3 Count by 'type'\n- High in order 3JHC, 2JHC, 1JHC, 3JHH, 2JHH, 3JHN, 2JHN, 1JHN."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Count by 'type'\ntype_index = train_df['type'].value_counts().index\ntype_cnt = train_df['type'].value_counts()\n\nplt.figure(figsize=(11, 4))\nsns.barplot(x=type_index, y=type_cnt)\nplt.xlabel('type'); plt.ylabel('Count')\nplt.title('Count by type')\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 2.4 Count by atom index 0, 1\n- Atom index 0 has the most number of distributions from 9 to 18.\n- Atom index 1 has the most number of distributions from 1 to 8."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Count by atom index 0, 1\nfor i in [0, 1]:\n    atom_index = train_df['atom_index_'+str(i)].value_counts().index\n    atom_cnt = train_df['atom_index_'+str(i)].value_counts()\n    \n    plt.figure(figsize=(11, 4))\n    sns.barplot(x=atom_index, y=atom_cnt)\n    plt.xlabel('atom index '+str(i)); plt.ylabel('Count')\n    plt.title('Count by atom index '+str(i))\n    plt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Load Structures Data\n**Columns**\n- molecule_name\n- atom_index\n- atom\n- x, y, z axis of atom"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_df = pd.read_csv(path_dir+'structures.csv')\n\nprint('Length of test set: {}'.format(len(structures_df)))\nstructures_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 3.1. 3Dimension plot by Molecule"},{"metadata":{"trusted":true},"cell_type":"code","source":"for name in structures_df['molecule_name'].unique()[:4]:\n    structures_molecule =structures_df[structures_df['molecule_name'] == name]\n\n    fig = plt.figure(figsize=(8, 5))\n    ax = fig.add_subplot(111, projection='3d')\n    ax.scatter(structures_molecule['x'], structures_molecule['y'], structures_molecule['z'], s=200, edgecolors='white')\n    ax.set_title(str(name)+ ' 3D plot')\n    ax.set_xlabel('x')\n    ax.set_ylabel('y')\n    ax.set_zlabel('z')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4. Preprocessing\n#### 4.1. Merge Train&Test - Structures Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"def mapping_atom_index(df, atom_idx):\n    atom_idx = str(atom_idx)\n    df = pd.merge(df, structures_df,\n                  left_on  = ['molecule_name', 'atom_index_'+atom_idx],\n                  right_on = ['molecule_name',  'atom_index'],\n                 how = 'left')\n    \n    df = df.drop('atom_index', axis=1)\n    df = df.rename(columns={'atom': 'atom_'+atom_idx,\n                            'x': 'x_'+atom_idx,\n                            'y': 'y_'+atom_idx,\n                            'z': 'z_'+atom_idx})\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_merge = mapping_atom_index(train_df, 0)\ntrain_merge = mapping_atom_index(train_merge, 1)\n\ntest_merge = mapping_atom_index(test_df, 0)\ntest_merge = mapping_atom_index(test_merge, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_tmp = train_merge[['id','molecule_name','type']]\ntest_tmp = test_merge[['id','molecule_name','type']]\n\ntrain_merge.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 4.2. Derived variables - 'Distance'\n- distance between *x axis* of atom index\n- distance between *y axis* of atom index\n- distance between *z axis* of atom index\n- distance between *atom*"},{"metadata":{"trusted":true},"cell_type":"code","source":"def dist_between_atom(df):\n    # distance between axis of atom\n    df['x_dist'] = (df['x_0'] - df['x_1'])**2\n    df['y_dist'] = (df['y_0'] - df['y_1'])**2\n    df['z_dist'] = (df['z_0'] - df['z_1'])**2\n    \n    # distance between atom\n    df['atom_dist'] = (df['x_dist']+df['y_dist']+df['z_dist'])**0.5\n    \n    return df\n    \ntrain_dist = dist_between_atom(train_merge)\ntest_dist = dist_between_atom(test_merge)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dist.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 4.3. Label encoding \n- type, atom_0, atom_1"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Label encoding\ncategorical_features = ['type', 'atom_0', 'atom_1']\nfor col in categorical_features:\n    le = LabelEncoder()\n    le.fit(list(train_dist[col].values) + list(test_dist[col].values))\n    train_dist[col] = le.transform(list(train_dist[col].values))\n    test_dist[col] = le.transform(list(test_dist[col].values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_le = train_dist.copy()\ntest_le = test_dist.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_le.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 4.4. Standardization\n- z = (x - u) / s"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train\ntrain_data = train_le.drop(['id','molecule_name','scalar_coupling_constant'], axis=1)\ntrain_target = train_le['scalar_coupling_constant']\n# test\ntest_data = test_le.drop(['id','molecule_name',], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# z-score standardization\ntrain_scale = (train_data - train_data.mean()) / train_data.mean()\ntrain_scale = train_scale.fillna(0)\ntest_scale = (test_data - train_data.mean()) / train_data.mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 4.5. Variable Correlations"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_corr = train_scale.copy()\ntrain_corr['scalar_coupling_constant'] = train_target\ncorrmat = train_corr.corr()\ntop_corr_features = corrmat.index[abs(corrmat['scalar_coupling_constant']) >= 0.1]\n\nplt.figure(figsize=(10,7))\nsns.heatmap(train_corr[top_corr_features].corr(), annot=True, cmap=\"RdYlGn\")\nplt.title('Variable Correlations')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5. Training by LightGBM "},{"metadata":{},"cell_type":"markdown","source":"#### 5.1. Training by 'type' through LightGBM"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_scale = train_scale.drop('type', axis=1)\ntrain_scale['type'] = train_tmp['type']\ntrain_scale['scalar_coupling_constant'] = train_target\n\ntest_scale = test_scale.drop('type', axis=1)\ntest_scale[['id', 'type']] = test_tmp[['id', 'type']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"score_by_type = []    # List of Validation score by type \nfeature_importance_df = []\ntest_pred_df = pd.DataFrame(columns=['id', 'scalar_coupling_constant'])   # Dataframe for submission\n\n# Extract data by type\ntypes = train_tmp['type'].unique()\nfor typ in types:\n    print('---Type of '+str(typ)+'---')\n    \n    train = train_scale[train_scale['type'] == typ]\n    target = train['scalar_coupling_constant']\n    train = train.drop(['type','scalar_coupling_constant'], axis=1)\n    \n    # Split train set / valid set\n    x_train, x_val, y_train, y_val = train_test_split(train, target, random_state=42)\n    \n    # LightGBM\n    categorical_features = ['atom_0','atom_1']\n    lgb_train = lgb.Dataset(x_train, y_train, categorical_feature=categorical_features)\n    lgb_val = lgb.Dataset(x_val, y_val, categorical_feature=categorical_features)\n\n    # Parameters of LightGBM\n    params = {'num_leaves': 128,\n              'min_child_samples': 79,\n              'objective': 'regression',\n              'max_depth': 9,\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.13,\n              'reg_lambda': 0.36,\n              'colsample_bytree': 1.0\n             }\n    # Training\n    lgb_model = lgb.train(params, lgb_train, valid_sets=[lgb_train, lgb_val], \n                          num_boost_round=15000,    # Number of boosting iterations.\n                          early_stopping_rounds=500,    # early stopping for valid set\n                          verbose_eval=2500)    # eval metric on the valid set is printed at 2500 each boosting\n    \n    # Feature Importances\n    feature_importance = lgb_model.feature_importance()\n    df_fi = pd.DataFrame({'columns':x_train.columns, 'importances':feature_importance})\n    df_fi = df_fi[df_fi['importances'] > 0].sort_values(by=['importances'], ascending=False)\n    feature_importance_df.append(df_fi)\n    \n    # Predict Validation set\n    score_by_type.append(list(lgb_model.best_score['valid_1'].values()))\n    \n    # Predict Test set\n    test = test_scale[test_scale['type'] == typ]\n    test_id = test['id']\n    test = test.drop(['id','type'], axis=1)\n    \n    test_preds = lgb_model.predict(test)\n    test_pred_df = pd.concat([test_pred_df, pd.DataFrame({'id':test_id, 'scalar_coupling_constant':test_preds})], axis=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 5.2. Validation MAE by type"},{"metadata":{"trusted":true},"cell_type":"code","source":"for typ, score in zip(types, score_by_type):\n    print('Type {} valid MAE  : {}'.format(str(typ), score))\n\nprint('\\nAverage of valid MAE  : {}'.format(np.mean(score_by_type)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 5.3. Feature Importances Plot by Type"},{"metadata":{"trusted":true},"cell_type":"code","source":"for typ, df_fi in zip(types, feature_importance_df):\n    fig = plt.figure(figsize=(12, 6))\n    ax = sns.barplot(df_fi['columns'], df_fi['importances'])\n    ax.set_xticklabels(df_fi['columns'], rotation=80, fontsize=13)\n    plt.title('Type '+str(typ)+' feature importance')\n    plt.tight_layout()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 5.4. Save prediction of test set to *.csv "},{"metadata":{"trusted":true},"cell_type":"code","source":"test_pred_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_pred_df.to_csv('lgb_submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"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.9"}},"nbformat":4,"nbformat_minor":1}