{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\n\nVisit my script of this kernel on: https://www.kaggle.com/kabure/lightgbm-full-pipeline-model and if you like it, please <b>upvote the kernel</b> =)\n\n## Description\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":{},"cell_type":"markdown","source":"<b>Disclaimer:</b> I don't have great knowledge about molecules and atoms. This is a absolutelly new world to me.  \n\nI am sure that it will be very challenging and fun. \n\nTo start, I will need to learn a lot but at this moment I have some questions, like:\n- What is this data?\n- What this columns means?\n- What is the distribution of the data?\n- How it works and correlated with other columns? \n- The extra data have some important or interesting contribution to the competition? \n\n\n## NOTE: This kernel are under construction. \n> Votes up and stay tuned. If you want the full code, fork this kernel. \n"},{"metadata":{},"cell_type":"markdown","source":"# Importing the libraries"},{"metadata":{"trusted":true},"cell_type":"code","source":"# To manipulate data\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport random\n\n# Standard plotly imports\nimport plotly.plotly as py\nimport plotly.graph_objs as go\nimport plotly.tools as tls\nfrom plotly.offline import iplot, init_notebook_mode\nimport cufflinks\n\n# Using plotly + cufflinks in offline mode\ninit_notebook_mode(connected=True)\ncufflinks.go_offline(connected=True)\n\nfrom functools import partial\nfrom hyperopt import fmin, hp, tpe, Trials, space_eval, STATUS_OK, STATUS_RUNNING\nfrom lightgbm import LGBMRegressor\nimport lightgbm\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder\n\nfrom tqdm import tqdm_notebook\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Importing data sets"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/champs-scalar-coupling/train.csv')\n# df_pot_energy = pd.read_csv('../input/potential_energy.csv')\n# df_mul_charges = pd.read_csv('../input/mulliken_charges.csv')\n# df_scal_coup_contrib = pd.read_csv('../input/scalar_coupling_contributions.csv')\n# df_magn_shield_tensor = pd.read_csv('../input/magnetic_shielding_tensors.csv')\n# df_dipole_moment = pd.read_csv('../input/dipole_moments.csv')\ndf_structure = pd.read_csv('../input/champs-scalar-coupling/structures.csv')\ndf_test = pd.read_csv('../input/champs-scalar-coupling/test.csv')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Looking some informations of all datasets\n- shape\n- first rows\n- nunique values"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"print(df_train.shape)\nprint(\"\")\nprint(df_train.head())\nprint(\"\")\nprint(df_train.nunique())\nprint(\"\")\nprint(df_train.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(df_structure.shape)\nprint(\"\")\nprint(df_structure.head())\nprint(\"\")\nprint(df_structure.nunique())\nprint(\"\")\nprint(df_structure.columns)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Exploring"},{"metadata":{},"cell_type":"markdown","source":"## Understanding the Target Distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['scalar_coupling_constant'].sample(200000).iplot(kind='hist', title='Scalar Coupling Constant Distribuition',\n                                                          xTitle='Scalar Coupling value', yTitle='Probability', histnorm='percent' )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"We have a clear distribution"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Looking the different Types"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\n\ng = plt.subplot(211)\ng = sns.countplot(x='type', data=df_train, )\ng.set_title(\"Count of Different Molecule Types\", fontsize=22)\ng.set_xlabel(\"Molecular Type Name\", fontsize=18)\ng.set_ylabel(\"Count Molecules in each Type\", fontsize=18)\n\ng1 = plt.subplot(212)\ng1 = sns.boxplot(x='type', y='scalar_coupling_constant', data=df_train )\ng1.set_title(\"Count of Different Molecule Types\", fontsize=22)\ng1.set_xlabel(\"Molecular Type Name\", fontsize=18)\ng1.set_ylabel(\"Scalar Coupling distribution\", fontsize=18)\n\nplt.subplots_adjust(wspace = 0.5, hspace = 0.5,top = 0.9)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Now, I will create a interactive button to set all chart in on chunk of s"},{"metadata":{},"cell_type":"markdown","source":"## Atom index 0 and Atom index 1 Counting distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\n\ng = plt.subplot(211)\ng = sns.countplot(x='atom_index_0', data=df_train, color='darkblue' )\ng.set_title(\"Count of Atom index 0\", fontsize=22)\ng.set_xlabel(\"index 0 Number\", fontsize=18)\ng.set_ylabel(\"Count\", fontsize=18)\n\ng1 = plt.subplot(212)\ng1 = sns.countplot(x='atom_index_1',data=df_train, color='darkblue' )\ng1.set_title(\"Count of Atom index 1\", fontsize=22)\ng1.set_xlabel(\"index 1 Number\", fontsize=18)\ng1.set_ylabel(\"Count\", fontsize=18)\n\nplt.subplots_adjust(wspace = 0.5, hspace = 0.5,top = 0.9)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Interesting. \n\nWe can see that in the index 0 the highet number of atoms has index 9 until 19 ~ 20; <br>\nIn the Index 1 the highest part of atoms has values between 0 and eight; \n\n"},{"metadata":{},"cell_type":"markdown","source":"## Let's cross the index and see if we have any clear pattern on data distribution of atoms\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"cross_index = ['atom_index_0','atom_index_1'] #seting the desired \n\ncm = sns.light_palette(\"green\", as_cmap=True)\npd.crosstab(df_train[cross_index[0]], df_train[cross_index[1]]).style.background_gradient(cmap = cm)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"## Scale Coupling Distribution by Atom Index 0 and Index 1"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(17,12))\n\ng = plt.subplot(211)\ng = sns.boxenplot(x='atom_index_0', y='scalar_coupling_constant', data=df_train, color='darkred' )\ng.set_title(\"Count of Atom index 0\", fontsize=22)\ng.set_xlabel(\"index 0 Number\", fontsize=18)\ng.set_ylabel(\"Count\", fontsize=18)\n\ng1 = plt.subplot(212)\ng1 = sns.boxenplot(x='atom_index_1', y='scalar_coupling_constant', data=df_train, color='darkblue' )\ng1.set_title(\"Count of Atom index 1\", fontsize=22)\ng1.set_xlabel(\"index 1 Number\", fontsize=18)\ng1.set_ylabel(\"Scalar Coupling distribution\", fontsize=18)\n\nplt.subplots_adjust(wspace = 0.5, hspace = 0.5,top = 0.9)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Cool."},{"metadata":{},"cell_type":"markdown","source":"## Now I will use index cross to get the mean scalar coupling by "},{"metadata":{"trusted":true},"cell_type":"code","source":"scalar_index_cross = ['atom_index_0', 'atom_index_1'] #seting the desired \n\ncm = sns.light_palette(\"green\", as_cmap=True)\npd.crosstab(df_train[scalar_index_cross[0]], df_train[scalar_index_cross[1]], \n            values=df_train['scalar_coupling_constant'], aggfunc=['mean']).style.background_gradient(cmap = cm)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that exists a well defined pattern of scallar coupling and indexes"},{"metadata":{},"cell_type":"markdown","source":"### Maping the atom structure"},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n    df = pd.merge(df, df_structure, 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\ndf_train = map_atom_info(df_train, 0)\ndf_train = map_atom_info(df_train, 1)\n\ndf_test = map_atom_info(df_test, 0)\ndf_test = map_atom_info(df_test, 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Calculating the distance "},{"metadata":{"trusted":true},"cell_type":"code","source":"## This is a very performative way to compute the distances\ntrain_p_0 = df_train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = df_train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = df_test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = df_test[['x_1', 'y_1', 'z_1']].values\n\n## linalg.norm, explanation:\n## This function is able to return one of eight different matrix norms, \n## or one of an infinite number of vector norms (described below),\n## depending on the value of the ord parameter.\ndf_train['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ndf_test['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\n\ndf_train['dist_x'] = (df_train['x_0'] - df_train['x_1']) ** 2\ndf_test['dist_x'] = (df_test['x_0'] - df_test['x_1']) ** 2\ndf_train['dist_y'] = (df_train['y_0'] - df_train['y_1']) ** 2\ndf_test['dist_y'] = (df_test['y_0'] - df_test['y_1']) ** 2\ndf_train['dist_z'] = (df_train['z_0'] - df_train['z_1']) ** 2\ndf_test['dist_z'] = (df_test['z_0'] - df_test['z_1']) ** 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['dist'].sample(200000).iplot(kind='hist', title='Scalar Coupling Constant Distribuition',\n                                                          xTitle='Scalar Coupling value', yTitle='Probability', histnorm='percent' )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Interesting. "},{"metadata":{},"cell_type":"markdown","source":"### Now, lets see distribution by:\n - Type\n - Index 0\n - Index 1\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"g = sns.FacetGrid(df_train, col=\"type\", col_wrap=2, height=4, aspect=1.5)\ng.map(sns.violinplot, \"dist\");\nplt.subplots_adjust(top=0.9)\ng.fig.suptitle('Violin Distribution of Dist by Each Type', fontsize=25)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Looking Boxplot of Index 0 and 1"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(14,12))\n\ng1 = plt.subplot(211)\ng1 = sns.boxenplot(x='atom_index_0', y='dist', data=df_train, color='blue' )\ng1.set_title(\"Distance Distribution by Atom index 0\", fontsize=22)\ng1.set_xlabel(\"Index 0 Number\", fontsize=18)\ng1.set_ylabel(\"Distance Distribution\", fontsize=18)\n\ng2 = plt.subplot(212)\ng2 = sns.boxenplot(x='atom_index_1', y='dist', data=df_train, color='green' )\ng2.set_title(\"Distance Distribution by Atom index 1\", fontsize=22)\ng2.set_xlabel(\"Index 1 Number\", fontsize=18)\ng2.set_ylabel(\"Distance Distribution\", fontsize=18)\n\nplt.subplots_adjust(wspace = 0.5, hspace = 0.3,top = 0.9)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"### Networks using the calculated distances\nWe have molecules, atom pairs, so this means data, which is interconnected. <br>\nNetwork graphs should be useful to visualize such data!<br>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"## graph from:\n### https://www.kaggle.com/artgor/molecular-properties-eda-and-models\n\nimport networkx as nx\n\nfig, ax = plt.subplots(figsize = (20, 12))\nfor i, t in enumerate(df_train['type'].unique()):\n    train_type = df_train.loc[df_train['type'] == t]\n    G = nx.from_pandas_edgelist(train_type, 'atom_index_0', 'atom_index_1', ['dist'])\n    plt.subplot(2, 4, i + 1);\n    nx.draw(G, with_labels=True);\n    plt.title(f'Graph for type {t}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"But there is a little problem: as we saw earlier, there are atoms which are very rare, as a result graphs will be skewed due to them. Now I'll drop atoms for each type which are present in less then 1% of connections"},{"metadata":{},"cell_type":"markdown","source":"## Removing rare atoms for each type "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 12))\nfor i, t in enumerate(df_train['type'].unique()):\n    train_type = df_train.loc[df_train['type'] == t]\n    bad_atoms_0 = list(train_type['atom_index_0'].value_counts(normalize=True)[train_type['atom_index_0'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms_1 = list(train_type['atom_index_1'].value_counts(normalize=True)[train_type['atom_index_1'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms = list(set(bad_atoms_0 + bad_atoms_1))\n    train_type = train_type.loc[(train_type['atom_index_0'].isin(bad_atoms_0) == False) & (train_type['atom_index_1'].isin(bad_atoms_0) == False)]\n    G = nx.from_pandas_edgelist(train_type, 'atom_index_0', 'atom_index_1', ['scalar_coupling_constant'])\n    plt.subplot(2, 4, i + 1);\n    nx.draw(G, with_labels=True);\n    plt.title(f'Graph for type {t}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's a very interesting graph that show the distance"},{"metadata":{},"cell_type":"markdown","source":"### Now, I will plot some molecules in a 3D graph to we see the difference in their distances"},{"metadata":{"trusted":true},"cell_type":"code","source":"number_of_colors = df_structure.atom.value_counts().count() # total number of different collors that we will use\n\n# Here I will generate a bunch of hexadecimal colors \ncolor = [\"#\"+''.join([random.choice('0123456789ABCDEF') for j in range(6)])\n             for i in range(number_of_colors)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_structure['color'] = np.nan\n\nfor idx, col in enumerate(df_structure.atom.value_counts().index):\n    listcol = ['#C15477', '#7ECF7B', '#4BDBBD', '#338340', '#F9E951']\n    df_structure.loc[df_structure['atom'] == col, 'color'] = listcol[idx]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"## Building a \ntrace1 = go.Scatter3d(\n    x=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_000003'].x,\n    y=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_000003'].y,\n    z=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_000003'].z, \n    hovertext=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_000003'].atom,\n    mode='markers', name=\"3 Atoms Molecule\",visible=False, \n    marker=dict(\n        size=10,\n        color=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_000003'].color,                \n    )\n) \n\ntrace2 = go.Scatter3d(\n    x=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].x,\n    y=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].y,\n    z=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].z, \n    hovertext=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].atom,\n    mode='markers', name=\"8 Atoms Molecule\",visible=True, \n    marker=dict(symbol='circle',\n                size=6,\n                color=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].color,\n                colorscale='Viridis',\n                line=dict(color='rgb(50,50,50)', width=0.5)\n               ),\n    hoverinfo='text'\n)\n\ntrace3 = go.Scatter3d(\n    x=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_051136'].x,\n    y=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_051136'].y,\n    z=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_051136'].z, \n    hovertext=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_051136'].atom,\n    mode='markers', name=\"17 Atoms Molecule\",visible=False, \n    marker=dict(\n        size=10,\n        color=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].color,             \n    )\n) \n\ntrace4 = go.Scatter3d(\n    x=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_088951'].x,\n    y=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_088951'].y,\n    z=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_088951'].z, \n    hovertext=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_088951'].atom,\n    mode='markers', name=\"25 Atoms Molecule\",visible=False, \n    marker=dict(\n        size=10,\n        color=df_structure[df_structure['molecule_name'] == 'dsgdb9nsd_002116'].color,             \n    )\n) \ndata = [trace1, trace2, trace3, trace4]\n\nupdatemenus = list([\n    dict(active=-1,\n         showactive=True,\n         buttons=list([  \n            dict(\n                label = '3 Atoms',\n                 method = 'update',\n                 args = [{'visible': [True, False, False, False]}, \n                     {'title': 'Molecule with 3 Atoms'}]),\n             \n             dict(\n                  label = '8 Atoms',\n                 method = 'update',\n                 args = [{'visible': [False, True, False, False]},\n                     {'title': 'Molecule with 8 Atoms'}]),\n\n            dict(\n                 label = '17 Atoms',\n                 method = 'update',\n                 args = [{'visible': [False, False, True, False]},\n                     {'title': 'Molecule with 17 Atoms'}]),\n\n            dict(\n                 label = '25 Atoms',\n                 method = 'update',\n                 args = [{'visible': [False, False, False, True]},\n                     {'title': 'Molecule with 25 Atoms'}])\n        ]),\n    )\n])\n\n\nlayout = dict(title=\"The distance between atoms in some molecules <br>(Select from Dropdown)<br> Molecule of 8 Atoms\", \n              showlegend=False,\n              updatemenus=updatemenus)\n\nfig = dict(data=data, layout=layout)\n\niplot(fig)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Wow, cool. <br>\nTry to zoom in !!!! "},{"metadata":{},"cell_type":"markdown","source":"## Model \n- Feature Engineering\n- Preprocessing\n- Modeeling\n- Model comparison\n- Feature importances "},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"## Function to reduce the DF size\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)    \n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['type_0'] = df_train['type'].apply(lambda x: x[0])\ndf_test['type_0'] = df_test['type'].apply(lambda x: x[0])\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"## All feature engineering references I got on Artgor's Kernel:\n## https://www.kaggle.com/artgor/brute-force-feature-engineering\n\ndef create_features(df):\n    df['molecule_couples'] = df.groupby('molecule_name')['id'].transform('count')\n    df['molecule_dist_mean'] = df.groupby('molecule_name')['dist'].transform('mean')\n    df['molecule_dist_min'] = df.groupby('molecule_name')['dist'].transform('min')\n    df['molecule_dist_max'] = df.groupby('molecule_name')['dist'].transform('max')\n    df['atom_0_couples_count'] = df.groupby(['molecule_name', 'atom_index_0'])['id'].transform('count')\n    df['atom_1_couples_count'] = df.groupby(['molecule_name', 'atom_index_1'])['id'].transform('count')\n    df[f'molecule_atom_index_0_x_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['x_1'].transform('std')\n    df[f'molecule_atom_index_0_y_1_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('mean')\n    df[f'molecule_atom_index_0_y_1_mean_diff'] = df[f'molecule_atom_index_0_y_1_mean'] - df['y_1']\n    df[f'molecule_atom_index_0_y_1_mean_div'] = df[f'molecule_atom_index_0_y_1_mean'] / df['y_1']\n    df[f'molecule_atom_index_0_y_1_max'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('max')\n    df[f'molecule_atom_index_0_y_1_max_diff'] = df[f'molecule_atom_index_0_y_1_max'] - df['y_1']\n    df[f'molecule_atom_index_0_y_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('std')\n    df[f'molecule_atom_index_0_z_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['z_1'].transform('std')\n    df[f'molecule_atom_index_0_dist_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('mean')\n    df[f'molecule_atom_index_0_dist_mean_diff'] = df[f'molecule_atom_index_0_dist_mean'] - df['dist']\n    df[f'molecule_atom_index_0_dist_mean_div'] = df[f'molecule_atom_index_0_dist_mean'] / df['dist']\n    df[f'molecule_atom_index_0_dist_max'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('max')\n    df[f'molecule_atom_index_0_dist_max_diff'] = df[f'molecule_atom_index_0_dist_max'] - df['dist']\n    df[f'molecule_atom_index_0_dist_max_div'] = df[f'molecule_atom_index_0_dist_max'] / df['dist']\n    df[f'molecule_atom_index_0_dist_min'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('min')\n    df[f'molecule_atom_index_0_dist_min_diff'] = df[f'molecule_atom_index_0_dist_min'] - df['dist']\n    df[f'molecule_atom_index_0_dist_min_div'] = df[f'molecule_atom_index_0_dist_min'] / df['dist']\n    df[f'molecule_atom_index_0_dist_std'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('std')\n    df[f'molecule_atom_index_0_dist_std_diff'] = df[f'molecule_atom_index_0_dist_std'] - df['dist']\n    df[f'molecule_atom_index_0_dist_std_div'] = df[f'molecule_atom_index_0_dist_std'] / df['dist']\n    df[f'molecule_atom_index_1_dist_mean'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('mean')\n    df[f'molecule_atom_index_1_dist_mean_diff'] = df[f'molecule_atom_index_1_dist_mean'] - df['dist']\n    df[f'molecule_atom_index_1_dist_mean_div'] = df[f'molecule_atom_index_1_dist_mean'] / df['dist']\n    df[f'molecule_atom_index_1_dist_max'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('max')\n    df[f'molecule_atom_index_1_dist_max_diff'] = df[f'molecule_atom_index_1_dist_max'] - df['dist']\n    df[f'molecule_atom_index_1_dist_max_div'] = df[f'molecule_atom_index_1_dist_max'] / df['dist']\n    df[f'molecule_atom_index_1_dist_min'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('min')\n    df[f'molecule_atom_index_1_dist_min_diff'] = df[f'molecule_atom_index_1_dist_min'] - df['dist']\n    df[f'molecule_atom_index_1_dist_min_div'] = df[f'molecule_atom_index_1_dist_min'] / df['dist']\n    df[f'molecule_atom_index_1_dist_std'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('std')\n    df[f'molecule_atom_index_1_dist_std_diff'] = df[f'molecule_atom_index_1_dist_std'] - df['dist']\n    df[f'molecule_atom_index_1_dist_std_div'] = df[f'molecule_atom_index_1_dist_std'] / df['dist']\n    df[f'molecule_atom_1_dist_mean'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('mean')\n    df[f'molecule_atom_1_dist_min'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('min')\n    df[f'molecule_atom_1_dist_min_diff'] = df[f'molecule_atom_1_dist_min'] - df['dist']\n    df[f'molecule_atom_1_dist_min_div'] = df[f'molecule_atom_1_dist_min'] / df['dist']\n    df[f'molecule_atom_1_dist_std'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('std')\n    df[f'molecule_atom_1_dist_std_diff'] = df[f'molecule_atom_1_dist_std'] - df['dist']\n    df[f'molecule_type_0_dist_std'] = df.groupby(['molecule_name', 'type_0'])['dist'].transform('std')\n    df[f'molecule_type_0_dist_std_diff'] = df[f'molecule_type_0_dist_std'] - df['dist']\n    df[f'molecule_type_dist_mean'] = df.groupby(['molecule_name', 'type'])['dist'].transform('mean')\n    df[f'molecule_type_dist_mean_diff'] = df[f'molecule_type_dist_mean'] - df['dist']\n    df[f'molecule_type_dist_mean_div'] = df[f'molecule_type_dist_mean'] / df['dist']\n    df[f'molecule_type_dist_max'] = df.groupby(['molecule_name', 'type'])['dist'].transform('max')\n    df[f'molecule_type_dist_min'] = df.groupby(['molecule_name', 'type'])['dist'].transform('min')\n    df[f'molecule_type_dist_std'] = df.groupby(['molecule_name', 'type'])['dist'].transform('std')\n    df[f'molecule_type_dist_std_diff'] = df[f'molecule_type_dist_std'] - df['dist']\n\n    df = reduce_mem_usage(df)\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndf_train = create_features(df_train)\ndf_test = create_features(df_test)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"good_columns = [\n'molecule_atom_index_0_dist_min',\n'molecule_atom_index_0_dist_max',\n'molecule_atom_index_1_dist_min',\n'molecule_atom_index_0_dist_mean',\n'molecule_atom_index_0_dist_std',\n'dist',\n'molecule_atom_index_1_dist_std',\n'molecule_atom_index_1_dist_max',\n'molecule_atom_index_1_dist_mean',\n'molecule_atom_index_0_dist_max_diff',\n'molecule_atom_index_0_dist_max_div',\n'molecule_atom_index_0_dist_std_diff',\n'molecule_atom_index_0_dist_std_div',\n'atom_0_couples_count',\n'molecule_atom_index_0_dist_min_div',\n'molecule_atom_index_1_dist_std_diff',\n'molecule_atom_index_0_dist_mean_div',\n'atom_1_couples_count',\n'molecule_atom_index_0_dist_mean_diff',\n'molecule_couples',\n'atom_index_1',\n'molecule_dist_mean',\n'molecule_atom_index_1_dist_max_diff',\n'molecule_atom_index_0_y_1_std',\n'molecule_atom_index_1_dist_mean_diff',\n'molecule_atom_index_1_dist_std_div',\n'molecule_atom_index_1_dist_mean_div',\n'molecule_atom_index_1_dist_min_diff',\n'molecule_atom_index_1_dist_min_div',\n'molecule_atom_index_1_dist_max_div',\n'molecule_atom_index_0_z_1_std',\n'y_0',\n'molecule_type_dist_std_diff',\n'molecule_atom_1_dist_min_diff',\n'molecule_atom_index_0_x_1_std',\n'molecule_dist_min',\n'molecule_atom_index_0_dist_min_diff',\n'molecule_atom_index_0_y_1_mean_diff',\n'molecule_type_dist_min',\n'molecule_atom_1_dist_min_div',\n'atom_index_0',\n'molecule_dist_max',\n'molecule_atom_1_dist_std_diff',\n'molecule_type_dist_max',\n'molecule_atom_index_0_y_1_max_diff',\n'molecule_type_0_dist_std_diff',\n'molecule_type_dist_mean_diff',\n'molecule_atom_1_dist_mean',\n'molecule_atom_index_0_y_1_mean_div',\n'molecule_type_dist_mean_div',\n'type']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in ['atom_index_0', 'atom_index_1', 'atom_1', 'type_0', 'type']:\n    if f in good_columns:\n        lbl = LabelEncoder()\n        lbl.fit(list(df_train[f].values) + list(df_test[f].values))\n        df_train[f] = lbl.transform(list(df_train[f].values))\n        df_test[f] = lbl.transform(list(df_test[f].values))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Selection"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Threshold for removing correlated variables\nthreshold = 0.95\n\n# Absolute value correlation matrix\ncorr_matrix = df_train.corr().abs()\n\n# Getting the upper triangle of correlations\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Select columns with correlations above threshold\nto_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n\nprint('There are %d columns to remove.' % (len(to_drop)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = df_train.drop(columns = to_drop)\ndf_test = df_test.drop(columns = to_drop)\n\nprint('Training shape: ', df_train.shape)\nprint('Testing shape: ', df_test.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Preprocessing "},{"metadata":{},"cell_type":"markdown","source":"Now, we will define our validation dataset and drop the unnecessary features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Split the 'features' and 'income' data into training and testing sets\nX_train, X_val, y_train, y_val = train_test_split(df_train.drop('scalar_coupling_constant',\n                                                                axis=1), \n                                                  df_train['scalar_coupling_constant'], \n                                                  test_size = 0.10, \n                                                  random_state = 0)\n\ndf_val = pd.DataFrame({\"type\":X_val[\"type\"]})\ndf_val['scalar_coupling_constant'] = y_val\n\nX_train = X_train.drop(['id', 'atom_1','atom_0', 'type','molecule_name'], axis=1).values\ny_train = y_train.values\n\nX_val = X_val.drop(['id', 'atom_1','atom_0', 'type','molecule_name'], axis=1).values\ny_val = y_val.values\n\nX_test = df_test.drop(['id', 'atom_1','atom_0', 'type','molecule_name'], axis=1).values\n\nprint(\"Training set has {} samples.\".format(X_train.shape[0]))\nprint(\"Validation set has {} samples.\".format(X_val.shape[0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Hyperopt \nHyperopt's job is to find the best value of a scalar-valued, possibly-stochastic function over a set of possible arguments to that function. It's a library for serial and parallel optimization over awkward search spaces, which may include real-valued, discrete, and conditional dimensions.\n- What the above means is that it is a optimizer that could minimize/maximize the loss function/accuracy(or whatever metric) for you.\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define searched space\nhyper_space = {'objective': 'regression',\n               'metric':'mae',\n               'boosting':'gbdt',\n               #'n_estimators': hp.choice('n_estimators', [25, 40, 50, 75, 100, 250, 500]),\n               'max_depth':  hp.choice('max_depth', [5, 8, 10, 12, 15]),\n               'num_leaves': hp.choice('num_leaves', [100, 250, 500, 650, 750, 1000,1300]),\n               'subsample': hp.choice('subsample', [.3, .5, .7, .8, 1]),\n               'colsample_bytree': hp.choice('colsample_bytree', [ .6, .7, .8, .9, 1]),\n               'learning_rate': hp.choice('learning_rate', [.1, .2, .3]),\n               'reg_alpha': hp.choice('reg_alpha', [.1, .2, .3, .4, .5, .6]),\n               'reg_lambda':  hp.choice('reg_lambda', [.1, .2, .3, .4, .5, .6]),               \n               'min_child_samples': hp.choice('min_child_samples', [20, 45, 70, 100])}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Defining the Metric to score our optimizer\ndef metric(df, preds):\n    df['diff'] = (df['scalar_coupling_constant'] - preds).abs()\n    return np.log(df.groupby('type')['diff'].mean().map(lambda x: max(x, 1e-9))).mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Creating the function that we will use to optimize the hyper parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"lgtrain = lightgbm.Dataset(X_train, label=y_train)\nlgval = lightgbm.Dataset(X_val, label=y_val)\n\ndef evaluate_metric(params):\n    \n    model_lgb = lightgbm.train(params, lgtrain, 500, \n                          valid_sets=[lgtrain, lgval], early_stopping_rounds=20, \n                          verbose_eval=500)\n\n    pred = model_lgb.predict(X_val)\n\n    score = metric(df_val, pred)\n    \n    print(score)\n \n    return {\n        'loss': score,\n        'status': STATUS_OK,\n        'stats_running': STATUS_RUNNING\n    }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Initiating the optimizer"},{"metadata":{"trusted":true,"_kg_hide-input":false,"_kg_hide-output":true},"cell_type":"code","source":"# Trail\ntrials = Trials()\n\n# Set algoritm parameters\nalgo = partial(tpe.suggest, \n               n_startup_jobs=-1)\n\n# Seting the number of Evals\nMAX_EVALS= 15\n\n# Fit Tree Parzen Estimator\nbest_vals = fmin(evaluate_metric, space=hyper_space, verbose=1,\n                 algo=algo, max_evals=MAX_EVALS, trials=trials)\n\n# Print best parameters\nbest_params = space_eval(hyper_space, best_vals)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hyper-parameters that we got in Hyperopt"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"BEST PARAMETERS: \" + str(best_params))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's use these parameters to train and predict "},{"metadata":{},"cell_type":"markdown","source":"## Predicting with the optimized params"},{"metadata":{"trusted":true},"cell_type":"code","source":"model_lgb = lightgbm.train(best_params, lgtrain, 4000, \n                      valid_sets=[lgtrain, lgval], early_stopping_rounds=30, \n                      verbose_eval=500)\n\nlgb_pred = model_lgb.predict(X_test)\n\ndf_test['scalar_coupling_constant'] = lgb_pred\n\ndf_test[['id', 'scalar_coupling_constant']].to_csv(\"molecular_struct_sub.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# I will keep working on this dataset exploration\n## If you liked votesup the kernel "},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"Some references:<br>\nhttps://www.kaggle.com/inversion/atomic-distance-benchmark<br>\nhttps://www.kaggle.com/tunguz/atomic-distances-with-h2o-automl<br>\nhttps://www.kaggle.com/artgor/molecular-properties-eda-and-models<br>\nhttps://www.kaggle.com/hrmello/is-type-related-to-scalar-coupling<br>"}],"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}