{"cells":[{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"},"cell_type":"markdown","source":"# Graph Based feature generation - Walk along the atomic bonds, get data\n## My first public kernal ever - upvote if you like it and surely leave a comment to improve it\n### Features generated\n* Atoms along the path\n* Atomic bonds for all the atoms along the path\n* Count of C, N, H, F and Os bonded to the atoms along the path\n\n### Next steps\n* Very slow at calculations\n* 90k unique molecules will take at least 10 hours for calculation\n* Need to test on a sample of 10-20% data and see about features importances\n* Wonder if anymore parallel processing speed can be improved\n\n### Inspirations\nhttps://www.kaggle.com/mykolazotko/3d-visualization-of-molecules-with-plotly"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport datetime\nimport time\n\nimport plotly.graph_objs as go\nfrom plotly.offline import init_notebook_mode, iplot\nfrom sympy.geometry import Point3D\n\nimport random\nimport os\n\nimport networkx as nx\n\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\n\nstructures = reduce_mem_usage(pd.read_csv('../input/structures.csv'))\ntrain_df = reduce_mem_usage(pd.read_csv('../input/train.csv'))\ntrain_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"molecule_list = [mol for mol in train_df['molecule_name'].drop_duplicates()]\nnew_list = random.sample(molecule_list, 9000)\nlen(new_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = structures[structures['molecule_name'].isin(new_list)]\ntrain_df = train_df[train_df['molecule_name'].isin(new_list)]\ntrain_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Add Fudge factor\natomic_radii = dict(C=0.77, F=0.71, H=0.38, N=0.75, O=0.73)\nfudge_factor = 0.05\natomic_radii = {k:v + fudge_factor for k,v in atomic_radii.items()}\n\ndef Build_Graph(molecule_name):\n    ## Get molecule structure\n    molecule = structures[structures.molecule_name == molecule_name]\n    \n    ## List of atom names\n    elements = molecule.atom.tolist()\n    \n    ## Get all xyz cordinates\n    coordinates = molecule[['x', 'y', 'z']].values\n    x_coordinates = coordinates[:, 0]\n    y_coordinates = coordinates[:, 1]\n    z_coordinates = coordinates[:, 2]\n    \n    ## Add radius to each atom to the list of atoms\n    radii = [atomic_radii[element] for element in elements]\n\n    ids = np.arange(coordinates.shape[0])\n    bonds = pd.DataFrame(columns = [\"Atom_0\", \"Atom_1\", \"Dist\"])\n    coordinates_compare, radii_compare, ids_compare = coordinates, radii, ids\n\n    ## For each atom in the atom list, \n        ## calcualte distances to all the other atoms,    \n    for _ in range(len(ids)):\n        coordinates_compare = np.roll(coordinates_compare, -1, axis=0)\n        radii_compare = np.roll(radii_compare, -1, axis=0)\n        ids_compare = np.roll(ids_compare, -1, axis=0)\n        distances = np.linalg.norm(coordinates - coordinates_compare, axis=1)\n        bond_distances = (radii + radii_compare) * 1.3\n        mask = np.logical_and(distances > 0.1, distances <  bond_distances)\n        distances = distances.round(2)\n        \n        ## Build tupple for bonds\n        for i, j, dist in zip(ids[mask], ids_compare[mask], distances[mask]):\n            if i < j:\n                bonds = bonds.append(pd.DataFrame([[i, j , dist]], columns = [\"Atom_0\", \"Atom_1\", \"Dist\"]))\n            else:\n                bonds = bonds.append(pd.DataFrame([[j, i , dist]], columns = [\"Atom_0\", \"Atom_1\", \"Dist\"]))\n        bonds = bonds.drop_duplicates()\n        subset = bonds[['Atom_0', 'Atom_1']]\n        tuples = [tuple(x) for x in subset.values]\n    \n    ## Use the tupple list of atoms and bonds to build graph\n    G=nx.Graph()\n    G.add_edges_from(tuples)\n\n    ## Add Node properties - Will be usefull in getting coordinates and atom name for each node\n    Node_Properties = (molecule\n                           .reset_index(drop=True)\n                           .groupby('atom_index')[['atom', 'x', 'y', 'z']]\n                           .apply(lambda x: x.to_dict()).to_dict()\n                      )\n    nx.set_node_attributes(G, Node_Properties)\n    return G","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def touch(path):\n    with open(path, 'a'):\n        os.utime(path, None)\n\ndef Crawl_Path(Full_Graph, start, end):\n    temp = list()\n    for node in nx.shortest_path(Full_Graph, start, end):\n        atom = Full_Graph.nodes[node].get(\"atom\").get(node)\n        x_Cord = Full_Graph.nodes[node].get(\"x\").get(node)\n        y_Cord = Full_Graph.nodes[node].get(\"y\").get(node)\n        z_Cord = Full_Graph.nodes[node].get(\"z\").get(node)\n\n        ### Get list of all node connected to the node along the path\n        Neighbors = [Full_Graph.nodes[n].get(\"atom\").get(n) for n in Full_Graph.neighbors(node)]\n\n        ### Get counts of C, H, N, O, F connected to the node along the path\n        C_Conns = Neighbors.count('C')\n        H_Conns = Neighbors.count('H')\n        N_Conns = Neighbors.count('N')\n        O_Conns = Neighbors.count('O')\n        F_Conns = Neighbors.count('F')\n        temp.append({'atom': atom, 'x_Coord': x_Cord, 'y_Coord': y_Cord, 'z_Coord': z_Cord, 'F_Conns': F_Conns,\n                        'C_Conns': C_Conns, 'H_Conns': H_Conns, 'N_Conns': N_Conns, 'O_Conns': O_Conns})\n    return temp\n    \ndef Get_Path_Details(temp_DF):\n    molecule_name = str(temp_DF['molecule_name'].iloc[0])\n    Full_Graph = Build_Graph(molecule_name)\n    temp_DF['Path'] = temp_DF.apply(lambda row: Crawl_Path(Full_Graph, row['atom_index_0'], row['atom_index_1']), axis=1)\n    return temp_DF","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from joblib import Parallel, delayed\nimport multiprocessing\n\nprint(datetime.datetime.now())\ndef applyParallel(dfGrouped, func):\n    retLst = Parallel(n_jobs=4)(delayed(func)(group) for name, group in dfGrouped)\n    return pd.concat(retLst)\n\nout = pd.DataFrame()\nout = out.append(applyParallel(train_df.groupby('molecule_name'), Get_Path_Details))\n\nprint(datetime.datetime.now())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def Map_PathInfo(t, Path):\n    t = len(Path)\n    if t == 2:\n        return (Path[0]['atom'], Path[0]['x_Coord'], Path[0]['y_Coord'], Path[0]['z_Coord']\n                    , Path[0]['F_Conns'], Path[0]['C_Conns'], Path[0]['H_Conns'], Path[0]['N_Conns'], Path[0]['O_Conns']\n                , Path[1]['atom'], Path[1]['x_Coord'], Path[1]['y_Coord'], Path[1]['z_Coord']\n                    , Path[1]['F_Conns'], Path[1]['C_Conns'], Path[1]['H_Conns'], Path[1]['N_Conns'], Path[1]['O_Conns']\n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0\n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)\n\n    elif t == 3:\n        return (Path[0]['atom'], Path[0]['x_Coord'], Path[0]['y_Coord'], Path[0]['z_Coord']\n                    , Path[0]['F_Conns'], Path[0]['C_Conns'], Path[0]['H_Conns'], Path[0]['N_Conns'], Path[0]['O_Conns']\n                \n                , Path[1]['atom'], Path[1]['x_Coord'], Path[1]['y_Coord'], Path[1]['z_Coord']\n                    , Path[1]['F_Conns'], Path[1]['C_Conns'], Path[1]['H_Conns'], Path[1]['N_Conns'], Path[1]['O_Conns']\n                \n                , Path[2]['atom'], Path[2]['x_Coord'], Path[2]['y_Coord'], Path[2]['z_Coord']\n                    , Path[2]['F_Conns'], Path[2]['C_Conns'], Path[2]['H_Conns'], Path[2]['N_Conns'], Path[2]['O_Conns']\n                \n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)\n\n    elif t == 4:\n        return (Path[0]['atom'], Path[0]['x_Coord'], Path[0]['y_Coord'], Path[0]['z_Coord']\n                    , Path[0]['F_Conns'], Path[0]['C_Conns'], Path[0]['H_Conns'], Path[0]['N_Conns'], Path[0]['O_Conns']\n                    , Path[1]['atom'], Path[1]['x_Coord'], Path[1]['y_Coord'], Path[1]['z_Coord']\n                    , Path[1]['F_Conns'], Path[1]['C_Conns'], Path[1]['H_Conns'], Path[1]['N_Conns'], Path[1]['O_Conns']\n                    , Path[2]['atom'], Path[2]['x_Coord'], Path[2]['y_Coord'], Path[2]['z_Coord']\n                    , Path[2]['F_Conns'], Path[2]['C_Conns'], Path[2]['H_Conns'], Path[2]['N_Conns'], Path[2]['O_Conns']\n                    , Path[3]['atom'], Path[3]['x_Coord'], Path[3]['y_Coord'], Path[3]['z_Coord']\n                    , Path[3]['F_Conns'], Path[3]['C_Conns'], Path[3]['H_Conns'], Path[3]['N_Conns'], Path[3]['O_Conns'])\n    \n    else:\n        return ('', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0\n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0\n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0\n                , '', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)\n\nprint(datetime.datetime.now())\n\n(out['Atom_0'], out['x_Cord_0'], out['y_Cord_0'], out['Z_Cord_0']\n        , out['F_Conns_0'], out['C_Conns_0'], out['H_Conns_0']\n        , out['N_Conns_0'], out['O_Conns_0']\n    , out['Atom_1'], out['x_Cord_1'], out['y_Cord_1'], out['Z_Cord_1']\n        , out['F_Conns_1'], out['C_Conns_1'], out['H_Conns_1']\n        , out['N_Conns_1'], out['O_Conns_1']\n    , out['Atom_2'], out['x_Cord_2'], out['y_Cord_2'], out['Z_Cord_2']\n        , out['F_Conns_2'], out['C_Conns_2'], out['H_Conns_2']\n        , out['N_Conns_2'], out['O_Conns_2']\n    , out['Atom_3'], out['x_Cord_3'], out['y_Cord_3'], out['Z_Cord_3']\n        , out['F_Conns_3'], out['C_Conns_3'], out['H_Conns_3']\n        , out['N_Conns_3'], out['O_Conns_3']) = (\n\n                zip(*out.apply(lambda row: Map_PathInfo(row['type'], row['Path']), axis = 1))\n                                                                                        )\n\nprint(datetime.datetime.now())\n\nout.to_pickle(r'../Graph based features')","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.8"}},"nbformat":4,"nbformat_minor":1}