{"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 os\nimport math\nimport tqdm\nimport numpy as np  # linear algebra\nimport pandas as pd  # data processing, CSV file I/O (e.g. pd.read_csv)\n\nos.listdir('/kaggle')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class Coupling:\n    def __init__(self, index, atom_index_0, atom_index_1, coupling_type, scalar_coupling_constant=None):\n        self.index = index\n        self.atom_index_0 = atom_index_0\n        self.atom_index_1 = atom_index_1\n        self.coupling_type = coupling_type\n        self.scalar_coupling_constant = scalar_coupling_constant\n\nclass Atom:\n    def __init__(self, index, name, x, y, z):\n        self.index = index\n        self.name = name\n        self.x = x\n        self.y = y\n        self.z = z\n\n        self.coupling_ids = []\n\nclass Molecule:\n    def __init__(self, name):\n        self.name = name\n        self.atoms = {}\n        self.couplings = []\n\n    def add_coupling(self, coupling):\n        self.couplings.append(coupling)\n        self.atoms[coupling.atom_index_0].coupling_ids.append(len(self.couplings) - 1)\n        self.atoms[coupling.atom_index_1].coupling_ids.append(len(self.couplings) - 1)\n\ndef get_molecule_index_from_name(name):\n    return int(name.split('_')[-1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_structures(filename):\n    molecules = {}\n    atom_names = []\n\n    with open(filename, 'r') as f:\n        structures_df = pd.read_csv(f)\n        atom_names = list(structures_df['atom'].unique())\n\n        with tqdm.tqdm(total=structures_df.shape[0]) as tqdm_bar:\n            for i, row in structures_df.iterrows():\n                molecule_index = get_molecule_index_from_name(row.values[0])\n                atom_index = row.values[1]\n                atom_name = row.values[2]\n                x = row.values[3]\n                y = row.values[4]\n                z = row.values[5]\n\n                if molecule_index not in molecules.keys():\n                    molecules[molecule_index] = Molecule(molecule_index)\n\n                molecules[molecule_index].atoms[atom_index] = Atom(atom_index, atom_name, x, y, z)\n                tqdm_bar.update(1)\n\n    return molecules, atom_names\n\nprint('Reading structures...')\nmolecules, atom_names = get_structures('/kaggle/input/champs-scalar-coupling/structures.csv')\nmolecule_count = len(molecules)\nprint('Molecules: {}'.format(molecule_count))\nprint('Atom names: {}'.format(atom_names))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_train_data(filename):\n    train_molecule_indexes = []\n\n    with open(filename, 'r') as f:\n        train_df = pd.read_csv(f)\n        coupling_types = list(train_df['type'].unique())\n\n        with tqdm.tqdm(total=train_df.shape[0]) as tqdm_bar:\n            for i, row in train_df.iterrows():\n                coupling_index = row.values[0]\n                molecule_index = get_molecule_index_from_name(row.values[1])\n                atom_index_0 = row.values[2]\n                atom_index_1 = row.values[3]\n                coupling_type = row.values[4]\n                scalar_coupling_constant = row.values[5]\n\n                coupling = Coupling(coupling_index, atom_index_0, atom_index_1, coupling_type,\n                    scalar_coupling_constant)\n                molecules[molecule_index].add_coupling(coupling)\n                train_molecule_indexes.append(molecule_index)\n                tqdm_bar.update(1)\n\n    train_molecule_indexes = sorted(list(set(train_molecule_indexes)))\n    return train_molecule_indexes, coupling_types\n\nprint('Reading train data...')\ntrain_molecule_indexes, coupling_types = get_train_data('/kaggle/input/champs-scalar-coupling/train.csv')\ntrain_molecule_count = len(train_molecule_indexes)\nprint('Train molecules: {}'.format(train_molecule_count))\nprint('Coupling types: {}'.format(coupling_types))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_test_data(filename):\n    test_molecule_indexes = []\n\n    with open(filename, 'r') as f:\n        test_df = pd.read_csv(f)\n        coupling_types = list(test_df['type'].unique())\n\n        with tqdm.tqdm(total=test_df.shape[0]) as tqdm_bar:\n            for i, row in test_df.iterrows():\n                coupling_index = row.values[0]\n                molecule_index = get_molecule_index_from_name(row.values[1])\n                atom_index_0 = row.values[2]\n                atom_index_1 = row.values[3]\n                coupling_type = row.values[4]\n\n                coupling = Coupling(coupling_index, atom_index_0, atom_index_1, coupling_type)\n                molecules[molecule_index].add_coupling(coupling)\n                test_molecule_indexes.append(molecule_index)\n                tqdm_bar.update(1)\n\n    test_molecule_indexes = sorted(list(set(test_molecule_indexes)))\n    return test_molecule_indexes\n\nprint('Reading test data...')\ntest_molecule_indexes = get_test_data('/kaggle/input/champs-scalar-coupling/test.csv')\ntest_molecule_count = len(test_molecule_indexes)\nprint('Test molecules: {}'.format(test_molecule_count))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max_adjacent_couplings = 9\ntrain_coupling_count = 0\ntest_coupling_count = 0\n\nfor molecule_index, molecule in molecules.items():\n    if molecule_index in train_molecule_indexes:\n        train_coupling_count += len(molecule.couplings)\n    else:\n        test_coupling_count += len(molecule.couplings)\n\nprint('Max adjacent couplings: {}'.format(max_adjacent_couplings))\nprint('Train couplings count: {}'.format(train_coupling_count))\nprint('Test couplings count: {}'.format(test_coupling_count))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"one_input_dim = 7\nside_input_count = max_adjacent_couplings + 1\n\nX_train = []\nfor i in range(2):\n    for j in range(side_input_count):\n        X_train.append(np.zeros(shape=(train_coupling_count, one_input_dim), dtype=np.float))\n\nY_train = np.zeros(shape=(train_coupling_count), dtype=np.float)\n\nX_test = []\nfor i in range(2):\n    for j in range(side_input_count):\n        X_test.append(np.zeros(shape=(test_coupling_count, one_input_dim), dtype=np.float))\n\ndef get_distance(x0, y0, z0, x1, y1, z1):\n    return math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2 + (z1 - z0) ** 2)\n\ndef get_angle_cos(x0, y0, z0, x1, y1, z1):\n    return (x0 * x1 + y0 * y1 + z0 * z1) / math.sqrt((x0 ** 2 + y0 ** 2 + z0 ** 2) * (x1 ** 2 + y1 ** 2 + z1 ** 2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_dataset(X_arr, molecule_indexes, total_coupling_count, with_output):\n    i = 0\n\n    with tqdm.tqdm(total=total_coupling_count) as tqdm_bar:\n        for molecule_index in molecule_indexes:\n            molecule = molecules[molecule_index]\n            coupling_count = len(molecule.couplings)\n\n            for coupling_index in range(coupling_count):\n                atom_index_0 = molecule.couplings[coupling_index].atom_index_0\n                atom_index_1 = molecule.couplings[coupling_index].atom_index_1\n\n                c_x = molecule.atoms[atom_index_1].x - molecule.atoms[atom_index_0].x\n                c_y = molecule.atoms[atom_index_1].y - molecule.atoms[atom_index_0].y\n                c_z = molecule.atoms[atom_index_1].z - molecule.atoms[atom_index_0].z\n\n                if with_output:\n                    Y_train[i] = molecule.couplings[coupling_index].scalar_coupling_constant\n\n                for side_index, atom_index in enumerate([ atom_index_0, atom_index_1 ]):\n                    coupling_type = molecule.couplings[coupling_index].coupling_type\n                    atom_name = molecule.atoms[atom_index].name\n\n                    x0 = molecule.atoms[atom_index].x\n                    y0 = molecule.atoms[atom_index].y\n                    z0 = molecule.atoms[atom_index].z\n\n                    j = 0\n                    X_arr[side_index * side_input_count + j][i, 0] = coupling_types.index(coupling_type)\n                    X_arr[side_index * side_input_count + j][i, 1] = atom_names.index(atom_name)\n                    X_arr[side_index * side_input_count + j][i, 2] = x0 - c_x\n                    X_arr[side_index * side_input_count + j][i, 3] = y0 - c_y\n                    X_arr[side_index * side_input_count + j][i, 4] = z0 - c_z\n                    X_arr[side_index * side_input_count + j][i, 5] = get_distance(c_x, c_y, c_z, x0, y0, z0)\n                    X_arr[side_index * side_input_count + j][i, 6] = get_angle_cos(c_x, c_y, c_z, x0, y0, z0)\n\n                    distances = []\n\n                    for coupling_index_to in molecule.atoms[atom_index].coupling_ids:\n                        if coupling_index != coupling_index_to:\n                            _atom_index_0 = molecule.couplings[coupling_index_to].atom_index_0\n                            _atom_index_1 = molecule.couplings[coupling_index_to].atom_index_1\n\n                            atom_index_to = 0\n\n                            if atom_index != _atom_index_0:\n                                atom_index_to = _atom_index_0\n                            else:\n                                atom_index_to = _atom_index_1\n\n                            x1 = molecule.atoms[atom_index_to].x\n                            y1 = molecule.atoms[atom_index_to].y\n                            z1 = molecule.atoms[atom_index_to].z\n\n                            distance = get_distance(x0, y0, z0, x1, y1, z1)\n                            angle = get_angle_cos(x0, y0, z0, x1, y1, z1)\n                            distances.append(tuple((coupling_index_to, atom_index_to, distance, angle)))\n\n                    distances.sort(key=lambda x: x[2])\n\n                    j += 1\n                    for coupling_index_to, atom_index_to, distance, angle in distances[:max_adjacent_couplings]:\n                        atom_name = molecule.atoms[atom_index_to].name\n                        coupling_type = molecule.couplings[coupling_index_to].coupling_type\n\n                        x1 = molecule.atoms[atom_index_to].x\n                        y1 = molecule.atoms[atom_index_to].y\n                        z1 = molecule.atoms[atom_index_to].z\n\n                        X_arr[side_index * side_input_count + j][i, 0] = coupling_types.index(coupling_type)\n                        X_arr[side_index * side_input_count + j][i, 1] = atom_names.index(atom_name)\n                        X_arr[side_index * side_input_count + j][i, 2] = x1 - c_x\n                        X_arr[side_index * side_input_count + j][i, 3] = y1 - c_y\n                        X_arr[side_index * side_input_count + j][i, 4] = z1 - c_z\n                        X_arr[side_index * side_input_count + j][i, 5] = distance\n                        X_arr[side_index * side_input_count + j][i, 6] = angle\n                        j += 1\n\n                i += 1\n                tqdm_bar.update(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Creating train dataset...')\nmake_dataset(X_train, train_molecule_indexes, train_coupling_count, with_output=True)\nprint('X_train: {} x {}'.format(X_train[0].shape, len(X_train)))\nprint('Y_train: {}'.format(Y_train.shape))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Creating test dataset...')\nmake_dataset(X_test, test_molecule_indexes, test_coupling_count, with_output=False)\nprint('X_test: {} x {}'.format(X_test[0].shape, len(X_test)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.models import Model\nfrom keras.layers import *\nfrom keras.callbacks import EarlyStopping, ReduceLROnPlateau\n\ndef get_model(input_dims):\n    input_now = Input(shape=(input_dims[0],))\n    x = Dense(48)(input_now)\n    x = LeakyReLU(alpha=0.05)(x)\n    x = Dense(16)(x)\n    x = LeakyReLU(alpha=0.05)(x)\n    model_one = Model(inputs=input_now, outputs=x, name='model_one')\n\n    input_now = Input(shape=(input_dims[1],))\n    x = Dense(48)(input_now)\n    x = LeakyReLU(alpha=0.05)(x)\n    x = Dense(16)(x)\n    x = LeakyReLU(alpha=0.05)(x)\n    model_two = Model(inputs=input_now, outputs=x, name='model_two')\n\n    input_now = Input(shape=(160,))\n    x = Dense(1024)(input_now)\n    x = LeakyReLU(alpha=0.05)(x)\n    x = Dense(512)(x)\n    x = LeakyReLU(alpha=0.05)(x)\n    x = Dense(256)(x)\n    x = LeakyReLU(alpha=0.05)(x)\n    model_side = Model(inputs=input_now, outputs=x, name='model_side')\n\n    inputs = []\n    sides = []\n\n    for i, side_dims in enumerate([ input_dims[:len(input_dims) // 2], input_dims[len(input_dims) // 2:] ]):\n        nodes = []\n\n        input_now = Input(shape=(side_dims[0],), name='input_0_{}'.format(i))\n        x = model_one(input_now)\n        x = BatchNormalization(name='batchnorm_0_{}'.format(i))(x)\n        inputs.append(input_now)\n        nodes.append(x)\n\n        for j, input_dim in enumerate(side_dims[1:]):\n            input_now = Input(shape=(input_dim,), name='input_1_{}_{}'.format(i, j))\n            x = model_two(input_now)\n            x = BatchNormalization(name='batchnorm_1_{}_{}'.format(i, j))(x)\n            inputs.append(input_now)\n            nodes.append(x)\n\n        x = concatenate(nodes)\n        x = model_side(x)\n        x = BatchNormalization(name='batchnorm_2_{}'.format(i))(x)\n        sides.append(x)\n\n    x = concatenate(sides)\n    x = Dense(256, name='dense_0')(x)\n    x = BatchNormalization(name='batchnorm_3')(x)\n    x = LeakyReLU(alpha=0.05, name='leaky_relu_0')(x)\n    x = Dense(128, name='dense_1')(x)\n    x = BatchNormalization(name='batchnorm_4')(x)\n    x = LeakyReLU(alpha=0.05, name='leaky_relu_1')(x)\n    x = Dense(32, name='dense_2')(x)\n    x = BatchNormalization(name='batchnorm_5')(x)\n    x = LeakyReLU(alpha=0.05, name='leaky_relu_2')(x)\n    output = Dense(1, activation='linear', name='output')(x)\n    model = Model(inputs=inputs, outputs=output)\n    model.compile(optimizer='adam', loss='mean_absolute_error')\n    return model\n\ninput_dims = [ X_val.shape[1] for X_val in X_train ]\nmodel = get_model(input_dims)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"epochs = 250\nbatch_size = 2048\nval_split = 0.1\n\nearly_stopping = EarlyStopping(monitor='val_loss', min_delta=0.0001, patience=40, verbose=1, mode='auto',\n    restore_best_weights=True)\nreduce_on_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.1, patience=25, min_lr=1e-6, mode='auto',\n    verbose=1)\n\nhist = model.fit(X_train, Y_train, validation_split=val_split, batch_size=batch_size, epochs=epochs,\n    callbacks=[ early_stopping, reduce_on_lr ], verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndef plot_hist(hist, filename):\n    plt.plot(hist.history['loss'])\n    plt.plot(hist.history['val_loss'])\n    plt.title('Mean Absolute Error')\n    plt.ylabel('Loss')\n    plt.xlabel('Epochs')\n    plt.legend(['Train', 'Val'], loc='upper left')\n    plt.savefig(filename)\n\nplot_hist(hist, '/kaggle/working/plot.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_results = model.predict(X_test, batch_size=batch_size, verbose=1)[:,0].tolist()\n\nwith open('/kaggle/input/champs-scalar-coupling/sample_submission.csv', 'r') as f:\n    submit = pd.read_csv(f)\n    submit[\"scalar_coupling_constant\"] = submission_results\n    submit.to_csv(\"/kaggle/working/submission.csv\", index=False)\n\nsubmit.head()","execution_count":null,"outputs":[]}],"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}