{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Synopsis\n\nThe number of bonds between two atoms is an interresting feature to play with but is not provided with the competition data.\n\nThis kernel builds a dataset with the reconstructed bonding scheme for each molecule:  \n* The algorithm used to do so is to build a spanning tree conecting each atom of the molecule. \n* Connection are made based on the expecting bonding distance for each pair of nearby atoms.  \n* A partial ionized group detection is performed.  \n* Shortest cycles are detected.\n\nNote: There is another alternative to build bonding data, which I've not tested yet, see post here:  \n[Helpful Information: How to obtain bonding data](Helpful Information: How to obtain bonding data)\n\n## Changelog\n\n* V4 : manual inspection of one molecule failing bond matching : dsgdb9nsd_000271 = alanine isomere\n* V5 : added Ase 3D mouse visualization for alanine \"bug\"\n* V6 : added support for ionized molecules and `4N` bonds after discovering the [zwitterionic](https://en.wikipedia.org/wiki/Zwitterion) form of alanine in `dsgdb9nsd_000202`. Changed type from `1CH` to `1.0CH`\n* V9 : added cycle detection code with statistics and examples\n* V11: save cycles in .csv output files\n\n## Table `train_bonds` and `test_bonds`\nThe [Output Dataset](https://www.kaggle.com/asauve/predicting-molecular-properties-bonds) will be updated following bond detction improvements. \nIt is built with the same keys as for the competition tables `structure` and `train`/`test`.\n\nTables `train_bonds` and `test_bonds` have one row per bond :\n```\nmolecule_name | atom_index_0 | atom_index_1 | nbond | error | bond_type\n```\n* `nbond` is the number of covalent bonds between `atom_index_0` (relative to `structure`) and `atom_index_1`. It can take values `1.5` for `COO-` groups.\n* `error` is `0` if the molecule has a  consistent reconstruction (meaning that each atom has the expected number ofcovalent bonds). \nOtherwise it is set to `1` for each bond of the molecule\n* `bond_type` is <nbond><atom0><atom1>. Atom appear in lexicographic order to ensure unicity for 2 atoms combination.   \n    Example: `1.0CH`, `2.0CC`,  `1.5CO`\n\n## Table `train_charges` and `test_charges`\nIt has been found that some molecules are ionized after discovering the zwiterionic form of *beta-alanine* in molecule `dsgdb9nsd_000202`.  \nHence support for atom charge has been included in tables `train_charges` and `test_charges`, plus an optional 4th bond for `N`. \nThese tables have the very same key as `structure` and one row is:\n```\nmolecule_name | atom_index | charge\n```\n* charge can be `0`, `-1`, `+0.5` (for `O` in carboxyle) or `+1` (may evolve if needed)\n\n## Table `train_cycles` and `test_cycles`\n\nCycles are detected and stored in a separated table because one atom can be a member of several cycles.  \nThese tables have the very same key as `structure` and one row is:\n```\nmolecule_name | cycle_index | cycle_seq | atom_index\n```\n* `cycle_index` starts from zero and allow to select one cycle\n* `cycle_seq` is the index of atom in the cycle sequence.\nIt is guaranteed that two consecutive atoms in *cycle sequence* have a covalent bond.\n* `atom_index` is the atom idex from table `structure` which can be member of several cycle.   \n\n# General information\n\n![](https://storage.googleapis.com/kaggle-competitions/kaggle/14313/logos/thumb76_76.png?t=2019-05-16-16-56-19)\n\nThis kernel uses data from [Predicting Molecular Properties](https://www.kaggle.com/c/champs-scalar-coupling) which is intended to predict interactions between atoms in the domain of Nuclar Magnetic Resonnance (NMR). More precisely it is the *scalar coupling constant* between atoms which is to be predicted.\n\nAs this challenge is based uppon molecular topological properties, it can be useful to have an appropriate way of representing molecules.\nThese visualizations can then be used to infer useful hints to understand the coupling properties, engineer appropriate features and debug\nprediction failures.\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Others Work\n\n* [Is type related to scalar coupling? EDA](https://www.kaggle.com/hrmello/is-type-related-to-scalar-coupling)  \n\"ALL 1JHC couplings are above 60.\"\n* [Competition metric](https://www.kaggle.com/abhishek/competition-metric)\n* [Molecular Properties EDA and models](https://www.kaggle.com/artgor/molecular-properties-eda-and-models)   \nthis kernel scores 0.626 and uses features `dist_speedup` which is atoms relative distance calculated with L2 norm\n* [Exploring Molecular Properties Data](https://www.kaggle.com/robikscube/exploring-molecular-properties-data)   \npropose a nice exploratory work on extraneous features in the train dataset\n* [How To: Easy Visualization of Molecules.](https://www.kaggle.com/borisdee/how-to-easy-visualization-of-molecules/comments)  \nBoris shows in  this kernel hos to use the Ase library for mouse movable molecular display (no bond type support)"},{"metadata":{},"cell_type":"markdown","source":"# Load data"},{"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)\nfrom ase import Atoms  # 3D mouse clickable funky visualization\nimport ase.visualize\nimport networkx as nx # for cycle detection\n\nfrom mpl_toolkits.mplot3d import Axes3D\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set()\n\nimport os\nprint(\"Modules loading Done!\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"csv_files = [ \"train.csv\", \"test.csv\", \"structures.csv\" ]\ncsv_vars  = [ filename[:-4] for filename in csv_files ]\nfor filename, var in zip( csv_files, csv_vars ):\n    eval_str = f\"{var:32s} = pd.read_csv('../input/{filename}')\"\n    print(eval_str)\n    exec(eval_str)\n    print(f\"{'nb of cols ':32s}= \"+str(eval(f\"len({var})\")))\n    exec(f\"display({var}.head())\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Molecule display code"},{"metadata":{},"cell_type":"markdown","source":"The following code snipet (unhide to see) is taken from [Training set molecule visualization](https://www.kaggle.com/asauve/training-set-molecule-visualization). \nIt builds a 3D molecule visualization with \n* scalar_coupling_constant as transparent connections: thickness is a factor of scalar value\n  * plain for 1J spin-spin coupling\n  * dashed for 2J\n  * dotted for 3J\n* covalent bonds are drawn as black lines: thickness increase with number of bonds"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def view3d_molecule(name):\n    m = structures[structures.molecule_name == name]\n    positions = m[['x','y','z']].values\n    # pip install ase\n    return ase.visualize.view(Atoms(positions=positions, symbols=m.atom.values), \n                              viewer=\"x3d\") \ncpk = { \n    'C': (\"black\", 2),\n    'H': (\"white\", 1),\n    'O': (\"red\",   2),\n    'N': (\"blue\", 2),\n    'F': (\"green\", 2) }\n\nbond_colors = {'1.0':'black', '1.5':'darkgreen', '2.0':'green', '3.0':'red'}\n\ndef bond_type_to_pair(bond_type):\n    return bond_type[3:]\ndef bond_type_to_n(bond_type):\n    return bond_type[0:3]\n\ndef plot_molecule(name, ax=None, bonds=None, charges=None, elev=0, azim=-60):\n    \"\"\"bonds = if provided add bonds display from the bond table dataset in https://www.kaggle.com/asauve/predicting-molecular-properties-bonds\n    elev = 3D elevation angle [degree] for the molecule view\n    azim = 3D azimut angle [degree]\n    \"\"\"\n    if not ax:\n        fig = plt.figure()\n        ax = fig.add_subplot(111, projection='3d')\n    if (elev != 0) or (azim != -60):\n        ax.view_init(elev=elev, azim=azim)\n    \n    # atoms location\n    m = structures[structures.molecule_name == name].sort_values(by='atom_index')\n    if (charges is not None):\n        charges = charges[charges.molecule_name == name].sort_values(by='atom_index')\n        if len(charges) != len(m):\n            print(f\"Warning bad charges data for molecule {name}\")\n    \n    # formula\n    acount = {a : 0 for a in cpk}\n    for a in m.atom:\n        acount[a] += 1\n    formula = \"\"\n    for a in acount:\n        if acount[a] == 1:\n            formula += a\n        elif acount[a] > 1:\n            formula += \"%s_{%d}\" % (a, acount[a])\n\n    ax.set_title(f'{name} ${formula}$')\n    \n    # display couplings (coupling is not molecular bonds!)\n    couples = train[train.molecule_name == name][['atom_index_0', 'atom_index_1', 'type', 'scalar_coupling_constant']]\n    for c in couples.itertuples():\n        m0 = m[m.atom_index == c.atom_index_0]\n        m1 = m[m.atom_index == c.atom_index_1]\n        ax.plot([float(m0.x), float(m1.x)],[float(m0.y), float(m1.y)],[float(m0.z), float(m1.z)],\n               linestyle = ['', '-', '--', 'dotted'][int(c.type[0])],\n               color     = ['', 'black', 'green', 'red' ][int(c.type[0])],\n               linewidth = abs(float(c.scalar_coupling_constant))/5,\n               alpha=0.2), \n    if bonds is not None:\n        for b in bonds[bonds.molecule_name == name].itertuples():\n            m0 = m[m.atom_index == b.atom_index_0]\n            m1 = m[m.atom_index == b.atom_index_1]\n            ax.plot([float(m0.x), float(m1.x)],[float(m0.y), float(m1.y)],[float(m0.z), float(m1.z)], 'black', \n                    linewidth=2*float(b.nbond),\n                    color=bond_colors[bond_type_to_n(b.bond_type)])\n            \n    # display atoms\n    ax.scatter(m.x, m.y, m.z, c=[cpk[a][0] for a in m.atom], s=[100*cpk[a][1] for a in m.atom], edgecolor='black')\n        \n    # display atom index and charges\n    for row in m.itertuples():\n        x = float(row.x) - 0.15 if row.x > ax.get_xlim()[0] + 0.15 else row.x\n        y = float(row.y) - 0.15 if row.y > ax.get_ylim()[0] + 0.15 else row.y\n        z = float(row.z) - 0.15 if row.z > ax.get_zlim()[0] + 0.15 else row.z\n        ax.text(x, y, z, str(row.atom_index), color='darkviolet')\n        if charges is not None:\n            ch = float(charges[charges.atom_index == row.atom_index].charge)\n            if ch != 0:\n                x = float(row.x) + 0.15 if row.x < ax.get_xlim()[1] - 0.15 else row.x\n                y = float(row.y) + 0.15 if row.y > ax.get_ylim()[1] - 0.15 else row.y\n                z = float(row.z) + 0.15 if row.z < ax.get_zlim()[1] - 0.15 else row.z\n                ax.text(x, y, z, f\"{ch:+.1f}\", color='orangered' if ch > 0 else 'blue',\n                        bbox=dict(boxstyle='round', facecolor='white', alpha=0.5, \n                                  edgecolor='black'))\n#plot_molecule(\"dsgdb9nsd_000140\", bonds=train_bonds, charges=train_charges)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This example shows the three kind of coupling types. The 1 bond coupling having the largest coupling constant."},{"metadata":{},"cell_type":"markdown","source":"# Bonds reconstruction algorithm\n\nUseful links :\n* <a href=\"https://en.wikipedia.org/wiki/Atomic_radii_of_the_elements_(data_page)\">Atomic radii of the elements</a>\n\nThe algorithm for bond selection is\n* Select a molecule\n* **STEP 1** : 1-bond connect   \n For each atom\n    * For each #VALENCE closest match:\n        * does the distance correspond to the expected bond distance?\n            * does both atoms have remaining bonds?\n              * decease remaining valence and create bond\n              * mark both atoms as leaf if one of them has zero remaining valence\n* **STEP 2** : n-bond greedy connect   \nWhile there are remaining leafs atoms with available valence\n    * For each connected neighbour \n      * add as many bonds as possible\n      * mark both connected atoms as leaf\n* **STEP 3** : ionized radical search\nIf there are remaining valence on some atoms\n    * Look for disconnected NH3  => convert to connected NH3+\n    * Look for O with one available bond connected to CO => convert to COO-\n    * Other cases may be handled later\n* **STEP 4** : cycle search\n    * Use `networkx` pip library for shortest cycle detection\n    * Index atoms by order in their respective cycle\n* **CHECK**   \nHas any atom remaining valence > 0?\n    * mark the molecule with error=1\n  "},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.spatial import cKDTree as KDTree\nfrom tqdm import tqdm_notebook as tqdm\n\n\n# \"usual\" valence of atoms\nVALENCE_MAX = { 'C': 4, 'H' : 1, 'N' : 4, 'O' : 2, 'F': 1 }\nVALENCE_STD = { 'C': 4, 'H' : 1, 'N' : 3, 'O' : 2, 'F': 1 }\n\n# expected distances in [A] for covalence 1 bond\nBOND_DIST_C1 = { 'C' : 0.77, 'H' : 0.38, 'N' : 0.75, 'O' : 0.73, 'F' : 0.71 }\n\n# order used for finding bonds by atom type\nBOND_ORDER = { 'H' : 0, 'F' : 0, 'O' : 1, 'N' : 2, 'C' : 3}\n\ndef add_bond(n_avail, nbond, a0, a1, d1=None):\n    key = tuple(sorted((a0, a1)))\n    if key in nbond:\n        nbond[key][0] += 1.0\n    elif d1 is not None:\n        nbond[key] = [1.0, d1]\n    else:\n        raise Exception(f\"{a0},{a1} added after phase 1\")\n    n_avail[a0] -= 1\n    n_avail[a1] -= 1\n    #display(f\"{a0},{a1} => {n_avail}, {nbond}\")\n\ndef get_bonded_atoms(atoms, nbond, i):\n    \"\"\"returns: [sorted atoms list], [sorted atom index] )\"\"\"\n    bonded = []\n    for (a0, a1), (n, _) in nbond.items():\n        if a0 == i:\n            bonded.append((a1, atoms[a1]))\n        elif a1 == i:\n            bonded.append((a0, atoms[a0]))\n    bonded = sorted(bonded, key=lambda b: b[1])\n    return \"\".join([b[1] for b in bonded]), [b[0] for b in bonded]\n\n\ndef search_bonds(kdt, n_avail, nbond, connected, isleaf, coords, atoms, atoms_idx, a0, \n                 connect_once=True, VALENCE=VALENCE_STD):\n    atom0 = atoms[a0]\n    if n_avail[a0] == 0:\n        return\n\n    # select closest atoms ORDERED BY DISTANCE: closest first\n    # note: the first answer is the atom itself and must be removed\n    next_dist, next_i = kdt.query(coords[a0], min(1+VALENCE[atom0], len(atoms)))\n    next_dist = next_dist[1:]; # remove a0 from list\n    next_i = next_i[1:]\n\n    # for each #VALENCE closest atoms\n    found = False\n    for d1,a1 in zip(next_dist, next_i):\n        if connect_once and (a1 in connected[a0]):\n            continue  # enforce 1-bond only in STEP 1\n        atom1 = atoms[a1]\n        predicted_bond = BOND_DIST_C1[atom0] + BOND_DIST_C1[atom1]\n        if abs(d1 / predicted_bond) < 1.2: # keep only atoms in the 20% expected distance or closer\n            if n_avail[a1] > 0:\n                add_bond(n_avail, nbond, a0, a1, d1)\n                connected[a0][a1] = 1\n                connected[a1][a0] = 1\n                if (n_avail[a0] == 0) or (n_avail[a1] == 0):\n                    isleaf[a0] = 1\n                    isleaf[a1] = 1\n                found = True\n                #print(\"leaF/Trunk & avail: \"+ \", \".join([f\"{i}:{atoms[i]}={leaflabel[isleaf[i]]}{n_avail[i]}\"\n                #                 for i in ordered_atoms_index]))\n\n        else:\n            #print(f\"-- match failure in molecule_name={name} {a0}:{atom0}-{a1}:{atoms[a1]}={d1} predicted={predicted_bond}\")\n            pass\n    return found\n    \ndef compute_bonds(structures, molecules):\n    out_name = []; out_a0 = []; out_a1 = []; out_n = []; out_dist = []; out_error = []; out_type = []\n    cycle_name = []; cycle_index = []; cycle_seq = [];  cycle_atom_index = []\n    charge_name = []; charge_atom_index = []; charge_value = [];\n    for imol,name in tqdm(list(enumerate(molecules))):\n        molecule = structures.loc[name ]\n        error    = 0\n        atoms    = molecule.atom.values\n        atoms_idx= molecule.atom_index.values\n        \n        n_avail  = np.asarray([VALENCE_STD[a] for a in atoms])\n        n_charge = np.zeros(len(atoms), dtype=np.float16)\n        isleaf   = np.zeros(len(atoms), dtype=np.bool) # is the atom in the leafs of connection tree?\n        coords   = molecule[['x', 'y', 'z']].values\n        kdt = KDTree(coords) # use an optimized structure for closest match query\n        nbond = {}\n        connected = { i:{} for i in atoms_idx}\n        \n        # select Hydrogen first to avoid butadyne-like ordering failures (molecule_name=dsgdb9nsd_000023)\n        ordered_atoms_index = list(atoms_idx)\n        ordered_atoms_index.sort(key = lambda i: BOND_ORDER[atoms[i]])\n        ordered_atoms_index = np.asarray(ordered_atoms_index)\n        \n        # STEP 1: 1-bond connect each atom with closest match\n        #         only one bond for each atom pair is done in step 1\n        for a0 in ordered_atoms_index:\n            search_bonds(kdt, n_avail, nbond, connected, isleaf, coords, atoms, atoms_idx, a0,\n                         connect_once=True, VALENCE=VALENCE_STD)\n     \n        # STEP 2: greedy connect n-bonds, progressing from leafs of connection tree\n        while (((n_avail > 0).sum() > 0) and isleaf).sum() > 0:\n            progress = False\n            for a0 in ordered_atoms_index:\n                #print(\"leaF/Trunk & avail: \" + \", \".join([f\"{i}:{atoms[i]}={leaflabel[leaf[i]]}{n_avail[i]}\"\n                #                                          for i in ordered_atoms_index]))\n                if (n_avail[a0] > 0) and isleaf[a0]:\n                    for a1 in connected[a0]:\n                        if (n_avail[a0] > 0) and (n_avail[a1] > 0):\n                            add_bond(n_avail, nbond, a0, a1) \n                            progress = True\n                            if (n_avail[a0] == 0) or (n_avail[a1] == 0):\n                                isleaf[a0] = 1\n                                isleaf[a1] = 1\n            if not progress:\n                break\n        \n        # gather remaining multiple bonds\n        if n_avail.sum() > 0:\n            for key in nbond.keys():\n                a0, a1 = key\n                while (n_avail[a0] > 0) and (n_avail[a1] > 0):\n                    add_bond(n_avail, nbond, a0, a1)\n    \n        # STEP 3: search for known ionized radicals\n        if n_avail.sum() > 0:\n            for (i,a) in zip(atoms_idx, atoms):\n                if a == 'N':\n                    # NH3+\n                    bonded_str, bonded_idx = get_bonded_atoms(atoms, nbond, i)\n                    if (bonded_str == \"HHH\") and (n_avail[i] == 0):\n                        # add a valence unit and search a dangling bond nearby\n                        n_avail[i] += 1; n_charge[i] += 1\n                        if search_bonds(kdt, n_avail, nbond, connected, isleaf, coords, atoms, atoms_idx,\n                                        i, connect_once=False, VALENCE=VALENCE_MAX):\n                            print(f\"++ NH3+ found for {name} atom_index={i}\")\n                        else:\n                            print(f\"** NH3+ bonding failure for {name} atom_index={i}\")\n                            \n                elif (a == 'O') and (n_avail[i] == 1):\n                    # COO-\n                    bonded_str, bonded_idx = get_bonded_atoms(atoms, nbond, i)\n                    if bonded_str == \"C\":\n                        C_i = bonded_idx[0]\n                        C_bonded_str, C_bonded_idx = get_bonded_atoms(atoms, nbond, C_i)\n                        if (\"OO\" in C_bonded_str):\n                            has_2CO = False \n                            #print (C_bonded_str, C_bonded_idx, nbond, name)\n                            for a1, i1 in zip(C_bonded_str, C_bonded_idx):\n                                key = tuple(sorted((C_i, i1)))\n                                if (a1 == 'O') and  (nbond[key][0] == 2):\n                                    has_2CO = True\n                            if (len(C_bonded_idx) == 3) and has_2CO:\n                                # found carboxyle!\n                                n_avail[i] -= 1\n                                print(f\"**  COO- found for {name} C_atom_index={C_i}\")\n                                for a1, i1 in zip(C_bonded_str, C_bonded_idx):\n                                    if a1 == 'O':\n                                        n_charge[i1] = -0.5\n                                        key = tuple(sorted((C_i, i1)))\n                                        nbond[key][0] = 1.5                    \n                                    \n        # detect cycles : algo complexity in O(m^2 * n)\n        #    paper : https://link.springer.com/article/10.1007/s00453-007-9064-z\n        #    nx doc: https://networkx.github.io/documentation/latest/reference/algorithms/generated/networkx.algorithms.cycles.minimum_cycle_basis.html\n        graph = nx.Graph([bond for bond in nbond.keys()])\n        unordered_cycles = nx.minimum_cycle_basis(graph)\n        \n        # index atoms by their sequential order in the cycle: i.e follow bonds\n        # Note: this code can be written in a much cleaner way!\n        if len(unordered_cycles) > 0:\n            for icycle, c in enumerate(unordered_cycles):\n                available = {i:1 for i in c}\n                a0 = c[0]\n                cycle = [ a0 ]\n                del(available[a0])\n                for index in range(1, len(c)):\n                    # get atoms bonded to a0\n                    bonded = [ b for b in nbond.keys() if a0 in b ]\n                    bonded = list(map(lambda b: b[0] if b[1] == a0 else b[1], bonded))\n                    \n                    # get next atom and remove it from cycle\n                    assert(len(bonded) > 0)\n                    found = False\n                    for a1 in bonded:\n                        if (a1 in bonded) and (a1 in available):\n                            cycle.append(a1)\n                            del(available[a1])\n                            a0 = a1\n                            found = True\n                            break\n                    assert(found)\n                    \n                # and add cycles found to the cycle dataframe lists\n                cycle_name.extend([name] * len(cycle))\n                cycle_index.extend([icycle] * len(cycle))\n                cycle_seq.extend(np.arange(len(cycle)))\n                cycle_atom_index.extend(cycle)\n        \n        # display info on failed molecules\n        if n_avail.sum() > 0:\n            error = 1\n            print(f\"   Remaining bondings={n_avail.sum()} for molecule_name={name}, atoms: \" +\n                  \", \".join([f\"{i}:{atoms[i]}\" for i in atoms_idx if n_avail[i] > 0]))\n\n        # inputs for DataFrame bonds\n        for (a0, a1), (n, dist) in nbond.items():\n            # append to python lists which is 7x faster than toa pd.DataFrame\n            out_name.append(name); out_a0.append(a0); out_a1.append(a1); out_n.append(n)\n            out_dist.append(dist); out_error.append(error)\n            out_type.append(f\"{n:0.1f}\" + \"\".join(sorted(f\"{atoms[a0]}{atoms[a1]}\")) )\n    \n        # inputs for DataFrame charges\n        charge_name.extend(        [name] * len(atoms) )\n        charge_atom_index.extend(  molecule.atom_index.values )\n        charge_value.extend(       n_charge )\n        \n    bonds = pd.DataFrame({'molecule_name':out_name, 'atom_index_0':out_a0,'atom_index_1':out_a1, 'nbond':out_n, \n                          'L2dist':out_dist, 'error':out_error, 'bond_type':out_type})\n    charges = pd.DataFrame({'molecule_name':charge_name, 'atom_index':charge_atom_index, \n                            'charge': charge_value})\n    cycles = pd.DataFrame({'molecule_name' : cycle_name, 'cycle_index' : cycle_index, \n                           'cycle_seq' : cycle_seq, 'atom_index' : cycle_atom_index})\n    return bonds, charges, cycles\n            \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build bonds for train / test sets"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"train_bonds, train_charges , train_cycles = \\\n    compute_bonds(structures.set_index('molecule_name'), \n                  train.molecule_name.unique())\n#debug_name = \"dsgdb9nsd_000202\"\n#train_bonds, train_charges = compute_bonds(structures.set_index('molecule_name'), \n#                                           [ debug_name ])\n#plot_molecule(debug_name, bonds=train_bonds, charges=train_charges)\n#view3d_molecule(debug_name)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"test_bonds, test_charges, test_cycles = \\\n    compute_bonds(structures.set_index('molecule_name'), \n                  test.molecule_name.unique())\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"Failed molecules ratio in test set = {test_bonds.groupby('molecule_name')['error'].max().sum() / len(test_bonds.molecule_name.unique()):0.3f}\")\nprint(f\"Failed molecules ratio in train set = {train_bonds.groupby('molecule_name')['error'].max().sum() / len(train_bonds.molecule_name.unique()):0.3f}\")\n\ntest_nionized = test_charges.groupby('molecule_name')['charge'].apply(lambda arr: (arr.values != 0).max()).sum()\nprint(f\"Ionized molecules ratio in test set = {test_nionized/len(test_bonds.molecule_name.unique()):0.3f}\")\n\ntrain_nionized = train_charges.groupby('molecule_name')['charge'].apply(lambda arr: (arr.values != 0).max()).sum()\nprint(f\"Ionized molecules ratio in train set = {train_nionized/len(train_bonds.molecule_name.unique()):0.3f}\")\n\nprint(f\"{len(test_bonds.molecule_name.unique())} unique molecules in test set\")\nprint(f\"{len(train_bonds.molecule_name.unique())} unique molecules in train set\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check bond type consistency\n\n## Distribution of bond length"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(15,15))\nbond_types = train_bonds.bond_type.unique()\nunique_pairs = sorted({bond_type_to_pair(b):1 for b in bond_types}.keys())\nprint(\"Bond types: \", bond_types)\nfor (i, pair) in enumerate(unique_pairs):\n    ax = fig.add_subplot(3, 3, i + 1)\n    for bond_type in sorted([b for b in bond_types if pair == bond_type_to_pair(b)]):\n        for (df, name) in [(train_bonds, \"train\"), (test_bonds, \"test\")]:\n            #errors = df[(df.bond_type == bond_type) & (df.error == 1)].L2dist.values\n            #if len(errors) > 0:\n            #    sns.distplot(errors, color=\"red\", kde=True)\n            selection = df[(df.bond_type == bond_type) & (df.error == 0)].L2dist.values\n            sns.distplot(selection, label=f\"{bond_type} {name} {len(selection)/len(df)*100:.2f}%\", kde=True, ax=ax)\n    ax.set_ylim(1, ax.get_ylim()[1])\n    ax.set_yscale('log')\n    ax.set_xlabel(\"L2dist $[\\mathring{A}]$\")\n    ratio_tr  = (train_bonds.bond_type.apply(bond_type_to_pair) == pair).sum() / len(train_bonds)\n    ratio_tst = (test_bonds.bond_type.apply(bond_type_to_pair) == pair).sum() / len(test_bonds)\n    ax.set_title(f\"{pair} in {ratio_tr*100:.1f}% of train, {ratio_tst*100:.1f}% of test\")\n    ax.legend()\nfig.suptitle(\"Bond length by atom pair and number of bonds\")\nfig.tight_layout(rect=[0, 0.03, 1, 0.95])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Train and test set does match very well for the relative distribution of bonds length. This is good news as the train and test should behave the same for the prediction of the *scalar coupling constant*.\n\nThe number of bond does also build distributions peaking at different bond distance, which is consistent with the expected behaviour: the number of bonds decrease the distance between atoms. The distribution of the `1.0CC`, `2.0CC` and `3.0CC` is a perfect example of very well separated distributions.\n\nAdding the `COO-` handling improved a lot things by separating the previous `1CO` bimodal distribution into two well defined peaks, one for `1.5CO` and the second fonr `1.0C0`.  \n=> some other bimodal distributions like `1.0CN`, which is contained in `6.5%` of bonds, can be expected to be resolved the same way if needed.\n"},{"metadata":{},"cell_type":"markdown","source":"## Close analysis of `1.0CN` double peak in train set\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure( figsize = (14,2) )\nselection = train_bonds[(train_bonds.bond_type == '1.0CN') & (df.error == 0)].L2dist\nax = sns.distplot(selection, label=f\"'1.0CN' {name} {len(selection)/len(df)*100:.2f}%\", kde=True)\nxticks = plt.xticks(np.linspace(1.3, 1.6, 16))\ntext = plt.title(\"distribution of 1.0CO bond length for train set\")\nax.annotate(\"Peak1\", (1.37, 8))\nax.annotate(\"Peak2\", (1.46, 12))\nxl = plt.xlim(1.3, 1.6)\nyl = plt.ylim(0, 14)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def select_range(df, l0, l1):\n    return df[(df.bond_type == '1.0CN') & (df.error.values == 0) & (df.L2dist.values > l0) & (df.L2dist.values < l1)].molecule_name.unique()\npeak1_molecules = select_range(train_bonds, 1.37, 1.375)\npeak2_molecules = select_range(train_bonds, 1.465, 1.470)\n\nncols = 3\nnrows = 3\nfig1 = plt.figure(figsize = (18, 12))\nfor i in range(ncols*nrows):\n    ax = fig1.add_subplot(nrows, ncols, i + 1, projection=\"3d\")\n    plot_molecule(peak1_molecules[i], bonds=train_bonds, charges=train_charges, ax=ax)\ntext = fig1.suptitle(\"molecules with 1CN from Peak 1\")\nfig1.tight_layout(rect=[0, 0.03, 1, 0.95])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig2 = plt.figure(figsize = (18, 12))\nfor i in range(ncols*nrows):\n    ax = fig2.add_subplot(nrows, ncols, i + 1, projection=\"3d\")\n    plot_molecule(peak2_molecules[i], bonds=train_bonds, charges=train_charges, ax=ax)\ntext = fig2.suptitle(\"molecules with 1.0CN from Peak 2\")\nfig2.tight_layout(rect=[0, 0.03, 1, 0.95])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Double bonds are in green, and triple in red, for further details on the visualization see [Training set molecule visualization](https://www.kaggle.com/asauve/training-set-molecule-visualization?scriptVersionId=15133086).\n\nThe produced bonding schemes seems healthy for both set of molecules in the selected examples.\n\nObservable differences in the selected set:\n* For Peak 1 (shorter bonds): the C contained in the `1.0CN` is double bonded with another molecule\n* For peak 2 (larger bonds): NC bonds are part of a 3 atom CCN cycle.\n\n=> There are strong structural differences in the bonding of the two sets of molecules that may explain the `1.0CN` bimodal distribution.\n\nThese observation could be confirmed on large scale with :\n* either a neighbour configuration (easy) selection: how many bonds have the atoms from the pair with which other atoms?\n* either a cycle detection code (long)\n"},{"metadata":{},"cell_type":"markdown","source":"# Example of molecule with 4 bonds on `N` and ionization\n\nThe following molecule failed to link all atoms if VALENCE[N]=3 and left two atoms with a dangling covalent bond.\n\nThe solution has been to add support for additional optional bonding (thanks to @mrhermann comments) and support for charged atoms."},{"metadata":{"trusted":true},"cell_type":"code","source":"alanine271 = \"dsgdb9nsd_000271\"\ndisplay(structures[structures.molecule_name == alanine271])\ndisplay(train_bonds[train_bonds.molecule_name == alanine271])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_molecule(alanine271, bonds=train_bonds, charges=train_charges, azim=20, elev=30)\nplot_molecule(alanine271, bonds=train_bonds, charges=train_charges, azim=110, elev=30)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"view3d_molecule(alanine271)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This molecule correspond zwitterionic form of [beta-alanine](https://en.wikipedia.org/wiki/%CE%92-Alanine)  $H_2N-CH_2-CH_2-COOH$ \n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def pairlen(A,B):\n    return np.linalg.norm([A.x.values-B.x.values, A.y.values-B.y.values, A.z.values-B.z.values])\nm = structures[structures.molecule_name == alanine271]\nN  = m[m.atom_index==0]\nC4 = m[m.atom_index==4];   C5 = m[m.atom_index==5]; C6 = m[m.atom_index==6]\nH1 = m[m.atom_index==1]; H2 = m[m.atom_index==2]; H3 = m[m.atom_index==3]; H10 = m[m.atom_index==10]; H11 = m[m.atom_index==11]\nO7 = m[m.atom_index==7]; O8 = m[m.atom_index==8]\nprint(\"N-C4=\", pairlen(N, C4))\nprint(\"N-C5=\", pairlen(N, C5))\nprint(\"C4-C5\", pairlen(C4, C5))\nprint(\"N-H1=\", pairlen(N, H1))\nprint(\"N-H2=\", pairlen(N, H2))\nprint(\"N-H3=\", pairlen(N, H3))\nprint(\"C5-H1\", pairlen(C5, H1))\nprint(\"C4-H3\", pairlen(C4, H3))\nprint(\"O8-H11\", pairlen(O8, H11))\nprint(\"O7-H10\", pairlen(O7, H10))\nprint(\"C6-O7\", pairlen(C6, O7))\nprint(\"C6-O8\", pairlen(C6, O8))\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The NH3 bonding seem consistent with the distance to its nearby `H`, but:\n* C4 has only 3 covalent bonds\n  * N-C4 = 1.56 and is in the upper range of realistic NC bonds : this may be the correct bond\n* O8 has only one covalent bond \n  * none of the nearby atom is a candidate for bonding (expected O-H ~ 0.97)\n\nI can't figure out how this could work to build an Alanine isomere from both `O` location.  \nIs there a bug in the `structure` for this molecule?"},{"metadata":{},"cell_type":"markdown","source":"## Distribution of `molecule_name` index\n\nThe index contained in the molecule name contains no information about molecule structure AFAIK, but looking at the distribution can provide hints about the train / test split."},{"metadata":{"trusted":true},"cell_type":"code","source":"# select molecule index\nfor df, name in [(train_bonds, \"train\"), (test_bonds, \"test\")]:\n    index_list = [float(name[10:] ) for name in df.molecule_name.unique() ]\n    ax = sns.distplot(index_list, label=f\"{name}\")\n    ax.set_ylabel(\"molecule_name index distribution\")\nax.set_xlim(0, ax.get_xlim()[1])\nlegend = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Very well: molecule name index is equally distributed among both datasets!"},{"metadata":{},"cell_type":"markdown","source":"# Molecules with cycles"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_cycles.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Number of cycles per molecule\n\nA high percentage of molecules in `train` and `test` dataset have 1 to 4 cycles. One again, train and test set are very homogenous for the cycle repartition."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_count_ratio(df, x, y, hue, ax=None, add_text=False,\n                     text_dy=0.05, text_rot=90, text_fmt=\"%.2f\"):\n    # https://github.com/mwaskom/seaborn/issues/1027      # normalize countplot\n    prop_df = (df[x]\n           .groupby(df[hue])\n           .value_counts(normalize=True)\n           .rename(y)\n           .reset_index())\n    ax = sns.barplot(x=x, y=y, hue=hue, data=prop_df, ax=ax)\n    if add_text:\n        # https://stackoverflow.com/questions/31749448/how-to-add-percentages-on-top-of-bars-in-seaborn\n        for p in ax.patches:\n            height = p.get_height()\n            ax.text(p.get_x() + p.get_width()/2.,\n                    height + text_dy,\n                    text_fmt % (height,),\n                    ha=\"center\", rotation = text_rot)\n    return ax\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_nbcycles():\n    df1 = train_cycles.copy()\n    df1[\"dataset\"] = [\"train\"] * len(df1)\n    df2 = test_cycles.copy()\n    df2[\"dataset\"] = [\"test\"] * len(df2)\n    df = pd.concat([df1, df2])\n    grp = (df.groupby(by=[\"molecule_name\", \"dataset\"])\n        .max()\n        .reset_index())\n    plot_count_ratio(grp, \"cycle_index\", \"ratio\", \"dataset\", add_text=True, text_fmt=\"%.3f\")\n\nplot_nbcycles()\ntxt = plt.title(\"Number of cycles per molecule\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Number of atoms per cycle"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_cycle_size():\n    df1 = train_cycles.copy()\n    df1[\"dataset\"] = [\"train\"] * len(df1)\n    df2 = test_cycles.copy()\n    df2[\"dataset\"] = [\"test\"] * len(df2)\n    df = pd.concat([df1, df2])\n    df.cycle_seq += 1\n    grp = df.groupby(by=[\"molecule_name\", \"cycle_index\", \"dataset\"]).max().reset_index()\n    plot_count_ratio(grp, \"cycle_seq\", \"ratio\", \"dataset\", add_text=True, text_fmt=\"%.4f\")\n    \nplot_cycle_size()\ntxt = plt.title(\"          Number of atoms per cycle\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here again test and train are very consistent.   \nMost of cycles count from 3 to 5 atoms while aromatic cycles account only for 8.5% of the total number of cycles.  \n\nCuriously there are up to 9 atoms per cycles. The 9-cycle case is real as can be seen bellow (you can rotate the molecule in the bottom view to check by yourself). This one is really twisted."},{"metadata":{"trusted":true},"cell_type":"code","source":"#grp = train_cycles.groupby(by=[\"molecule_name\",\"cycle_seq\"]).max().reset_index()\n#grp[grp.cycle_seq == 8].molecule_name\nplot_molecule(\"dsgdb9nsd_024553\", bonds=train_bonds, charges=train_charges)\nview3d_molecule(\"dsgdb9nsd_024553\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cycle composition"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_cycle_composition():\n    df1 = train_cycles.copy()\n    df1[\"dataset\"] = [\"train\"] * len(df1)\n    df2 = test_cycles.copy()\n    df2[\"dataset\"] = [\"test\"] * len(df2)\n    df = pd.concat([df1, df2])\n    \n    # add atom name\n    df = df.merge(how=\"left\", on = [\"molecule_name\", \"atom_index\"], \n                   right = structures[[\"molecule_name\", \"atom_index\", \"atom\"]])\n\n    plot_count_ratio(df, \"atom\", \"ratio\", \"dataset\", add_text=True, text_dy=0.1, text_fmt=\"%.2f\")\nplot_cycle_composition()\ntxt = plt.title(\"Cycle composition\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"No surprise here: cycles a made of 82% `C` atoms."},{"metadata":{},"cell_type":"markdown","source":"## A curiosity: cage-like molecule"},{"metadata":{},"cell_type":"markdown","source":"There is one example of molecule with 4 cycles:  in pratice there are 5 cycles in this one but the last one remains undetected by the code as all its bonds are already included in other cycles.\n\nThis funky molecule name is *[prismane](https://fr.wikipedia.org/wiki/3-Prismane)* and has an interresting story as it was a candidate for benzene structures before it was discovered benzene was a cycle."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_molecule(\"dsgdb9nsd_000719\", bonds=train_bonds, charges=train_charges)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Write outputs\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_bonds.to_csv('train_bonds.csv', index=False)\ntrain_charges.to_csv('train_charges.csv', index=False)\ntrain_cycles.to_csv('train_cycles.csv', index=False)\n\ntest_bonds.to_csv('test_bonds.csv', index=False)\ntest_charges.to_csv('test_charges.csv', index=False)\ntest_cycles.to_csv('test_cycles.csv', index=False)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What to do next?\n\n* Compare with output from the [Open Babel](http://openbabel.org/wiki/Babel) software \n* Fix molecules with failed bonding (more than 2%)\n  * some of them are cycles which build double bonds in the wrong order\n  * some ionized groups are not well handled\n* (**DONE**) Cycles plays an important role in deshielding of atoms.\nThe next step would be to build a cycle detection code in order to test new features related to cycles.\n* (**DONE** by adding charge support) There are issues with some molecules: can a seasonned chemist guess how `dsgdb9nsd_000271` could be bonded correctly?\n"},{"metadata":{},"cell_type":"markdown","source":"# Credits\n\nThanks to [M. Hermann](https://www.kaggle.com/mrhermann) for his advices about bond properties and the Zwitterionic ionized form.\n\nThanks to "}],"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}