{"cells":[{"metadata":{},"cell_type":"markdown","source":"Forked from https://www.kaggle.com/fnands/makegraphinput\n\nintroducing a validation set for use in MPNN\n\n# **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\nfrom tqdm import tqdm_notebook\nfrom sklearn import preprocessing\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import train_test_split\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":{"trusted":true},"cell_type":"code","source":"train_bonds.head()","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\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"angs.head()","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":"#why not standard scaler??\n#scale_min  = train['scalar_coupling_constant'].min()\n#scale_max  = train['scalar_coupling_constant'].max()\n#scale_mid = (scale_max + scale_min)/2\n#scale_norm = scale_max - scale_mid\n\n#train['scalar_coupling_constant'] = (train['scalar_coupling_constant'] - scale_mid)/scale_norm\ntrain['scalar_coupling_constant']=(train['scalar_coupling_constant']-train['scalar_coupling_constant'].min())/(train['scalar_coupling_constant'].max()-train['scalar_coupling_constant'].min())\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":{"trusted":true},"cell_type":"code","source":"train['scalar_coupling_constant'].hist()","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'])\n#why not standard scaler??\n#normalized_df=(df-df.min())/(df.max()-df.min())\n#structures[['x', 'y', 'z']] = structures[['x', 'y', 'z']]/10.\nstructures[['x', 'y', 'z']]=(structures[['x', 'y', 'z']]-structures[['x', 'y', 'z']].min())/(structures[['x', 'y', 'z']].max()-structures[['x', 'y', 'z']].min())\nnuclear_charge = {'H':1.0, 'C':6.0, 'N':7.0, 'O':8.8, 'F':9.0}\nstructures['nuclear_charge'] = [nuclear_charge[x] for x in structures['atom'].values]\nstructures['nuclear_charge'] = structures['nuclear_charge'] / 9.0\n#structures['nuclear_charge']=(structures['nuclear_charge']-structures['nuclear_charge'].min())/(structures['nuclear_charge'].max()-structures['nuclear_charge'].min())\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures[['x', 'y', 'z']].hist()","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()\ntrain_mol_names, valid_mol_names = train_test_split(train_mol_names, test_size=0.2, random_state=42)\nprint(train_mol_names.shape)\nprint(valid_mol_names.shape)\n\nvalid = train.loc[train['molecule_name'].isin(valid_mol_names)]\ntrain = train.loc[train['molecule_name'].isin(train_mol_names)]\n\nvalid_bonds = train_bonds.loc[train_bonds['molecule_name'].isin(valid_mol_names)]\ntrain_bonds = train_bonds.loc[train_bonds['molecule_name'].isin(train_mol_names)]\n\nprint(train.shape)\nprint(valid.shape)\n\ntest_mol_names  = test['molecule_name'].unique()\n\ntrain_structures = structures.loc[structures['molecule_name'].isin(train_mol_names)]\nvalid_structures = structures.loc[structures['molecule_name'].isin(valid_mol_names)]\ntest_structures = structures.loc[structures['molecule_name'].isin(test_mol_names)]\n\ntrain_struct_group = train_structures.groupby('molecule_name')\nvalid_struct_group = valid_structures.groupby('molecule_name')\ntest_struct_group  = test_structures.groupby('molecule_name')\n\ntrain_group = train.groupby('molecule_name')\nvalid_group = valid.groupby('molecule_name')\ntest_group  = test.groupby('molecule_name')\n\ntrain_bond_group = train_bonds.groupby('molecule_name')\nvalid_bond_group = valid_bonds.groupby('molecule_name')\ntest_bond_group  = test_bonds.groupby('molecule_name')\n\ntrain_angs = angs.loc[angs['molecule_name'].isin(train_mol_names)]\nvalid_angs = angs.loc[angs['molecule_name'].isin(valid_mol_names)]\ntest_angs = angs.loc[angs['molecule_name'].isin(test_mol_names)]\n\ntrain_angs_group = train_angs.groupby('molecule_name')\nvalid_angs_group = valid_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','nuclear_charge'] \n#Values our edges will have (minus distance, for now)\nbond_vals = ['nbond_1', 'nbond_1.5', 'nbond_2', 'nbond_3']\nj_coup_vals = ['1JHC', '1JHN', '2JHC', '2JHH', '2JHN', '3JHC', '3JHH', '3JHN']\n#j_coup_vals = ['2JHC', '2JHN','2JHH']\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_valid_mols = len(valid_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\n\nprint(node_dim)\nprint(bond_dim)\nprint(ang_dim)\nprint(edge_dim)\nprint(j_coup_dim)\nprint(j_coup_vals)","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\nvalid_nodes_array     = np.zeros((n_valid_mols, max_size, node_dim), dtype=np.float32) \nvalid_in_edges_array  = np.zeros((n_valid_mols, max_size, max_size, edge_dim + add_edge_dim),dtype=np.float32) \nvalid_out_edges_array = np.zeros((n_valid_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) \nprint(len(valid_group))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#how can we parallize this?\ndef make_arrs(val_group, struct_group, bond_group, ang_group, mode='train'):\n    debug=False\n    i = 0\n    maxit = 2\n    for values, structs, bonds, angles in tqdm_notebook(zip(val_group, struct_group, bond_group, ang_group),total=len(val_group)):\n        if i>maxit and debug: break\n        # Calculate distances\n        distances = np.zeros((max_size, max_size, add_edge_dim))\n        coords = structs[1][['x','y','z']].values\n        \n        dists  = distance_matrix(coords, coords)\n        #can we bin distances here?\n        \n        distances[:dists.shape[0],:dists.shape[1], 0] = dists\n        \n        # Create nodes\n        if debug:\n            print(structs)\n            print(structs[1])\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 is type descriptos one_hot_encoded -> use it to filter on type later on\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        if debug:\n            print(\"edges:\")\n            print(in_feats.shape)\n            print(in_bonds.shape)\n            print(ang_mat.shape)\n            print(distances.shape)\n        in_edges = np.concatenate((in_feats, in_bonds, ang_mat, distances),axis=2)\n  \n        if not mode=='test':           \n            out_edges = np.zeros((max_size, max_size, 1))\n            \n            # set irrelevant coupling values to zero\n            idx = values[1]['type'].isin(j_coup_vals)\n            values[1]['scalar_coupling_constant'].loc[~idx] = 0.0\n            \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            if debug:\n                print(idx)\n                print(values[1])\n                print(out_edges.shape)\n                print(out_edges)\n                input()\n            if mode == 'train':\n                train_nodes_array[i]      = nodes\n                train_in_edges_array[i]   = in_edges\n                train_out_edges_array[i]  = out_edges\n                \n            if mode == 'valid':\n                valid_nodes_array[i]      = nodes\n                valid_in_edges_array[i]   = in_edges\n                valid_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, mode = 'train')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#check distances\nprint(train_in_edges_array.shape)\n#distance matrix first molecule\nprint(train_in_edges_array[0,:,:,15])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_arrs(valid_group, valid_struct_group, valid_bond_group, valid_angs_group, mode = 'valid')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"make_arrs(test_group, test_struct_group, test_bond_group, test_angs_group, mode = 'test')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save as numpy arrays"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.savez_compressed(\"train_mol_names.npz\" , train_mol_names)\nnp.savez_compressed(\"valid_mol_names.npz\" , valid_mol_names)\n\nnp.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_valid.npz\" , valid_nodes_array)\nnp.savez_compressed(\"in_edges_valid.npz\" , valid_in_edges_array)\nnp.savez_compressed(\"out_edges_valid.npz\" , valid_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":{"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}