{"cells":[{"metadata":{},"cell_type":"markdown","source":"> \n**Some features containing the atom_indices have high importance (but they should not be relevant!?)**\n\n* Usually only heavy atoms are stored by cheminformatics software. Explicit Hydrogens are added\nafter the last heavy atoms and have a larger atom index.\n\n* It seems they (i.e. atom_index_0) contain some info on minimum and average molecule size...\n"},{"metadata":{},"cell_type":"markdown","source":"Open babel code thanks to this kernel: https://www.kaggle.com/asauve/v7-estimation-of-mulliken-charges-with-open-babel"},{"metadata":{"trusted":true},"cell_type":"code","source":"!conda install -y -c openbabel openbabel \nimport openbabel as ob","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport openbabel as ob\nobConversion = ob.OBConversion()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain = pd.read_csv('../input/train.csv')\nstructures = pd.read_csv('../input/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain['type_0'] = Xtrain['type'].apply(lambda x: x[0]).astype(int)\nXtrain['type_1'] = Xtrain['type'].apply(lambda x: x[2:3])\nXtrain['type_2'] = Xtrain['type'].apply(lambda x: x[3:4])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"atom_index_1 (the heavy atoms only except for the HH coupling) have lower indices "},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x=\"type_2\", y=\"atom_index_1\",\n               split=True, data=Xtrain)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"atom_index_0 are connected to higher indices "},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x=\"type_1\", y=\"atom_index_0\",\n               split=True, data=Xtrain)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef read_ob_molecule(molecule_name, datadir=\"../input/structures\"):\n    mol = ob.OBMol()\n    path = f\"{datadir}/{molecule_name}.xyz\"\n    if not obConversion.ReadFile(mol, path):\n        raise FileNotFoundError(f\"Could not read molecule {path}\")\n    return mol","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now creating some open babel data on number of atoms"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Xtrain = Xtrain.iloc[:100000] # for speedup\nXtrain = Xtrain.sample(100000) # for speedup\ndf_structures = pd.read_csv('../input/structures.csv')\nstructures_idx = df_structures.set_index([\"molecule_name\"])\nmolecule_names = Xtrain.molecule_name.unique()\nob_natoms = []\nob_heavyatoms = []\nob_molweight = []\nob_molecule_name = []\nob_atom_index = []\nfor i,molecule_name in enumerate(molecule_names):\n    if i%10000 ==0:\n        print(\"OB molecule %d\"%(i))\n        \n    # fill data for output DF\n    ms = structures_idx.loc[molecule_name].sort_index()\n    natoms = len(ms)\n    ob_molecule_name.extend([molecule_name] * natoms)\n    ob_atom_index.extend(ms.atom_index.values)\n    \n    # calculate open babel charge for each method\n    mol = read_ob_molecule(molecule_name)\n    assert (mol.NumAtoms() == natoms)\n    mw = mol.GetMolWt()\n    ob_heavyatoms.extend([mol.NumHvyAtoms()]*natoms)\n    ob_natoms.extend([natoms]*natoms)\n    ob_molweight.extend([mw]*natoms)\n    \nob_data = pd.DataFrame({'molecule_name':ob_molecule_name,'natoms':ob_natoms,'nheavyatoms':  ob_heavyatoms, 'mw':ob_molweight,'atom_index':ob_atom_index})\nob_data.head(10)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob = pd.merge(Xtrain, ob_data, how='left',\n                        left_on=['molecule_name', 'atom_index_0'],\n                        right_on=['molecule_name', 'atom_index'],suffixes=('_at1', '_at2'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x=\"atom_index_0\", y=\"natoms\",\n               split=True, data=Xtrain_ob)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is some correlation with the average number of atoms at least for larger indices. "},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob['atom_index_sum'] = Xtrain_ob['atom_index_0']+Xtrain_ob['atom_index_1']\nXtrain_ob['atom_index_diff'] = Xtrain_ob['atom_index_0']-Xtrain_ob['atom_index_1']\nXtrain_ob['atom_index_abs'] = np.fabs(Xtrain_ob['atom_index_0']-Xtrain_ob['atom_index_1'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_atom_info(df, structures, 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob = Xtrain_ob.drop('atom_index', axis=1)\nXtrain_ob.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nXtrain_ob = map_atom_info(Xtrain_ob, structures, 0)\nXtrain_ob = map_atom_info(Xtrain_ob, structures, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_p_0 = Xtrain_ob[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = Xtrain_ob[['x_1', 'y_1', 'z_1']].values\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain_ob['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x=\"atom_index_sum\", y=\"natoms\",\n               split=True, data=Xtrain_ob)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For JHH there is some correlation with average natoms?"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.violinplot(x=\"atom_index_sum\", y=\"natoms\",\n               split=False, data=Xtrain_ob.loc[Xtrain_ob.type=='2JHH'])","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}