{"cells":[{"metadata":{},"cell_type":"markdown","source":"FORK!  \nThanks to Ryoji Nomura, https://www.kaggle.com/rio114  \nThanks to Alexandre Sauvé, for dramatic speed up!  \n"},{"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/'\n\n# kaggle cloud has no output folder\n# os.makedirs('output',exist_ok=True)\nOUTPUT = '.'\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\n\n","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\n","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. When computing 1000 molecules, it took 1000 sec in my home environment. That's why pre-computed csv is uploaded. Now its faster."},{"metadata":{"trusted":true},"cell_type":"code","source":"# define index for faster computations\ndf_structures_idx = df_structures.set_index('molecule_name')\ndf_train_idx = df_train.set_index('molecule_name')    \n    \n# only 5 hydrogen atoms are considered as inverse squared distance\nnum = 5\n\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])\nstart = time.time()\n\n# number of molecules to process\nmax_mol = 100000\nk = 0\nprint (\"Calculating \",max_mol,\" molecules.\")\n\n# from joblib import parallel_backend\n# with parallel_backend('threading', n_jobs=2):\n#    Parallel()(delayed)\n\nfor mol in mols[:max_mol]:\n    k += 1\n    if k and k % 1000 == 0:\n        print(k, ': ' , mol)\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: {:4.2f}\".format(elapsed_time) + \" [sec]\")\n","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"# kaggle cloud has no output folder\n# os.makedirs('output')\ndf_distance.to_csv(OUTPUT + 'distance-XXX.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# private dataset for speed purposes pre-calculted\n# df_dist = pd.read_csv(OUTPUT + 'distance1000.csv')\ndf_dist.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":"# key index error in pandas\ndef 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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train.head(15)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_dist.head(15)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"start = time.time()\n# df_train_dist = merge_atom(df_train, df_dist)\n# df_train_dist = pd.concat([df_train, df_dist], axis=1)\ndf_train_dist = merge_atom(df_train, df_distance) # corrected!: df_dist -> df_distance\nelapsed_time = time.time() - start\nprint (\"elapsed_time: {:4.2f}\".format(elapsed_time) + \" [sec]\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train_dist.head(15)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# takes a long time to write, unessecary if data is not processed further\n# df_train_dist.to_csv(OUTPUT + 'train_dist-XXX.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# private data included for speed purposes\n# df_train_dist = pd.read_csv(OUTPUT + 'train_dist1000.csv')\ndf_train_dist.head(15)\n","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":"# remove nan\nX = np.nan_to_num(X)\nX","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":"start = time.time()\nprint('start training regressor')\n\n# on the Kaggle cloud one should see 400% CPU use (quad thread)\n# use hidden_layer_sizes=(100,50) for better results\n# not sure this is best approach, depends on layer otpimization\n# simple regression maybe better see https://scikit-learn.org/stable/modules/classes.html\n\n## adam is a bit better but takes longer than 'sgd' solver\nmlp = MLPRegressor(activation='relu', solver='sgd', hidden_layer_sizes=(100,50))\nmlp.fit(X_train, y_train)\ny_pred = mlp.predict(X_val)\n\n# Ransac is much faster, but less accurate\n## from sklearn.linear_model import RANSACRegressor\n## from sklearn.datasets import make_regression\n## ransac = RANSACRegressor(random_state=1234)\n## ransac.fit(X_train, y_train)\n## y_pred = ransac.predict(X_val)\n\n## from sklearn.linear_model import TheilSenRegressor\n## from sklearn.datasets import make_regression\n## theilsen = TheilSenRegressor(random_state=1234)\n## theilsen.fit(X_train, y_train)\n## y_pred = theilsen.predict(X_val)\n\nelapsed_time = time.time() - start\nprint (\"elapsed_time: {:4.2f}\".format(elapsed_time) + \" [sec]\")","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":"print('len(y_val)  :',len(y_val))\nprint('len(y_pred) :', len(y_pred))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_val.view()[:10]\ny_val.view()[-10:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(y_val, y_pred, marker='.')\nplt.title('1JHC')\nplt.plot([60, 220], [60, 220])\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# not sure if correct, also axis labels crooked\nfrom sklearn.linear_model import Ridge\nfrom yellowbrick.regressor import ResidualsPlot\n\n# Instantiate the linear model and visualizer\nridge = Ridge()\nvisualizer = ResidualsPlot(ridge)\n\nvisualizer.fit(X_train, y_train)   # Fit the training data to the model\nvisualizer.score(X_val, y_pred)    # Evaluate the model on the test data\nvisualizer.poof()                  # Draw/show/poof the data\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred.view()[:10]\ny_pred.view()[-10:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import *\nfrom math import sqrt\n\n# current sklearn on kaggle version has no max_error\n# https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/metrics/regression.py\ndef max_error(y_true, y_pred):\n    return np.max(np.abs(y_true - y_pred))\n\nprint(\"Mean squared error     : %.6f\" %    mean_squared_error(y_val, y_pred))\nprint(\"Median absolute error  : %.6f\" % median_absolute_error(y_val, y_pred))\nprint(\"Mean absolute error    : %.6f\" %   mean_absolute_error(y_val, y_pred))\nprint(\"Maximum residual error : %.6f\" %             max_error(y_val, y_pred))\nprint(\"                  RMSE : %.6f\" % sqrt(mean_squared_error(y_val, y_pred)))\nprint(\"                    R2 : %.6f\" %                 r2_score(y_val, y_pred))          \n\n","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}