{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport os\nimport pandas as pd\n\nDATA_DIR = '../input/'\nprint(os.listdir(DATA_DIR))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(DATA_DIR + 'train.csv', index_col=0)\ntest_df = pd.read_csv(DATA_DIR + 'test.csv')\nstructures_df = pd.read_csv(DATA_DIR + 'structures.csv')\nsc_contri_df = pd.read_csv(DATA_DIR + 'scalar_coupling_contributions.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mol = train_df.molecule_name.unique()\ntest_mol =  test_df.molecule_name.unique()\nstructure_mol =  structures_df.molecule_name.unique()\nprint('Train unique molecules', len(train_mol))\nprint('Test unique molecules', len(test_mol))\nprint('Structures unique molecules', len(structure_mol))\nprint('Common unique molecules', len(set(train_mol).intersection(set(test_mol))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('1st atom index types', len(train_df.atom_index_0.unique()))\nprint('2nd atom index types', len(train_df.atom_index_1.unique()))\nprint('All atom indices', len(set(train_df.atom_index_0.unique()).union(train_df.atom_index_1.unique())))\nprint(set(train_df.atom_index_0.unique()).union(train_df.atom_index_1.unique()))\nprint('Atom types', structures_df.atom.unique())\n\nprint('Molecule type ', test_df['type'].unique())\nprint('Molecule types are same in train,test', set(train_df['type'].unique()) == set(test_df['type'].unique()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for contri_col in ['fc','sd','pso','dso']:\n    train_df[contri_col] = sc_contri_df[contri_col]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Count of different molecule type"},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(figsize=(10,10))\nvc_train = train_df['type'].value_counts().to_frame('train')\nvc_test =test_df['type'].value_counts().to_frame('test')\npd.concat([vc_train, vc_test],axis=1).plot(ax=ax)\nax.set_xticklabels(['0'] + vc_train.index.tolist(), fontsize=10)\nax.legend()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## All 4 contributions sum up to the target\n Maximum difference is of the order of 1e-3 which is much lower than std of the columns. We can therefore train 4 models instead of training just one model. For tree based models it will be helpful. For linear regression it will be useless as the weights of the four models can be summed to get the weights of one model which gives same performance."},{"metadata":{"trusted":true},"cell_type":"code","source":"(sc_contri_df[['fc','sd','pso','dso']].sum(axis=1) - train_df.scalar_coupling_constant).describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sc_contri_df[['fc','sd','pso','dso']].describe().loc[['mean','std','min','max']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Correlation plot between different target components"},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(figsize=(10,10))\nsns.heatmap(sc_contri_df[['fc','sd','pso','dso']].corr(), annot=True, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Relationship between atom and atom index\nOne information is that all indices >= 14 are mostly H. Not much useful though"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_df = structures_df.groupby('atom')['atom_index'].apply(lambda x: x.unique())\ndf = pd.DataFrame(index=['C','F','H','N','O'], columns=list(range(29)))\ndf.index.name = 'atom'\ndf.columns.name ='atom_index'\ndf[:] = 0\ndf = df.stack().to_frame('present').reset_index(level=1)\n\nfor atom in set(df.index.values):\n    df.loc[(df.index ==atom) &(df.atom_index.isin(temp_df.loc[atom])), 'present'] = 1\n    \n_, ax= plt.subplots(figsize=(10,5))\nsns.heatmap(df.set_index('atom_index',append=True)['present'].unstack())\nax.set_title('Relationship between atom and atom_index')\ndel df, temp_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Target variation"},{"metadata":{"trusted":true},"cell_type":"code","source":"Y_COL ='scalar_coupling_constant'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"###  wrt to atom indices\nThere is no particular difference in target variable wrt atom_index_0, \n#### wrt atom_index_1, there is some pattern"},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(nrows=2, figsize=(10, 20))\nsns.violinplot(x='atom_index_0', y=Y_COL, data=train_df, ax=ax[0])\nsns.violinplot(x='atom_index_1', y=Y_COL, data=train_df, ax=ax[1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## wrt molecule type\nThere is some information here"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x='type', y=Y_COL, data=train_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### wrt molecule type and atom_index_1"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder\ntrain_df['type_enc'] = LabelEncoder().fit_transform(train_df['type'])\ntrain_df['idx1_type'] = train_df['atom_index_1'] + train_df['type_enc']*train_df['atom_index_1'].max()\n_, ax = plt.subplots(figsize=(20,5))\nsns.violinplot(x='idx1_type', y=Y_COL, data=train_df,ax=ax)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## In train df, for which bonds do we have to estimate scalar_coupling_constant?\nIt forms pretty nice graphs to visualize. I'm however yet to get something out from them into my model :D "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Taken from https://www.kaggle.com/artgor/molecular-properties-eda-and-models\nimport networkx as nx\nfig, ax = plt.subplots(figsize = (20, 12))\nfor i, t in enumerate(train_df['type'].unique()):\n    train_df_type = train_df.loc[train_df['type'] == t]\n    bad_atoms_0 = list(train_df_type['atom_index_0'].value_counts(normalize=True)[train_df_type['atom_index_0'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms_1 = list(train_df_type['atom_index_1'].value_counts(normalize=True)[train_df_type['atom_index_1'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms = list(set(bad_atoms_0 + bad_atoms_1))\n    train_df_type = train_df_type.loc[(train_df_type['atom_index_0'].isin(bad_atoms_0) == False) & (train_df_type['atom_index_1'].isin(bad_atoms_1) == False)]\n    G = nx.from_pandas_edgelist(train_df_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":"## How atom index varies on train/test data.\nTrain and test are very similar in this. In general, across many different features, train data and test data are very similar. It is a good thing for us. It is no wonder that performance on validation closely matches that of leaderboard :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(ncols=2, nrows=2, figsize=(10,10))\nsns.violinplot(x='atom_index_0', data=train_df, ax=ax[0,0])\nax[0,0].set_title('Train')\nsns.violinplot(x='atom_index_1', data=train_df, ax=ax[1,0])\nax[1,0].set_title('Train')\nsns.violinplot(x='atom_index_0', data=test_df, ax=ax[0,1])\nax[0,1].set_title('Test')\nsns.violinplot(x='atom_index_1', data=test_df, ax=ax[1,1])\nax[1,1].set_title('Test')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## What is the meaning of atom index. Is it correct to say that for every molecule, atom index is just a numbering of its atoms?\nWe see below that for each molecule, atom_index is just a way to number the atoms. However, it is not random. \n#### C atoms are sequentially numbered according to their position in the carbon chain. H atoms numbering starts after backbone structure have been numbered."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_molecule(mol_name, figsize=(10,10)):\n    from mpl_toolkits.mplot3d import Axes3D\n    temp_df = structures_df[structures_df.molecule_name == mol]\n\n    marker_size = {'C':120,'H':10,'N':70,'F':90,'O':80}\n    marker_color ={'C':'blue','H':'orange','N':'green','F':'black','O':'violet'}\n    fig = plt.figure(figsize=figsize)\n    ax = fig.add_subplot(111, projection='3d')\n    for atom in temp_df.atom.unique():\n        dta = temp_df[temp_df.atom == atom]\n        ax.scatter(dta.x, dta.y, dta.z, s=marker_size[atom], c=marker_color[atom],label=atom)\n        for _, row in dta.iterrows():\n            ax.text(row.x, row.y, row.z,row.atom_index,fontsize=20)\n    ax.legend()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that for each molecule atom_index starts from 0 and goes to some value covering all integers in between. It is not that some integer is missing."},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_df = structures_df.groupby('molecule_name')['atom_index'].agg({\n    'sum':'sum',\n    'count':'count'})\ntemp_df['all_present_sum'] = temp_df['count']*(temp_df['count'] -1)/2\nassert(temp_df['all_present_sum'] == temp_df['sum']).all()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mol = 'dsgdb9nsd_121674'\nplot_molecule(mol,figsize=(20,20))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## How large are the molecules.\nMost molecules comprise of about 10-25 atoms"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_df = structures_df.molecule_name.value_counts()\nprint(temp_df.describe())\ntemp_df.hist()\ndel temp_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Some information can be extracted if we just target sign and target magnitude separately"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\ntrain_df['target_sign'] = (train_df['scalar_coupling_constant'] > 0).astype(np.int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_df =pd.merge(train_df, structures_df, how='left',left_on=['molecule_name','atom_index_0'],\n         right_on=['molecule_name','atom_index'])\n\ntemp_df =pd.merge(temp_df, structures_df, how='left',left_on=['molecule_name','atom_index_1'],\n         right_on=['molecule_name','atom_index'])\nprint(temp_df.groupby('atom_y')['target_sign'].describe()[['mean','std','min','25%']])\ndel temp_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## For one molecule, how many rows are in train/test data.\nIt might make sense to stratify wrt molecules as is suggested in some kernels as we see lot of variation. In my case however, I haven't got much benefit till now doing it.\nHere again, to our good luck, train is very similar to test. "},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(figsize=(20,5))\nsns.distplot(train_df.molecule_name.value_counts(),ax=ax,label='train')\nsns.distplot(test_df.molecule_name.value_counts(),ax=ax, label='test')\nax.set_xlabel('Number of rows')\nax.legend()","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":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"nbformat":4,"nbformat_minor":1}