{"cells":[{"metadata":{},"cell_type":"markdown","source":"# **Graph creator**\n\nA graph is a relativly natural way of representing molecules, and many method make use of structuring the data in this way.\n\nThis kernel shows a basic example of one can structure our data as a graph.   \n\nHere we will create an array for our node values, and an adjacency matrix for our edge values. \n\n\n(**Note**: One can argue whether an adjacency matrix is really the best way to go here as we have an undirected graph, and it is therefore a bit innefficienct ( n^2 as opposed to n(n-1)/2 ), but it is easy to work with)\n\nThis is by no means the fastest way of doing this, but it is straightforward and only has to be run once, and the output can then be used. "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom scipy.spatial import distance_matrix\n\nfrom sklearn import preprocessing\nimport os\nprint(os.listdir(\"../input\"))\n\ndatadir = \"../input/\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read in train, test and structures files"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv(datadir + 'champs-scalar-coupling/train.csv')\ntest = pd.read_csv(datadir + 'champs-scalar-coupling/test.csv')\nstructures = pd.read_csv(datadir + 'champs-scalar-coupling/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read in bonds files\nTaken from:   https://www.kaggle.com/asauve/dataset-with-number-of-bonds-between-atoms  \n(thanks Alexandre Sauvé!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_bonds = pd.read_csv(datadir + 'predicting-molecular-properties-bonds/train_bonds.csv')\ntest_bonds = pd.read_csv(datadir + 'predicting-molecular-properties-bonds/test_bonds.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Reads from angles file\nTaken from: https://www.kaggle.com/soerendip/calculate-angles-and-dihedrals-with-networkx\n(thanks Rakete!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"angs = pd.read_csv(datadir + \"angle-and-dihedral-for-the-champs-structures/angles.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Normalize targets so they have are centered around 0 and have max of 1, and one-hot encode coupling types"},{"metadata":{"trusted":true},"cell_type":"code","source":"scale_min  = train['scalar_coupling_constant'].min()\nscale_max  = train['scalar_coupling_constant'].max()\nscale_mid = (scale_max + scale_min)/2\nscale_norm = scale_max - scale_mid\n\ntrain['scalar_coupling_constant'] = (train['scalar_coupling_constant'] - scale_mid)/scale_norm\n\n# One hot encoding gets  too big for Kaggle, let's try label\n# use npz now, back to OH\ntrain[['1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN']] =  pd.get_dummies(train['type'])\ntest[['1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN']]  =  pd.get_dummies(test['type'])\n\n#le = preprocessing.LabelEncoder()\n#le.fit(['1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN'])\n#train['l_type'] = (le.transform(train['type']) + 1)/8.\n#test['l_type'] = (le.transform(test['type']) + 1)/8.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pre-process the structures by one-hot encoding the atom types, and normalize distances to have around max of 1"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures[['C', 'F' ,'H', 'N', 'O']] = pd.get_dummies(structures['atom'])\nstructures[['x', 'y', 'z']] = structures[['x', 'y', 'z']]/10.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Process bonds"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_bonds[['nbond_1', 'nbond_1.5', 'nbond_2', 'nbond_3']] = pd.get_dummies(test_bonds['nbond'])#test_bonds['nbond']/3\ntrain_bonds[['nbond_1', 'nbond_1.5', 'nbond_2', 'nbond_3']] = pd.get_dummies(train_bonds['nbond'])#train_bonds['nbond']/3\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Process angles"},{"metadata":{"trusted":true},"cell_type":"code","source":"angs['dihedral'] = angs['dihedral']/np.pi\n# Should I rather one-hot this?\nangs['shortest_path_n_bonds'] = angs['shortest_path_n_bonds']/6.0\nangs = angs.fillna(0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Find training and testing molecules, and split structrues into test and train. Then group by molecule\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mol_names = train['molecule_name'].unique()\ntest_mol_names  = test['molecule_name'].unique()\n\ntrain_structures = structures.loc[structures['molecule_name'].isin(train_mol_names)]\ntest_structures = structures.loc[structures['molecule_name'].isin(test_mol_names)]\n\ntrain_struct_group = train_structures.groupby('molecule_name')\ntest_struct_group  = test_structures.groupby('molecule_name')\n\ntrain_group = train.groupby('molecule_name')\ntest_group  = test.groupby('molecule_name')\n\ntrain_bond_group = train_bonds.groupby('molecule_name')\ntest_bond_group  = test_bonds.groupby('molecule_name')\n\ntrain_angs = angs.loc[angs['molecule_name'].isin(train_mol_names)]\ntest_angs = angs.loc[angs['molecule_name'].isin(test_mol_names)]\n\ntrain_angs_group = train_angs.groupby('molecule_name')\ntest_angs_group  = test_angs.groupby('molecule_name')\n\n# Find max nodes in graph:\nmax_size = train_struct_group.size().max()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Define node and edge values"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Values our nodes will have\nnode_vals = ['C', 'F' ,'H', 'N', 'O']#, 'x', 'y', 'z']\n#Values our edges will have (minus distance, for now)\nbond_vals = ['nbond_1', 'nbond_1.5', 'nbond_2', 'nbond_3']#['nbond']\nj_coup_vals = ['1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN']#'l_type']\nang_vals = ['shortest_path_n_bonds','cosinus','dihedral']\nedge_vals = j_coup_vals + bond_vals + ang_vals\n\n# Find amount of training molecules\nn_train_mols = len(train_mol_names)\nn_test_mols = len(test_mol_names)\n\n# Find dim of edges and nodes\nbond_dim  = len(bond_vals)\nj_coup_dim= len(j_coup_vals)\nang_dim   = len(ang_vals)\nnode_dim  = len(node_vals)\nedge_dim  = len(edge_vals) \n\n# Additional edge dims for distances \nadd_edge_dim = 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pre-allocate arrays that we will fill later\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_nodes_array     = np.zeros((n_train_mols, max_size, node_dim), dtype=np.float32) \ntrain_in_edges_array  = np.zeros((n_train_mols, max_size, max_size, edge_dim + add_edge_dim),dtype=np.float32) \ntrain_out_edges_array = np.zeros((n_train_mols, max_size, max_size, 1),dtype=np.float32) \n\ntest_nodes_array     = np.zeros((n_test_mols, max_size, node_dim), dtype=np.float32) \ntest_in_edges_array  = np.zeros((n_test_mols, max_size, max_size, edge_dim + add_edge_dim),dtype=np.float32) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_arrs(val_group, struct_group, bond_group, ang_group, test):\n    i = 0\n    for values, structs, bonds, angles in zip(val_group, struct_group, bond_group, ang_group):\n        if (not i%1000):\n            print(i)\n\n        # Calculate distances\n        distances = np.zeros((max_size, max_size, 1))\n        coords = structs[1][['x','y','z']].values\n        dists  = distance_matrix(coords, coords)\n        distances[:dists.shape[0],:dists.shape[1], 0] = dists \n        \n        # Create nodes\n        mol_info = structs[1][node_vals].values\n        nodes = np.zeros((max_size, node_dim))\n        nodes[:mol_info.shape[0], :mol_info.shape[1]] = mol_info\n\n        # Create edges\n        in_feats = np.zeros((max_size, max_size, j_coup_dim))\n        ind = values[1][['atom_index_0', 'atom_index_1' ]].values\n        in_feats[ind[:,0], ind[:,1], 0:j_coup_dim] = values[1][j_coup_vals].values\n        in_feats[ind[:,1], ind[:,0], 0:j_coup_dim] = in_feats[ind[:,0], ind[:,1], 0:j_coup_dim]\n\n        # Create bonds\n        in_bonds = np.zeros((max_size, max_size, bond_dim))\n        ind_bonds = bonds[1][['atom_index_0', 'atom_index_1' ]].values\n        in_bonds[ind_bonds[:,0], ind_bonds[:,1]] = bonds[1][bond_vals].values\n        in_bonds[ind_bonds[:,1], ind_bonds[:,0]] = in_bonds[ind_bonds[:,0], ind_bonds[:,1]]\n        \n        # Create angles\n        ind_angs = angles[1][['atom_index_0', 'atom_index_1' ]].values\n        ang_mat  = np.zeros((max_size, max_size, ang_dim))\n        ang_mat[ind_angs[:,0], ind_angs[:,1]]  = angles[1][ang_vals]\n        ang_mat[ind_angs[:,1], ind_angs[:,0]]  = ang_mat[ind_angs[:,0], ind_angs[:,1]]\n        \n        # concat all edge values \n        in_edges = np.concatenate((in_feats, in_bonds, ang_mat, distances),axis=2)\n\n\n\n        \n        if not test:           \n            out_edges = np.zeros((max_size, max_size, 1))\n            out_edges[ind[:,0], ind[:,1], 0] = values[1]['scalar_coupling_constant' ].values\n            out_edges[ind[:,1], ind[:,0], 0] = out_edges[ind[:,0], ind[:,1], 0]\n        \n\n            train_nodes_array[i]      = nodes\n            train_in_edges_array[i]   = in_edges\n            train_out_edges_array[i]  = out_edges\n        else:\n            test_nodes_array[i]      = nodes\n            test_in_edges_array[i]   = in_edges\n        i = i + 1\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_arrs(train_group, train_struct_group, train_bond_group, train_angs_group, test = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_arrs(test_group, test_struct_group, test_bond_group, test_angs_group, test = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save as numpy arrays"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.savez_compressed(\"nodes_train.npz\" , train_nodes_array)\nnp.savez_compressed(\"in_edges_train.npz\" , train_in_edges_array)\nnp.savez_compressed(\"out_edges_train.npz\" , train_out_edges_array)\n\nnp.savez_compressed(\"nodes_test.npz\" , test_nodes_array)\nnp.savez_compressed(\"in_edges_test.npz\" , test_in_edges_array)","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}