{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":" Code based on  :https://www.kaggle.com/todnewman/keras-neural-net-for-champs\n>  and https://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481"},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\")\nwarnings.filterwarnings(action=\"ignore\",category=DeprecationWarning)\nwarnings.filterwarnings(action=\"ignore\",category=FutureWarning)\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nimport os\nimport math\nimport gc\nimport copy\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler, MinMaxScaler\n\n\nfrom keras.layers import Dense, Input, Activation\nfrom keras.layers import BatchNormalization,Add,Dropout\nfrom keras.optimizers import Adam\nfrom keras.models import Model, load_model\nfrom keras import callbacks\nfrom keras import backend as K\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_PATH = '../input/champs-scalar-coupling'\nSUBMISSIONS_PATH = './'\n# use atomic numbers to recode atomic names\nATOMIC_NUMBERS = {\n    'H': 1,\n    'C': 6,\n    'N': 7,\n    'O': 8,\n    'F': 9\n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_dtypes = {\n    'molecule_name': 'category',\n    'atom_index_0': 'int8',\n    'atom_index_1': 'int8',\n    'type': 'category',\n    'scalar_coupling_constant': 'float32'\n}\ntrain_csv = pd.read_csv(f'{DATA_PATH}/train.csv', index_col='id', dtype=train_dtypes)\ntrain_csv['molecule_index'] = train_csv.molecule_name.str.replace('dsgdb9nsd_', '').astype('int32')\ntrain_csv = train_csv[['molecule_index', 'atom_index_0', 'atom_index_1', 'type', 'scalar_coupling_constant']]\ntrain_csv.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submit = pd.read_csv(f'{DATA_PATH}/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_csv = pd.read_csv(f'{DATA_PATH}/test.csv', index_col='id', dtype=train_dtypes)\ntest_csv['molecule_index'] = test_csv['molecule_name'].str.replace('dsgdb9nsd_', '').astype('int32')\ntest_csv = test_csv[['molecule_index', 'atom_index_0', 'atom_index_1', 'type']]\ntest_csv.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_dtypes = {\n    'molecule_name': 'category',\n    'atom_index': 'int8',\n    'atom': 'category',\n    'x': 'float32',\n    'y': 'float32',\n    'z': 'float32'\n}\nstructures_csv = pd.read_csv(f'{DATA_PATH}/structures.csv', dtype=structures_dtypes)\nstructures_csv['molecule_index'] = structures_csv.molecule_name.str.replace('dsgdb9nsd_', '').astype('int32')\nstructures_csv = structures_csv[['molecule_index', 'atom_index', 'atom', 'x', 'y', 'z']]\nstructures_csv['atom'] = structures_csv['atom'].replace(ATOMIC_NUMBERS).astype('int8')\nstructures_csv.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_type_dataframes(base, structures, coupling_type):\n    base = base[base['type'] == coupling_type].drop('type', axis=1).copy()\n    base = base.reset_index()\n    base['id'] = base['id'].astype('int32')\n    structures = structures[structures['molecule_index'].isin(base['molecule_index'])]\n    return base, structures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def add_coordinates(base, structures, index):\n    df = pd.merge(base, structures, how='inner',\n                  left_on=['molecule_index', f'atom_index_{index}'],\n                  right_on=['molecule_index', 'atom_index']).drop(['atom_index'], axis=1)\n    df = df.rename(columns={\n        'atom': f'atom_{index}',\n        'x': f'x_{index}',\n        'y': f'y_{index}',\n        'z': f'z_{index}'\n    })\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def add_atoms(base, atoms):\n    df = pd.merge(base, atoms, how='inner',\n                  on=['molecule_index', 'atom_index_0', 'atom_index_1'])\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def merge_all_atoms(base, structures):\n    df = pd.merge(base, structures, how='left',\n                  left_on=['molecule_index'],\n                  right_on=['molecule_index'])\n    df = df[(df.atom_index_0 != df.atom_index) & (df.atom_index_1 != df.atom_index)]\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def add_center(df):\n    df['x_c'] = ((df['x_1'] + df['x_0']) * np.float32(0.5))\n    df['y_c'] = ((df['y_1'] + df['y_0']) * np.float32(0.5))\n    df['z_c'] = ((df['z_1'] + df['z_0']) * np.float32(0.5))\n\ndef add_distance_to_center(df):\n    df['d_c'] = ((\n        (df['x_c'] - df['x'])**np.float32(2) +\n        (df['y_c'] - df['y'])**np.float32(2) + \n        (df['z_c'] - df['z'])**np.float32(2)\n    )**np.float32(0.5))\n\ndef add_distance_between(df, suffix1, suffix2):\n    df[f'd_{suffix1}_{suffix2}'] = ((\n        (df[f'x_{suffix1}'] - df[f'x_{suffix2}'])**np.float32(2) +\n        (df[f'y_{suffix1}'] - df[f'y_{suffix2}'])**np.float32(2) + \n        (df[f'z_{suffix1}'] - df[f'z_{suffix2}'])**np.float32(2)\n    )**np.float32(0.5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef add_distances(df):\n    n_atoms = 1 + max([int(c.split('_')[1]) for c in df.columns if c.startswith('x_')])\n    \n    for i in range(1, n_atoms):\n        for vi in range(min(4, i)):\n            add_distance_between(df, i, vi)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def add_n_atoms(base, structures):\n    dfs = structures['molecule_index'].value_counts().rename('n_atoms').to_frame()\n    return pd.merge(base, dfs, left_on='molecule_index', right_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_couple_dataframe(some_csv, structures_csv, coupling_type, n_atoms=15):\n    base, structures = build_type_dataframes(some_csv, structures_csv, coupling_type)\n    base = add_coordinates(base, structures, 0)\n    base = add_coordinates(base, structures, 1)\n    \n    base = base.drop(['atom_0', 'atom_1'], axis=1)\n    atoms = base.drop('id', axis=1).copy()\n    if 'scalar_coupling_constant' in some_csv:\n        atoms = atoms.drop(['scalar_coupling_constant'], axis=1)\n        \n    add_center(atoms)\n    atoms = atoms.drop(['x_0', 'y_0', 'z_0', 'x_1', 'y_1', 'z_1'], axis=1)\n\n    atoms = merge_all_atoms(atoms, structures)\n    \n    add_distance_to_center(atoms)\n    \n    atoms = atoms.drop(['x_c', 'y_c', 'z_c', 'atom_index'], axis=1)\n    atoms.sort_values(['molecule_index', 'atom_index_0', 'atom_index_1', 'd_c'], inplace=True)\n    atom_groups = atoms.groupby(['molecule_index', 'atom_index_0', 'atom_index_1'])\n    atoms['num'] = atom_groups.cumcount() + 2\n    atoms = atoms.drop(['d_c'], axis=1)\n    atoms = atoms[atoms['num'] < n_atoms]\n\n    atoms = atoms.set_index(['molecule_index', 'atom_index_0', 'atom_index_1', 'num']).unstack() \n    atoms.columns = [f'{col[0]}_{col[1]}' for col in atoms.columns]\n    atoms = atoms.reset_index()\n    \n    # downcast back to int8\n    for col in atoms.columns:\n        if col.startswith('atom_'):\n            atoms[col] = atoms[col].fillna(0).astype('int8')\n            \n    atoms['molecule_index'] = atoms['molecule_index'].astype('int32')\n    \n    full = add_atoms(base, atoms)\n    add_distances(full)\n    full.sort_values('id', inplace=True)\n    \n    return full","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def take_n_atoms(df, n_atoms, four_start=4):\n    labels = []\n    for i in range(2, n_atoms):\n        label = f'atom_{i}'\n        labels.append(label)\n\n    for i in range(n_atoms):\n        num = min(i, 4) if i < four_start else 4\n        for j in range(num):\n            labels.append(f'd_{i}_{j}')\n    if 'scalar_coupling_constant' in df:\n        labels.append('scalar_coupling_constant')\n    return df[labels]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_nn_model(input_shape):\n    inp = Input(shape=(input_shape,))\n    x = Dense(2048, activation=\"relu\")(inp)\n    x = BatchNormalization()(x)\n    #x = Dropout(0.4)(x)\n    x = Dense(1024, activation=\"relu\")(x)\n    x = BatchNormalization()(x)\n    x = Dense(1024, activation=\"relu\")(x)\n    x = BatchNormalization()(x)\n    x = Dense(512, activation=\"relu\")(x)\n    x = BatchNormalization()(x)\n    x = Dense(512, activation=\"relu\")(x)\n    x = BatchNormalization()(x)\n    out = Dense(1, activation=\"linear\")(x)  \n   # out1 = Dense(2, activation=\"linear\")(x)#mulliken charge 2\n   # out2 = Dense(6, activation=\"linear\")(x)#tensor 6(xx,yy,zz)\n   # out3 = Dense(12, activation=\"linear\")(x)#tensor 12(others) \n   # out4 = Dense(1, activation=\"linear\")(x)#scalar_coupling_constant \n    #model = Model(inputs=inp, outputs=[out,out1,out2,out3,out4])\n    model = Model(inputs=inp, outputs=[out])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_history(history, label):\n    plt.plot(history.history['loss'])\n    plt.plot(history.history['val_loss'])\n    plt.title('Loss for %s' % label)\n    plt.ylabel('Loss')\n    plt.xlabel('Epoch')\n    _= plt.legend(['Train','Validation'], loc='upper left')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Set up GPU preferences\nconfig = tf.ConfigProto( device_count = {'GPU': 1 , 'CPU': 2} ) \nconfig.gpu_options.allow_growth = True\nconfig.gpu_options.per_process_gpu_memory_fraction = 0.6\nsess = tf.Session(config=config) \nK.set_session(sess)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"from datetime import datetime\n\nmol_types=train_csv[\"type\"].unique()\ncv_score=[]\ncv_score_total=0\nepoch_n = 1000\nverbose = 1\nbatch_size = 2048\n    \n# Set to True if we want to train from scratch.  False will reuse saved models as a starting point.\nretrain =True\n\nstart_time=datetime.now()\ntest_prediction=np.zeros(len(test_csv))\ninput_features = ['atom_2', 'atom_3', 'atom_4', 'atom_5', 'atom_6', 'atom_7',\n       'atom_8','atom_9', 'atom_10','d_1_0', 'd_2_0', 'd_2_1', 'd_3_0',\n       'd_3_1', 'd_3_2', 'd_4_0', 'd_4_1', 'd_4_2', 'd_4_3', 'd_5_0',\n       'd_5_1', 'd_5_2', 'd_5_3', 'd_6_0', 'd_6_1', 'd_6_2', 'd_6_3',\n       'd_7_0', 'd_7_1', 'd_7_2', 'd_7_3', 'd_8_0', 'd_8_1', 'd_8_2',\n       'd_8_3', 'd_9_0', 'd_9_1', 'd_9_2', 'd_9_3', 'd_10_0', 'd_10_1', 'd_10_2',\n       'd_10_3']\n\n\n\n# Loop through each molecule type\nfor mol_type in mol_types:\n\n    model_name_wrt = ('/kaggle/working/molecule_model_%s.hdf5' % mol_type)\n    print('Training %s' % mol_type, 'out of', mol_types, '\\n')\n\n    full = build_couple_dataframe(train_csv, structures_csv, mol_type, n_atoms=11)\n    full2 = build_couple_dataframe(test_csv, structures_csv, mol_type, n_atoms=11)\n    df_train_ = take_n_atoms(full, 11)\n    df_test_ = take_n_atoms(full2, 11)\n    df_train_  = df_train_.fillna(0)\n    df_test_  = df_test_.fillna(0)\n    \n    # Standard Scaler from sklearn does seem to work better here than other Scalers\n    input_data=StandardScaler().fit_transform(pd.concat([df_train_.loc[:,input_features],df_test_.loc[:,input_features]]))   \n    #input_data=StandardScaler().fit_transform(df_train_.loc[:,input_features])\n    target_data=df_train_.loc[:,\"scalar_coupling_constant\"].values\n\n    # Simple split to provide us a validation set to do our CV checks with\n    train_index, cv_index = train_test_split(np.arange(len(df_train_)),random_state=111, test_size=0.1)\n    # Split all our input and targets by train and cv indexes\n    train_target=target_data[train_index]\n    cv_target=target_data[cv_index]\n    train_input=input_data[train_index]\n    cv_input=input_data[cv_index]\n    test_input=input_data[len(df_train_):,:]\n\n    # Build the Neural Net\n    nn_model=create_nn_model(train_input.shape[1])\n    \n    # If retrain==False, then we load a previous saved model as a starting point.\n    if not retrain:\n        nn_model = load_model(model_name_rd)\n        \n    nn_model.compile(loss='mae', optimizer=Adam())#, metrics=[auc])\n    \n    # Callback for Early Stopping... May want to raise the min_delta for small numbers of epochs\n    es = callbacks.EarlyStopping(monitor='val_loss', min_delta=0.0001, patience=40,verbose=1, mode='auto', restore_best_weights=True)\n    # Callback for Reducing the Learning Rate... when the monitor levels out for 'patience' epochs, then the LR is reduced\n    rlr = callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.08,patience=28, min_lr=1e-6, mode='auto', verbose=1)\n    # Save the best value of the model for future use\n    sv_mod = callbacks.ModelCheckpoint(model_name_wrt, monitor='val_loss', save_best_only=True, period=1)\n    history = nn_model.fit(train_input,[train_target], \n            validation_data=(cv_input,[cv_target]), \n            callbacks=[es, rlr, sv_mod], epochs=epoch_n, batch_size=batch_size, verbose=verbose)\n    \n    cv_predict=nn_model.predict(cv_input)\n    plot_history(history, mol_type)\n    accuracy=np.mean(np.abs(cv_target-cv_predict[:,0]))\n    print(np.log(accuracy))\n    cv_score.append(np.log(accuracy))\n    cv_score_total+=np.log(accuracy)\n    \n    # Predict on the test data set using our trained model\n    test_predict=nn_model.predict(test_input)\n    \n    # for each molecule type we'll grab the predicted values\n    test_prediction[test_csv[\"type\"]==mol_type]=test_predict[:,0]\n    K.clear_session()\n\ncv_score_total/=len(mol_types)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print ('Total training time: ', datetime.now() - start_time)\n\ni=0\nfor mol_type in mol_types: \n    print(mol_type,\": cv score is \",cv_score[i])\n    i+=1\nprint(\"total cv score is\",cv_score_total)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def submits(predictions):\n \n    submit[\"scalar_coupling_constant\"] = predictions\n    submit.to_csv(\"/kaggle/working/submission.csv\", index=False)\nsubmits(test_prediction)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add more layers to get a better score! However,maybe,features are really more important than algorithms..."}],"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}