{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# 1JHN Molecule Animations\nIn this kernel we display EDA to help us understand how molecular structure predicts `scalar_coupling_constant`. In Kaggle's \"Prediction Molecular Properties\" competition we are given the structure of molecules as `x,y,z` coordinates in 3D. Then we must predict the `scalar_coupling_constant` between two specific atoms within a molecule.\n\nWe will only be exploring 1JHN couplings. So each plot below shows 1 hydrogen atom at (0,0,0) and 1 nitrogen atom at (1,0,0). We learn from EDA that the additional atoms that connect to the nitrogen affect the value of the `sc_constant`. Typically nitrogen connects with 3 atoms. So after connecting with (the original) hydrogen it most often pairs with 2 more. If it only pairs with one more, it has a small `sc_constant`. If it pairs with another hydrogen it has a small `sc_constant`. If it pairs with larger molecules like carbon or more nitrogen, it has a large `sc_constant`.\n\nThe distance of the bond between the original hydrogen and nitrogen correlates with the `sc_constant`. Also the distances to the additional two atoms that nitrogen connects with help predict `sc_constant`. This is seen in EDA below and confirmed when we build a model using these facts with LGBM."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# LOAD LIBRARIES\nimport numpy as np, pandas as pd, sys, time\nfrom matplotlib.animation import ArtistAnimation\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nfrom IPython.display import Image, display\nimport random, os, gc\nimport lightgbm as lgb ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"source_file = '../input/champs-scalar-coupling/structures.csv'\nstructures_df = pd.read_csv(source_file)\nstructures_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_file = '../input/champs-scalar-coupling/train.csv'\ntrain_df = pd.read_csv(train_file)\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Display Examples\nBelow is the histogram of `scalar_coupling_constant` for bond type `1JHN`. The `sc_constant` ranges from 30 to 70. We will display 10 examples from small `sc_constant` to large `sc_constant` and see if we can detect any patterns. In each animation, we position the hydrogen atom at coordinates (0,0,0). Then we scale the molecule so distance between hydrogen and nitrogen is 1. Next we position the molecule so that nitrogen is at coordinates (1,0,0). Finally we animate the molecule by rotating it around the x-axis. By positioning and scaling all molecules into the same standard position, we can now observe patterns of similaries and differences."},{"metadata":{"trusted":true},"cell_type":"code","source":"BTYPE = '1JHN'\nDISPLAY_BONDS = False\ntrain_df = train_df[train_df.type==BTYPE].copy()\nplt.figure(figsize=(20,4))\nsns.distplot(train_df.scalar_coupling_constant)\nplt.title('Bond Type '+BTYPE,fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# DIVIDE SCC RANGE INTO NUM DIVS\nDIVS=10\n# DISPLAY NUM SAMPLES FROM EACH DIV\nSAMPLES=1\n\nmn = train_df.scalar_coupling_constant.min()\nmx = train_df.scalar_coupling_constant.max()\nstep = (mx-mn)/DIVS\n#print(BTYPE,'sc_constant min =',mn,'max =',mx)\n#print('We will display',SAMPLES,'sample(s) from each of',DIVS,'divisions')\n\n# CHOOSE SAMPLES\nsamples_idx = []; \nfor k in range(DIVS):\n    ln = len(train_df[ (train_df.scalar_coupling_constant>mn+k*step)&(train_df.scalar_coupling_constant<mn+(k+1)*step)])\n    pick = np.min((ln,SAMPLES))\n    if pick!=0:\n        samples_idx += train_df[ (train_df.scalar_coupling_constant>mn+k*step)&(train_df.scalar_coupling_constant<mn+(k+1)*step)].sample(pick).index.values.tolist()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# PUT STRUCTURES IN NUMPY HASH TABLE FOR FAST LOOKOUT\n# https://www.kaggle.com/cpmpml/ultra-fast-distance-matrix-computation\nnp_xyz = structures_df[['x','y','z']].values.astype('float32')\n#np_atoms = structures_df['atom'].values\nmp = {'H':0,'C':1,'O':2,'N':3,'F':4}\nnp_atoms = structures_df['atom'].map(mp).values.astype('int8')\nss = structures_df.groupby('molecule_name').size()\nss = ss.cumsum()\nssx = np.zeros(len(ss) + 1, 'int')\nssx[1:] = ss; ssx[0] = 0\nnp_dict = {}\nfor k in range(len(ss)):\n    np_dict[ss.index[k]] = k\n\n# STANDARDIZE ATOM 1 AND 2 INTO (X,Y) = (0,0) AND (1,0)\ndef standardizer(mol_name, t1, t2, phi2=0, a_count=0):\n    # USES GLOBALS NP_ATOMS, NP_XYZ, SSX, NP_DICT\n    atom_df = np_atoms[ ssx[np_dict[mol_name]]:ssx[np_dict[mol_name]+1] ]\n    mol_df = np_xyz[ ssx[np_dict[mol_name]]:ssx[np_dict[mol_name]+1], ].copy()\n    mol_df -= mol_df[t1,]\n    bond_dist = np.sqrt(mol_df[t2,].dot(mol_df[t2,])).astype('float32')\n    mol_df /= bond_dist\n    \n    if mol_df[t2,0]==0: theta = np.sign(mol_df[t2,2])*np.pi/2\n    else: theta = np.arctan(mol_df[t2,2]/mol_df[t2,0])\n    M = np.array([[np.cos(-theta),-np.sin(-theta)],[np.sin(-theta),np.cos(-theta)]])\n    mol_df[:,[0,2]] = (M.dot(mol_df[:,[0,2]].T)).T    \n    \n    if mol_df[t2,0]==0: phi = np.sign(mol_df[t2,1])*np.pi/2\n    else: phi = np.arctan(mol_df[t2,1]/mol_df[t2,0])\n    if mol_df[t2,0] < 0: phi += np.pi\n    M = np.array([[np.cos(-phi),-np.sin(-phi)],[np.sin(-phi),np.cos(-phi)]])\n    mol_df[:,[0,1]] = (M.dot(mol_df[:,[0,1]].T)).T         \n      \n    M = np.array([[np.cos(-phi2),-np.sin(-phi2)],[np.sin(-phi2),np.cos(-phi2)]])\n    mol_df[:,[1,2]] = (M.dot(mol_df[:,[1,2]].T)).T \n            \n    # CLOSEST NEIGHBORS TO ATOM_INDEX_1\n    neighbor_types = [np.nan]*a_count\n    neighbor_dists = [np.nan]*a_count\n    if a_count>0:\n        dst = np.sum( (mol_df-[1,0,0])*(mol_df-[1,0,0]), axis=1 )\n        dst[t1] = 1e3; dst[t2] = 1e3\n        idst = np.argsort(dst)\n        ln = np.min((len(idst)-2,a_count))\n        neighbor_types[:ln] = atom_df[idst[:ln]]\n        neighbor_dists[:ln] = dst[idst[:ln]]\n        \n    return (atom_df, mol_df, bond_dist, neighbor_types, neighbor_dists)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def animate(j,title='',display_bonds=DISPLAY_BONDS):\n    # CONSTANTS\n    atoms = ['H','C','O','N','F']\n    colors = ['white','lightgray','red','blue','green']\n    sizes = [25, 70, 60, 65, 50]\n    lkeys = []\n    for k in range(5):\n        lkeys.append(mpatches.Patch(facecolor=colors[k], label=atoms[k], edgecolor='black'))\n    \n    # EDA PLOTS\n    mol_name = train_df.loc[j,'molecule_name']\n    target_1 = train_df.loc[j,'atom_index_0']\n    target_2 = train_df.loc[j,'atom_index_1']\n    constant = train_df.loc[j,'scalar_coupling_constant']\n    btype = train_df.loc[j,'type']\n    frames = []; scl = 1.0\n    fig = plt.figure(figsize=(8,5))\n    for jj,alpha in enumerate(np.linspace(0,2*np.pi,16)):\n        atom_df, mol_df, dis, _, _ = standardizer(mol_name, target_1, target_2, alpha)\n        if jj==0: # SCALE CIRCLES IN PLOT\n            mxx = np.sqrt(np.max( mol_df[:,1]**2 + mol_df[:,2]**2 ))\n            if btype[0]=='2': mxx *= 1.5\n            elif btype[0]=='3': mxx *= 2.0\n            if mxx<0.1: scl = 200\n            else: scl = 87.*(1+0.2*(4-mxx))/mxx\n            if display_bonds: scl /= 3.\n        # DISPLAY BONDS\n        all = []\n        if display_bonds:\n            ds = mol_df[:,None,:]-mol_df\n            dist = np.sqrt(np.einsum('ijk,ijk->ij',ds,ds))\n            dist = (dist<1.65/dis).astype(np.int8)\n            for k in range(len(mol_df)):\n                for kk in range(k+1,len(mol_df)):\n                    if dist[k,kk]==1:\n                        pp, = plt.plot([mol_df[k,0],mol_df[kk,0]],[mol_df[k,1],mol_df[kk,1]],color='black', zorder=-1 )\n                        all.append( pp )\n        # ARRANGE ATOMS FROM BACK TO FRONT\n        mol_df = np.concatenate([mol_df,atom_df.reshape((-1,1))],axis=1)\n        midx = np.argsort(mol_df[:,2],axis=0)\n        mol_df = mol_df[midx,]\n        atom_df = mol_df[:,3].astype('int8')\n        # DISPLAY ATOMS FROM BACK TO FRONT\n        for k in range(len(mol_df)):\n            all.append( plt.scatter(mol_df[k,0],mol_df[k,1],color=colors[atom_df[k]],edgecolor='black',s=scl*sizes[atom_df[k]]) )\n        xlm = plt.xlim()\n        plt.legend(handles=lkeys)\n        pre = btype\n        if title!='': pre = title\n        plt.title(pre+', bond '+str(j)+', sc_constant = '+str(np.round(constant,2)),fontsize=16)\n        plt.axis('equal')\n        frames.append(all)\n    # MAKE ANIMATION\n    ani = ArtistAnimation(fig, frames)\n    ani.save('data'+str(j)+'.gif', writer='imagemagick', fps=5)\n    plt.close()\n    #print('Bond',j,',Molecule',mol_name)\n    with open('data'+str(j)+'.gif','rb') as file:\n        display(Image(file.read()))\n        \n# DISPLAY OUR 10 SAMPLES\nfor j in samples_idx: animate(j)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We observe that the atoms connected to nitrogen are associated with `sc_constant`. For the largest `sc_constant`, nitrogen is connected to more nitrogen (either 2N or 1N,1C). For middle `sc_constant`, nitrogen is connected to 2 carbons. For small `sc_constant`, nitrogen is connected to only 1 atom or 2 atoms that include hydrogen. (There is an exception where nitrogen connects to 2 carbons and those carbons connect mainly to hydrogen. That case has low `sc_constant`).\n\n# Build LGBM Model\nWe will build an LGBM model using the types of the closest atoms to the original nitrogen and their distances to nitrogen. From the plots above, we see that these features are predictive."},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# INITIAL ARRAYS\nneighbors = 6\nbond_dists = np.zeros((len(train_df),1))\nneighbor_types = np.zeros((len(train_df),neighbors))\nneighbor_dists = np.zeros((len(train_df),neighbors))\nbond_names = []\n\n# COMPUTE FEATURES\nfor i,k in enumerate(train_df.index):\n    bond_names.append(k)\n    _, _, a, b, c  = standardizer(train_df.molecule_name[k], \n        train_df.atom_index_0[k], train_df.atom_index_1[k], phi2=0, a_count=neighbors)\n    bond_dists[i]=a; neighbor_types[i,]=b; neighbor_dists[i,]=c\n    #if i%10000==0: print(i)\n    \n# TRAINING DATA\nX_train = np.concatenate([bond_dists,neighbor_types,neighbor_dists],axis=1)\ny_train = train_df['scalar_coupling_constant'].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# COMPETITION METRIC\ndef mean_log_mae(y_true, y_pred, floor=1e-9):\n    mae = np.mean(np.abs(y_true-y_pred))\n    return \"log_mae\", np.log( np.max(mae, floor) ), False\n\n# CREATE 80% TRAIN AND 20% VALIDATION USING UNIQUE MOLECULES\nmol_names = train_df.molecule_name.unique()\nmols = np.random.choice(mol_names,int(0.8*len(mol_names)),replace=False)\ntrain_df['index'] = np.arange(len(train_df))\nidxT = train_df[train_df.molecule_name.isin(mols)]['index'].values\nidxV = np.setdiff1d( np.arange(len(X_train)), idxT)\n\n# TRAIN LGBM MODEL\nlgbm = lgb.LGBMRegressor(n_estimators=1000, colsample_bytree=0.9, objective='regression', num_leaves=256,\n            max_depth=-1, learning_rate=0.05, metric='mae')  \nh = lgbm.fit(X_train[idxT,], y_train[idxT], eval_set=[(X_train[idxV,], y_train[idxV])], \n             eval_metric=mean_log_mae, verbose=250)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Analysis\nWe notice that there is a linear correlation between hydrogen and nitrogen bond distance and `sc_constant` from the plot below. Using just this bond distance, we can predict `sc_constant` with `mae = 3`. That's pretty good. Now using the types of atoms that nitrogen is additionally connected to and their bond distances we see above from LGBM that we can predict `sc_constant` within `mae = 0.5`. That's very good."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.scatter(bond_dists,y_train)\nplt.title('Distance between N and original H',fontsize=16)\nplt.xlabel('Bond distance')\nplt.ylabel('Bond coupling')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def display_plot(rank,annot=[]):\n    atoms = ['H','C','O','N','F']\n    colors2 = ['lightgray','darkgray','red','blue','green']\n    lkeys = []\n    for k in range(5):\n        lkeys.append(mpatches.Patch(facecolor=colors2[k], label=atoms[k]))\n    postfix=['st','nd','rd','th','th','th','th','th','th','th']\n    plt.figure(figsize=(20,5))\n    for k in range(5):\n        plt.scatter(neighbor_dists[neighbor_types[:,rank]==k,rank],y_train[neighbor_types[:,rank]==k],color=colors2[k])\n    plt.title('Distance between N and '+str(rank+1)+postfix[rank]+ \n              ' closest atom with type identified by color',fontsize=16)\n    plt.xlabel('Bond distance')\n    plt.ylabel('Bond coupling')\n    plt.legend(handles=lkeys)\n    for k in annot:\n        plt.text(k[0],k[1],k[2],fontsize=20)\n    plt.show()\n    \nneigh = 0\ndisplay_plot(neigh,[[0.98,50,'A'],[1.55,30,'B'],[1.75,70,'C'],[2,40,'D']])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The plot above is the distance to the 1st neighbor of nitrogen (excluding the original hydrogen). The color in the plot indicates what type of atom it is. We see that if the closest atom is carbon with bond distance under 1.6, then `sc_constant` is between 30 and 40 (class B). We see that if closest atom is nitrogen, then `sc_constant` is above 60 (class C). If the closest atom is not those 2 cases nor hydrogen with distance 1 (class A above), then it is carbon and it linearly correlates with `sc_constant`. This is class D. From eyeball it looks like `sc_constant = m*x+b` where slope `m= -66` and `b=182`. Below we display classes A and B. We see that the 1st neighbor is hydrogen and carbon respectively."},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE A FROM PICTURE ABOVE\nneigh = 0 \nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==0)&(neighbor_dists[:,neigh]<1.2)).flatten()[5]]\nanimate(bn,'Class A')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE B FROM PICTURE ABOVE\nneigh = 0 \nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==1)&(neighbor_dists[:,neigh]<1.6)).flatten()[5]]\nanimate(bn,'Class B')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The plot below is the distance to the second closest neighbor to nitrogen (excluding the original hydrogen). We observe that if the second neighbor is hydrogen then `sc_constant` is around 40 (class E). If the second neighbor is nitrogen with distance less than 2 (then this new nitrogen is connected to original nitrogen), then `sc_constant` is around 70 (class F). If the second neighbor is nitrogen with distance greater than 5 then that means that the original nitrogen has only 1 connection and this new nitrogen is 2 bonds away (i.e. not connect directly with first nitrogen) and `sc_constant` is around 30 (class J). You get the idea. Below we plot molecules from classes E, F, I. "},{"metadata":{"trusted":true},"cell_type":"code","source":"neigh = 1\ndisplay_plot(neigh,[[0.98,40,'E'],[1.75,70,'F'],[2,50,'G'],[4,33,'H'],[4.8,40,'I'],[5.3,33,'J']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE E FROM PICTURE ABOVE\nneigh = 1\nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==0)&(neighbor_dists[:,neigh]<1.5)).flatten()[5]]\nanimate(bn,'Class E')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE F FROM PICTURE ABOVE\nneigh = 1\nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==3)&(neighbor_dists[:,neigh]<2)).flatten()[5]]\nanimate(bn,'Class F')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE G FROM PICTURE ABOVE\nneigh = 1\nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==2)&(neighbor_dists[:,neigh]>4.5)).flatten()[5]]\nanimate(bn,'Class I')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below are the 3rd closest neighbors to nitrogen. Class K is where carbon is the third closest atom and the distance is less than 2.5. This is only possible if carbon is connected to the nitrogen. Therefore the original nitrogen is connected to at least 4 atoms. This is rare since nitrogen usually connects to 3 atoms. Class K is most likely nitrogen connecting with 3 hydrogens (one is the original) and 1 carbon. The animated plot below confirms this."},{"metadata":{"trusted":true},"cell_type":"code","source":"neigh = 2\ndisplay_plot(neigh,[[2.2,42,'K']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE K FROM PICTURE ABOVE\nneigh = 2\nbn = bond_names[np.argwhere((neighbor_types[:,neigh]==1)&(neighbor_dists[:,neigh]<2.5)).flatten()[2]]\nanimate(bn,'Class K')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below are plots of the 4th and 5th closest neighbors to nitrogen. Classes L and M are interesting. For a molecule to achieve a distance of 10+ and 15+ with the 4th and 5th closest atom respectively means that the molecule must be long and narrow. These classes of molecules have a small `sc_constant` near 30."},{"metadata":{"trusted":true},"cell_type":"code","source":"display_plot(3,[[10,40,'L']])\ndisplay_plot(4,[[15,40,'M']])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Class M above are molecules whose 5th closest atom is distance 15 away from the original nitrogen. The animation below shows that these molecules are long and narrow."},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TYPE M FROM PICTURE ABOVE\nneigh = 4\nbn = bond_names[np.argwhere( neighbor_dists[:,neigh]>14 ).flatten()[5]]\nanimate(bn,'Class M')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\nViewing EDA helps us understand how models are predicting `sc_constant` in Kaggle's \"Prediction Molecular Properties\" competition. There are many patterns between the structure of the molecule, the distances between atoms, and `sc_constant`. If a model could understand all the different structures of molecules and have knowledge of distances, it could achieve a very accurate `sc_constant` prediction.\n\nWe also notice that using this naive approach of computing distances to the original nitrogen doesn't work as well for the 3rd, 4th, 5th closest atom etc. At that point, it's hard to tell whether the atoms are chemically connected or not. Using more features and/or GNN (graph neural network) can make better sense of deeper structure further away from the original nitrogen."}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}