{"cells":[{"metadata":{},"cell_type":"markdown","source":"Importing Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set()\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"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.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder\nfrom sklearn import preprocessing","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading Data Into Memory"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\nmulliken_charge = pd.read_csv('../input/mulliken_charges.csv')\nstructures = pd.read_csv('../input/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mapping Molecule Name With Structures"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Map the atom structure data into train and test files\n\ndef map_atom_info(df, atom_idx):\n    df = pd.merge(df, structures, how = 'left',\n                  left_on  = ['molecule_name', f'atom_index_{atom_idx}'],\n                  right_on = ['molecule_name',  'atom_index'])\n    \n    df = df.drop('atom_index', axis=1)\n    df = df.rename(columns={'atom': f'atom_{atom_idx}',\n                            'x': f'x_{atom_idx}',\n                            'y': f'y_{atom_idx}',\n                            'z': f'z_{atom_idx}'})\n    return df\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mapping the Mulliken Charge With Molecule Name"},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_mulliken_charge(df,atom_idx) :\n    df = pd.merge(df,mulliken_charge,how = 'left',\n                 left_on = ['molecule_name',f'atom_index_{atom_idx}'],\n                 right_on = ['molecule_name','atom_index']\n                 )\n    df = df.rename(columns={'mulliken_charge': f'mulliken_charge_{atom_idx}'}\n                  )\n    df = df.drop('atom_index',axis = 1)\n    return df\n\ntrain = map_mulliken_charge(train,0)\ntrain = map_mulliken_charge(train,1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mapping the Potential Energiers With Molecule Name"},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy = pd.read_csv('../input/potential_energy.csv')\ntrain = train.merge(potential_energy, on=\"molecule_name\", how = 'inner')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Computing And Merging Dipole Moments "},{"metadata":{"trusted":true},"cell_type":"code","source":"dipole_moments = pd.read_csv('../input/dipole_moments.csv')\ndipole_moment = np.sqrt(dipole_moments.X ** 2 + dipole_moments.Y ** 2 + dipole_moments.Z ** 2)\ndipole_moments['dipole_moment'] = dipole_moment\ndipole_moments = dipole_moments.drop(['X','Y','Z'],axis = 1)\ntrain = train.merge(dipole_moments,on='molecule_name',how = 'inner')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = reduce_mem_usage(train)\ntest = reduce_mem_usage(test)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm_notebook as tqdm\natomic_radius = {'H':0.38, 'C':0.77, 'N':0.75, 'O':0.73, 'F':0.71} # Without fudge factor\n\nfudge_factor = 0.05\natomic_radius = {k:v + fudge_factor for k,v in atomic_radius.items()}\nprint(atomic_radius)\n\nelectronegativity = {'H':2.2, 'C':2.55, 'N':3.04, 'O':3.44, 'F':3.98}\n\n#structures = pd.read_csv(structures, dtype={'atom_index':np.int8})\n\natoms = structures['atom'].values\natoms_en = [electronegativity[x] for x in tqdm(atoms)]\natoms_rad = [atomic_radius[x] for x in tqdm(atoms)]\n\nstructures['EN'] = atoms_en\nstructures['rad'] = atoms_rad\n\ndisplay(structures.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"i_atom = structures['atom_index'].values\np = structures[['x', 'y', 'z']].values\np_compare = p\nm = structures['molecule_name'].values\nm_compare = m\nr = structures['rad'].values\nr_compare = r\n\nsource_row = np.arange(len(structures))\nmax_atoms = 28\n\nbonds = np.zeros((len(structures)+1, max_atoms+1), dtype=np.int8)\nbond_dists = np.zeros((len(structures)+1, max_atoms+1), dtype=np.float32)\n\nprint('Calculating the bonds')\n\nfor i in tqdm(range(max_atoms-1)):\n    p_compare = np.roll(p_compare, -1, axis=0)\n    m_compare = np.roll(m_compare, -1, axis=0)\n    r_compare = np.roll(r_compare, -1, axis=0)\n    \n    mask = np.where(m == m_compare, 1, 0) #Are we still comparing atoms in the same molecule?\n    dists = np.linalg.norm(p - p_compare, axis=1) * mask\n    r_bond = r + r_compare\n    \n    bond = np.where(np.logical_and(dists > 0.0001, dists < r_bond), 1, 0)\n    \n    source_row = source_row\n    target_row = source_row + i + 1 #Note: Will be out of bounds of bonds array for some values of i\n    target_row = np.where(np.logical_or(target_row > len(structures), mask==0), len(structures), target_row) #If invalid target, write to dummy row\n    \n    source_atom = i_atom\n    target_atom = i_atom + i + 1 #Note: Will be out of bounds of bonds array for some values of i\n    target_atom = np.where(np.logical_or(target_atom > max_atoms, mask==0), max_atoms, target_atom) #If invalid target, write to dummy col\n    \n    bonds[(source_row, target_atom)] = bond\n    bonds[(target_row, source_atom)] = bond\n    bond_dists[(source_row, target_atom)] = dists\n    bond_dists[(target_row, source_atom)] = dists\n\nbonds = np.delete(bonds, axis=0, obj=-1) #Delete dummy row\nbonds = np.delete(bonds, axis=1, obj=-1) #Delete dummy col\nbond_dists = np.delete(bond_dists, axis=0, obj=-1) #Delete dummy row\nbond_dists = np.delete(bond_dists, axis=1, obj=-1) #Delete dummy col\n\nprint('Counting and condensing bonds')\n\nbonds_numeric = [[i for i,x in enumerate(row) if x] for row in tqdm(bonds)]\nbond_lengths = [[dist for i,dist in enumerate(row) if i in bonds_numeric[j]] for j,row in enumerate(tqdm(bond_dists))]\nbond_lengths_mean = [ np.mean(x) for x in bond_lengths]\nn_bonds = [len(x) for x in bonds_numeric]\n\n\nbond_data = {'n_bonds':n_bonds, 'bond_lengths_mean': bond_lengths_mean }\nbond_df = pd.DataFrame(bond_data)\nstructures = structures.join(bond_df)\ndisplay(structures.head(20))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\ntrain['dist_x'] = (train['x_0'] - train['x_1']) ** 2\ntest['dist_x'] = (test['x_0'] - test['x_1']) ** 2\ntrain['dist_y'] = (train['y_0'] - train['y_1']) ** 2\ntest['dist_y'] = (test['y_0'] - test['y_1']) ** 2\ntrain['dist_z'] = (train['z_0'] - train['z_1']) ** 2\ntest['dist_z'] = (test['z_0'] - test['z_1']) ** 2\n\ntrain['type_0'] = train['type'].apply(lambda x: x[0])\ntest['type_0'] = test['type'].apply(lambda x: x[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import gc\ndel structures,mulliken_charge,dipole_moments,potential_energy,dipole_moment\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"molecules = train.pop('molecule_name')\ntest = test.drop('molecule_name', axis=1)\n\nscalar_coupling_constant = train.pop('scalar_coupling_constant')\npotential_energy = train.pop('potential_energy')\nmulliken_charge_0 = train.pop('mulliken_charge_0')\nmulliken_charge_1 = train.pop('mulliken_charge_1')\ndipole_moment = train.pop('dipole_moment')\n\n\nfor f in ['atom_1', 'type_0', 'type','atom_0']:\n        lbl = LabelEncoder()\n        lbl.fit(list(train[f].values) + list(test[f].values))\n        train[f] = lbl.transform(list(train[f].values))\n        test[f] = lbl.transform(list(test[f].values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = train.drop('id',axis = 1)\ntest = test.drop('id',axis = 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = preprocessing.StandardScaler()\ntrain = scaler.fit_transform(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Importing Neural Network Libraries"},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nimport keras.backend as K\nfrom keras import metrics\n\nimport keras\nfrom keras.engine.input_layer import Input\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport random, os, sys\nfrom keras.models import *\nfrom keras.layers import *\nfrom keras.callbacks import *\nfrom keras.initializers import *\nimport tensorflow as tf\nfrom keras.engine.topology import Layer\nfrom keras import callbacks\n\npd.set_option('precision', 30)\nnp.set_printoptions(precision = 30)\n\nnp.random.seed(368)\ntf.set_random_seed(368)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Model : The Architecture I'm using here is sort of a hybrid architecture where I'm Taking outputs at various points and then passing those outputs to other layer.\nThese Outputs are the mulliken charges and the potential energy.The final output is the scaler coupling constant.\nUsing the test set my neural network is trying to find out mulliken charges,potential energy in previous layers in the last layer Im using mulliken charges,potential energy and features computed at previous layers to find out scaler coupling constant."},{"metadata":{"trusted":true},"cell_type":"code","source":"def nn_model() :\n    i  = Input(shape = (24,))\n    \n    # Initial Block\n    x  = Dense(64,activation = 'relu')(i)\n    x  = BatchNormalization()(x)\n    x  = Dense(32,activation = 'relu')(x)\n    x  = BatchNormalization()(x)\n    x  = Dense(16,activation = 'relu')(x)\n    x  = BatchNormalization()(x)\n    \n\n    \n    # Mulliken Charge 0 Block\n    x1 = Dense(64,activation = 'relu')(i)\n    x1 = BatchNormalization()(x1)\n    x1 = Dense(32,activation = 'relu')(x1)\n    x1 = BatchNormalization()(x1)\n    x1 = Dense(16,activation = 'relu')(x1)\n    x1 = BatchNormalization()(x1)\n\n    x1_output = Dense(1,activation = 'linear',name = 'mulliken_charge_0')(x1)\n    \n    \n    # Mulliken Charge 1 Block\n    x2 = Dense(64,activation = 'relu')(i)\n    x2 = BatchNormalization()(x2)\n    x2 = Dense(32,activation = 'relu')(x2)\n    x2 = BatchNormalization()(x2)\n    x2 = Dense(16,activation = 'relu')(x2)\n    x2 = BatchNormalization()(x2)  \n    \n    x2_output = Dense(1,activation = 'linear',name = 'mulliken_charge_1')(x2)\n    \n    # Dipole Moment Block\n    x3 = Dense(128,activation = 'relu')(i)\n    x3 = BatchNormalization()(x3)\n    x3 = Dense(64,activation = 'relu')(x3)\n    x3 = BatchNormalization()(x3)\n    x3 = Dense(32,activation = 'relu')(x3)\n    x3 = BatchNormalization()(x3)\n    x3 = Dense(16,activation = 'relu')(x3)\n    x3 = BatchNormalization()(x3)\n    \n    x3_output = Dense(1,activation = 'linear',name = 'dipole_moment')(x3)\n    \n    concat = concatenate([x,x1_output,x2_output,x3_output])\n    \n    # Scalar Coupling Constant Block\n    x4 = Dense(64,activation = 'relu')(concat)\n    x4 = BatchNormalization()(x4)\n    x4 = Dense(32,activation = 'relu')(x4)\n    x4 = BatchNormalization()(x4)\n    x4 = Dense(16,activation = 'relu')(x4)\n    x4 = BatchNormalization()(x4)\n    \n    x4_output = Dense(1,activation = 'linear',name = 'scaler_coupling_constant')(x4)\n    \n    \n    return Model(inputs = [i] , outputs = [x4_output,x3_output,x2_output,x1_output])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = nn_model()\nmodel.compile(loss='mean_absolute_error', optimizer='adam')\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Model Visualization"},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import SVG\nfrom keras.utils.vis_utils import model_to_dot\nSVG(model_to_dot(model,show_shapes = True).create(prog='dot', format='svg'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = model.fit(x = train,y = [scalar_coupling_constant.values,dipole_moment.values,mulliken_charge_1.values,mulliken_charge_0.values],\n                    validation_split=0.1,epochs=100,verbose=1,batch_size = 1024)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Training History Plots :"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['loss'])\nplt.plot(history.history['val_loss'])\nplt.title('Loss')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n_= plt.legend(['Train','Validation'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['scaler_coupling_constant_loss'])\nplt.plot(history.history['val_scaler_coupling_constant_loss'])\nplt.title('scaler_coupling_constant_loss')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n_= plt.legend(['Train','Validation'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['dipole_moment_loss'])\nplt.plot(history.history['val_dipole_moment_loss'])\nplt.title('dipole_moment_loss')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n_= plt.legend(['Train','Validation'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['mulliken_charge_1_loss'])\nplt.plot(history.history['val_mulliken_charge_1_loss'])\nplt.title('mulliken_charge_1_loss')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n_= plt.legend(['Train','Validation'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['mulliken_charge_0_loss'])\nplt.plot(history.history['val_mulliken_charge_0_loss'])\nplt.title('mulliken_charge_0_loss')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\n_= plt.legend(['Train','Validation'], loc='upper left')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Predicting Using Test Set"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_preds = model.predict(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv('../input/sample_submission.csv', index_col='id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = sample_submission.copy()\npredictions['scalar_coupling_constant'] = y_preds[0]\npredictions.to_csv('submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a href=\"submission.csv\"> Download File </a>\n\n"}],"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}