{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"%%capture\n\nimport os\nfrom os.path import join\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom scipy.spatial.distance import euclidean\nfrom scipy import stats\nfrom tqdm import tqdm_notebook\nimport time\n\n\ntry:\n    import ase\nexcept:\n    !pip install ase\n    import ase\nfrom ase import Atoms\nimport ase.visualize","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = '../input'\n\nstructure_dir = join(DATA_DIR, 'structures')\nos.listdir(DATA_DIR)","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data quantity\n\nTotal size: ~1GB  \n(excluding the redundant 500MB from the structures folder (the same info is in structures.csv)","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv(join(DATA_DIR, 'train.csv'))\ntest = pd.read_csv(join(DATA_DIR, 'test.csv'))\nstructures = pd.read_csv(join(DATA_DIR, 'structures.csv'))\n\nprint(f'number of structures: {len(os.listdir(structure_dir))}')\nprint(f'training interactions: {len(train)}')\nprint(f'test interactions: {len(test)}')\nprint(f'fraction of interactions used for testing: {round(len(test) / (len(test) + len(train)), 2)}')\ntrain.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Structures","metadata":{}},{"cell_type":"code","source":"assert len(structures.molecule_name.unique()) == len(os.listdir(structure_dir))\nstructures.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"atom_count = structures.groupby('atom').size()\nplt.bar(atom_count.index, atom_count.values)\nplt.title('atom count');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3D visualization\n\nMost importantly, let's look at a few random molecules to get a rough idea about what we're dealing with.\n\nThe code for the 3D-visualization is from this cool kernel:  \nhttps://www.kaggle.com/borisdee/how-to-easy-visualization-of-molecules  \nMany thanks to the creator for sharing!","metadata":{}},{"cell_type":"code","source":"def view_molecule(structures=structures, name=None):\n    \"\"\"\n    Modified from:\n    https://www.kaggle.com/borisdee/how-to-easy-visualization-of-molecules\n    \"\"\"\n    if name is None:\n        name = np.random.choice(structures.molecule_name.unique())\n    print(f'molecule-name: {name}')\n    \n    molecule_df = structures.query('molecule_name == @name')\n    molecule    = Atoms(positions=molecule_df[['x', 'y', 'z']].values,\n                        symbols=molecule_df.atom.values)\n    \n    return ase.visualize.view(molecule, viewer=\"x3d\")\n\n\nview_molecule(name=None)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Number of atoms per molecule","metadata":{}},{"cell_type":"code","source":"atom_counts = structures.groupby('molecule_name').size()\nprint(f'minimum number of atoms per molecule: {atom_counts.min()}')\nprint(f'maximum number of atoms per molecule: {atom_counts.max()}')\natom_counts.hist(bins=atom_counts.max())\nplt.title('number of atoms per molecule');\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_structures = structures[structures.molecule_name.isin(set(train.molecule_name))]\nprint(f'number of traininig-structures: {len(train_structures.molecule_name.unique())}')\ntrain_counts = train_structures.groupby('molecule_name').size()\ntrain_counts.hist()\nplt.title('training: number of atoms per molecule')\nplt.show()\n\ntest_structures = structures[structures.molecule_name.isin(set(test.molecule_name))]\nprint(f'number of test-structures: {len(test_structures.molecule_name.unique())}')\ntest_structures.groupby('molecule_name').size().hist()\nplt.title('test: number of atoms per molecule');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Number of labelled interactions per molecule","metadata":{}},{"cell_type":"code","source":"train_j_counts = train.groupby('molecule_name').size()\nprint(f'min and max number of j-couplings: {min(train_j_counts.values)}, {max(train_j_counts.values)}')\ntrain_j_counts.hist()\nplt.title('training-set: number of labelled interactions per molecule')\nplt.show()\n\ntest.groupby('molecule_name').size().hist()\nplt.title('test-set: number of labelled interactions per molecule');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ...per coupling type","metadata":{}},{"cell_type":"code","source":"for sc_type, type_df in train.groupby('type'):\n    counts = type_df.groupby('molecule_name').size()\n    counts.hist()\n    plt.title(f'{sc_type}: number of labelled interactions')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Target: scalar coupling constant\n\nA measure of magnetic interaction between two atoms. Determined by the 3D environment of the atoms and measured with NMR.\n\nThe relationshipb between 3D structure and SCC is not trivial and requires expensive quantum mechanical calculations. Let's see if the relationship can be learned from data with a (graph) neural net.\n","metadata":{}},{"cell_type":"code","source":"sc_types = train.type.unique()\nprint(f'{len(sc_types)} scalar coupling types:\\n{sorted(sc_types)}')\nprint(f'\\nbonds-apart: {set(int(x[0]) for x in sc_types)}')\nprint(f'\\nsc-interactions: {set(x[2:4] for x in sc_types)}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Number of interactions per scalar-coupling type","metadata":{}},{"cell_type":"code","source":"type_count = train.groupby('type').size() / len(train)\nplt.bar(type_count.index, type_count.values)\nplt.title('training: number of interactions per type')\nplt.ylabel('fraction of interactions')\nplt.xlabel('scalar coupling type')\nplt.show()\n\ntype_count = test.groupby('type').size() / len(test)\nplt.bar(type_count.index, type_count.values)\nplt.ylabel('fraction of interactions')\nplt.xlabel('scalar coupling type')\nplt.title('test: number of interactions per type');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_sc_type_features(df):\n    df['bonds']    = df.type.map(lambda x: int(x[0]))\n    df['atom_pair'] = df.type.map(lambda x: x[2:4])\n    return df\n\n\ntrain = add_sc_type_features(train)\ntest  = add_sc_type_features(test)\n\n\nbonds_counts = train.bonds.value_counts()\nplt.bar(bonds_counts.index, bonds_counts.values)\nplt.xlabel('bonds apart')\nplt.ylabel('number of interactions')\nplt.title('train: bonds between sc-atoms')\nplt.show()\n\norder_counts = test.bonds.value_counts()\nplt.bar(order_counts.index, order_counts.values)\nplt.xlabel('bonds apart')\nplt.ylabel('number of interactions')\nplt.title('test: bonds between sc-atoms')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"atom_pair_types = train.atom_pair.value_counts()\nplt.bar(atom_pair_types.index, atom_pair_types.values)\nplt.xlabel('atom-pair')\nplt.ylabel('number of interactions')\nplt.title('train: sc atom-pairs')\nplt.show()\n\natom_pair_types = train.atom_pair.value_counts()\nplt.bar(atom_pair_types.index, atom_pair_types.values)\nplt.xlabel('atom-pair')\nplt.ylabel('number of interactions')\nplt.title('test: bonds sc atom-pairs')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Critical information:\n\nThe 8 scalar coupling types in this dataset have very different ranges of sc-constant!\n=> The most straightforward guess would be to predict the median for each type. This would already be much better than just predicting the global median.\n\nAlso, that's the reason, why training a separate model for each type can pay off in this competition.","metadata":{}},{"cell_type":"code","source":"sns.boxplot(data=train, x='type', y='scalar_coupling_constant', fliersize=0.5)\nplt.show();\n\nsns.boxplot(data=train, x='bonds', y='scalar_coupling_constant', fliersize=0.5)\nplt.show();\n\nsns.boxplot(data=train, x='atom_pair', y='scalar_coupling_constant', fliersize=0.5)\nplt.show();","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.scalar_coupling_constant.hist()\nplt.title('(training) target-variable distribution');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training-set Features\n\nThese features are provided only for the training set. I had a close look at them but ended up not using any of them for this competition.","metadata":{}},{"cell_type":"markdown","source":"### Molecule (graph-level) features","metadata":{}},{"cell_type":"markdown","source":"molecule dipole moment: https://en.wikipedia.org/wiki/Dipole#Molecular_dipoles","metadata":{}},{"cell_type":"code","source":"dipole = pd.read_csv(join(DATA_DIR, 'dipole_moments.csv'))\n\ndisplay(dipole.head())\nplt.hist(dipole[['X', 'Y', 'Z']].values.ravel())\nplt.title('diple moments distribution');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"potential energy in molecules: https://www.quora.com/Where-is-potential-energy-stored-in-molecules","metadata":{}},{"cell_type":"code","source":"pot_energy = pd.read_csv(join(DATA_DIR, 'potential_energy.csv'))\n\ndisplay(pot_energy.head())\npot_energy.potential_energy.hist()\nplt.title('potential energy distribution');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Atom (node-leve) features","metadata":{}},{"cell_type":"markdown","source":"magnetic shield tensors: basically structural information obtaind from NMR  \nHence, the information should be more readily available from the structures directly.","metadata":{}},{"cell_type":"code","source":"mag_shield = pd.read_csv(join(DATA_DIR, 'magnetic_shielding_tensors.csv'))\n\ndisplay(mag_shield.head())\nplt.hist(mag_shield.drop(['molecule_name', 'atom_index'], axis=1).values.ravel())\nplt.title('magnetic shield distribution');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Mulliken charge: https://en.wikipedia.org/wiki/Mulliken_population_analysis  \n... basically parital charges of atoms. But dependent on the calculation method.","metadata":{}},{"cell_type":"code","source":"mulliken = pd.read_csv(join(DATA_DIR, 'mulliken_charges.csv'))\n\ndisplay(mulliken.head())\nmulliken.mulliken_charge.hist()\nplt.title('mulliken charge distribution');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Interaction (edge-level) features","metadata":{}},{"cell_type":"code","source":"sc_contributions = pd.read_csv(join(DATA_DIR, 'scalar_coupling_contributions.csv'))\ndisplay(sc_contributions.head())\n\nfor contribution in ('fc', 'sd', 'pso', 'dso'):\n    sc_contributions[contribution].hist()\n    plt.title(f'scalar coupling contribution: {contribution}')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"check that scalar coupling contributions add up to target variable:","metadata":{}},{"cell_type":"code","source":"sc = pd.merge(sc_contributions, train, on=('molecule_name', 'atom_index_0', 'atom_index_1', 'type'))\nassert len(sc) == len(train) == len(sc_contributions)\n\nsc['contribution_sum'] = sc.fc + sc.sd + sc.pso + sc.dso\nsc['difference'] = (sc.scalar_coupling_constant - sc.contribution_sum).map(abs)\nprint('maximal differnce between sum of scalar coupling contributions and\\nthe target variable '\n      f'scalar coupling: {sc.difference.max()}')\nsc.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distance distributions\n\nCode from these 3 kernel is used for distance computations:  \nhttps://www.kaggle.com/cpmpml/ultra-fast-distance-matrix-computation  \nhttps://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark  \nhttps://www.kaggle.com/artgor/brute-force-feature-engineering","metadata":{}},{"cell_type":"code","source":"#Code in this cell is from this kernel:\n#https://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark\n\ndef map_atom_info(df, atom_idx, structures_df=structures):\n    df = pd.merge(df, structures_df, 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\ntrain = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 1)\n\ntrain.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# code in this cell is from:\n# https://www.kaggle.com/artgor/brute-force-feature-engineering\n\ntrain_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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Distances between atom-pairs with scalar coupling (target)","metadata":{}},{"cell_type":"code","source":"train.dist.hist()\nplt.title('train: distances of target sc-interactions')\nplt.xlabel('distance [Å]')\nplt.ylabel('number of interactions (atom-pairs)');\nplt.show()\n\ntest.dist.hist()\nplt.title('test: distances of target sc-interactions')\nplt.xlabel('distance [Å]')\nplt.ylabel('number of interactions (atom-pairs)')\nplt.show()\n\nprint(f'maximal sc-distance: {max(*train.dist, *test.dist):.2f}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('training-set:')\nfor n_bonds, df in train.groupby('bonds'):\n    df.dist.hist()\n    plt.title(f'train: distances of {n_bonds}-bond sc')\n    plt.xlabel('distance [Å]')\n    plt.ylabel('number of interactions (atom-pairs)');\n    plt.show()\n    \nprint('test-set:')\nfor n_bonds, df in test.groupby('bonds'):\n    df.dist.hist()\n    plt.title(f'test: distances of {n_bonds}-bond sc')\n    plt.xlabel('distance [Å]')\n    plt.ylabel('number of interactions (atom-pairs)');\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see a clear separation of one-bond j-coupled atoms and two-bond j-coupled atoms. This means that from euclidean distance alone, we could tell, whether the coupled atoms are one or two bonds apart.\n\nBetween 2-bond and 3-bond j-coupled atoms, there is a littel overlap. I.e. in the vast majority, the euclidean distance indicates whether two coupled atoms are two or three bonds apart. In a small number of cases, 2-bond and 3-bond coupled atoms can have similar euclidean distances. This is owed mostly to the different bond-angles resulting in different euclidean distnaces. Different bond lengths of different bond types (C-H vs. C-C, vs. N-H, etc.) also effect the euclidean distance. However, because bonds in organic molecules do not have vastly different lengths (~1-1.5 Angström), this effect is probably smaller than the effect of bond-angles.\n\nAll in all, euclidian distance pretty much includes the information of the number of bonds between two atoms (only for 1-3 bonds, of course) and vice versa.","metadata":{}},{"cell_type":"markdown","source":"### Maximum atom-atom distance within molecules\n#### how much space do the molecules take up?\n\nThis information is only relevant if a 3D-conv-net trained directly on the 3D structure would be used.","metadata":{}},{"cell_type":"code","source":"def get_dist_matrix(df_structures_idx, molecule, bug_fixed=True):\n    \"\"\"\n    Function from:\n    https://www.kaggle.com/cpmpml/ultra-fast-distance-matrix-computation\n    \"\"\"\n    df_temp = df_structures_idx.loc[molecule]\n    locs = df_temp[['x','y','z']].values\n    num_atoms = len(locs)\n    loc_tile = np.tile(locs.T, (num_atoms,1,1))\n    \n    # bug fix: first: sum up the squares, THEN take the square root\n    if bug_fixed:\n        return np.sqrt( ((loc_tile - loc_tile.T)**2).sum(axis=1) )\n    else:\n        return np.sqrt( (loc_tile - loc_tile.T)**2 ).sum(axis=1)\n\n\nmol_name2dist_matrix = {}\nstructures.index = structures.molecule_name\n\nfor mol_name in tqdm_notebook(structures.molecule_name.unique()):\n    mol_name2dist_matrix[mol_name] = get_dist_matrix(structures, mol_name)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max_dists = [matrix.max() for matrix in mol_name2dist_matrix.values()]\nplt.hist(max_dists)\nplt.title('Maximum atom-atom distance within molecule')\nplt.ylabel('number of molecules')\nplt.xlabel('distance [Å]')\nplt.show()\n\nmin_dist = min(matrix[matrix != 0].min() for matrix in mol_name2dist_matrix.values())\nprint(f'minimum atom-to-atom position (over all molecules): {min_dist:.4f}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Check distances\n\nJust double- and triple- checking the distance computations... :-D","metadata":{}},{"cell_type":"code","source":"# distance unit-test:\nstructures_idx = structures.set_index('molecule_name')\nmolname = 'dsgdb9nsd_000001'\n\nmdf = structures_idx.loc[molname]\n\n#dist_matrix_bug   = get_dist_matrix(structures_idx, molname, bug_fixed=False)\ndist_matrix_fixed = get_dist_matrix(structures_idx, molname)\n\na = mdf.iloc[0][['x', 'y', 'z']].values\nb = mdf.iloc[1][['x', 'y', 'z']].values\ndisplay(mdf.head(2))\n\nprint('scipy-euclidian: ', euclidean(a, b))\nprint('numpy-norm: ', np.linalg.norm(a - b))\nprint('generic: ', np.sqrt(((a - b)**2).sum()))\n#print(\"CPMP's calculation (bug): \", dist_matrix_bug[0, 1])\nprint(\"CPMP's calculation (fixed): \", dist_matrix_fixed[0, 1])","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}