{"cells":[{"metadata":{"_uuid":"20be50e6-eafa-400a-b662-be25cdec10f0","_cell_guid":"92778900-a124-417f-af33-dbded344b482","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!conda install -y -c rdkit rdkit;\nfrom rdkit import Chem\n!conda install -y -c openbabel openbabel;\n#import openbabel as ob\n#import pybel as pb","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35f4d157-a002-4f0b-b904-c8613bd36e99","_cell_guid":"b8464036-2d98-44bc-ba23-1c077183c199","trusted":true},"cell_type":"code","source":"# import data\ntrain = pd.read_csv(\"../input/train.csv\")\ntest = pd.read_csv(\"../input/test.csv\")\n#structures = pd.read_csv(\"../input/structures.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame().to_csv('magic.csv',index=False)\nimport pybel as pb","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1abc68e1-a252-4055-b9f2-11b998b315d7","_cell_guid":"cd9371d6-1658-4c30-b012-9410e7f2a1a3","trusted":true},"cell_type":"markdown","source":"We can extract features from the coupling type. See my previous kernel on [here](https://www.kaggle.com/edeanf/190621-predmolprop-exploration)."},{"metadata":{"_uuid":"ccdac487-0f73-4fa4-98d8-e84ab13a4eb4","_cell_guid":"77bd36fa-41f4-4b77-95c7-14e1679b8b5e","trusted":true},"cell_type":"code","source":"def PrintDataframe(df):\n    print(df.to_string(),'\\n')\n    \ndef DecomposeType(df):\n    df['num_bonds'] = df.type.map(lambda s: int(s[0]))\n    df['atom_end_type'] = df.type.map(lambda s: s[-1])\n\nDecomposeType(train)\nDecomposeType(test)\nPrintDataframe(train.head())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"31086061-74a8-4330-9743-3fe0244080a6","_cell_guid":"ba84013f-cb6a-4444-9d4c-f6e28dd5cf5e","trusted":true},"cell_type":"markdown","source":"Although there are some features that can be extracted from the structures.csv, using the .xyz files is more powerful. I tried using the following api from jensengroup [github](https://github.com/jensengroup/xyz2mol) in order to import .xyz files into RDkit's Mol object, but there are errors. Jensen also [mentions](https://proteinsandwavefunctions.blogspot.com/2018/01/xyz2mol-converting-xyz-file-to-rdkit.html) that one can use openbabel to convert to smile and then convert to rdkit Mol object, but the ordering of the atoms can change, which is bad since we need to be able to reference in original indexing. So I wrote my own code which uses open babel"},{"metadata":{"trusted":true},"cell_type":"code","source":"from collections import Counter\n\ndef ProcessXYZ(filepath):\n    '''\n    reads a filepath and returns RDKit Mol Object, 3D Distance Matrix, and Adjacency Matrix\n    '''\n    \n    # load as openbabel molecule\n    pbMol = pb.readfile(\"xyz\", filepath).__next__()\n    #obMol = pbMol.OBMol\n    num_atoms = len(pbMol.atoms)#obMol.NumAtoms()\n    \n    # extract atomic properties\n    rwMol = Chem.RWMol()\n    conf = Chem.Conformer(num_atoms)\n    AtomCounts = Counter()\n    HybOrb = []\n    atoms =[]\n    PCList = []\n    atom_coord = []\n    atom_spins = []\n    atom_valence = []\n    atom_fc = []\n    for i,atom in enumerate(pbMol.atoms):#ob.OBMolAtomIter(obMol):\n        rwMol.AddAtom(Chem.Atom(atom.atomicnum))\n        conf.SetAtomPosition(i, atom.coords)\n        atom_coord.append(np.array(atom.coords))\n        AtomCounts = AtomCounts+Counter(atom.type)\n        atoms.append(atom.type)\n        HybOrb.append(int(atom.hyb))\n        PCList.append(atom.partialcharge)\n        atom_spins.append(int(atom.spin))\n        atom_valence.append(int(atom.valence))\n        atom_fc.append(atom.formalcharge)\n    # make RDKit Mol Object\n    mol = rwMol.GetMol()\n    mol.AddConformer(conf)\n    \n    # Get 3D Distance Matrix\n    DistMat = np.array(Chem.Get3DDistanceMatrix(mol))\n    \n    # Now make adjacency matrix\n    # note that Atom.GetIdx() appears to use one-based indexing\n    AdjMat = np.zeros((num_atoms, num_atoms)).astype(int)\n    \n    def PBMolBondIter(pbMol):\n        for i in range(pbMol.OBMol.NumBonds()):\n            yield pbMol.OBMol.GetBond(i)\n    \n    for bond in PBMolBondIter(pbMol):#ob.OBMolBondIter(obMol):\n        i,j = bond.GetBeginAtomIdx()-1, bond.GetEndAtomIdx()-1\n        BO = bond.GetBO()\n        AdjMat[i][j]=BO\n        AdjMat[j][i]=BO\n    return DistMat, AdjMat, atoms, HybOrb, pbMol.OBMol.NumBonds(),PCList, atom_coord, atom_spins, atom_valence, atom_fc","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"753587c7-d5eb-445d-afe6-516bc789b348","_cell_guid":"af034e9e-6a47-4a44-9052-9e488b317d47","trusted":true},"cell_type":"markdown","source":"Because a single molecule is used multiple times, it is more straightforward if not more efficienct to iterate through the Dataframe instead of using assign. I imagine the efficiency depends on how expensive it is to load the .xyz file repeatedly and how much space it would take to store all the molecules in memory. In any case, we can optimize code later. Let's just get it to work."},{"metadata":{"trusted":true},"cell_type":"code","source":"def unit_vector(vector):\n    \"\"\" Returns the unit vector of the vector \"\"\"\n    return vector / np.linalg.norm(vector)\n\ndef angle_between(v1, v2):\n    \"\"\" Returns the angle in radians between vectors 'v1' and 'v2' \"\"\"\n    v1_u = unit_vector(v1)\n    v2_u = unit_vector(v2)\n    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) * 180/np.pi\n\ndef angle_2bond(atom1,atom2,atom3):\n    vect12 = atom1-atom2\n    vect32 = atom3-atom2\n    return angle_between(vect32,vect12)\n\ndef dihedral_angle(atom1, atom2, atom3, atom4):\n    \"\"\"\n    Returns the dihedral angle between atoms 1 and 4\n    given the x, y, z coordinates of all four atoms\n    \"\"\"\n    # First calculate the vectors between different atom combinations\n    vec_12 = atom1-atom2\n    vec_23 = atom2-atom3\n    vec_34 = atom3-atom4\n\n    # Then calculate the normal vector for the plane that contains bonds 2-3 and 3-4\n    norm = np.cross(vec_23, vec_34)\n\n    # Now calculate the angle between the normal vector and vec_12\n    # Subtract 90 deg to get the angle between the plane and vec_12\n    return  angle_between(vec_12,norm) - 90\n\n\n# Dictionary of distance features for each atom\ndef GetFeatureDict(DistMat,atoms):\n    feature_dict ={}\n    temp=DistMat[np.nonzero(DistMat)]\n    feature_dict['min']=temp.min()\n    feature_dict['mean']=temp.mean()\n    feature_dict['max']=temp.max()\n    temp=DistMat[np.nonzero(DistMat*[atom[0]=='C' for atom in atoms])]\n    if temp.size>0:\n        feature_dict['Cmin']=temp.min()\n        feature_dict['Cmean']=temp.mean()\n        feature_dict['Cmax']=temp.max()\n    else:\n        feature_dict['Cmin']=-1\n        feature_dict['Cmean']=-1\n        feature_dict['Cmax']=-1\n    temp=DistMat[np.nonzero(DistMat*[atom[0]=='O' for atom in atoms])]\n    if temp.size>0:\n        feature_dict['Omin']=temp.min()\n        feature_dict['Omean']=temp.mean()\n        feature_dict['Omax']=temp.max()\n    else:\n        feature_dict['Omin']=-1\n        feature_dict['Omean']=-1\n        feature_dict['Omax']=-1\n    temp=DistMat[np.nonzero(DistMat*[atom[0]=='N' for atom in atoms])]\n    if temp.size>0:\n        feature_dict['Nmin']=temp.min()\n        feature_dict['Nmean']=temp.mean()\n        feature_dict['Nmax']=temp.max()\n    else:\n        feature_dict['Nmin']=-1\n        feature_dict['Nmean']=-1\n        feature_dict['Nmax']=-1\n    temp=DistMat[np.nonzero(DistMat*[atom[0]=='F' for atom in atoms])]\n    if temp.size>0:\n        feature_dict['Fmin']=temp.min()\n        feature_dict['Fmean']=temp.mean()\n        feature_dict['Fmax']=temp.max()\n    else:\n        feature_dict['Fmin']=-1\n        feature_dict['Fmean']=-1\n        feature_dict['Fmax']=-1\n    return feature_dict","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"02d02c5f-b94e-40cb-a408-370996286e0f","_cell_guid":"2eef820d-6641-4b86-bc01-97f41e5d1855","trusted":true},"cell_type":"code","source":"import scipy.sparse as sparse\ndef ExtendDataframe(df):\n    # initialize the columns\n    num_mol_bonds =[]\n    min_d = []\n    mean_d = []\n    max_d = []\n    space_dr=[]\n    bond_dr=[]\n    bond_1=[]\n    bond_2=[]\n    bond_3=[]\n    \n    atom_0_pc=[]\n    atom_end_pc=[]\n    atom_0_fc=[]\n    atom_end_fc=[]\n    atom_0_val=[]\n    atom_end_val=[]\n    atom_0_sm=[]\n    atom_end_sm=[]\n    \n    atom_0_type2=[]\n    atom_2_type=[]\n    atom_3_type=[]\n    atom_end_type2=[]\n    atom_2_hyb=[]\n    atom_3_hyb=[]\n    atom_end_hyb=[]\n    path_count = []\n    \n    atom_0_min = []\n    atom_0_mean = []\n    atom_0_max = []\n    atom_0_Cmin = []\n    atom_0_Cmean = []\n    atom_0_Cmax = []\n    atom_0_Omin = []\n    atom_0_Omean = []\n    atom_0_Omax = []\n    atom_0_Nmin = []\n    atom_0_Nmean = []\n    atom_0_Nmax = []\n    atom_0_Fmin = []\n    atom_0_Fmean = []\n    atom_0_Fmax = []\n\n    atom_end_min = []\n    atom_end_mean = []\n    atom_end_max = []\n    atom_end_Cmin = []\n    atom_end_Cmean = []\n    atom_end_Cmax = []\n    atom_end_Omin = []\n    atom_end_Omean = []\n    atom_end_Omax = []\n    atom_end_Nmin = []\n    atom_end_Nmean = []\n    atom_end_Nmax = []\n    atom_end_Fmin = []\n    atom_end_Fmean = []\n    atom_end_Fmax = []\n\n    bond2_angle = []\n    bond3_angle = []\n    \n    # initialize the molecule currently loaded\n    loaded_name = ''\n    \n    # Depth First Search for looking up bonds between atom_0 and atom_1\n    def DFS(path):\n        # a path is a list of atom indices starting from atom_0 which is at index i\n        # only if the (num_bonds+1)th atom has index j is the path saved into paths\n        # where j is the index of atom_1\n        m = path[-1]\n        if len(path) == num_bonds+1:\n            if m == j:\n                paths.append(path.copy())\n        else:\n            for _,n in zip(*SparseAdjMat[m].nonzero()):\n                if n not in path:\n                    path.append(n)\n                    DFS(path)\n                    path.pop()\n    \n    for row in df.itertuples(index=False):\n        mol_name = row.molecule_name\n        i = row.atom_index_0\n        j = row.atom_index_1\n        num_bonds = row.num_bonds\n        \n        # construct molecule properties\n        if loaded_name != mol_name:\n            # load new molecule\n            loaded_name = mol_name\n            filepath = '../input/structures/'+loaded_name+'.xyz'\n            #mol,DistMat,AdjMat,AtomCounts,NumAtoms,atoms,HybOrb = ProcessXYZ(filepath)\n            DistMat,AdjMat,atoms,HybOrb,bond_count,PCList,atom_coord,atom_spins, atom_val, atom_fc = ProcessXYZ(filepath)\n            SparseAdjMat=sparse.csr_matrix(AdjMat)\n            \n            num_mol_bonds.append(bond_count)\n            min_d.append(DistMat[np.nonzero(DistMat)].min())\n            mean_d.append(DistMat[np.nonzero(DistMat)].mean())\n            max_d.append(DistMat.max())\n            atom_dist_dict={}\n        else:\n            num_mol_bonds.append(num_mol_bonds[-1])\n            min_d.append(min_d[-1])\n            mean_d.append(mean_d[-1])\n            max_d.append(max_d[-1])\n        \n        # through-space distance and atom properties\n        space_dr.append(DistMat[i][j])\n        atom_0_type2.append(atoms[i])\n        atom_end_type2.append(atoms[j])\n        atom_0_pc.append(PCList[i])\n        atom_end_pc.append(PCList[j])\n        \n        atom_0_fc.append(atom_fc[i])\n        atom_end_fc.append(atom_fc[j])\n        atom_0_val.append(atom_val[i])\n        atom_end_val.append(atom_val[j])\n        atom_0_sm.append(atom_spins[i])\n        atom_end_sm.append(atom_spins[j])\n        \n        # find path between the atom_0 and atom_1\n        paths = []\n        if num_bonds>1:\n            DFS([i])\n            path_count.append(len(paths))\n        \n        if num_bonds==1:\n            bond_dr.append(DistMat[i][j])\n            bond_1.append(AdjMat[i][j])\n            bond_2.append(-1)\n            bond_3.append(-1)\n            atom_2_type.append('')\n            atom_3_type.append('')\n            atom_2_hyb.append(-1)\n            atom_3_hyb.append(-1)\n            atom_end_hyb.append(HybOrb[j])\n            path_count.append(1)\n            bond2_angle.append(-1)\n            bond3_angle.append(-1)\n        elif num_bonds==2:\n            path = paths[0]\n            bond_dr.append(DistMat[path[0]][path[1]]+DistMat[path[1]][path[2]])\n            bond_1.append(AdjMat[path[0]][path[1]])\n            bond_2.append(AdjMat[path[1]][path[2]])\n            bond_3.append(-1)\n            atom_2_type.append(atoms[path[1]])\n            atom_3_type.append('')\n            atom_2_hyb.append(HybOrb[path[1]])\n            atom_3_hyb.append(-1)\n            atom_end_hyb.append(HybOrb[j])\n            bond2_angle.append(abs(angle_2bond(atom_coord[path[0]],atom_coord[path[1]],atom_coord[path[2]])))\n            bond3_angle.append(-1)\n        else:\n            path = paths[0]\n            bond_dr.append(DistMat[path[0]][path[1]]+DistMat[path[1]][path[2]]+DistMat[path[2]][path[3]])\n            bond_1.append(AdjMat[path[0]][path[1]])\n            bond_2.append(AdjMat[path[1]][path[2]])\n            bond_3.append(AdjMat[path[2]][path[3]])\n            atom_2_type.append(atoms[path[1]])\n            atom_3_type.append(atoms[path[2]])\n            atom_2_hyb.append(HybOrb[path[1]])\n            atom_3_hyb.append(HybOrb[path[2]])\n            atom_end_hyb.append(HybOrb[j])\n            bond2_angle.append(-1)\n            bond3_angle.append(float(abs(dihedral_angle(atom_coord[path[0]],atom_coord[path[1]],atom_coord[path[2]],atom_coord[path[3]]))))\n        \n        # distance features for atom end\n        if(i not in atom_dist_dict):\n            atom_dist_dict[i]=GetFeatureDict(DistMat[i][:],atoms)\n        feature_list = [\n            ('min',atom_0_min),('mean',atom_0_mean),('max',atom_0_max),\n            ('Cmin',atom_0_Cmin),('Cmean',atom_0_Cmean),('Cmax',atom_0_Cmax),\n            ('Omin',atom_0_Omin),('Omean',atom_0_Omean),('Omax',atom_0_Omax),\n            ('Nmin',atom_0_Nmin),('Nmean',atom_0_Nmean),('Nmax',atom_0_Nmax),\n            ('Fmin',atom_0_Fmin),('Fmean',atom_0_Fmean),('Fmax',atom_0_Fmax)\n        ]\n        for name, feature in feature_list:\n            feature.append(atom_dist_dict[i][name])\n        \n        if(j not in atom_dist_dict):\n            atom_dist_dict[j]=GetFeatureDict(DistMat[j][:],atoms)    \n        feature_list = [\n            ('min',atom_end_min),('mean',atom_end_mean),('max',atom_end_max),\n            ('Cmin',atom_end_Cmin),('Cmean',atom_end_Cmean),('Cmax',atom_end_Cmax),\n            ('Omin',atom_end_Omin),('Omean',atom_end_Omean),('Omax',atom_end_Omax),\n            ('Nmin',atom_end_Nmin),('Nmean',atom_end_Nmean),('Nmax',atom_end_Nmax),\n            ('Fmin',atom_end_Fmin),('Fmean',atom_end_Fmean),('Fmax',atom_end_Fmax)\n        ]\n        for name, feature in feature_list:\n            feature.append(atom_dist_dict[j][name])\n        \n    cols = {'num_mol_bonds': num_mol_bonds,\n            'min_d':min_d,\n            'mean_d':mean_d,\n            'max_d':max_d,\n            'space_dr': space_dr,\n            'bond_dr': bond_dr,\n            'bond_1': bond_1,\n            'bond_2': bond_2,\n            'bond_3': bond_3,\n            'atom_0_pc': atom_0_pc,\n            'atom_end_pc': atom_end_pc,\n            'atom_0_fc': atom_0_fc,\n            'atom_end_fc': atom_end_fc,\n            'atom_0_val': atom_0_val,\n            'atom_end_val': atom_end_val,\n            'atom_0_sm': atom_0_sm,\n            'atom_end_sm':atom_end_sm,\n            'atom_0_type2': atom_0_type2,\n            'atom_2_type': atom_2_type,\n            'atom_3_type': atom_3_type,\n            'atom_end_type2': atom_end_type2,\n            'atom_2_hyb': atom_2_hyb,\n            'atom_3_hyb': atom_3_hyb,\n            'atom_end_hyb': atom_end_hyb,\n            'path_count': path_count,\n            'atom_0_min': atom_0_min,\n            'atom_0_mean': atom_0_mean,\n            'atom_0_max': atom_0_max,\n            'atom_0_Cmin': atom_0_Cmin,\n            'atom_0_Cmean': atom_0_Cmean,\n            'atom_0_Cmax': atom_0_Cmax,\n            'atom_0_Omin': atom_0_Omin,\n            'atom_0_Omean': atom_0_Omean,\n            'atom_0_Omax': atom_0_Omax,\n            'atom_0_Nmin': atom_0_Nmin,\n            'atom_0_Nmean': atom_0_Nmean,\n            'atom_0_Nmax': atom_0_Nmax,\n            'atom_0_Fmin': atom_0_Fmin,\n            'atom_0_Fmean': atom_0_Fmean,\n            'atom_0_Fmax': atom_0_Fmax,\n            'atom_end_min': atom_end_min,\n            'atom_end_mean': atom_end_mean,\n            'atom_end_max': atom_end_max,\n            'atom_end_Cmin': atom_end_Cmin,\n            'atom_end_Cmean': atom_end_Cmean,\n            'atom_end_Cmax': atom_end_Cmax,\n            'atom_end_Omin': atom_end_Omin,\n            'atom_end_Omean': atom_end_Omean,\n            'atom_end_Omax': atom_end_Omax,\n            'atom_end_Nmin': atom_end_Nmin,\n            'atom_end_Nmean': atom_end_Nmean,\n            'atom_end_Nmax': atom_end_Nmax,\n            'atom_end_Fmin': atom_end_Fmin,\n            'atom_end_Fmean': atom_end_Fmean,\n            'atom_end_Fmax': atom_end_Fmax,\n            'bond2_angle': bond2_angle,\n            'bond3_angle': bond3_angle\n            }\n    return df.assign(**cols)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ab1e58fb-cbc6-4a9a-a7bf-fb94b345627f","_cell_guid":"6a6927cc-9c18-4e14-a667-512c451101d7","trusted":true},"cell_type":"code","source":"#train_sample = train.sample(10000)\nimport time\nstart = time.time()\n#train_extend = ExtendDataframe(train_sample)\ntrain_extend = ExtendDataframe(train)\ntest_extend = ExtendDataframe(test)\nend = time.time()\nprint((end-start)/60)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend.atom_0_type2.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend.atom_end_type2.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend[train_extend.columns[:20]].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d0993e5b-2ac6-4fb9-a2f1-f09602d98373","_cell_guid":"395e92de-8120-4a8b-94e2-22c84cd3a3fc","trusted":true},"cell_type":"code","source":"train_extend[train_extend.columns[6:21]].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend[train_extend.columns[21:36]].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend.atom_end_sm.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend[train_extend.columns[36:51]].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_extend[train_extend.columns[51:]].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_extend[train_extend.bond2_angle!=-1].bond2_angle)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_extend[train_extend.bond3_angle!=-1].bond3_angle)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e0691e9b-9029-4e17-951d-e695236f9fbc","_cell_guid":"46033990-0764-4c4c-855f-58e7bbc3239a","trusted":true},"cell_type":"code","source":"# save for model building\ntrain_extend.to_csv('train_extend.csv', index=False)\ntest_extend.to_csv('test_extend.csv', index=False)","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}