{"cells":[{"metadata":{},"cell_type":"markdown","source":"Thanks to Alexandre Sauvé, I've made it speed up!"},{"metadata":{},"cell_type":"markdown","source":"Hi gens! My idea is applying **Coulomb Interaction** which force is propotional to inverse squared distance (1/r^2). I guess inverse distance (1/r) can be also applicable when focusing on potential. Anyway, i've considered inversed squared distance, here. If we want to use inverse distance, preprocessed data can be converted easily.\n\n1. get assigned atoms from train data which are included 'atom_index_0' or 'atom_index_1' of molecule\n2. get distances from each atom belonging to the molecule and pickup 'num = 5' nearest regarding to each atom [H, C, N, O, F]. Though in this competition we focus on bondings of H-H, H-C, H-N, properties of bondings are strongly affected by O, F atoms. That's why I'd like to consider interaction as I mentioned.\n3. mearge distance array according to atom_index_0 and atom_index_1 then dimension of feature of bonding is 50 = num x atoms x 2.\n4. feed the feature into model. model can be built for each bonding type, 1JHH, 1JHC, 2JHC etc."},{"metadata":{},"cell_type":"markdown","source":"## contents\n\n* [Preparations](#Preparations)\n* [Compute dictances](#Compute-distances)\n* [Merge DataFrames](#Merge-DataFrames)\n* [Train MLP regression](#Train-MLP-regression)\n* [Visualize prediction](#Visualize-prediction)\n"},{"metadata":{},"cell_type":"markdown","source":"## Preparations"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport time\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.neural_network import MLPRegressor","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FOLDER = '../input/champs-scalar-coupling/'\nOUTPUT = '../input/preprocessed/'\nos.listdir(FOLDER)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"My idea may be nice because only structure data is used for train. In other words, other property data such as moment and potential can be ignored."},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_mulliken_charges = pd.read_csv(FOLDER + 'mulliken_charges.csv')\n# df_sample =  pd.read_csv(FOLDER + 'sample_submission.csv')\n# df_magnetic_shielding_tensors = pd.read_csv(FOLDER + 'magnetic_shielding_tensors.csv')\ndf_train = pd.read_csv(FOLDER + 'train.csv')\n# df_test = pd.read_csv(FOLDER + 'test.csv')\n# df_dipole_moments = pd.read_csv(FOLDER + 'dipole_moments.csv')\n# df_potential_energy = pd.read_csv(FOLDER + 'potential_energy.csv')\ndf_structures = pd.read_csv(FOLDER + 'structures.csv')\n# df_scalar_coupling_contributions = pd.read_csv(FOLDER + 'scalar_coupling_contributions.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Compute distances\n\nInverse squared distances are computed by functions below."},{"metadata":{},"cell_type":"markdown","source":"This function is to get distances each other in a molecule. The output is (n, n) matrix. \"n\" is the number of atoms in molecule."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_dist_matrix(df_structures_idx, molecule):\n    df_temp = df_structures_idx.loc[molecule]\n    locs = df_temp[['x','y','z']].values\n    num_atoms = len(locs)\n    loc_tile = np.tile(locs.T, (num_atoms,1,1))\n    dist_mat = ((loc_tile - loc_tile.T)**2).sum(axis=1)\n    return dist_mat","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This function is to get assigned atoms which we are interested in for the bondings. Assigned atoms seems to be only H, C, N. You know, O and F are not in scope for our task."},{"metadata":{"trusted":true},"cell_type":"code","source":"def assign_atoms_index(df_idx, molecule):\n    se_0 = df_idx.loc[molecule]['atom_index_0']\n    se_1 = df_idx.loc[molecule]['atom_index_1']\n    if type(se_0) == np.int64:\n        se_0 = pd.Series(se_0)\n    if type(se_1) == np.int64:\n        se_1 = pd.Series(se_1)\n    assign_idx = pd.concat([se_0, se_1]).unique()\n    assign_idx.sort()\n    return assign_idx","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is to get distances which origins are assigned atoms. \n\nOrigins are atom_index_0 in df_train. Distances are called from distance matrix generated by function defined above, but not all. Only \"num_pickup (default 5)\" nearest for each atoms H, C, N, O, F are called. For example, if there are 10 H in a molecule, only 5 H are considered as inverse squared distance. Other 4 H are ignored."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pickup_dist_matrix(df_idx, df_structures_idx, molecule, num_pickup=5, atoms=['H', 'C', 'N', 'O', 'F']):\n    pickup_dist_matrix = np.zeros([0, len(atoms)*num_pickup])\n    assigned_idxs = assign_atoms_index(df_idx, molecule) # [0, 1, 2, 3, 4, 5, 6] -> [1, 2, 3, 4, 5, 6]\n    dist_mat = get_dist_matrix(df_structures_idx, molecule)\n    for idx in assigned_idxs: # [1, 2, 3, 4, 5, 6] -> [2]\n        df_temp = df_structures_idx.loc[molecule]\n        locs = df_temp[['x','y','z']].values\n\n        dist_arr = dist_mat[idx] # (7, 7) -> (7, )\n\n        atoms_mole = df_structures_idx.loc[molecule]['atom'].values # ['O', 'C', 'C', 'N', 'H', 'H', 'H']\n        atoms_mole_idx = df_structures_idx.loc[molecule]['atom_index'].values # [0, 1, 2, 3, 4, 5, 6]\n\n        mask_atoms_mole_idx = atoms_mole_idx != idx # [ True,  True, False,  True,  True,  True,  True]\n        masked_atoms = atoms_mole[mask_atoms_mole_idx] # ['O', 'C', 'N', 'H', 'H', 'H']\n        masked_atoms_idx = atoms_mole_idx[mask_atoms_mole_idx]  # [0, 1, 3, 4, 5, 6]\n        masked_dist_arr = dist_arr[mask_atoms_mole_idx]  # [ 5.48387003, 2.15181049, 1.33269675, 10.0578779, 4.34733927, 4.34727838]\n        masked_locs = locs[masked_atoms_idx]\n\n        sorting_idx = np.argsort(masked_dist_arr) # [2, 1, 5, 4, 0, 3]\n        sorted_atoms_idx = masked_atoms_idx[sorting_idx] # [3, 1, 6, 5, 0, 4]\n        sorted_atoms = masked_atoms[sorting_idx] # ['N', 'C', 'H', 'H', 'O', 'H']\n        sorted_dist_arr = 1/masked_dist_arr[sorting_idx] #[0.75035825,0.46472494,0.23002898,0.23002576,0.18235297,0.09942455]\n\n        target_matrix = np.zeros([len(atoms), num_pickup])\n        for a, atom in enumerate(atoms):\n            pickup_atom = sorted_atoms == atom # [False, False,  True,  True, False,  True]\n            pickup_dist = sorted_dist_arr[pickup_atom] # [0.23002898, 0.23002576, 0.09942455]\n\n            num_atom = len(pickup_dist)\n            if num_atom > num_pickup:\n                target_matrix[a, :num_pickup] = pickup_dist[:num_pickup]\n            else:\n                target_matrix[a, :num_atom] = pickup_dist\n        \n        pickup_dist_matrix = np.vstack([pickup_dist_matrix, target_matrix.reshape(-1)])\n    return pickup_dist_matrix #(num_atoms, num_pickup*5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below is execution however it takes long time even if it's been made speed up. Please try this kernel in your own environment!"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_structures_idx = df_structures.set_index('molecule_name')\ndf_train_idx = df_train.set_index('molecule_name')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num = 5\nmols = df_train['molecule_name'].unique()\nnum_div = len(mols) // 5\ndist_mat = np.zeros([0, num*5])\natoms_idx = np.zeros([0], dtype=np.int32)\nmolecule_names = np.empty([0])\n\nstart = time.time()\n\nfor mol in mols[:100]:\n    \n    assigned_idxs = assign_atoms_index(df_train_idx, mol)\n    dist_mat_mole = get_pickup_dist_matrix(df_train_idx, df_structures_idx, mol, num_pickup=num)\n    mol_name_arr = [mol] * len(assigned_idxs) \n    \n    molecule_names = np.hstack([molecule_names, mol_name_arr])\n    atoms_idx = np.hstack([atoms_idx, assigned_idxs])\n    dist_mat = np.vstack([dist_mat, dist_mat_mole])\n    \ncol_name_list = []\natoms = ['H', 'C', 'N', 'O', 'F']\nfor a in atoms:\n    for n in range(num):\n        col_name_list.append('dist_{}_{}'.format(a, n))\n        \nse_mole = pd.Series(molecule_names, name='molecule_name')\nse_atom_idx = pd.Series(atoms_idx, name='atom_index')\ndf_dist = pd.DataFrame(dist_mat, columns=col_name_list)\ndf_distance = pd.concat([se_mole, se_atom_idx,df_dist], axis=1)\n\nelapsed_time = time.time() - start\nprint (\"elapsed_time:{0}\".format(elapsed_time) + \"[sec]\")","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"# df_distance.to_csv(OUTPUT + 'distance1000.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_dist = pd.read_csv(OUTPUT + 'distance1000.csv')\ndf_distance.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge DataFrames\n\nBelow is picking up atoms that are assigned for each target bonding by keys of atom_index."},{"metadata":{"trusted":true},"cell_type":"code","source":"def merge_atom(df, df_distance):\n    df_merge_0 = pd.merge(df, df_distance, left_on=['molecule_name', 'atom_index_0'], right_on=['molecule_name', 'atom_index'])\n    df_merge_0_1 = pd.merge(df_merge_0, df_distance, left_on=['molecule_name', 'atom_index_1'], right_on=['molecule_name', 'atom_index'])\n    del df_merge_0_1['atom_index_x'], df_merge_0_1['atom_index_y']\n    return df_merge_0_1","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"start = time.time()\ndf_train_dist = merge_atom(df_train, df_distance) # corrected!: df_dist -> df_distance\nelapsed_time = time.time() - start\nprint (\"elapsed_time:{0}\".format(elapsed_time) + \"[sec]\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_train_dist.to_csv(OUTPUT + 'train_dist1000.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_train_dist = pd.read_csv(OUTPUT + 'train_dist1000.csv')\ndf_train_dist.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train MLP regression\nFinaly, I feed data into model. Here I use simple MLP. I expect better model can be found.\nHere, I've chosen multi layer perceptron regression for checking my preprocessings. We would find better model, such as lightGBM."},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"df_1JHC = df_train_dist.query('type == \"1JHC\"')\ny = df_1JHC['scalar_coupling_constant'].values\nX = df_1JHC[df_1JHC.columns[6:]].values\nprint(X.shape)\nprint(y.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mlp = MLPRegressor(hidden_layer_sizes=(100,50))\nmlp.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualize prediction\n\nBelow is checking scatter of validation and its prediction.\n\nLooks good! \n\nValidation data points (y_val, y_pred) are almost on line! I expect models for other bonding (2JHH, 2JHC,,) can be built the same way. And, accuracy can be better."},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred = mlp.predict(X_val)\nplt.scatter(y_val, y_pred)\nplt.title('1JHC')\nplt.plot([80, 200], [80, 200])\nplt.show()","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":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"nbformat":4,"nbformat_minor":1}