{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import gc\nimport math\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport random\npd.set_option('mode.chained_assignment', None)\nimport keras\nimport tensorflow as tf\nfrom keras import backend as K\n\nK.set_image_data_format('channels_last')\nfrom keras.models import Model, Input, load_model\nfrom keras.layers import Input\nfrom keras.layers.core import Dropout, Lambda\nfrom keras.layers.convolutional import Conv2D, Conv2DTranspose\nfrom keras.layers.pooling import MaxPooling2D\nfrom keras.layers.merge import concatenate\nfrom keras.optimizers import Adam\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom keras.layers.normalization import BatchNormalization\nfrom keras.layers import Activation, Dense\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let us convert the problem into a standard Neural Network Image problem whereby we take an input image and produce an output image - have a look at the many kaggle comptetions that use RLE or run length encoding!\n\nThe idea here is to create \"images\" for every molecule\n\nWe will then create a target output \"image\" for each molecule and you should then train your favourite NN to reproduce the output target images\n\nOnce you have done that, you can grab the actual scalar coupling constants by just looking up there values from the 2d matrix output from your neural network\n\n(As these matrices are symmetric you should average over (x,y) and (y,x) indices)\n\nNote I am producing the smallest possible matrices to save memory - it will be up to you to create a generator to expand them to 29 by 29.\n\nGood Luck"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_absolute_error\ndef group_mean_log_mae(y_true, y_pred, types, floor=1e-9):\n    \"\"\"\n    Fast metric computation for this competition: https://www.kaggle.com/c/champs-scalar-coupling\n    Code is from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n    \"\"\"\n    maes = (y_true-y_pred).abs().groupby(types).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\nmolecules = train.molecule_name.unique()\ntrain=train[train.molecule_name.isin(molecules[:10000])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')\nstructures = structures[structures.molecule_name.isin(train.molecule_name)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"typedict = train[['molecule_name','atom_index_0','atom_index_1','type']]\ntypedict = typedict.set_index(['molecule_name','atom_index_0','atom_index_1'])\ntypedict = typedict['type'].map({'2JHH': 1,\n                                 '3JHH': 2,\n                                 '1JHC': 3,\n                                 '2JHC': 4,\n                                 '3JHC': 5,\n                                 '1JHN': 6,\n                                 '2JHN': 7,\n                                 '3JHN': 8 })\n\ntypedict = typedict.to_dict()\ntypedict[('dsgdb9nsd_000001',1,0)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def TypeMunger(x):\n    if (x.molecule_name,x.atom_index_x,x.atom_index_y) in typedict:\n        return typedict[(x.molecule_name,x.atom_index_x,x.atom_index_y)]\n    elif (x.molecule_name,x.atom_index_y,x.atom_index_x) in typedict:\n        return typedict[(x.molecule_name,x.atom_index_y,x.atom_index_x)]\n    return 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rawimagedata = {}\nsizesofmatrices = {}\n\nfor k,groupdf in tqdm((structures.groupby('molecule_name'))):\n    # I am just mapping the atom types to numerics as an example feel free to one hot encode them\n    groupdf.atom =  groupdf.atom.map({'H': 1, 'C': 2, 'N':3,'O':4,'F':5})\n    inputimage = groupdf.merge(groupdf,on=['molecule_name'],how='outer')\n    #Fermi Contact seems to love r^-3!\n    inputimage['recipdistancecubed'] = 1/np.sqrt((inputimage.x_x-inputimage.x_y)**2+\n                                                 (inputimage.y_x-inputimage.y_y)**2+\n                                                 (inputimage.z_x-inputimage.z_y)**2)**3\n    inputimage.recipdistancecubed = inputimage.recipdistancecubed.replace(np.inf,0)\n    inputimage['H1'] = (inputimage.atom_x==1).astype(int)\n    inputimage['C1'] = (inputimage.atom_x==2).astype(int)\n    inputimage['N1'] = (inputimage.atom_x==3).astype(int)\n    inputimage['O1'] = (inputimage.atom_x==4).astype(int)\n    inputimage['F1'] = (inputimage.atom_x==5).astype(int)\n    inputimage['H2'] = (inputimage.atom_y==1).astype(int)\n    inputimage['C2'] = (inputimage.atom_y==2).astype(int)\n    inputimage['N2'] = (inputimage.atom_y==3).astype(int)\n    inputimage['O2'] = (inputimage.atom_y==4).astype(int)\n    inputimage['F2'] = (inputimage.atom_y==5).astype(int)\n    inputimage['bondtype'] = inputimage.apply(lambda x: TypeMunger(x),axis=1)\n    inputimage['2JHH'] = (inputimage.atom_y==1).astype(int)\n    inputimage['3JHH'] = (inputimage.atom_y==2).astype(int)\n    inputimage['1JHC'] = (inputimage.atom_y==3).astype(int)\n    inputimage['2JHC'] = (inputimage.atom_y==4).astype(int)\n    inputimage['3JHC'] = (inputimage.atom_y==5).astype(int)\n    inputimage['1JHN'] = (inputimage.atom_y==6).astype(int)\n    inputimage['2JHN'] = (inputimage.atom_y==7).astype(int)\n    inputimage['3JHN'] = (inputimage.atom_y==8).astype(int)\n    \n    \n    sizesofmatrices[k] = int(math.sqrt(inputimage.shape[0]))\n    rawimagedata[k] = inputimage[['H1','C1','N1','O1','F1',\n                                  'H2','C2','N2','O2','F2',\n                                  '2JHH','3JHH','1JHC','2JHC','3JHC','1JHN','2JHN','3JHN',\n                                  'recipdistancecubed']].values.reshape(sizesofmatrices[k],sizesofmatrices[k],19)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"targetimages = {}\nfor k,groupdf in tqdm((train.groupby('molecule_name'))):\n\n    outputimage = pd.DataFrame({'molecule_name':k,'atom_index':np.arange(sizesofmatrices[k])})\n    outputimage = outputimage.merge(outputimage,on=['molecule_name'],how='outer')\n    outputimage = outputimage.merge(groupdf,\n                                    left_on=['molecule_name','atom_index_x','atom_index_y'],\n                                    right_on=['molecule_name','atom_index_0','atom_index_1'],how='left')\n    outputimage = outputimage.merge(groupdf,\n                                    left_on=['molecule_name','atom_index_x','atom_index_y'],\n                                    right_on=['molecule_name','atom_index_1','atom_index_0'],how='left')\n    outputimage['sc'] = outputimage.scalar_coupling_constant_x.fillna(0)+outputimage.scalar_coupling_constant_y.fillna(0)\n    targetimages[k] = outputimage[['sc']].values.reshape(sizesofmatrices[k],sizesofmatrices[k],1)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note the output target matrix is symmetric so you will get better results averaging (x,y) with (y,x)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_generator(mymoleculearr, batch_size=16, number_of_batches=None, shuffle=False):\n    if number_of_batches is None:\n        number_of_batches = np.ceil(len(mymoleculearr) / batch_size)\n\n    counter = 0\n    while True:\n        idx_start = batch_size * counter\n        idx_end = batch_size * (counter + 1)\n        x_batch = []\n        y_batch = []\n        for key in mymoleculearr[idx_start:idx_end]:   \n            localimage = rawimagedata[key]\n            a = np.zeros((32,32,19))\n            a[:localimage.shape[0],:localimage.shape[1]] = localimage\n            x_batch.append(a)\n            localmask = targetimages[key]\n            b = np.zeros((32,32,1))\n            b[:localmask.shape[0],:localmask.shape[1]] = localmask\n            y_batch.append(b)\n        \n        yield (np.array(x_batch),np.array(y_batch))\n        counter += 1\n        if (counter == number_of_batches):\n            counter = 0\n            if shuffle:\n                np.random.shuffle(mymoleculearr)\n           \ndef test_generator(testmoleculearr, batch_size=16, number_of_batches=None):\n    if number_of_batches is None:\n        number_of_batches = np.ceil(len(testmoleculearr) / batch_size)\n        print(number_of_batches)\n    counter = 0\n    while True:\n        idx_start = batch_size * counter\n        idx_end = batch_size * (counter + 1)\n        x_batch = []\n        for key in testmoleculearr[idx_start:idx_end]:    \n            localimage = rawimagedata[key]\n            a = np.zeros((32,32,19))\n            a[:localimage.shape[0],:localimage.shape[1]] = localimage\n            x_batch.append(a)\n                \n        yield (np.array(x_batch))\n        counter += 1\n        if (counter == number_of_batches):\n            counter = 0\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"allmolecules = [*targetimages]\nnp.random.shuffle(allmolecules)\nsplitlen = int(len(allmolecules)*.8)\nprint(splitlen)\ntrainmolarr = allmolecules[:splitlen]\nvalmolarr = allmolecules[splitlen:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 16\nNO_OF_TRAINING_IMAGES = len(trainmolarr)\nNO_OF_VAL_IMAGES = len(valmolarr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_sample():\n    batch = next(train_generator(allmolecules, 16))\n    x = batch[0][0]\n    y = batch[1][0]\n    size = (5, 5)\n    plt.figure(figsize=size)\n    plt.imshow(x[:, :, 0], cmap='gray')\n    plt.show()\n    plt.figure(figsize=size)\n    plt.imshow(y[:, :, 0], cmap='gray')\n    plt.show();\n    return","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_sample()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def conv2d_block(input_tensor, n_filters, kernel_size=3, batchnorm=True):\n    # first layer\n    x = Conv2D(filters=n_filters, kernel_size=(kernel_size, kernel_size), kernel_initializer=\"he_normal\",\n               padding=\"same\")(input_tensor)\n    if batchnorm:\n        x = BatchNormalization()(x)\n    x = Activation(\"elu\")(x)\n    # second layer\n    x = Conv2D(filters=n_filters, kernel_size=(kernel_size, kernel_size), kernel_initializer=\"he_normal\",\n               padding=\"same\")(x)\n    if batchnorm:\n        x = BatchNormalization()(x)\n    x = Activation(\"elu\")(x)\n    return x\n\ndef get_unet(input_img, n_filters=16, dropout=0.5, batchnorm=False):\n    # contracting path\n    c1 = conv2d_block(input_img, n_filters=n_filters*1, kernel_size=3, batchnorm=batchnorm)\n    p1 = MaxPooling2D((2, 2)) (c1)\n    p1 = Dropout(dropout*0.5)(p1)\n\n    c2 = conv2d_block(p1, n_filters=n_filters*2, kernel_size=3, batchnorm=batchnorm)\n    p2 = MaxPooling2D((2, 2)) (c2)\n    p2 = Dropout(dropout)(p2)\n\n    c3 = conv2d_block(p2, n_filters=n_filters*4, kernel_size=3, batchnorm=batchnorm)\n    p3 = MaxPooling2D((2, 2)) (c3)\n    p3 = Dropout(dropout)(p3)\n\n    c4 = conv2d_block(p3, n_filters=n_filters*8, kernel_size=3, batchnorm=batchnorm)\n    p4 = MaxPooling2D(pool_size=(2, 2)) (c4)\n    p4 = Dropout(dropout)(p4)\n    \n    c5 = conv2d_block(p4, n_filters=n_filters*16, kernel_size=3, batchnorm=batchnorm)\n    \n    # expansive path\n    u6 = Conv2DTranspose(n_filters*8, (3, 3), strides=(2, 2), padding='same') (c5)\n    u6 = concatenate([u6, c4])\n    u6 = Dropout(dropout)(u6)\n    c6 = conv2d_block(u6, n_filters=n_filters*8, kernel_size=3, batchnorm=batchnorm)\n\n    u7 = Conv2DTranspose(n_filters*4, (3, 3), strides=(2, 2), padding='same') (c6)\n    u7 = concatenate([u7, c3])\n    u7 = Dropout(dropout)(u7)\n    c7 = conv2d_block(u7, n_filters=n_filters*4, kernel_size=3, batchnorm=batchnorm)\n\n    u8 = Conv2DTranspose(n_filters*2, (3, 3), strides=(2, 2), padding='same') (c7)\n    u8 = concatenate([u8, c2])\n    u8 = Dropout(dropout)(u8)\n    c8 = conv2d_block(u8, n_filters=n_filters*2, kernel_size=3, batchnorm=batchnorm)\n\n    u9 = Conv2DTranspose(n_filters*1, (3, 3), strides=(2, 2), padding='same') (c8)\n    u9 = concatenate([u9, c1], axis=3)\n    u9 = Dropout(dropout)(u9)\n    c9 = conv2d_block(u9, n_filters=n_filters*1, kernel_size=3, batchnorm=batchnorm)\n    \n    outputs = Conv2D(1, (1, 1), activation='linear') (c9)\n    model = Model(inputs=[input_img], outputs=[outputs])\n    return model\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def mymetric(y_true,y_pred):\n    y_true = K.flatten(y_true)\n    y_pred = K.flatten(y_pred)\n    a = K.gather(y_true,tf.where((y_true>0)|(y_true<0)))\n    b = K.gather(y_pred,tf.where((y_true>0)|(y_true<0)))\n    return K.mean(K.abs(a-b))\n \ndef myloss(y_true,y_pred):\n    y_true = K.flatten(y_true)\n    y_pred = K.flatten(y_pred)\n    a = K.gather(y_true,tf.where((y_true>0)|(y_true<0)))\n    b = K.gather(y_pred,tf.where((y_true>0)|(y_true<0)))\n    return K.mean(K.square(a-b))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = get_unet(Input((32,32,19)), n_filters=16, dropout=0.0, batchnorm=False)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.compile(optimizer=Adam(lr=0.001),loss=[myloss],metrics=[mymetric] )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 16","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"K.tensorflow_backend._get_available_gpus()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_gen = train_generator(trainmolarr,  batch_size=BATCH_SIZE, number_of_batches=None, shuffle=True)\nval_gen = train_generator(valmolarr, batch_size=BATCH_SIZE, number_of_batches=None, shuffle=False)\n\ncallbacks = [\n            EarlyStopping(monitor='val_loss', patience=10, verbose=0),\n       ]\n\nresults = model.fit_generator(train_gen, epochs=100, \n                          steps_per_epoch = len(trainmolarr)//BATCH_SIZE,\n                          validation_data=val_gen, \n                          validation_steps =  len(valmolarr)//BATCH_SIZE,   \n                          callbacks=callbacks)\nmodel.save('Model.h5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#model.load_weights('Model.h5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"val_gen = test_generator(valmolarr, batch_size=BATCH_SIZE, number_of_batches = np.ceil(len(valmolarr) / BATCH_SIZE))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"out = model.predict_generator(val_gen, np.ceil(len(valmolarr) / BATCH_SIZE))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"out.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lookuptable = {}\nfor i,m in enumerate(valmolarr[:out.shape[0]]):\n    lookuptable[m] = out[i]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"o =train[train.molecule_name.isin(valmolarr[:out.shape[0]])]\no.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def LookUp(a):\n\n    return ((lookuptable[a.molecule_name][a.atom_index_0][a.atom_index_1]+lookuptable[a.molecule_name][a.atom_index_1][a.atom_index_0])/2)[0]\n\no['predictions'] = o.apply(lambda x: LookUp(x),axis=1)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"o[o.molecule_name==valmolarr[0]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for typ in o.type.unique():\n    print(typ,group_mean_log_mae(o[o.type==typ].scalar_coupling_constant,\n                               o[o.type==typ].predictions,o[o.type==typ].type))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"group_mean_log_mae(o.scalar_coupling_constant,o.predictions,o.type)","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}