{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\nfrom sklearn import linear_model\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.model_selection import train_test_split\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.preprocessing import OneHotEncoder\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load all the data into a single dataframe. Use one hot encoding for he atom types"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"le = LabelEncoder()\noh = OneHotEncoder(sparse = False)\ndf = pd.read_csv('../input/train.csv',sep=',', header=0, usecols=['molecule_name','atom_index_0','atom_index_1','type','scalar_coupling_constant'] )\ndf_m_coupling_contributions = pd.read_csv('../input/scalar_coupling_contributions.csv',\n                       sep=',', header=0, usecols=['molecule_name','atom_index_0','atom_index_1','fc','sd','pso', 'dso'])\ndf = pd.merge(df, df_m_coupling_contributions, on=['molecule_name','atom_index_0','atom_index_1'])\ndf_m_diple_moments = pd.read_csv('../input/dipole_moments.csv',sep=',', header=0, usecols=['molecule_name','X','Y','Z'])\ndf = pd.merge(df, df_m_diple_moments, on='molecule_name')\ndf_m_pot_engy = pd.read_csv('../input/potential_energy.csv',sep=',', header=0, usecols=['molecule_name','potential_energy'])\ndf = pd.merge(df, df_m_pot_engy, on='molecule_name')\n\ndf_a_str = pd.read_csv('../input/structures.csv',\n                       sep=',', header=0, usecols=['molecule_name','atom_index','atom','x','y','z'])\n#df_a_str['atom'] = le.fit_transform(df_a_str['atom'])\nf = df_a_str['atom'].values\nf = np.reshape(f, (-1,1))\nf = oh.fit_transform(f)\nohdf = pd.DataFrame(f)\ndf_a_str = pd.concat([df_a_str, ohdf], axis=1)\ndf_a_str = df_a_str.rename(index=str, columns={0: 'A0',1:'A1',2:'A2',3:'A3',4:'A4'})\ndf_a_str.drop(columns=['atom'], inplace = True)\ndf_a_mag_sh_tensor = pd.read_csv('../input/magnetic_shielding_tensors.csv',\n                       sep=',', header=0, usecols=['molecule_name','atom_index','XX','YX','ZX','XY','YY','ZY','XZ','YZ','ZZ'])\n\ndf_a_mlkn_charges = pd.read_csv('../input/mulliken_charges.csv',sep=',', header=0, usecols=['molecule_name','atom_index','mulliken_charge'])\n\ndf_a_str = pd.merge(df_a_str, df_a_mag_sh_tensor, on=['molecule_name','atom_index'])\n\ndf_a_str = pd.merge(df_a_str, df_a_mlkn_charges, on=['molecule_name','atom_index'])\n\ndf_atom_1_prop = df_a_str.rename(index=str, columns={'atom_index': 'atom_index_0','A0':'A0_0','A1':'A1_0','A2':'A2_0','A3':'A3_0','A4':'A4_0','x':'x_0','y':'y_0','z':'z_0', 'XX':'XX_0', 'YX':'YX_0', 'ZX':'ZX_0', 'XY':'XY_0', 'YY':'YY_0', 'ZY':'ZY_0', 'XZ':'XZ_0', 'YZ':'YZ_0', 'ZZ':'ZZ_0', 'mulliken_charge':'mulliken_charge_0'})\ndf = pd.merge(df, df_atom_1_prop, on=['molecule_name','atom_index_0'])\n\ndf_atom_2_prop = df_a_str.rename(index=str, columns={'atom_index': 'atom_index_1','A0':'A0_1','A1':'A1_1','A2':'A2_1','A3':'A3_1','A4':'A4_1','x':'x_1','y':'y_1','z':'z_1', 'XX':'XX_1', 'YX':'YX_1', 'ZX':'ZX_1', 'XY':'XY_1', 'YY':'YY_1', 'ZY':'ZY_1', 'XZ':'XZ_1', 'YZ':'YZ_1', 'ZZ':'ZZ_1', 'mulliken_charge':'mulliken_charge_1'})\ndf = pd.merge(df, df_atom_2_prop, on=['molecule_name','atom_index_1'])\nss = StandardScaler()\ndf[['scalar_coupling_constant', 'fc', 'sd', 'pso', 'dso', 'X', 'Y', 'Z',\n       'potential_energy', 'x_0', 'y_0', 'z_0', 'A0_0', 'A1_0', 'A2_0', 'A3_0',\n       'A4_0', 'XX_0', 'YX_0', 'ZX_0', 'XY_0', 'YY_0', 'ZY_0', 'XZ_0', 'YZ_0',\n       'ZZ_0', 'mulliken_charge_0', 'x_1', 'y_1', 'z_1', 'A0_1', 'A1_1',\n       'A2_1', 'A3_1', 'A4_1', 'XX_1', 'YX_1', 'ZX_1', 'XY_1', 'YY_1', 'ZY_1',\n       'XZ_1', 'YZ_1', 'ZZ_1', 'mulliken_charge_1']] = ss.fit_transform(df[['scalar_coupling_constant', 'fc', 'sd', 'pso', 'dso', 'X', 'Y', 'Z',\n       'potential_energy', 'x_0', 'y_0', 'z_0', 'A0_0', 'A1_0', 'A2_0', 'A3_0',\n       'A4_0', 'XX_0', 'YX_0', 'ZX_0', 'XY_0', 'YY_0', 'ZY_0', 'XZ_0', 'YZ_0',\n       'ZZ_0', 'mulliken_charge_0', 'x_1', 'y_1', 'z_1', 'A0_1', 'A1_1',\n       'A2_1', 'A3_1', 'A4_1', 'XX_1', 'YX_1', 'ZX_1', 'XY_1', 'YY_1', 'ZY_1',\n       'XZ_1', 'YZ_1', 'ZZ_1', 'mulliken_charge_1']])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A molecule has variable number of atom pairs. Each pair has a scalar coupling constant value. This function defines a (m x n x p) matrix for a molecule where m is the number of atom pair, n is 2, since its a pair, and p is the number of feature for each atom. \nA more sophisticated feature generation could help this model to be more accurate."},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_atom_pairs(name, molecule):\n    df = molecule.apply(list)\n    atom_pair_y = np.zeros((df.shape[0], 8))\n    \n    atom_pair = np.zeros((df.shape[0], 2, 18))\n    atom_pair[:,0,:] = df.as_matrix(columns=['x_0','y_0','z_0','XX_0','YX_0','ZX_0','XY_0',\n                                            'YY_0','ZY_0','XZ_0','YZ_0','ZZ_0','mulliken_charge_0',\n                                             'A0_0','A1_0','A2_0','A3_0','A4_0'])\n    atom_pair[:,1,:] = df.as_matrix(columns=['x_1','y_1','z_1','XX_1','YX_1','ZX_1','XY_1',\n                                            'YY_1','ZY_1','XZ_1','YZ_1','ZZ_1','mulliken_charge_1',\n                                             'A0_1','A1_1','A2_1','A3_1','A4_1'])\n   \n    atom_pair_y = df.as_matrix(columns=['potential_energy','X','Y','Z','fc','sd','pso','dso'])\n    return atom_pair, atom_pair_y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"moleculelist = []\nmolecule_ylist = []","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are unable to use the entire dataset, as kaggle has the memory limitation, so 10000 molecules are used to train the model. Bidirectional does not accept variable size input sequence, therefore each molecure is appened with 0 for upto 650 pairs of atoms.\nUsing the entire training dataset will improve the performance of this model."},{"metadata":{"trusted":true},"cell_type":"code","source":"molecules = df.groupby('molecule_name')\nc = 0\nfor name, molecule in molecules:\n    atoms, molecule_y = build_atom_pairs(name, molecule)\n    amolecule = np.zeros((650,atoms.shape[1], atoms.shape[2]))\n    amolecule[:atoms.shape[0],:atoms.shape[1],:atoms.shape[2]] = atoms\n    amolecule = amolecule.transpose([0,2,1]).reshape(amolecule.shape[0], -1)\n    amolecule_y = np.zeros((650,molecule_y.shape[1]))\n    amolecule_y[:molecule_y.shape[0],:molecule_y.shape[1]] = molecule_y\n    moleculelist.append(amolecule)\n    molecule_ylist.append(amolecule_y)\n    c = c + 1\n    if c > 10000:\n        break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras import layers\nfrom keras.optimizers import Adam\nfrom keras.layers import Input, Dense, Activation, BatchNormalization, Flatten\nfrom keras.layers import Dropout\nfrom keras.models import Model\nfrom keras.preprocessing import image\nfrom keras.utils import layer_utils\nfrom keras.utils.data_utils import get_file\nfrom keras.applications.imagenet_utils import preprocess_input\nfrom keras.utils.vis_utils import model_to_dot\nfrom keras.utils import plot_model\nfrom keras.models import Sequential\nfrom keras.layers import Bidirectional\nfrom keras.layers import LSTM","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The model intend to use bidirectional RNN. This is because, if we consider a molecule M, containing p atoms in it, then each of pair combination has a scalar coupling constant and each pair within the molecule has an affect on each other's scalar couplingg constant value. So a bidirectional RNN, ensure that the affect of a pair is distributed among other pairs."},{"metadata":{"trusted":true},"cell_type":"code","source":"def BRNNModel(inputdim):\n    model = Sequential()\n    model.add(Bidirectional(LSTM(100, return_sequences=True, input_dim= inputdim )))\n    #model.add(Bidirectional(LSTM(10)))\n    model.add(Dense(8))\n    model.add(Activation('relu'))\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This code is not used"},{"metadata":{"trusted":true},"cell_type":"code","source":"def batch_generator(X, y, batch_size):\n    number_of_batches = X.shape[0]/batch_size\n    counter=0\n    shuffle_index = np.arange(np.shape(y)[0])\n#    X =  X[shuffle_index, :]\n#    y =  y[shuffle_index]\n    while 1:\n        index_batch = shuffle_index[batch_size*counter:batch_size*(counter+1)]\n        X_batch = X[index_batch,:].todense()\n        y_batch = y[index_batch]\n        counter += 1\n        yield(np.array(X_batch),y_batch)\n        if (counter > number_of_batches):\n#            np.random.shuffle(shuffle_index)\n            counter=0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = np.asarray(moleculelist)\ny = np.asarray(molecule_ylist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test, Y_train, Y_test = train_test_split(X, y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = BRNNModel(X.shape[1])\nmodel.compile(optimizer = 'adam', loss = \"mean_squared_error\", metrics = [\"accuracy\"])#\"adam\"\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.fit(X_train, Y_train, epochs=30, batch_size=16, verbose=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#model.fit_generator(generator=batch_generator(X_train, Y_train, batch_size=32), \n#                    nb_epoch=10, samples_per_epoch=X_train.shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = model.evaluate(x = X_test, y = Y_test)\nprint()\nprint (\"Loss = \" + str(preds[0]))\nprint (\"Test Accuracy = \" + str(preds[1]))","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":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}