{"cells":[{"metadata":{},"cell_type":"markdown","source":"### 1. Introduction\n\nIn this kernel I share the features that made me able to get -LB 2.0 with NN and LGB approaches and helped me place on top 5%.\n\n<b>They include:</b>\n- Distance between atoms;\n- Dihedral angles\n- Eigen Values\n- Forces like Yukawa and Van der Waals\n- ACSF descriptors\n- Mulliken charges\n- bonds properties\n- a few others\n\nThe most <u>creative features</u> that you will find here that are not in other public kernels are forces resultants. Basically I calculated the axis between Atom0 and Atom1 (called it X), the axis between atom0 and the molecule center (called it Y) and the axis orthonal to X and Y (called it Z). For each one of these axis, I calculated the resultant forces like Yukawa and Van der Waals. According to my LGB model, those were the most important features split and gain wise.\n\nYou can also check my NN and LGB implementation and tricks in the kernel [NN and LGB tricks and pipeline for top 5% LB ](https://www.kaggle.com/felipemello/nn-and-lgb-tricks-and-pipeline-for-top-5-lb). Hope you guys like it.\n\n<b>Please, if you find the content here interesting, consider upvoting the kernel to reward my time editing and sharing it. Thank you very much :)</b>\n\nWe will only calculte features here for '1JHN'. for speed and RAM reasons, but you can easily change that by altering the variable: 'mol_types' on section 3"},{"metadata":{},"cell_type":"markdown","source":"### 2. Load libs and utils"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom operator import itemgetter\nfrom scipy.sparse.csgraph import shortest_path\nfrom scipy.stats import kurtosis, skew\nimport gc\nfrom itertools import combinations, permutations\nimport time\nimport gc\nimport warnings\nimport copy\n\n\npd.options.display.precision = 15\nnp.set_printoptions(suppress=True)\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def 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.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\ndef get_dihedrals(dihedral_paths, vec_mat):\n    \n    x0, x1, x2, x3 = dihedral_paths.T\n    \n    b0 = -1.0 * vec_mat[x1, x0]\n    b1 = vec_mat[x2, x1]\n    b2 = vec_mat[x3, x2]\n\n    b0xb1 = np.cross(b0, b1)\n    b1xb2 = np.cross(b1, b2)\n\n    b0xb1_x_b1xb2 = np.cross(b0xb1, b1xb2)\n\n    y = np.sum(b0xb1_x_b1xb2*b1, axis=1)/np.linalg.norm(b1, axis=1)\n    x = np.sum(b0xb1*b1xb2, axis=1)\n    \n    grad = np.arctan2(y, x)\n    return grad\n\ndef get_path_matrix(Pr, n_atoms, dist_mat):\n    \n    path_matrix = [0]*(n_atoms**2)\n    path_length_matrix = [0]*(n_atoms**2)\n    path_distance = [0]*(n_atoms**2)\n    for i in range(n_atoms):\n        for j in range(n_atoms):                 \n            path = [j]\n            k = j\n            while Pr[i, k] != -9999:\n                path.append(Pr[i, k])\n                k = Pr[i, k]\n            path = np.array(path[::-1])\n            path_matrix[i*n_atoms+j] = path\n            path_len = len(path)\n            path_length_matrix[i*n_atoms+j] = path_len\n            path_distance[i*n_atoms+j] = np.sum([dist_mat[path[i], path[i+1]] for i in range(path_len-1)])\n    \n    return np.array(path_length_matrix).reshape(n_atoms, n_atoms), np.array(path_matrix).reshape(n_atoms, n_atoms), np.array(path_distance).reshape(n_atoms, n_atoms)\n    \ndef get_fast_dist_matrix(xyz, df_start_indexes, struct_start_indexes, molecule_id):\n    struct_start_index, struct_end_index = struct_start_indexes[molecule_id], struct_start_indexes[molecule_id+1]\n    locs = xyz[struct_start_index:struct_end_index]    \n    num_atoms = struct_end_index - struct_start_index\n    loc_tile = np.tile(locs.T, (num_atoms,1,1))\n    dist_mat = np.sqrt(((loc_tile - loc_tile.T)**2).sum(axis=1))\n    vec_mat = loc_tile.T - loc_tile \n    vec_mat = np.transpose(vec_mat, axes=(1,2,0)).T\n    \n    df_start_index, df_end_index = df_start_indexes[molecule_id], df_start_indexes[molecule_id+1]\n    atom_index = atom_indexes[df_start_index:df_end_index]    \n    return dist_mat, vec_mat, atom_index, locs\n\ndef get_projection(u,v):\n    #v must be an unit-vector\n    return np.sum(u*v, axis=1).reshape(-1,1)*v\n\ndef calculate_force_resultant(unit_vec, force_mat, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat):\n    \n    forces = unit_vec*force_mat.reshape(n_atoms,n_atoms,1)\n    \n    #forces_resultant_per_axis has shape: [[X, Y, Z resultant_force on atom_num_0], [X, Y, Z resultant_force on atom_num_1], [...], [X, Y, Z resultant_force on atom_num_n]]\n    forces_resultant_per_axis = np.sum(forces, axis=1)\n    \n    u1 = forces_resultant_per_axis[atom_index[:,0]]\n    u2 = forces_resultant_per_axis[atom_index[:,1]]\n    v_X = main_X_unit_vec.reshape(n_combinations,3)\n    proj_X_1 = get_projection(u1,v_X)\n    proj_X_2 = get_projection(u2,v_X)\n    proj_X = proj_X_1 + proj_X_2\n    angle_X = np.sum(proj_X*v_X, axis=1)/(np.linalg.norm(proj_X, axis=1))\n    X_resultant = (np.linalg.norm(proj_X, axis = 1)*angle_X).reshape(-1,1)\n    X_resultant = np.nan_to_num(X_resultant)\n    \n    v_Y = main_Y_unit_vec.reshape(n_combinations,3)\n    proj_Y_1 = get_projection(u1,v_Y)\n    proj_Y_2 = get_projection(u2,v_Y)\n    proj_Y_1_norm = np.linalg.norm(proj_Y_1, axis=1)\n    proj_Y_2_norm = np.linalg.norm(proj_Y_2, axis=1)\n    angle_Y = np.sum(proj_Y_1*proj_Y_2, axis=1)/(proj_Y_1_norm*proj_Y_2_norm)\n    angle_Y = np.nan_to_num(angle_Y)\n    force_Y_momentum = (proj_Y_1_norm - proj_Y_2_norm*angle_Y)\n    \n    v_Z = main_Z_unit_vec.reshape(n_combinations,3)\n    proj_Z_1 = get_projection(u1,v_Z)\n    proj_Z_2 = get_projection(u2,v_Z)\n    proj_Z_1_norm = np.linalg.norm(proj_Z_1, axis=1)\n    proj_Z_2_norm = np.linalg.norm(proj_Z_2, axis=1)\n    angle_Z = np.sum(proj_Z_1*proj_Z_2, axis=1)/(proj_Z_1_norm*proj_Z_2_norm)\n    angle_Z = np.nan_to_num(angle_Z)\n    force_Z_momentum = (proj_Z_1_norm - proj_Z_2_norm*angle_Z)\n    \n    momentum_Y = (force_Y_momentum*dist_mat[atom_index[:,1], atom_index[:,0]]).reshape(n_combinations, 1)\n    momentum_Z = (force_Z_momentum*dist_mat[atom_index[:,1], atom_index[:,0]]).reshape(n_combinations, 1)\n    momentum = np.sqrt(momentum_Y**2 + momentum_Z**2)\n    \n    Y_Z_resultant = np.linalg.norm(proj_Y_1 + proj_Y_2 + proj_Z_1 + proj_Z_2, axis=1).reshape(n_combinations, 1)\n    Y_Z_resultant_0 = np.linalg.norm(proj_Y_1 + proj_Z_1, axis=1).reshape(n_combinations, 1)\n    Y_Z_resultant_1 = np.linalg.norm(proj_Y_2 + proj_Z_2, axis=1).reshape(n_combinations, 1)\n    Y_resultant =  np.linalg.norm(proj_Y_1 + proj_Y_2, axis=1).reshape(n_combinations, 1)\n    Z_resultant =  np.linalg.norm(proj_Z_1 + proj_Z_2, axis=1).reshape(n_combinations, 1)\n    proportion_between_X_and_Y_Z = 1 - Y_Z_resultant/(Y_Z_resultant+np.abs(X_resultant))\n    proportion_between_X_and_Y_Z = np.nan_to_num(proportion_between_X_and_Y_Z)\n    \n    return X_resultant, Y_Z_resultant, proportion_between_X_and_Y_Z, momentum, Y_Z_resultant_0, Y_Z_resultant_1, Y_resultant, Z_resultant\n\ndef f_cutoff(data, Rc):\n    \n    fc = 0.5*(np.cos(np.pi*data/Rc) + 1)\n    fc[data > Rc] = 0\n    \n    return fc\n\ndef fill_array(array, n_combinations, expected_length, fill):\n    \n    if expected_length > array.shape[1]:\n        array = np.concatenate([array, np.zeros((n_combinations,expected_length - array.shape[1])) + fill], axis=1)\n    \n    return array\n\ndef get_G2(Rij, eta, Rs, Rc):\n     return (np.exp(-eta*(Rij-Rs)**2) * f_cutoff(Rij, Rc)).sum(axis=-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Pick your configurations\n\nThese are the main configurations to run this kernel. Addapt it as you want."},{"metadata":{"trusted":true},"cell_type":"code","source":"#path to where the data competition is\noriginal_data_folder = '../input/champs-scalar-coupling'\n\n#'train' or 'test' depending on the dataset that you want to build features for\ndf_type = 'train' #or 'test'\n\n#verbose\nverbose = True\n\n#number o nearest atoms to consider when calculating features.\n#More atom = more features\n#After the 8th atom, it starts to become irrelevant to the model\nmax_n_nearest = 11\n\n#it will be a lot faster to generate features for all mol_types at once. \n#Just do one at a time if you have less than 16MB of RAM\nmol_types = ['1JHN'] \n\n#Number of CPU cores for parallel processing\nnumber_of_cores = 10","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4. Get the data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"################\n#Get all the data\n################\n\ndf = pd.read_csv(f'{original_data_folder}/{df_type}.csv')\ndf = df[df['type'].isin(mol_types)] #selecte only the rows with desired coupling types\ndf = reduce_mem_usage(df, verbose=True)\n\nunique_mol = df['molecule_name'].unique()\n\ndef get_data(data_path, verbose):\n    \n    data = pd.read_csv(data_path)\n    data = data.loc[data['molecule_name'].isin(unique_mol)]\n    data = reduce_mem_usage(data, verbose=True)\n    \n    return data\n\ndf_bonds = get_data(data_path = f'../input/predicting-molecular-properties-bonds/{df_type}_bonds.csv', verbose = verbose)\ndf_structures = get_data(data_path = '../input/champs-scalar-coupling/structures.csv', verbose = verbose)\ndf_structures_eigen = get_data(data_path = '../input/qm7-coulomb-matrix-eigenvalue-features/struct_eigen.csv', verbose = verbose)\nif df_type == 'test':\n    df_mulliken = get_data(data_path = '../input/predicting-mulliken-charges-with-acsf-descriptors/mulliken_charges_test_set.csv', verbose = verbose)\nelse:\n    df_mulliken = get_data(data_path = '../input/champs-scalar-coupling/mulliken_charges.csv', verbose = verbose)\n          \n\n################\n#Get indexes per molecule, so we can work with each molecule at a time.\n#start_indexes = index where the molecule start\n#end_indexes = index where the molecule end\n################\n\nxyz = df_structures[['x','y','z']].values\natom_indexes = df[['atom_index_0', 'atom_index_1']].values\n\nstruct_end_indexes = df_structures.groupby('molecule_name').size().cumsum()\ndf_end_indexes = df[['molecule_name', 'atom_index_0', 'atom_index_1']].groupby('molecule_name').size().cumsum()\ndf_bonds_end_indexes = df_bonds[['molecule_name', 'atom_index_0', 'atom_index_1']].groupby('molecule_name').size().cumsum()\n\nstruct_start_indexes = np.zeros(len(struct_end_indexes) + 1, 'int')\nstruct_start_indexes[1:] = struct_end_indexes\nstruct_start_indexes[0] = 0\n\ndf_start_indexes = np.zeros(len(df_end_indexes) + 1, 'int')\ndf_start_indexes[1:] = df_end_indexes\ndf_start_indexes[0] = 0\n\ndf_bonds_start_indexes = np.zeros(len(df_bonds_end_indexes) + 1, 'int')\ndf_bonds_start_indexes[1:] = df_bonds_end_indexes\ndf_bonds_start_indexes[0] = 0\n\n################\n#Data about atoms that we use in our calculations. The variable radius depends on the bonds, \n#while atom_radius is an acceptable average.\n#EN is the eletronegativity\n#Enconde is the atomic number\n################\n\natom_variable_radius = {'H': np.array([0.32, np.NaN, np.NaN]),\n                'C': np.array([.75, .67, .60]),\n                'O': np.array([.63, .57, .54]),\n                'N': np.array([.71, .60, .54]),\n                'F': np.array([.64, .59, .53])}\n\natoms_radius = {'H':0.38, 'C':0.77, 'N':0.75,'O':0.73, 'F':0.71}\natoms_EN = {'H':2.2, 'C':2.55, 'N':3.04,'O':3.44, 'F':3.98}\natoms_encoding = {0:0, '0':0, 'H':1, 'C':6, 'N':7,'O':8, 'F':9}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 5. Define function to calculate features and test it"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_mol_feats(molecule_id):\n    \n    ##########\n    #get some initial data\n    #########\n    \n    dist_mat, vec_mat, atom_index, locs = get_fast_dist_matrix(xyz, df_start_indexes, struct_start_indexes, molecule_id)  \n    \n    #Uncomment if you want to augument the data by changing atom_0 and atom_1 positions\n    #atom_index = np.array([[i,j] for j, i in atom_index])\n    \n    #n_nearest is the value of how many atoms you wanna analyze and it is defined by:\n    #mininum value between the chosen max_n_nearest and how many atoms are in this molecule\n    n_nearest = min(dist_mat.shape[1], max_n_nearest)\n    \n    #n_combinations is how many couplings are being calculated in this molecule\n    n_combinations = len(atom_index)\n    \n    #Number of atoms in the molecule\n    n_atoms = len(dist_mat)\n    \n    #unit vec between atoms\n    unit_vec = np.nan_to_num(vec_mat/np.linalg.norm(vec_mat, axis=2).reshape(-1, n_atoms, 1))\n    \n    #n_nearest atoms to atom_0 and atom_1 ordered by their distance to the center of atom_0 and atom_1\n    ordered_smallest_dist_indexes = np.argpartition(dist_mat, range(n_nearest))[:,1:n_nearest+1]\n    \n    \n    struct_start_index, struct_end_index = struct_start_indexes[molecule_id], struct_start_indexes[molecule_id+1]\n    df_start_index, df_end_index = df_start_indexes[molecule_id], df_start_indexes[molecule_id+1]\n    df_bonds_start_index, df_bonds_end_index = df_bonds_start_indexes[molecule_id], df_bonds_start_indexes[molecule_id+1]\n    \n    #taken from the structures_eigen df\n    eigen_cols = [ 'connectedness', 'coulomb_mean', 'eigv_max',\n                   'eigv_min', 'fiedler_eig', 'sv_0', 'sv_1', 'sv_2',\n                   'sv_3', 'sv_4', 'sv_min']\n    \n    eigen_info = df_structures_eigen.iloc[struct_start_index:struct_end_index]\n    eigen_info = eigen_info[eigen_cols].to_numpy()\n    \n    atom_names = np.array(df_structures[struct_start_index:struct_end_index]['atom'])\n    mol_name = df_structures.iloc[struct_start_index, 0]\n    atom_code = np.array(itemgetter(*atom_names)(atoms_encoding))\n    \n    ##########\n    #get distances from all n_nearest atoms to atom_0, atom_1 and other 2 closest atoms\n    #########\n    \n    locs0 = locs[atom_index[:,0]]\n    locs1 = locs[atom_index[:,1]]\n    \n    mid_point = ((locs0 + locs1)/2).reshape(-1, 1, 3)\n    mid_point_distances = locs-mid_point\n    mid_point_distances = np.sqrt(np.sum(mid_point_distances**2, axis=-1))\n    ordered_smallest_dist_indexes = np.argpartition(mid_point_distances, range(n_nearest))[:,:n_nearest]\n    \n    closest_atoms = []\n    for i, row in enumerate(ordered_smallest_dist_indexes):\n        atoms_0_and_1_index = atom_index[i,:]\n        nearest_ordered_atoms_excluding_0_and_1 = row[np.logical_and(row != atom_index[i,0], row != atom_index[i,1])][:n_nearest-2]\n        closest_atoms.append(np.concatenate([atoms_0_and_1_index, nearest_ordered_atoms_excluding_0_and_1]))\n\n    closest_atoms = np.array(closest_atoms)\n    \n    n_combs_atom_pairs = []\n    for c in range(n_combinations):\n        atom_pairs = []\n        for i in range(1, len(closest_atoms[0])):\n            for vi in range(min(4, i)):\n                atom_pairs.append([closest_atoms[c][i], closest_atoms[c][vi]])\n        \n        n_combs_atom_pairs.append(atom_pairs)\n        \n    n_combs_atom_pairs = np.array(n_combs_atom_pairs)     \n            \n    distances_atom_pairs = dist_mat[n_combs_atom_pairs[:,:,0], n_combs_atom_pairs[:,:,1]]       \n    \n    expected_number_of_pair_combinations = 6 + 4*(max_n_nearest-4)\n    distances_atom_pairs = fill_array(distances_atom_pairs, n_combinations, expected_length = expected_number_of_pair_combinations, fill = 0)\n    \n    ##########\n    #Calculate G4 angle descriptors\n    #########\n    \n    \n    all_Rij_Rik =  np.array(list(permutations(list(range(n_atoms)), 3)))\n    \n    \n    Rij = dist_mat[all_Rij_Rik[:,0], all_Rij_Rik[:,1]].reshape(n_atoms, -1)\n    Rik = dist_mat[all_Rij_Rik[:,0], all_Rij_Rik[:,2]].reshape(n_atoms, -1)\n    Rjk = dist_mat[all_Rij_Rik[:,2], all_Rij_Rik[:,1]].reshape(n_atoms, -1)\n    all_ij_unit_vec = unit_vec[all_Rij_Rik[:,1], all_Rij_Rik[:,0]].reshape(n_atoms, -1, 3)\n    all_ik_unit_vec = unit_vec[all_Rij_Rik[:,2], all_Rij_Rik[:,0]].reshape(n_atoms, -1, 3)\n    cosijk = np.sum(all_ij_unit_vec * all_ik_unit_vec, axis=-1)\n    \n\n    G4 = []\n    Rc = 10\n    \n    for eta, ksi, lamb, Rs in [[0.01, 4,  1, 2], [1, 4,  -1, 6], [0.1, 4,  -1, 6], [0.5, 2,  1, 6]]:\n        G4_values = (1+lamb*cosijk)**ksi\n        G4_values *= np.exp(-eta*(Rij**2 + Rik**2 + Rjk**2))\n        G4_values = G4_values * f_cutoff(Rij, Rc) * f_cutoff(Rjk, Rc) * f_cutoff(Rik, Rc)\n        G4_values = (2**(1-ksi))*np.sum(G4_values, axis=-1)\n        G4.append(G4_values)\n        \n    G4_selected = []\n\n    for row_G4 in G4:\n        G4_selected.append(fill_array(row_G4[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0))\n    \n    G4_selected = np.concatenate(G4_selected, axis=-1)\n\n    ##########\n    #get some bond based properties, like std, mean, kurtosis, etc\n    #########\n    \n    bond_info = df_bonds.iloc[df_bonds_start_index:df_bonds_end_index, [1, 2, 3, 4]].to_numpy() # 'atom_index_0', 'atom_index_1', 'n_bonds', 'L2dist'\n    \n    bonds_dist_mean = np.zeros(n_atoms)\n    bonds_dist_std = np.zeros(n_atoms)\n    bonds_dist_kurt = np.zeros(n_atoms)\n    bonds_dist_skew = np.zeros(n_atoms)\n    atom_cov_radius = np.zeros(n_atoms)\n    n_bonds = np.zeros(n_atoms)\n    \n    bond_matrix = np.zeros((n_atoms, n_atoms))\n    \n    column_atoms_0 = bond_info[:, 0].astype('int')\n    column_atoms_1 = bond_info[:, 1].astype('int')\n    \n    bond_matrix[column_atoms_0, column_atoms_1] = 1\n    bond_matrix[column_atoms_1, column_atoms_0] = 1\n    is_atom_bond = bond_matrix[n_combs_atom_pairs[:,:,0], n_combs_atom_pairs[:,:,1]]   \n    \n    for i in range(n_atoms):\n        atom_name = atom_names[i]\n        atom_column_0 = bond_info[:, 0] == i\n        atom_column_1 = bond_info[:, 1] == i\n        \n        bond_info_atom_i = bond_info[np.logical_or(atom_column_0, atom_column_1)]\n        bonds_lengths = bond_info_atom_i[:,3]\n        bonds_dist_mean[i] = np.mean(bonds_lengths)\n        bonds_dist_std[i] = np.std(bonds_lengths)\n        bonds_dist_kurt[i] = kurtosis(bonds_lengths)\n        bonds_dist_skew[i] = skew(bonds_lengths)\n        \n        highest_order_bond = max(bond_info_atom_i[:,2])\n        n_bonds[i] = int(np.sum(atom_column_0)+ np.sum(atom_column_1))\n        atom_cov_radius[i] = atom_variable_radius[atom_name][int(highest_order_bond-1)]\n        \n\n    bonds_dist_mean = fill_array(bonds_dist_mean[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0)    \n    bonds_dist_std = fill_array(bonds_dist_std[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0)    \n    bonds_dist_kurt = fill_array(bonds_dist_kurt[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0)    \n    bonds_dist_skew = fill_array(bonds_dist_skew[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0)    \n    n_bonds = fill_array(n_bonds[closest_atoms],n_combinations, expected_length = max_n_nearest, fill = 0)    \n    \n    expected_number_of_pair_combinations = 6 + 4*(max_n_nearest-4)\n    is_atom_bond = fill_array(is_atom_bond, n_combinations, expected_length = expected_number_of_pair_combinations, fill = 0)\n    \n    ##########\n    #Calculate dihedral angles and graph path distances between atoms\n    #########\n    \n    atom_radius = np.array(itemgetter(*list(range(n_atoms)))(atom_cov_radius))\n    atom_radius = np.tile(atom_radius, n_atoms).reshape(dist_mat.shape)\n    atom_radius_mult = atom_radius*atom_radius.T\n    atom_radius_sum = atom_radius + atom_radius.T\n    distance_minus_radius = dist_mat-atom_radius_sum\n    distance_minus_radius[distance_minus_radius<0] = 0\n\n    #GET DIHEDRALS AND PATHS\n    graph_matrix = np.zeros(distance_minus_radius.shape)\n    bonds_index = distance_minus_radius<=0.3\n    graph_matrix[bonds_index] = dist_mat[bonds_index]\n    \n    D, Pr = shortest_path(graph_matrix, directed=False, method='FW', return_predecessors=True)\n    \n    path_length_matrix, path_matrix, path_distance_matrix = get_path_matrix(Pr, n_atoms, dist_mat)\n    \n    dihedral_indices = np.where(path_length_matrix==4)\n    dihedral_matrix = np.zeros((n_atoms, n_atoms)) - 11\n    dihedral_paths = []\n    dihedral = []\n    if len(dihedral_indices[0]) > 0:\n        dihedral_paths = np.concatenate(path_matrix[dihedral_indices[0], dihedral_indices[1]]).reshape(-1, 4)\n        dihedral = get_dihedrals(dihedral_paths, vec_mat)\n        \n        #We use -11 to set apart atoms without dihedrals\n        #We use -10 to se apart cos(dihedrals) of atoms without dihedrals\n        dihedral_matrix[dihedral_indices[0], dihedral_indices[1]] = dihedral\n        \n        error = 0.001 #is used to compare floats\n        dihedral_main = dihedral_matrix[atom_index[:,0], atom_index[:,1]]\n        cos_dihedral_main = np.cos(dihedral_main)\n        cos_dihedral_main[(np.abs(cos_dihedral_main - np.cos(-11))) < error] = -10\n        \n        cos_2_times_dihedral_main = np.cos(2*dihedral_main)\n        cos_2_times_dihedral_main[(np.abs(cos_2_times_dihedral_main - np.cos(-2*11))) < error] = -10\n        \n    else:\n        dihedral_main = np.zeros((n_combinations, 1)) - 11\n        cos_dihedral_main = np.zeros((n_combinations, 1)) - 10\n        cos_2_times_dihedral_main = np.zeros((n_combinations, 1)) - 10\n    \n    cos_dihedral_main = cos_dihedral_main.reshape(n_combinations,1)\n    cos_2_times_dihedral_main = cos_2_times_dihedral_main.reshape(n_combinations,1)\n  \n    \n    ##########\n    #Calculate forces like coulomb, yukawa, harmonic, van der waals, and others\n    #########\n    \n\n    atoms_mulliken_charges = np.array(df_mulliken[struct_start_index:struct_end_index]['mulliken_charge'])\n    \n    atoms_mulliken = np.tile(atoms_mulliken_charges, n_atoms).reshape(dist_mat.shape)\n    atoms_mulliken_mult = atoms_mulliken*atoms_mulliken.T\n    \n    coulomb = atoms_mulliken_mult/dist_mat\n    coulomb[np.eye(n_atoms)==1] = 0\n    \n    yukawa = np.exp(-dist_mat)/dist_mat\n    yukawa[np.eye(n_atoms)==1] = 0\n    \n    harmonic_distance = atom_radius_mult/atom_radius_sum\n    harmonic_distance[np.eye(n_atoms)==1] = 0\n    \n    vander_numerator = harmonic_distance\n    vander_denominator = 6*(distance_minus_radius)**2\n    vander_denominator[vander_denominator < 1] = 1 #avoid division by 0\n    vander = vander_numerator/vander_denominator\n    vander[np.eye(n_atoms)==1] = 0\n\n    inv_dist_1 = 1/dist_mat\n    inv_dist_1[np.eye(n_atoms)==1] = 0\n    \n    inv_dist_2 = 1/dist_mat**2\n    inv_dist_2[np.eye(n_atoms)==1] = 0\n    \n    inv_dist_3 = 1/dist_mat**3\n    inv_dist_3[np.eye(n_atoms)==1] = 0\n    \n    atom_EN = np.array(itemgetter(*atom_names)(atoms_EN))\n    atom_EN = np.tile(atom_EN, n_atoms).reshape(dist_mat.shape)\n    atom_EN_sum = atom_EN + atom_EN.T\n    \n    inv_EN = 1/(dist_mat*0.5*atom_EN_sum)**2\n    inv_EN[np.eye(n_atoms)==1] = 0\n    \n    yukawa_resultant = np.linalg.norm(np.sum(unit_vec*yukawa.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    coulomb_resultant = np.linalg.norm(np.sum(unit_vec*coulomb.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    harmonic_resultant = np.linalg.norm(np.sum(unit_vec*harmonic_distance.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    vander_resultant = np.linalg.norm(np.sum(unit_vec*vander.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    inv_EN_resultant = np.linalg.norm(np.sum(unit_vec*inv_EN.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    inv_dist_2_resultant = np.linalg.norm(np.sum(unit_vec*inv_dist_2.reshape(n_atoms,n_atoms,1), axis=1), axis=-1)\n    \n    \n    yukawa_resultant = fill_array(yukawa_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    coulomb_resultant = fill_array(coulomb_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    harmonic_resultant = fill_array(harmonic_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    vander_resultant = fill_array(vander_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    inv_EN_resultant = fill_array(inv_EN_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    inv_dist_2_resultant = fill_array(inv_dist_2_resultant[closest_atoms], n_combinations, expected_length = max_n_nearest, fill = 0)\n    \n    ####################\n        \n    #Every pair of atom_0 and atom_1 form a plane with the center(0, 0, 0) of the atom.\n    #To get the plane we need cross(a,b) resulting in coeficients c1, c2, c3 and a point, in this case, we will use (0,0,0)\n    #The final place equation is c1(x -0) + c2(y - 0) + c3(z -0) = 0, which is just the cross product.\n    planes = np.cross(locs[atom_index[:,0]], locs[atom_index[:,1]])\n    \n    #lets get the unit orthogonal vec of this place by diving by the norm and reshape it for the next step\n    main_Z_unit_vec = ( planes/np.linalg.norm(planes,axis=1).reshape(-1,1) ).reshape(-1, 1, 3)\n    main_X_unit_vec = unit_vec[atom_index[:,1], atom_index[:,0]].reshape(-1, 1, 3)\n    main_Y_unit_vec = np.cross(main_Z_unit_vec, main_X_unit_vec)\n\n    ####################    \n    \n    #calculate resultant forces\n    yukawa_parallel, yukawa_orthogonal, yukawa_proportion_between_X_and_Y_Z, yukawa_momentum, yukawa_Y_Z_resultant_0, yukawa_Y_Z_resultant_1, yukawa_Y_resultant, yukawa_Z_resultant = calculate_force_resultant(unit_vec, yukawa, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    inv_EN_parallel, inv_EN_orthogonal, inv_EN_proportion_between_X_and_Y_Z, inv_EN_momentum, inv_EN_Y_Z_resultant_0, inv_EN_Y_Z_resultant_1, inv_EN_Y_resultant, inv_EN_Z_resultant = calculate_force_resultant(unit_vec, inv_EN, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    coulomb_parallel, coulomb_orthogonal, coulomb_proportion_between_X_and_Y_Z, coulomb_momentum, coulomb_Y_Z_resultant_0, coulomb_Y_Z_resultant_1, coulomb_Y_resultant, coulomb_Z_resultant = calculate_force_resultant(unit_vec, coulomb, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    harmonic_parallel, harmonic_orthogonal, harmonic_proportion_between_X_and_Y_Z, harmonic_momentum, harmonic_Y_Z_resultant_0, harmonic_Y_Z_resultant_1, harmonic_Y_resultant, harmonic_Z_resultant  = calculate_force_resultant(unit_vec, harmonic_distance, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    vander_parallel, vander_orthogonal, vander_proportion_between_X_and_Y_Z, vander_momentum, vander_Y_Z_resultant_0, vander_Y_Z_resultant_1, vander_Y_resultant, vander_Z_resultant  = calculate_force_resultant(unit_vec, vander, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    dist_1_parallel, dist_1_orthogonal, dist_1_proportion_between_X_and_Y_Z, dist_1_momentum, dist_1_Y_Z_resultant_0, dist_1_Y_Z_resultant_1, dist_1_Y_resultant, dist_1_Z_resultant  = calculate_force_resultant(unit_vec, dist_mat, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    inv_dist_1_parallel, inv_dist_1_orthogonal, inv_dist_1_proportion_between_X_and_Y_Z, inv_dist_1_momentum, inv_dist_1_Y_Z_resultant_0, inv_dist_1_Y_Z_resultant_1, inv_dist_1_Y_resultant, inv_dist_1_Z_resultant  = calculate_force_resultant(unit_vec, inv_dist_1, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    inv_dist_2_parallel, inv_dist_2_orthogonal, inv_dist_2_proportion_between_X_and_Y_Z, inv_dist_2_momentum, inv_dist_2_Y_Z_resultant_0, inv_dist_2_Y_Z_resultant_1, inv_dist_2_Y_resultant, inv_dist_2_Z_resultant  = calculate_force_resultant(unit_vec, inv_dist_2, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n    inv_dist_3_parallel, inv_dist_3_orthogonal, inv_dist_3_proportion_between_X_and_Y_Z, inv_dist_3_momentum, inv_dist_3_Y_Z_resultant_0, inv_dist_3_Y_Z_resultant_1, inv_dist_3_Y_resultant, inv_dist_3_Z_resultant  = calculate_force_resultant(unit_vec, inv_dist_3, n_atoms, atom_index, n_combinations, main_X_unit_vec, main_Y_unit_vec, main_Z_unit_vec, dist_mat)\n      \n    #Some final feactures based on angle, mulliken charge, and if there is a bond between atoms\n\n    unit_nearest_0 = unit_vec[closest_atoms[:, 2:],atom_index[:,0].reshape(n_combinations, 1)]\n    unit_nearest_1 = unit_vec[closest_atoms[:, 2:],atom_index[:,1].reshape(n_combinations, 1)]\n    angle_between_0_1_nearest_to_center = np.sum(unit_nearest_0*unit_nearest_1, axis=2)\n    \n    angle_between_0_1_nearest_to_center = fill_array(angle_between_0_1_nearest_to_center, n_combinations, expected_length = max_n_nearest, fill = 0)\n    \n    mulliken_charges = atoms_mulliken_charges[closest_atoms].reshape(n_combinations, -1)\n    mulliken_charges = fill_array(mulliken_charges, n_combinations, expected_length = max_n_nearest, fill = 0)\n   \n    atom_code = atom_code[closest_atoms].reshape(n_combinations, -1)\n    atom_code = fill_array(atom_code, n_combinations, expected_length = max_n_nearest, fill = 0)\n    \n    \n    ##########\n    #Put it all together:\n    #########\n    \n    one_column_data= np.concatenate([np.tile(np.array([sum(atom_names == 'C')]), n_combinations).reshape(n_combinations,1),\n                                     np.tile(np.array([sum(atom_names == 'H')]), n_combinations).reshape(n_combinations,1),\n                                     np.tile(np.array([sum(atom_names == 'F')]), n_combinations).reshape(n_combinations,1),\n                                     np.tile(np.array([sum(atom_names == 'N')]), n_combinations).reshape(n_combinations,1),\n                                     np.tile(np.array([sum(atom_names == 'O')]), n_combinations).reshape(n_combinations,1),\n                                     \n                                     cos_dihedral_main ,\n                                     cos_2_times_dihedral_main ,\n\n                                       yukawa_Y_Z_resultant_0, \n                                       yukawa_Y_Z_resultant_1 ,\n                                       yukawa_Y_resultant ,\n                                       yukawa_Z_resultant ,\n                                       yukawa_momentum,\n                                       \n                                       inv_EN_Y_Z_resultant_0, \n                                       inv_EN_Y_Z_resultant_1 ,\n                                       inv_EN_Y_resultant ,\n                                       inv_EN_Z_resultant ,\n                                       inv_EN_momentum,\n                                       \n                                       vander_Y_Z_resultant_0, \n                                       vander_Y_Z_resultant_1 ,\n                                       vander_Y_resultant ,\n                                       vander_Z_resultant ,\n                                       vander_momentum,\n                                       \n                                       coulomb_Y_Z_resultant_0, \n                                       coulomb_Y_Z_resultant_1 ,\n                                       coulomb_Y_resultant ,\n                                       coulomb_Z_resultant ,\n                                       coulomb_momentum,\n                                       \n                                       yukawa_parallel ,\n                                       coulomb_parallel ,\n                                       harmonic_parallel ,\n                                       vander_parallel ,\n                                       dist_1_parallel ,\n                                       inv_dist_1_parallel ,\n                                       inv_dist_2_parallel ,\n                                       inv_dist_3_parallel ,\n                                       \n                                       yukawa_orthogonal ,\n                                       coulomb_orthogonal ,\n                                       harmonic_orthogonal, \n                                       vander_orthogonal ,\n                                       dist_1_orthogonal,\n                                       inv_dist_1_orthogonal,\n                                       inv_dist_2_orthogonal,\n                                       inv_dist_3_orthogonal,\n                                       \n                                       yukawa_proportion_between_X_and_Y_Z, \n                                       coulomb_proportion_between_X_and_Y_Z, \n                                       harmonic_proportion_between_X_and_Y_Z, \n                                       vander_proportion_between_X_and_Y_Z ,\n                                       dist_1_proportion_between_X_and_Y_Z,\n                                       inv_dist_1_proportion_between_X_and_Y_Z,\n                                       inv_dist_2_proportion_between_X_and_Y_Z,\n                                       inv_dist_3_proportion_between_X_and_Y_Z,\n                                       \n                                       bonds_dist_kurt[:,0].reshape(n_combinations,1),\n                                       bonds_dist_kurt[:,1].reshape(n_combinations,1),\n                                       bonds_dist_kurt[:,2].reshape(n_combinations,1),\n                                       bonds_dist_skew[:,0].reshape(n_combinations,1),\n                                       bonds_dist_skew[:,1].reshape(n_combinations,1),\n                                       bonds_dist_skew[:,2].reshape(n_combinations,1)], axis=1) \n    \n    multiple_columns_data = [ distances_atom_pairs,\n                             is_atom_bond,\n                             n_bonds,\n                             G4_selected,\n                             bonds_dist_mean,\n                             bonds_dist_std,\n                             yukawa_resultant,\n                             coulomb_resultant,\n                             inv_EN_resultant,\n                             angle_between_0_1_nearest_to_center,\n                             mulliken_charges,\n                             atom_code,\n                             eigen_info[closest_atoms[:,0]],\n                             eigen_info[closest_atoms[:,1]],\n                             eigen_info[closest_atoms[:,2]]]\n    \n    results = np.concatenate([one_column_data, np.concatenate(multiple_columns_data, axis=1)], axis=1).astype('float32')\n    \n    return results","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's test our function and see how many features we have, considering the 11 nearest atoms"},{"metadata":{"trusted":true},"cell_type":"code","source":"original_max_n_nearest = copy.copy(max_n_nearest)\n\nmolecule_id = 0\nmax_n_nearest = 11\nfeatures = get_mol_feats(molecule_id)\nprint(features.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Good, we have 301 features. Let's check them."},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame(features)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is still kinda hard to understand what is what. After we add column names to the dataframe, it will make more sense.\n\nIf you think that 11 nearest atoms or 301 features is too much, you can decrease the number of nearest atoms to consider.\n\nExample: 8 nearest atoms, will generate 238 features, instead of 301"},{"metadata":{"trusted":true},"cell_type":"code","source":"max_n_nearest = 8\nmolecule_id = 0\nfeatures = get_mol_feats(molecule_id)\nprint(features.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 6. Run it in parallel\n\nNow let's run it in parallel for speed"},{"metadata":{"trusted":true},"cell_type":"code","source":"max_n_nearest = original_max_n_nearest\n\nimport multiprocessing\nimport traceback\np = multiprocessing.Pool(number_of_cores)\n\ntry:\n    t1 = time.time()\n    result = p.map(get_mol_feats, list(range(df_structures.molecule_name.nunique())))\n    t2 = time.time()\nexcept Exception as e:\n    print('Caught exception in worker thread (x = d):')\n    traceback.print_exc()\n    print()\n    raise e\n    \np.terminate()\np.join\n\nprint(f'Calculating features took {(t2 - t1)/60:0.2f} minutes')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 7. Put in dataframe format and save it\n\nTo save RAM memory, lets delete all that we don't need.\n\nLet's also create the names for the columns in our dataframe"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"del df_structures_eigen, df_bonds_end_indexes, df_bonds_start_indexes, xyz, unique_mol, df_start_indexes, df_end_indexes, df_bonds, df, df_structures, df_mulliken, atom_indexes\n\nimport gc\ngc.collect()\n\n#Let's create the names of the columns of our soon to be DataFrame\n\n#column names for distance atom pairs\ndistances_atom_pairs = []\nis_atom_bond = []\nfor i in range(1, max_n_nearest):\n    for vi in range(min(4, i)):\n        distances_atom_pairs.append(f'{i}_{vi}_distances_atom_pairs')\n        is_atom_bond.append(f'{i}_{vi}_atom_bond')\n\n\neigen_0 = [col + '_0' for col in [ 'connectedness', 'coulomb_mean', 'eigv_max', 'eigv_min', 'fiedler_eig', 'sv_0', 'sv_1', 'sv_2','sv_3', 'sv_4', 'sv_min']]\neigen_1 = [col + '_1' for col in [ 'connectedness', 'coulomb_mean', 'eigv_max', 'eigv_min', 'fiedler_eig', 'sv_0', 'sv_1', 'sv_2','sv_3', 'sv_4', 'sv_min']]\neigen_2 = [col + '_2' for col in [ 'connectedness', 'coulomb_mean', 'eigv_max', 'eigv_min', 'fiedler_eig', 'sv_0', 'sv_1', 'sv_2','sv_3', 'sv_4', 'sv_min']]\nn_neighbours = [f'{i}_n_neighbours' for i in range(1,max_n_nearest + 1)]\nG4_selected_1 = [f'{i}_G4_1_selected' for i in range(1,max_n_nearest + 1)]\nG4_selected_2 = [f'{i}_G4_2_selected' for i in range(1,max_n_nearest + 1)]\nG4_selected_3 = [f'{i}_G4_3_selected' for i in range(1,max_n_nearest + 1)]\nG4_selected_4 = [f'{i}_G4_5_selected' for i in range(1,max_n_nearest + 1)]\nG4_selected_5 = [f'{i}_G4_6_selected' for i in range(1,max_n_nearest + 1)]\nG4_selected_6 = [f'{i}_G4_7_selected' for i in range(1,max_n_nearest + 1)]\nbonds_dist_mean = [f'{i}_bonds_dist_mean' for i in range(1,max_n_nearest + 1)]\nbonds_dist_std = [f'{i}_bonds_dist_std' for i in range(1,max_n_nearest + 1)]\nyukawa_resultant = [f'{i}_yukawa_resultant' for i in range(1,max_n_nearest + 1)]\ncoulomb_resultant = [f'{i}_coulomb_resultant' for i in range(1,max_n_nearest + 1)]\nharmonic_resultant = [f'{i}_harmonic_resultant' for i in range(1,max_n_nearest + 1)]\nvander_resultant = [f'{i}_vander_resultant' for i in range(1,max_n_nearest + 1)]\ninv_EN_resultant = [f'{i}_inv_EN_resultant' for i in range(1,max_n_nearest + 1)]\ninv_dist_2_resultant = [f'{i}_inv_dist_2_resultant' for i in range(1,max_n_nearest + 1)]\nangle_between_0_1_nearest_to_center = [f'{i}_angle_between_0_1_nearest_to_center' for i in range(1,max_n_nearest + 1)]\nmulliken_charges = [f'{i}_mulliken_charges' for i in range(1,max_n_nearest + 1)]\natom_code = [f'{i}_atom_code' for i in range(1,max_n_nearest + 1)]\n             \nsingle_columns = ['N_C',\n                  'N_H',\n                  'N_F',\n                  'N_N',\n                  'N_O',\n                  'cos_dihedral_main',\n                  'cos_2_times_dihedral_main',\n                  \n                   'yukawa_Y_Z_resultant_0', \n                   'yukawa_Y_Z_resultant_1' ,\n                   'yukawa_Y_resultant',\n                   'yukawa_Z_resultant' ,\n                   'yukawa_momentum',\n                   \n                   'inv_EN_Y_Z_resultant_0', \n                   'inv_EN_Y_Z_resultant_1' ,\n                   'inv_EN_Y_resultant',\n                   'inv_EN_Z_resultant' ,\n                   'inv_EN_momentum',\n                   \n                   'vander_Y_Z_resultant_0', \n                   'vander_Y_Z_resultant_1' ,\n                   'vander_Y_resultant',\n                   'vander_Z_resultant' ,\n                   'vander_momentum',\n                   \n                   'coulomb_Y_Z_resultant_0', \n                   'coulomb_Y_Z_resultant_1' ,\n                   'coulomb_Y_resultant',\n                   'coulomb_Z_resultant' ,\n                   'coulomb_momentum',\n                   \n                   'yukawa_parallel',\n                   'coulomb_parallel',\n                   'harmonic_parallel',\n                   'vander_parallel',\n                   'dist_1_parallel',\n                   'inv_dist_1_parallel',\n                   'inv_dist_2_parallel',\n                   'inv_dist_3_parallel',\n                   \n                   'yukawa_orthogonal',\n                   'coulomb_orthogonal',\n                   'harmonic_orthogonal', \n                   'vander_orthogonal',\n                   'dist_1_orthogonal',\n                   'inv_dist_1_orthogonal',\n                   'inv_dist_2_orthogonal',\n                   'inv_dist_3_orthogonal',\n                   \n                   'yukawa_proportion_between_X_and_Y_Z', \n                   'coulomb_proportion_between_X_and_Y_Z', \n                   'harmonic_proportion_between_X_and_Y_Z', \n                   'vander_proportion_between_X_and_Y_Z' ,\n                   'dist_1_proportion_between_X_and_Y_Z',\n                   'inv_dist_1_proportion_between_X_and_Y_Z',\n                   'inv_dist_2_proportion_between_X_and_Y_Z',\n                   'inv_dist_3_proportion_between_X_and_Y_Z',\n                   \n                   '0_bonds_dist_kurt',\n                   '1_bonds_dist_kurt',\n                   '2_bonds_dist_kurt',\n                   '0_bonds_dist_skew',\n                   '1_bonds_dist_skew',\n                   '2_bonds_dist_skew']\n                             \n\n\nmul_columns = [  distances_atom_pairs,\n                   n_neighbours,\n                 is_atom_bond,\n                 G4_selected_1,\n                 G4_selected_2,\n                 G4_selected_3,\n                 G4_selected_4,\n                 bonds_dist_mean,\n                 bonds_dist_std,\n                 yukawa_resultant,\n                 coulomb_resultant,\n                 inv_EN_resultant,\n                 angle_between_0_1_nearest_to_center,\n                 mulliken_charges,\n                 atom_code,\n                 eigen_0,\n                 eigen_1,\n                 eigen_2]\n\n\nall_col = np.sum([single_columns, np.sum(mul_columns)])\n\n\ndtype_dict = dict()\nfor col in all_col:\n    \n    dtype_dict[col] = 'float32'\n    \n    if 'connectedness' in col:\n        dtype_dict[col] = 'int8'\n    if 'atom_code' in col:\n        dtype_dict[col] = 'int8'\n    if col.startswith('N_'):\n        dtype_dict[col] = 'int8'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And, finally, let's create a dataframe with all of our features and save it"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"data_concat = np.concatenate(result, axis=0)\ndel result\ngc.collect()\n\nfeatures_df = pd.DataFrame(data_concat, columns = dtype_dict).astype(dtype_dict)\ndel data_concat\ngc.collect()\n\ndf = pd.read_csv(f'{original_data_folder}/{df_type}.csv')\ndf = df[df['type'].isin(mol_types)] #selecte only the rows with desired coupling types\ndf = reduce_mem_usage(df, verbose=True)\ndf.reset_index(drop = True, inplace=True)\n\ndf = features_df.join(df)\ndel features_df\ngc.collect()\n\nfor t in mol_types:\n    df_only_selected_type = df[df['type'] == t]\n    df_only_selected_type.to_csv(f'{df_type}_{t}.csv')\n    del df_only_selected_type","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at our data"},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head(50)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see all the variables"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(list(df.columns))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check some graphs"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n%matplotlib inline\n\ndf_1JHN = df[df['type'] == '1JHN']\nselected_cols = ['yukawa_Y_Z_resultant_0', 'yukawa_Y_Z_resultant_1', 'yukawa_Y_resultant',\n                 'yukawa_Z_resultant', 'yukawa_momentum', 'inv_EN_Y_Z_resultant_0', 'inv_EN_Y_Z_resultant_1',\n                 'inv_EN_Y_resultant', 'inv_EN_Z_resultant', 'inv_EN_momentum', 'vander_Y_Z_resultant_0',\n                 'vander_Y_Z_resultant_1', 'vander_Y_resultant', 'vander_Z_resultant', 'vander_momentum',\n                 'coulomb_Y_Z_resultant_0', 'coulomb_Y_Z_resultant_1', 'coulomb_Y_resultant',\n                 'coulomb_Z_resultant', 'coulomb_momentum', 'yukawa_parallel', 'coulomb_parallel', \n                 'harmonic_parallel', 'vander_parallel', 'dist_1_parallel', 'inv_dist_1_parallel',\n                 'inv_dist_2_parallel', 'inv_dist_3_parallel', 'yukawa_orthogonal', 'coulomb_orthogonal',\n                 'harmonic_orthogonal', 'vander_orthogonal', 'dist_1_orthogonal', 'inv_dist_1_orthogonal',\n                 'inv_dist_2_orthogonal', 'inv_dist_3_orthogonal', 'yukawa_proportion_between_X_and_Y_Z',\n                 'coulomb_proportion_between_X_and_Y_Z', 'harmonic_proportion_between_X_and_Y_Z',\n                 'vander_proportion_between_X_and_Y_Z', 'dist_1_proportion_between_X_and_Y_Z',\n                 'inv_dist_1_proportion_between_X_and_Y_Z', 'inv_dist_2_proportion_between_X_and_Y_Z',\n                 'inv_dist_3_proportion_between_X_and_Y_Z', '1_yukawa_resultant', '2_yukawa_resultant',\n                '1_coulomb_resultant', '2_coulomb_resultant', '1_inv_EN_resultant', '2_inv_EN_resultant',\n                 '2_G4_3_selected']\n\nfor col in selected_cols:\n    print(col)\n    fig = plt.figure()\n    plt.subplots()\n    plt.scatter(df_1JHN.loc[0::20,'scalar_coupling_constant'], df_1JHN.loc[0::20,col], s=0.2)\n    plt.title(f'{col} VS coupling constant for 1JHN')\n    plt.ylabel(col)\n    plt.xlabel('scalar_coupling_constant')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 8. Calculate features for the test dataset"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df_type = 'test' #or 'train'\n\n################\n#Get all the data\n################\n\ndf = pd.read_csv(f'{original_data_folder}/{df_type}.csv')\ndf = df[df['type'].isin(mol_types)] #selecte only the rows with desired coupling types\ndf = reduce_mem_usage(df, verbose=True)\n\nunique_mol = df['molecule_name'].unique()\n\ndef get_data(data_path, verbose):\n    \n    data = pd.read_csv(data_path)\n    data = data.loc[data['molecule_name'].isin(unique_mol)]\n    data = reduce_mem_usage(data, verbose=True)\n    \n    return data\n\ndf_bonds = get_data(data_path = f'../input/predicting-molecular-properties-bonds/{df_type}_bonds.csv', verbose = verbose)\ndf_structures = get_data(data_path = '../input/champs-scalar-coupling/structures.csv', verbose = verbose)\ndf_structures_eigen = get_data(data_path = '../input/qm7-coulomb-matrix-eigenvalue-features/struct_eigen.csv', verbose = verbose)\nif df_type == 'test':\n    df_mulliken = get_data(data_path = '../input/predicting-mulliken-charges-with-acsf-descriptors/mulliken_charges_test_set.csv', verbose = verbose)\nelse:\n    df_mulliken = get_data(data_path = '../input/champs-scalar-coupling/mulliken_charges.csv', verbose = verbose)\n          \n\n################\n#Get indexes per molecule, so we can work with each molecule at a time.\n#start_indexes = index where the molecule start\n#end_indexes = index where the molecule end\n################\n\nxyz = df_structures[['x','y','z']].values\natom_indexes = df[['atom_index_0', 'atom_index_1']].values\n\nstruct_end_indexes = df_structures.groupby('molecule_name').size().cumsum()\ndf_end_indexes = df[['molecule_name', 'atom_index_0', 'atom_index_1']].groupby('molecule_name').size().cumsum()\ndf_bonds_end_indexes = df_bonds[['molecule_name', 'atom_index_0', 'atom_index_1']].groupby('molecule_name').size().cumsum()\n\nstruct_start_indexes = np.zeros(len(struct_end_indexes) + 1, 'int')\nstruct_start_indexes[1:] = struct_end_indexes\nstruct_start_indexes[0] = 0\n\ndf_start_indexes = np.zeros(len(df_end_indexes) + 1, 'int')\ndf_start_indexes[1:] = df_end_indexes\ndf_start_indexes[0] = 0\n\ndf_bonds_start_indexes = np.zeros(len(df_bonds_end_indexes) + 1, 'int')\ndf_bonds_start_indexes[1:] = df_bonds_end_indexes\ndf_bonds_start_indexes[0] = 0\n\n################\n#run in parallel\n################\n\np = multiprocessing.Pool(number_of_cores)\n\ntry:\n    t1 = time.time()\n    result = p.map(get_mol_feats, list(range(df_structures.molecule_name.nunique())))\n    t2 = time.time()\nexcept Exception as e:\n    print('Caught exception in worker thread (x = d):')\n    traceback.print_exc()\n    print()\n    raise e\n    \np.terminate()\np.join\n\nprint(f'Calculating features took {(t2 - t1)/60:0.2f} minutes')\n\n################\n#save to df\n################\n\ndel df_structures_eigen, df_bonds_end_indexes, df_bonds_start_indexes, xyz, unique_mol, df_start_indexes, df_end_indexes, df_bonds, df, df_structures, df_mulliken, atom_indexes\n\nimport gc\ngc.collect()\n\ndata_concat = np.concatenate(result, axis=0)\ndel result\ngc.collect()\n\nfeatures_df = pd.DataFrame(data_concat, columns = dtype_dict).astype(dtype_dict)\ndel data_concat\ngc.collect()\n\ndf = pd.read_csv(f'{original_data_folder}/{df_type}.csv')\ndf = df[df['type'].isin(mol_types)] #selecte only the rows with desired coupling types\ndf = reduce_mem_usage(df, verbose=True)\ndf.reset_index(drop = True, inplace=True)\n\ndf = features_df.join(df)\ndel features_df\ngc.collect()\n\nfor t in mol_types:\n    df_only_selected_type = df[df['type'] == t]\n    df_only_selected_type.to_csv(f'{df_type}_{t}.csv')\n    del df_only_selected_type","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Done! :)\n    \n<b>Please, if you found the content here interesting, consider upvoting the kernel to reward my time editing and sharing it. Thank you very much.</b>\n\nBelow is the final submission that we made usign the LGB and NN models that we built using the features calculated above."},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv('../input/preds-on-oof-and-test/final_model_submission.csv')\nsub.to_csv('final_sub.csv', index=False)","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.7.3"}},"nbformat":4,"nbformat_minor":1}