{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\n\n![](http://)\n\n# Table of Contents\n\n**1. [Introduction](#id1)**
\n**2. [Setup](#id2)**
\n**3. [EDA](#id3)**
\n**4. [Graph Neural Network modeling](#id4)**
\n**5. [Going further](#id5)**
"},{"metadata":{},"cell_type":"markdown","source":"
\n# Introduction\n\nIn this competition, you are given data of molecules with each atom position and need to predict `scalar_coupling_constant`, which is defined on pair of atoms.\n\nIn this kernel, I will introduce following.\n\n**EDA: Exploratory Data Analysis**\n - Understand competition data, visualize properties & molecules.\n - Explanation about QM9 data.\n\n**Baseline modeling using Graph Convolutional Neural Network based model**\n - I will use `WeaveNet` provided by `chainer-chemistry` library.\n\n**Going further**\n - Approach we can consider to get good score.\n - Related library","attachments":{}},{"metadata":{},"cell_type":"markdown","source":"
\n# Setup & Installing modules\n\nI will use following library which is not installed to Kaggle docker by default.\n\n - [RDKit](https://github.com/rdkit/rdkit): Chemistry preprocessing and visualization.\n - [cupy](https://github.com/cupy/cupy): To use GPU with Chainer. It supports numpy-like API for GPU array processing. \n - [chainer-chemistry](https://github.com/pfnet-research/chainer-chemistry): Many kinds of Graph-convolution based network and its data preprocessing is implemented.\n - [chaineripy](https://github.com/grafi-tt/chaineripy): To use GPU with Chainer. It supports numpy-like API for GPU array processing. \n\n`cupy`, `chainer-chemistry` and `chaineripy` can be install via pip.\n\nTo install `rdkit`, you can install via conda package, by the following command.\n\nAlso, please turn on \"GPU\" and \"Internet\" in the Settings tab on right side to run deep learning training on GPU and install library.\n"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"!pip install --quiet cupy-cuda100==5.4.0\n!pip install --quiet chainer-chemistry==0.5.0\n!pip install --quiet chaineripy\n!conda install -y --quiet -c rdkit rdkit","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Check correctly installed, and modules can be imported.\nimport chainer\nimport chainer_chemistry\nimport chaineripy\nimport cupy\nimport rdkit\n\nprint('chainer version: ', chainer.__version__)\nprint('cupy version: ', cupy.__version__)\nprint('chainer-chemistry version: ', chainer_chemistry.__version__)\nprint('rdkit version: ', rdkit.__version__)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## import library and Define util functions\n\nCode from [Visualize molecules with RDKit](https://www.kaggle.com/corochann/visualize-molecules-with-rdkit) kernel."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"from contextlib import contextmanager\nimport gc\nimport numpy as np # linear algebra\nimport numpy\nimport os\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom pathlib import Path\nfrom time import time, perf_counter\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nfrom chainer_chemistry.datasets.numpy_tuple_dataset import NumpyTupleDataset\n\nimport rdkit\nfrom rdkit import Chem","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# util functions...\n\n@contextmanager\ndef timer(name):\n t0 = perf_counter()\n yield\n t1 = perf_counter()\n print('[{}] done in {:.3f} s'.format(name, t1-t0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"\"\"\"\nCopied from\nhttps://github.com/jensengroup/xyz2mol/blob/master/xyz2mol.py\n\nModified `chiral_stereo_check` method for this task's purpose.\n\"\"\"\n##\n# Written by Jan H. Jensen based on this paper Yeonjoon Kim and Woo Youn Kim\n# \"Universal Structure Conversion Method for Organic Molecules: From Atomic Connectivity\n# to Three-Dimensional Geometry\" Bull. Korean Chem. Soc. 2015, Vol. 36, 1769-1777 DOI: 10.1002/bkcs.10334\n#\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nimport itertools\nfrom rdkit.Chem import rdmolops\nfrom collections import defaultdict\nimport copy\nimport networkx as nx # uncomment if you don't want to use \"quick\"/install networkx\n\nglobal __ATOM_LIST__\n__ATOM_LIST__ = [x.strip() for x in ['h ', 'he', \\\n 'li', 'be', 'b ', 'c ', 'n ', 'o ', 'f ', 'ne', \\\n 'na', 'mg', 'al', 'si', 'p ', 's ', 'cl', 'ar', \\\n 'k ', 'ca', 'sc', 'ti', 'v ', 'cr', 'mn', 'fe', 'co', 'ni', 'cu', \\\n 'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr', \\\n 'rb', 'sr', 'y ', 'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag', \\\n 'cd', 'in', 'sn', 'sb', 'te', 'i ', 'xe', \\\n 'cs', 'ba', 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy', \\\n 'ho', 'er', 'tm', 'yb', 'lu', 'hf', 'ta', 'w ', 're', 'os', 'ir', 'pt', \\\n 'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn', \\\n 'fr', 'ra', 'ac', 'th', 'pa', 'u ', 'np', 'pu']]\n\n\ndef get_atom(atom):\n global __ATOM_LIST__\n atom = atom.lower()\n return __ATOM_LIST__.index(atom) + 1\n\n\ndef getUA(maxValence_list, valence_list):\n UA = []\n DU = []\n for i, (maxValence, valence) in enumerate(zip(maxValence_list, valence_list)):\n if maxValence - valence > 0:\n UA.append(i)\n DU.append(maxValence - valence)\n return UA, DU\n\n\ndef get_BO(AC, UA, DU, valences, UA_pairs, quick):\n BO = AC.copy()\n DU_save = []\n\n while DU_save != DU:\n for i, j in UA_pairs:\n BO[i, j] += 1\n BO[j, i] += 1\n\n BO_valence = list(BO.sum(axis=1))\n DU_save = copy.copy(DU)\n UA, DU = getUA(valences, BO_valence)\n UA_pairs = get_UA_pairs(UA, AC, quick)[0]\n\n return BO\n\n\ndef valences_not_too_large(BO, valences):\n number_of_bonds_list = BO.sum(axis=1)\n for valence, number_of_bonds in zip(valences, number_of_bonds_list):\n if number_of_bonds > valence:\n return False\n\n return True\n\n\ndef BO_is_OK(BO, AC, charge, DU, atomic_valence_electrons, atomicNumList, charged_fragments):\n Q = 0 # total charge\n q_list = []\n if charged_fragments:\n BO_valences = list(BO.sum(axis=1))\n for i, atom in enumerate(atomicNumList):\n q = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n Q += q\n if atom == 6:\n number_of_single_bonds_to_C = list(BO[i, :]).count(1)\n if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:\n Q += 1\n q = 2\n if number_of_single_bonds_to_C == 3 and Q + 1 < charge:\n Q += 2\n q = 1\n\n if q != 0:\n q_list.append(q)\n\n if (BO - AC).sum() == sum(DU) and charge == Q and len(q_list) <= abs(charge):\n return True\n else:\n return False\n\n\ndef get_atomic_charge(atom, atomic_valence_electrons, BO_valence):\n if atom == 1:\n charge = 1 - BO_valence\n elif atom == 5:\n charge = 3 - BO_valence\n elif atom == 15 and BO_valence == 5:\n charge = 0\n elif atom == 16 and BO_valence == 6:\n charge = 0\n else:\n charge = atomic_valence_electrons - 8 + BO_valence\n\n return charge\n\n\ndef clean_charges(mol):\n # this hack should not be needed any more but is kept just in case\n #\n\n rxn_smarts = ['[N+:1]=[*:2]-[C-:3]>>[N+0:1]-[*:2]=[C-0:3]',\n '[N+:1]=[*:2]-[O-:3]>>[N+0:1]-[*:2]=[O-0:3]',\n '[N+:1]=[*:2]-[*:3]=[*:4]-[O-:5]>>[N+0:1]-[*:2]=[*:3]-[*:4]=[O-0:5]',\n '[#8:1]=[#6:2]([!-:6])[*:3]=[*:4][#6-:5]>>[*-:1][*:2]([*:6])=[*:3][*:4]=[*+0:5]',\n '[O:1]=[c:2][c-:3]>>[*-:1][*:2][*+0:3]',\n '[O:1]=[C:2][C-:3]>>[*-:1][*:2]=[*+0:3]']\n\n fragments = Chem.GetMolFrags(mol, asMols=True, sanitizeFrags=False)\n\n for i, fragment in enumerate(fragments):\n for smarts in rxn_smarts:\n patt = Chem.MolFromSmarts(smarts.split(\">>\")[0])\n while fragment.HasSubstructMatch(patt):\n rxn = AllChem.ReactionFromSmarts(smarts)\n ps = rxn.RunReactants((fragment,))\n fragment = ps[0][0]\n if i == 0:\n mol = fragment\n else:\n mol = Chem.CombineMols(mol, fragment)\n\n return mol\n\n\ndef BO2mol(mol, BO_matrix, atomicNumList, atomic_valence_electrons, mol_charge, charged_fragments):\n # based on code written by Paolo Toscani\n\n l = len(BO_matrix)\n l2 = len(atomicNumList)\n BO_valences = list(BO_matrix.sum(axis=1))\n\n if (l != l2):\n raise RuntimeError('sizes of adjMat ({0:d}) and atomicNumList '\n '{1:d} differ'.format(l, l2))\n\n rwMol = Chem.RWMol(mol)\n\n bondTypeDict = {\n 1: Chem.BondType.SINGLE,\n 2: Chem.BondType.DOUBLE,\n 3: Chem.BondType.TRIPLE\n }\n\n for i in range(l):\n for j in range(i + 1, l):\n bo = int(round(BO_matrix[i, j]))\n if (bo == 0):\n continue\n bt = bondTypeDict.get(bo, Chem.BondType.SINGLE)\n rwMol.AddBond(i, j, bt)\n mol = rwMol.GetMol()\n\n if charged_fragments:\n mol = set_atomic_charges(mol, atomicNumList, atomic_valence_electrons, BO_valences, BO_matrix, mol_charge)\n else:\n mol = set_atomic_radicals(mol, atomicNumList, atomic_valence_electrons, BO_valences)\n\n return mol\n\n\ndef set_atomic_charges(mol, atomicNumList, atomic_valence_electrons, BO_valences, BO_matrix, mol_charge):\n q = 0\n for i, atom in enumerate(atomicNumList):\n a = mol.GetAtomWithIdx(i)\n charge = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n q += charge\n if atom == 6:\n number_of_single_bonds_to_C = list(BO_matrix[i, :]).count(1)\n if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:\n q += 1\n charge = 0\n if number_of_single_bonds_to_C == 3 and q + 1 < mol_charge:\n q += 2\n charge = 1\n\n if (abs(charge) > 0):\n a.SetFormalCharge(int(charge))\n\n # shouldn't be needed anymore bit is kept just in case\n # mol = clean_charges(mol)\n\n return mol\n\n\ndef set_atomic_radicals(mol, atomicNumList, atomic_valence_electrons, BO_valences):\n # The number of radical electrons = absolute atomic charge\n for i, atom in enumerate(atomicNumList):\n a = mol.GetAtomWithIdx(i)\n charge = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n\n if (abs(charge) > 0):\n a.SetNumRadicalElectrons(abs(int(charge)))\n\n return mol\n\n\ndef get_bonds(UA, AC):\n bonds = []\n\n for k, i in enumerate(UA):\n for j in UA[k + 1:]:\n if AC[i, j] == 1:\n bonds.append(tuple(sorted([i, j])))\n\n return bonds\n\n\ndef get_UA_pairs(UA, AC, quick):\n bonds = get_bonds(UA, AC)\n if len(bonds) == 0:\n return [()]\n\n if quick:\n G = nx.Graph()\n G.add_edges_from(bonds)\n UA_pairs = [list(nx.max_weight_matching(G))]\n return UA_pairs\n\n max_atoms_in_combo = 0\n UA_pairs = [()]\n for combo in list(itertools.combinations(bonds, int(len(UA) / 2))):\n flat_list = [item for sublist in combo for item in sublist]\n atoms_in_combo = len(set(flat_list))\n if atoms_in_combo > max_atoms_in_combo:\n max_atoms_in_combo = atoms_in_combo\n UA_pairs = [combo]\n # if quick and max_atoms_in_combo == 2*int(len(UA)/2):\n # return UA_pairs\n elif atoms_in_combo == max_atoms_in_combo:\n UA_pairs.append(combo)\n\n return UA_pairs\n\n\ndef AC2BO(AC, atomicNumList, charge, charged_fragments, quick):\n # TODO\n atomic_valence = defaultdict(list)\n atomic_valence[1] = [1]\n atomic_valence[6] = [4]\n atomic_valence[7] = [4, 3]\n atomic_valence[8] = [2, 1]\n atomic_valence[9] = [1]\n atomic_valence[14] = [4]\n atomic_valence[15] = [5, 4, 3]\n atomic_valence[16] = [6, 4, 2]\n atomic_valence[17] = [1]\n atomic_valence[32] = [4]\n atomic_valence[35] = [1]\n atomic_valence[53] = [1]\n\n atomic_valence_electrons = {}\n atomic_valence_electrons[1] = 1\n atomic_valence_electrons[6] = 4\n atomic_valence_electrons[7] = 5\n atomic_valence_electrons[8] = 6\n atomic_valence_electrons[9] = 7\n atomic_valence_electrons[14] = 4\n atomic_valence_electrons[15] = 5\n atomic_valence_electrons[16] = 6\n atomic_valence_electrons[17] = 7\n atomic_valence_electrons[32] = 4\n atomic_valence_electrons[35] = 7\n atomic_valence_electrons[53] = 7\n\n # make a list of valences, e.g. for CO: [[4],[2,1]]\n valences_list_of_lists = []\n for atomicNum in atomicNumList:\n valences_list_of_lists.append(atomic_valence[atomicNum])\n\n # convert [[4],[2,1]] to [[4,2],[4,1]]\n valences_list = list(itertools.product(*valences_list_of_lists))\n\n best_BO = AC.copy()\n\n # implemenation of algorithm shown in Figure 2\n # UA: unsaturated atoms\n # DU: degree of unsaturation (u matrix in Figure)\n # best_BO: Bcurr in Figure\n #\n\n for valences in valences_list:\n AC_valence = list(AC.sum(axis=1))\n UA, DU_from_AC = getUA(valences, AC_valence)\n\n if len(UA) == 0 and BO_is_OK(AC, AC, charge, DU_from_AC, atomic_valence_electrons, atomicNumList,\n charged_fragments):\n return AC, atomic_valence_electrons\n\n UA_pairs_list = get_UA_pairs(UA, AC, quick)\n for UA_pairs in UA_pairs_list:\n BO = get_BO(AC, UA, DU_from_AC, valences, UA_pairs, quick)\n if BO_is_OK(BO, AC, charge, DU_from_AC, atomic_valence_electrons, atomicNumList, charged_fragments):\n return BO, atomic_valence_electrons\n\n elif BO.sum() >= best_BO.sum() and valences_not_too_large(BO, valences):\n best_BO = BO.copy()\n\n return best_BO, atomic_valence_electrons\n\n\ndef AC2mol(mol, AC, atomicNumList, charge, charged_fragments, quick):\n # convert AC matrix to bond order (BO) matrix\n BO, atomic_valence_electrons = AC2BO(AC, atomicNumList, charge, charged_fragments, quick)\n\n # add BO connectivity and charge info to mol object\n mol = BO2mol(mol, BO, atomicNumList, atomic_valence_electrons, charge, charged_fragments)\n\n return mol\n\n\ndef get_proto_mol(atomicNumList):\n mol = Chem.MolFromSmarts(\"[#\" + str(atomicNumList[0]) + \"]\")\n rwMol = Chem.RWMol(mol)\n for i in range(1, len(atomicNumList)):\n a = Chem.Atom(atomicNumList[i])\n rwMol.AddAtom(a)\n\n mol = rwMol.GetMol()\n\n return mol\n\n\ndef get_atomicNumList(atomic_symbols):\n atomicNumList = []\n for symbol in atomic_symbols:\n atomicNumList.append(get_atom(symbol))\n return atomicNumList\n\n\ndef read_xyz_file(filename):\n atomic_symbols = []\n xyz_coordinates = []\n\n with open(filename, \"r\") as file:\n for line_number, line in enumerate(file):\n if line_number == 0:\n num_atoms = int(line)\n elif line_number == 1:\n if \"charge=\" in line:\n charge = int(line.split(\"=\")[1])\n else:\n charge = 0\n else:\n atomic_symbol, x, y, z = line.split()\n atomic_symbols.append(atomic_symbol)\n xyz_coordinates.append([float(x), float(y), float(z)])\n\n atomicNumList = get_atomicNumList(atomic_symbols)\n\n return atomicNumList, charge, xyz_coordinates\n\n\ndef xyz2AC(atomicNumList, xyz):\n import numpy as np\n mol = get_proto_mol(atomicNumList)\n\n conf = Chem.Conformer(mol.GetNumAtoms())\n for i in range(mol.GetNumAtoms()):\n conf.SetAtomPosition(i, (xyz[i][0], xyz[i][1], xyz[i][2]))\n mol.AddConformer(conf)\n\n dMat = Chem.Get3DDistanceMatrix(mol)\n pt = Chem.GetPeriodicTable()\n\n num_atoms = len(atomicNumList)\n AC = np.zeros((num_atoms, num_atoms)).astype(int)\n\n for i in range(num_atoms):\n a_i = mol.GetAtomWithIdx(i)\n Rcov_i = pt.GetRcovalent(a_i.GetAtomicNum()) * 1.30\n for j in range(i + 1, num_atoms):\n a_j = mol.GetAtomWithIdx(j)\n Rcov_j = pt.GetRcovalent(a_j.GetAtomicNum()) * 1.30\n if dMat[i, j] <= Rcov_i + Rcov_j:\n AC[i, j] = 1\n AC[j, i] = 1\n\n return AC, mol\n\n\ndef chiral_stereo_check(mol):\n # Chem.SanitizeMol(mol)\n num_error = Chem.SanitizeMol(mol, Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES, catchErrors=True)\n if num_error != 0:\n print('error id', num_error)\n\n Chem.DetectBondStereochemistry(mol, -1)\n Chem.AssignStereochemistry(mol, flagPossibleStereoCenters=True, force=True)\n Chem.AssignAtomChiralTagsFromStructure(mol, -1)\n\n return mol\n\n\ndef xyz2mol(atomicNumList, charge, xyz_coordinates, charged_fragments, quick):\n # Get atom connectivity (AC) matrix, list of atomic numbers, molecular charge,\n # and mol object with no connectivity information\n AC, mol = xyz2AC(atomicNumList, xyz_coordinates)\n\n # Convert AC to bond order matrix and add connectivity and charge info to mol object\n new_mol = AC2mol(mol, AC, atomicNumList, charge, charged_fragments, quick)\n\n # Check for stereocenters and chiral centers\n new_mol = chiral_stereo_check(new_mol)\n\n return new_mol\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def mol_from_xyz(filepath, to_canonical=False):\n charged_fragments = True # alternatively radicals are made\n\n # quick is faster for large systems but requires networkx\n # if you don't want to install networkx set quick=False and\n # uncomment 'import networkx as nx' at the top of the file\n quick = True\n\n atomicNumList, charge, xyz_coordinates = read_xyz_file(filepath)\n # print('atomicNumList', atomicNumList, 'charge', charge)\n mol = xyz2mol(atomicNumList, charge, xyz_coordinates, charged_fragments, quick)\n\n # Canonical hack\n if to_canonical:\n smiles = Chem.MolToSmiles(mol, isomericSmiles=True)\n mol = Chem.MolFromSmiles(smiles)\n return mol","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# This script is referred from http://rdkit.blogspot.jp/2015/02/new-drawing-code.html\n# and http://cheminformist.itmol.com/TEST/wp-content/uploads/2015/07/rdkit_moldraw2d_2.html\nfrom __future__ import print_function\nfrom rdkit import Chem\nfrom rdkit.Chem.Draw import IPythonConsole\nfrom IPython.display import SVG\n\nfrom rdkit.Chem import rdDepictor\nfrom rdkit.Chem.Draw import rdMolDraw2D\ndef moltosvg(mol,molSize=(450,150),kekulize=True):\n mc = Chem.Mol(mol.ToBinary())\n if kekulize:\n try:\n Chem.Kekulize(mc)\n except:\n mc = Chem.Mol(mol.ToBinary())\n if not mc.GetNumConformers():\n rdDepictor.Compute2DCoords(mc)\n drawer = rdMolDraw2D.MolDraw2DSVG(molSize[0],molSize[1])\n drawer.DrawMolecule(mc)\n drawer.FinishDrawing()\n svg = drawer.GetDrawingText()\n return svg\n\ndef render_svg(svg):\n # It seems that the svg renderer used doesn't quite hit the spec.\n # Here are some fixes to make it work in the notebook, although I think\n # the underlying issue needs to be resolved at the generation step\n return SVG(svg.replace('svg:',''))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n\n# EDA: Exploratory Data Analysis\n\nLet's begin data analaysis!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Input data files are available in the \"../input/\" directory.\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Load each data...\n\ninput_dir = Path('../input/champs-scalar-coupling')\ntrain_df = pd.read_csv(input_dir/'train.csv')\ntest_df = pd.read_csv(input_dir/'test.csv')\nstructures_df = pd.read_csv(input_dir/'structures.csv')\nsample_submission_df = pd.read_csv(input_dir/'sample_submission.csv')\n\ndipole_moments_df = pd.read_csv(input_dir/'dipole_moments.csv')\nmagnetic_shielding_tensors_df = pd.read_csv(input_dir/'magnetic_shielding_tensors.csv')\nmulliken_charges_df = pd.read_csv(input_dir/'mulliken_charges.csv')\npotential_energy_df = pd.read_csv(input_dir/'potential_energy.csv')\nscalar_coupling_contributions_df = pd.read_csv(input_dir/'scalar_coupling_contributions.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**train_df** and **test_df** stores the main information of train/test data.\n\n`scalar_coupling_constant` is the value to predict in this competition.\nAs you can see, `molecule_name` is duprecated in the train data. Because `scalar_coupling_constant` is defined on **atom pair** basis.\nSo you need to predict several values in each molecule."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"train_df and test_df contains minimal information. \n**structures_df** contains more detailed information, the type & position of each atoms in the molecules."},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"sample_submission_df is the format to submit. When you predict values from test_df, we can assign predicted values in this format. The `id` column corresponds to the `test_df`."},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These 4 DataFrames are the main information.\n\nFollowing DataFrames are for additional information, if you can come up with the nice usage to combine these properties you can use it. However these information is **provided only on train data**.\nSo directly using these values to predict our target, `scalar_coupling_constant`, for test data is not possible."},{"metadata":{},"cell_type":"markdown","source":"dipole_moments.csv stores the molecular electric dipole moments for **each molecule**.\n\nDipole moment is a vector value, storing (x, y, z) information."},{"metadata":{"trusted":true},"cell_type":"code","source":"dipole_moments_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"potential_energy.csv stores the potential energy for **each molecule**.\n\nIt is scalar value, only one value is stored."},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy_df['potential_energy'].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"magnetic_shielding_tensors.csv stores the magnetic shielding tensors **for all atoms in the molecules each molecule**.\n\nIt is tensor value, (XX, YX, ZX, XY, YY, ZY, XZ, YZ, ZZ) information is stored for each atom specified by `atom_index`."},{"metadata":{"trusted":true},"cell_type":"code","source":"magnetic_shielding_tensors_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"mulliken_charges.csv stores the mulliken charges for **all atoms in the molecules**.\n\nIt is scalar value, one value is stored for each atom specified by `atom_index`.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_charges_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_charges_df['mulliken_charge'].hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"scalar_coupling_contributions.csv: \n> The scalar coupling constants in train.csv (or corresponding files) are a sum of four terms. scalar_coupling_contributions.csv contain all these terms.\n> the fifth column (fc) is the Fermi Contact contribution, the sixth column (sd) is the Spin-dipolar contribution, the seventh column (pso) is the Paramagnetic spin-orbit contribution and the eighth column (dso) is the Diamagnetic spin-orbit contribution.\n\nAs we can see, the first 3 column is the same format with train data, and there are 4 float values `fc`, `sd`, `pso`, and `dso`.\nThese sum is the actual `scalar_coupling_constant` to predict in this competition."},{"metadata":{"trusted":true},"cell_type":"code","source":"scalar_coupling_contributions_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So `scalar_coupling_constant` which we need to predict is defined on **atom pair** basis, but many information comes as **each atom** or **molecule** basis.\n\nHow to combine these information is the challenge in this competition."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Atom pair wise info\nprint('--- Atom pair wise info ---')\nprint(f'train_df {train_df.shape}')\nprint(f'test_df {test_df.shape}')\nprint(f'scalar_coupling_contributions_df {scalar_coupling_contributions_df.shape}') # only contain train info\nprint(f'sample_submission_df {sample_submission_df.shape}') # for test prediction\n\n# Molecule wise info\nprint('--- Molecule wise info ---')\nprint(f'potential_energy_df {potential_energy_df.shape}') # only contain train info\nprint(f'dipole_moments_df {dipole_moments_df.shape}') # only contain train info\n\n# Atom wise info\nprint('--- Atom wise info ---')\nprint(f'structures_df {structures_df.shape}') # contains both train/test info\nprint(f'magnetic_shielding_tensors_df {magnetic_shielding_tensors_df.shape}') # only contain train info\nprint(f'mulliken_charges_df {mulliken_charges_df.shape}') # only contain train info","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Number of information is molecule wise < atom wise < atom pair wise"},{"metadata":{},"cell_type":"markdown","source":"I don't use \"additional information\" in the following. Let's delete for memory efficiency."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"del dipole_moments_df\ndel magnetic_shielding_tensors_df\ndel mulliken_charges_df\ndel potential_energy_df\ndel scalar_coupling_contributions_df\n\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Number of atoms"},{"metadata":{},"cell_type":"markdown","source":"How many molecules or atoms are there in the dataset?"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mol_names = train_df['molecule_name'].unique()\ntest_mol_names = test_df['molecule_name'].unique()\n\nprint(f'Number of molecules: train {len(train_mol_names)}, test {len(test_mol_names)}, train+test {len(train_mol_names) + len(test_mol_names)}')\nprint(f'Number of total atoms: train+test {len(structures_df)}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So train/test dataset ratio is around 2:1."},{"metadata":{"trusted":true},"cell_type":"code","source":"print('name of molecues in train', train_mol_names)\nprint('name of molecues in test', test_mol_names)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Molecule name is assigned as increasing order from 000000 until 0133885, and randomly splitted into train and test."},{"metadata":{},"cell_type":"markdown","source":"**Atom type**: What kind of atom type are there in the train/test data?\n\nWe have only 5 atoms, H, C, N, O and F!"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures_df['atom'].unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How many number of atoms are there in each molecule?"},{"metadata":{"trusted":true},"cell_type":"code","source":"mol_name_counts = structures_df['molecule_name'].value_counts()\nmol_name_counts.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mol_name_counts.hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Number of atom in each molecule ranges **3~29**."},{"metadata":{},"cell_type":"markdown","source":"Many atoms are actually hydrogen H, how many \"heavy atom\" (C, N, O and F) are there?"},{"metadata":{"trusted":true},"cell_type":"code","source":"heavy_atoms_in_mol = structures_df[structures_df['atom'] != 'H']['molecule_name'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"heavy_atoms_in_mol.hist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Number of heavy atom in each molecule ranges **1~9**. And its histogram is very skewed, as the number of heavy atom increases, number of data increases."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Histogram of heavy atoms\nheavy_atoms_in_mol.value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems competition data actually comes from **QM9 dataset**!\n\nYou can refer QM9 dataset's explanation here.\n\n- http://quantum-machine.org/datasets/\n\n> Computational de novo design of new drugs and materials requires rigorous and unbiased exploration of chemical compound space. However, large uncharted territories persist due to its size scaling combinatorially with molecular size. We report computed geometric, energetic, electronic, and thermodynamic properties for 134k stable small organic molecules made up of CHONF. These molecules correspond to the subset of all 133,885 species with up to nine heavy atoms (CONF) out of the GDB-17 chemical universe of 166 billion organic molecules.\n\nAlso, you can download original QM9 dataset from below\n\n - https://figshare.com/collections/Quantum_chemistry_structures_and_properties_of_134_kilo_molecules/978904"},{"metadata":{},"cell_type":"markdown","source":"So this competition uses QM9 dataset as base molecule list, and calculated `scalar_coupling_constant` which itself is not calculated on original QM9 dataset. QM9 Dataset is actually **enumerated list of organic molecules with up to 9 heavy atoms (CHONF)**. Why the data contains much more 9 heavy atoms than 1 heavy atoms? Because there are much more combination to construct valid molecules when the number of heavy atoms increases!"},{"metadata":{},"cell_type":"markdown","source":"Also, as you can see below, molecule name is defined in alomost increasing order of heavy atoms."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"heavy_atoms_in_mol.sort_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nplt.scatter(np.arange(len(heavy_atoms_in_mol)), heavy_atoms_in_mol.sort_index().values)\nplt.xlabel('molecule_name')\nplt.ylabel('Number of heavy atoms')\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Molecule visualization"},{"metadata":{},"cell_type":"markdown","source":"1. Visualize molecules with RDKit (refer [Visualize molecules with RDKit](https://www.kaggle.com/corochann/visualize-molecules-with-rdkit) kernel)."},{"metadata":{"trusted":true},"cell_type":"code","source":"for mol_id in ['000001', '001000', '010000', '133885']:\n mol_name = f'dsgdb9nsd_{mol_id}'\n filepath = input_dir/f'structures/{mol_name}.xyz'\n mol = mol_from_xyz(filepath, to_canonical=True)\n print('mol_name:', mol_name)\n display(mol)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, molecule becomes larger and more complex structure with increasing order of mol_name."},{"metadata":{},"cell_type":"markdown","source":"### Scalar coupling constant"},{"metadata":{},"cell_type":"markdown","source":"Let's focus on `scalar_coupling_constant`, this competition's label to predict.\n\nAs we can see in `train_df` data, scalar coupling constant is defined on atom pair. \"coupling type\" is defined in each scalar coupling."},{"metadata":{"trusted":true},"cell_type":"code","source":"coupling_tyes = train_df['type'].unique()\nprint('coupling_tyes', coupling_tyes)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems we have 8 types of coupling type.\n\nIt is the format with number and \"J\" and some molecule types.\n\nBut what is this number 1, 2 or 3?\nI refer @geoman3 's answer in [this discussion](https://www.kaggle.com/c/champs-scalar-coupling/discussion/95204#latest-554649).\n\n![](https://i.imgur.com/1jGXP8L.png)\n\nIt seems the number indicates to how many bonds separated between these atom pair.\nWe have number up to 3, meaning that atom pair with too far each other is not in focus of this competition."},{"metadata":{},"cell_type":"markdown","source":"Let's see the histogram of `scalar_coupling_constant`."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nsns.distplot(train_df['scalar_coupling_constant'])\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How about histogram of `scalar_coupling_constant` for each coupling type?"},{"metadata":{"trusted":true},"cell_type":"code","source":"grid = sns.FacetGrid(train_df[['type', 'scalar_coupling_constant']], col='type', hue='type', col_wrap=4)\ngrid.map(sns.distplot, 'scalar_coupling_constant')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see the value ranges very different in each coupling type. \"1JHC\" coupling type is biggest, around 80. \"1JHN\" is the second biggest around 50. Other coupling types are around 0.\n\nIt may be better to handle each couping type explicitly (include coupling type feature, use different model etc)."},{"metadata":{},"cell_type":"markdown","source":"# Distance features\n\nMost simple feature we can come up with between atom pair is a distance between these 2 atoms. Let's see how distance related to `scalar_coupling_constant`."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Calculate distance feature...\n# --- dist hist ---\n\n# from https://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark\n\ndef map_atom_info(df, atom_idx):\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\ntrain_df = map_atom_info(train_df, 0)\ntrain_df = map_atom_info(train_df, 1)\n\ntest_df = map_atom_info(test_df, 0)\ntest_df = map_atom_info(test_df, 1)\n\ntrain_p_0 = train_df[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train_df[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test_df[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test_df[['x_1', 'y_1', 'z_1']].values\n\ntrain_df['dist_speedup'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest_df['dist_speedup'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# dist_speedup and atom positions are added to both train_df and test_df\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"grid = sns.FacetGrid(train_df[['type', 'dist_speedup']], col='type', hue='type', col_wrap=4)\ngrid.map(sns.distplot, 'dist_speedup')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"len(train_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df_sampled = train_df.sample(10000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.scatterplot(x=\"dist_speedup\", y=\"scalar_coupling_constant\", hue='type', data=train_df_sampled)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distance for 1J coupling is smaller since atom pair is neighbor to each other, and its scalar coupling constant is higher.\n\nDistance tend to be long from 2J to 3J coupling, but its scalar coupling constant value is similar.\n\nIt seems only distance feature is not enough, we need other information to differentiate scalar coupling."},{"metadata":{},"cell_type":"markdown","source":"
\n\n# Graph Neural Network modeling\n\nLet's use Graph convolutinal neural network for this task. I use [chainer-chemistry](https://github.com/pfnet-research/chainer-chemistry) which provides various kinds of graph convolution network and chemistry related data prerprocessing.\n\nConsider molecule as a \"graph\" and atom as a \"node\", we can use graph convolution on molecule.\nBenefit of graph convolutional neural network is that it can handle input in node order permutation way.\n\nHere I use WeaveNet implementation, original paper is the following.\n\n - [Molecular Graph Convolutions: Moving Beyond Fingerprints](https://arxiv.org/abs/1603.00856)\n \n \n\n\nUsually most of the graph convolution network deals with bond in molecule as \"edge\" information of graph, however we need to predict the value on \"atom pair\" in this competition.\n\nWeaveNet explicitly deals with both atom feature \"A\" and pair feature \"P\", so it is easy to use this network for baseline modeling.","attachments":{}},{"metadata":{},"cell_type":"markdown","source":"### Data preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_test_df = pd.concat([train_df, test_df], axis=0, sort=True)\ntrain_test_df['coupling_type_id'] = train_test_df['type'].map(\n {'1JHC': 0, '1JHN': 1, '2JHH': 2, '2JHN': 3, '2JHC': 4, '3JHH': 5, '3JHC': 6, '3JHN': 7})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import chainer\n\nfrom chainer_chemistry.dataset.preprocessors.common import construct_atomic_number_array\nfrom chainer_chemistry.dataset.preprocessors.weavenet_preprocessor import construct_atom_type_vec, construct_formal_charge_vec, \\\n construct_partial_charge_vec, construct_atom_ring_vec, construct_hybridization_vec, construct_hydrogen_bonding, \\\n construct_aromaticity_vec, construct_num_hydrogens_vec\nfrom chainer_chemistry.dataset.preprocessors.weavenet_preprocessor import construct_pair_feature\nfrom chainer_chemistry.datasets.numpy_tuple_dataset import NumpyTupleDataset\nfrom chainer_chemistry.links.scaler.standard_scaler import StandardScaler\n\nimport os\nimport numpy\nfrom tqdm import tqdm_notebook as tqdm\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Atom feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_max_atoms = 29\natom_list = ['H', 'C', 'N', 'O', 'F']\ninclude_unknown_atom = False\n\n# Each feature can be calculated using functions provided by chainer-chemistry using RDKit `mol` instance\nconstruct_atomic_number_array(mol)\natom_type_vec = construct_atom_type_vec(\n mol, num_max_atoms, atom_list=atom_list,\n include_unknown_atom=include_unknown_atom)\nformal_charge_vec = construct_formal_charge_vec(mol, num_max_atoms)\npartial_charge_vec = construct_partial_charge_vec(mol, num_max_atoms)\natom_ring_vec = construct_atom_ring_vec(mol, num_max_atoms)\nhybridization_vec = construct_hybridization_vec(mol, num_max_atoms)\nhydrogen_bonding = construct_hydrogen_bonding(mol, num_max_atoms)\naromaticity_vec = construct_aromaticity_vec(mol, num_max_atoms)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Each feature is in shape (atom_id, feature_dim), where atom_id is padded to \"num_max_atoms\" size when atom is not exist.\natom_type_vec.shape, formal_charge_vec.shape, partial_charge_vec.shape, atom_ring_vec.shape, hybridization_vec.shape, hydrogen_bonding.shape, aromaticity_vec.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"WEAVE_DEFAULT_NUM_MAX_ATOMS=20\ndef construct_atom_feature(mol, add_Hs,\n num_max_atoms=WEAVE_DEFAULT_NUM_MAX_ATOMS,\n atom_list=None, include_unknown_atom=False):\n \"\"\"construct atom feature\n\n Args:\n mol (Mol): mol instance\n add_Hs (bool): if the `mol` instance was added Hs, set True.\n num_max_atoms (int): number of max atoms\n atom_list (list): list of atoms to extract feature. If None, default\n `ATOM` is used as `atom_list`\n include_unknown_atom (bool): If False, when the `mol` includes atom\n which is not in `atom_list`, it will raise\n `MolFeatureExtractionError`.\n If True, even the atom is not in `atom_list`, `atom_type` is set\n as \"unknown\" atom.\n\n Returns (numpy.ndarray): 2 dimensional array. First axis size is\n `num_max_atoms`, representing each atom index.\n Second axis for feature.\n\n \"\"\"\n atom_type_vec = construct_atom_type_vec(\n mol, num_max_atoms, atom_list=atom_list,\n include_unknown_atom=include_unknown_atom)\n # TODO(nakago): Chilarity\n formal_charge_vec = construct_formal_charge_vec(\n mol, num_max_atoms=num_max_atoms)\n partial_charge_vec = construct_partial_charge_vec(\n mol, num_max_atoms=num_max_atoms)\n atom_ring_vec = construct_atom_ring_vec(\n mol, num_max_atoms=num_max_atoms)\n hybridization_vec = construct_hybridization_vec(\n mol, num_max_atoms=num_max_atoms)\n hydrogen_bonding = construct_hydrogen_bonding(\n mol, num_max_atoms=num_max_atoms)\n aromaticity_vec = construct_aromaticity_vec(\n mol, num_max_atoms=num_max_atoms)\n if add_Hs:\n num_hydrogens_vec = construct_num_hydrogens_vec(\n mol, num_max_atoms=num_max_atoms)\n feature = numpy.hstack((atom_type_vec, formal_charge_vec,\n partial_charge_vec, atom_ring_vec,\n hybridization_vec, hydrogen_bonding,\n aromaticity_vec, num_hydrogens_vec))\n else:\n feature = numpy.hstack((atom_type_vec, formal_charge_vec,\n partial_charge_vec, atom_ring_vec,\n hybridization_vec, hydrogen_bonding,\n aromaticity_vec))\n return feature","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### pair feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"pair_feature = construct_pair_feature(mol, num_max_atoms=num_max_atoms)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Each feature is in shape (atom_id*atom_id, feature_dim).\npair_feature.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef to_ndarray(feature_list):\n feat_array = np.empty(len(feature_list), dtype=np.ndarray)\n feat_array[:] = feature_list[:]\n return feat_array\n\ndef create_features(mol_names):\n atom_array_list = []\n adj_array_list = []\n to_array_list = []\n from_array_list = []\n label_array_list = []\n coupling_type_id_array_list = []\n\n for mol_name in tqdm(mol_names):\n mol = mol_from_xyz(input_dir/f'structures/{mol_name}.xyz')\n\n atom_array = construct_atom_feature(mol, add_Hs=True,\n num_max_atoms = num_max_atoms, atom_list=atom_list,\n include_unknown_atom=False)\n pair_feature = construct_pair_feature(\n mol, num_max_atoms=num_max_atoms) \n# atom_ids = construct_atomic_number_array(mol)\n# df_mol = structure_df.loc[structure_df['molecule_name'] == mol_name, :]\n# atom_id_array = df_mol['atom_id'].values.astype(np.int32)\n# assert np.alltrue(atom_ids == atom_id_array), print(atom_ids, atom_id_array)\n\n# xyz_array = df_mol[['x', 'y', 'z']].values\n# dist_array = np.linalg.norm(xyz_array[:, None, :] - xyz_array[None, :, :], axis=2).astype(np.float32)\n\n df_mol2 = train_test_df.loc[train_test_df['molecule_name'] == mol_name, :]\n from_array = df_mol2['atom_index_0'].values.astype(np.int32)\n to_array = df_mol2['atom_index_1'].values.astype(np.int32)\n label_array = df_mol2['scalar_coupling_constant'].values.astype(np.float32)\n #label_array = df_mol2['scalar_coupling_constant_scaled'].values.astype(np.float32)\n coupling_type_id_array = df_mol2['coupling_type_id'].values.astype(np.int32)\n\n atom_array_list.append(atom_array)\n adj_array_list.append(pair_feature)\n from_array_list.append(from_array)\n to_array_list.append(to_array)\n coupling_type_id_array_list.append(coupling_type_id_array)\n label_array_list.append(label_array)\n return tuple(map(to_ndarray, [atom_array_list, adj_array_list, from_array_list, to_array_list, coupling_type_id_array_list, label_array_list]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# For debugging\n\nk = create_features(train_mol_names[:2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Calculating all features takes time (around half day), so I already uploaded calculated data as \"dataset\" [here](https://www.kaggle.com/corochann/champs-weave-dataset). We can use this cache for now."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_dataset(filepath):\n if not os.path.exists(filepath):\n return None\n load_data = numpy.load(filepath, allow_pickle=True)\n result = []\n i = 0\n while True:\n key = 'arr_{}'.format(i)\n if key in load_data.keys():\n result.append(load_data[key])\n i += 1\n else:\n break\n return NumpyTupleDataset(*result)\n\nwith timer('train'):\n train = load_dataset('../input/champs-weave-dataset/train_data_v2.npz')\nwith timer('test'):\n test = load_dataset('../input/champs-weave-dataset/test_data_v2.npz')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let me explain each data.\n\n`train[0]` is to get 0-th molecule's data."},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_array, dist_array, from_array, to_array, coupling_type_id_array, label_array = train[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- atom_array: stores feature for each atom, shape (num_atom, feature_dim)\n- dist_array: stores feature for each atom pair, shape (num_atom * num_atom, feature_dim)\n\n- from_array: stores index for atom pair to be predicted\n- to_array : stores index for atom pair to be predicted\n- coupling_type_id_array: stores coupling type for each atom pair to be predicted\n- label_array: stores `scalar_couling_constant` for each atom pair to be predicted"},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_array.shape, dist_array.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from_array, to_array, label_array, coupling_type_id_array","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### converter\n\nEach data (`train[0], train[1]...`) have different length for the label.\n`converter` defines how to make batch from data for neural network training."},{"metadata":{"trusted":true},"cell_type":"code","source":"from chainer.dataset.convert import _concat_arrays, to_device\nfrom itertools import chain\n\ndef concat_with_padding(batch, device, i, padding):\n return to_device(device, _concat_arrays([example[i] for example in batch], padding))\n\n\ndef flatten_array(list_of_list, dtype=np.float32):\n flatten_list = list(chain.from_iterable(list_of_list))\n return np.array(flatten_list, dtype=dtype)\n\nclass EdgeConverter:\n\n def __init__(self, return_label=True):\n self.return_label = return_label\n\n def __call__(self, batch, device=None):\n # atom_array, dist_array, from_array, to_array, coupling_type_id_array, label_array, \n atom = concat_with_padding(batch, device, 0, padding=None)\n adj = concat_with_padding(batch, device, 1, padding=None)\n from_array = to_device(device, flatten_array([example[2] for example in batch], dtype=np.int32))\n to_array = to_device(device,flatten_array([example[3] for example in batch], dtype=np.int32))\n coupling_type = to_device(device, flatten_array([example[4] for example in batch], dtype=np.int32))\n label = to_device(device, flatten_array([example[5] for example in batch], dtype=np.float32))\n batch_indices = to_device(device, flatten_array([np.ones(len(example[2]), dtype=np.int32) * i for i, example in enumerate(batch)], dtype=np.int32))\n\n if self.return_label:\n return atom, adj, from_array, to_array, coupling_type, batch_indices, label[:, None]\n else:\n return atom, adj, from_array, to_array, coupling_type, batch_indices\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"converter = EdgeConverter()\n#print('0', train[0])\n#print('1', train[1])\nbatch = [train[0], train[1]]\n\ninput_data = converter(batch, device=-1)\natom, adj, from_array, to_array, coupling_type, batch_indices, label = input_data\nprint('input_data:', atom.shape, adj.shape, from_array.shape, to_array.shape, coupling_type.shape, batch_indices.shape, label.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`train[0]` and `train[1]` are stacked to create 1-minibatch in the above example."},{"metadata":{},"cell_type":"markdown","source":"### Scaler\n\nEach coupling constant value is different, so I will define scaler for each coupling constant."},{"metadata":{"trusted":true},"cell_type":"code","source":"class CouplingTypeScaler(chainer.Link):\n\n def __init__(self):\n super(CouplingTypeScaler, self).__init__()\n self.mean = None\n self.register_persistent('mean')\n self.std = None\n self.register_persistent('std')\n\n \n def set_param(self, mean, std):\n self.mean = mean\n self.std = std\n \n def inverse_transform(self, x, coupling_indices):\n assert x.ndim == 2, x.shape\n assert x.shape[1] == 1, x.shape\n assert x.shape[0] == len(coupling_indices)\n x = x * self.std[coupling_indices][:, None] + self.mean[coupling_indices][:, None]\n return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# scale each scalar couplint constant\nscaler = StandardScaler()\n\nmean_list = []\nstd_list = []\n\ntypes = ['1JHC', '1JHN', '2JHH', '2JHN', '2JHC', '3JHH', '3JHC', '3JHN']\nfor target_type in types:\n print('scaling type', target_type)\n label = train_df.loc[train_df['type'] == target_type, 'scalar_coupling_constant'].values\n train_df.loc[train_df['type'] == target_type, 'scalar_coupling_constant_scaled'] = scaler.fit_transform(label)\n print('mean', scaler.mean[0], 'std', scaler.std[0])\n mean_list.append(scaler.mean[0])\n std_list.append(scaler.std[0])\n\ncoupling_type_scaler = CouplingTypeScaler()\ncoupling_type_scaler.set_param(\n np.array(mean_list, dtype=np.float32), np.array(std_list, dtype=np.float32))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model definition"},{"metadata":{"trusted":true},"cell_type":"code","source":"import chainer\nfrom chainer import functions\nfrom chainer import links\n\nfrom chainer_chemistry.config import MAX_ATOMIC_NUM\nfrom chainer_chemistry.config import WEAVE_DEFAULT_NUM_MAX_ATOMS\nfrom chainer_chemistry.links.readout.general_readout import GeneralReadout\nfrom chainer_chemistry.links.connection.embed_atom_id import EmbedAtomID\n\nfrom chainer_chemistry.models.weavenet import LinearLayer, WeaveModule\n\nclass WeaveNet(chainer.Chain):\n \"\"\"WeaveNet implementation\n\n Args:\n weave_channels (list): list of int, output dimension for each weave\n module\n hidden_dim (int): hidden dim\n n_atom (int): number of atom of input array\n n_sub_layer (int): number of layer for each `AtomToPair`, `PairToAtom`\n layer\n n_atom_types (int): number of atom id\n readout_mode (str): 'sum' or 'max' or 'summax'\n \"\"\"\n\n def __init__(self, weave_channels=None, hidden_dim=16,\n n_atom=WEAVE_DEFAULT_NUM_MAX_ATOMS,\n n_sub_layer=1, n_atom_types=MAX_ATOMIC_NUM,\n readout_mode='sum'):\n weave_channels = weave_channels or WEAVENET_DEFAULT_WEAVE_CHANNELS\n weave_module = [\n WeaveModule(n_atom, c, n_sub_layer, readout_mode=readout_mode)\n for c in weave_channels\n ]\n\n super(WeaveNet, self).__init__()\n with self.init_scope():\n self.embed = EmbedAtomID(out_size=hidden_dim, in_size=n_atom_types)\n self.weave_module = chainer.ChainList(*weave_module)\n # self.readout = GeneralReadout(mode=readout_mode)\n self.readout_mode = readout_mode\n\n def __call__(self, atom_x, pair_x, train=True):\n if atom_x.dtype == self.xp.int32:\n # atom_array: (minibatch, atom)\n atom_x = self.embed(atom_x)\n\n for i in range(len(self.weave_module)):\n atom_x, pair_x = self.weave_module[i].forward(atom_x, pair_x)\n# if i == len(self.weave_module) - 1:\n# # last layer, only `atom_x` is needed.\n# atom_x = self.weave_module[i].forward(atom_x, pair_x,\n# atom_only=True)\n# else:\n# # not last layer, both `atom_x` and `pair_x` are needed\n# atom_x, pair_x = self.weave_module[i].forward(atom_x, pair_x)\n \n # x = self.readout(atom_x, axis=1)\n # return x\n return pair_x\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import chainer\nimport chainer.links as L\nimport chainer.functions as F\n\nfrom chainer_chemistry.models.mlp import MLP\n\n\ndef calc_pair_feature(h1, h2):\n # Any permutation invariant is okay\n # calc_pair_feature(h1, h2) == calc_pair_feature(h2, h1)\n return F.concat([h1 + h2, h1 * h2], axis=1)\n\n\nclass EdgePredictor(chainer.Chain):\n\n def __init__(self, graph_conv, coupling_type_scaler=None, num_coupling_type=8, edge_embed_dim=16):\n super(EdgePredictor, self).__init__()\n with self.init_scope():\n self.graph_conv = graph_conv\n self.coupling_type_scaler = coupling_type_scaler\n self.embed = L.EmbedID(num_coupling_type, edge_embed_dim)\n self.mlp = MLP(out_dim=1, hidden_dim=16, n_layers=2)\n\n def forward(self, atom, adj, from_array, to_array, coupling_array, batch_indices):\n batchsize, num_node, ch = atom.shape\n # (batch, node*node, ch)\n h = self.graph_conv(atom, adj)\n assert h.shape[:2] == (batchsize, num_node*num_node)\n assert h.ndim == 3\n\n indices = num_node * from_array + to_array\n h = h[batch_indices, indices]\n h_coupling = self.embed(coupling_array)\n\n # (edge, ch)\n h = F.concat([h, h_coupling], axis=1)\n h = self.mlp(h)\n # (edge, 1)\n h = self.coupling_type_scaler.inverse_transform(h, coupling_array)\n return h\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import copy\nfrom logging import getLogger\n\nimport numpy\n\nimport chainer\nfrom chainer import cuda\nfrom chainer.dataset import convert\nfrom chainer import reporter\nfrom chainer.training.extensions import Evaluator\nfrom chainer.training.extensions import Evaluator\n\n\ndef _get_1d_numpy_array(v):\n \"\"\"Convert array or Variable to 1d numpy array\n\n Args:\n v (numpy.ndarray or cupy.ndarray or chainer.Variable): array to be\n converted to 1d numpy array\n\n Returns (numpy.ndarray): Raveled 1d numpy array\n\n \"\"\"\n if isinstance(v, chainer.Variable):\n v = v.data\n return cuda.to_cpu(v).ravel()\n\n\n\nclass LogMAEEvaluator(Evaluator):\n\n def __init__(self, iterator, target, converter=convert.concat_examples,\n device=None, eval_hook=None, eval_func=None, \n name=None, logger=None):\n super(LogMAEEvaluator, self).__init__(\n iterator, target, converter=converter, device=device,\n eval_hook=eval_hook, eval_func=eval_func)\n self.name = name\n self.logger = logger or getLogger()\n\n def evaluate(self):\n iterator = self._iterators['main']\n eval_func = self.eval_func or self._targets['main']\n\n if self.eval_hook:\n self.eval_hook(self)\n\n if hasattr(iterator, 'reset'):\n iterator.reset()\n it = iterator\n else:\n it = copy.copy(iterator)\n\n y_total = []\n t_total = []\n coupling_total = []\n for batch in it:\n in_arrays = self.converter(batch, self.device)\n with chainer.no_backprop_mode(), chainer.using_config('train',\n False):\n y = eval_func(*in_arrays[:-1])\n t = in_arrays[-1]\n coupling = in_arrays[-2]\n y_data = _get_1d_numpy_array(y)\n t_data = _get_1d_numpy_array(t)\n coupling_data = _get_1d_numpy_array(coupling)\n y_total.append(y_data)\n t_total.append(t_data)\n coupling_total.append(coupling_data)\n\n y_total = numpy.concatenate(y_total).ravel()\n t_total = numpy.concatenate(t_total).ravel()\n coupling_total = numpy.concatenate(coupling_total).ravel()\n # metrics_value = self.metrics_fun(y_total, t_total)\n\n metrics_dict = {}\n num_coupling = 8\n num_valid_coupling = 0\n log_mae_mean = 0\n for i in range(num_coupling): # 8 is num_coupling type\n this_coupling = coupling_total == i\n if np.sum(this_coupling) == 0:\n continue\n \n log_mae = np.log(np.mean(np.abs(y_total[this_coupling] - t_total[this_coupling])))\n log_mae_mean += log_mae\n num_valid_coupling += 1\n metrics_dict[f'logmae_{i}'] = log_mae\n #print('log_mae', i, log_mae, np.sum(this_coupling == i))\n log_mae_mean /= num_valid_coupling\n metrics_dict['logmae_mean'] = log_mae_mean\n \n observation = {}\n with reporter.report_scope(observation):\n reporter.report(metrics_dict, self._targets['main'])\n return observation\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train script"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy\nimport os\n\nfrom chainer.datasets import split_dataset_random\nfrom chainer import functions as F\nfrom chainer import iterators\nfrom chainer import optimizers\nfrom chainer import training\nfrom chainer.training import extensions as E\n\nfrom chainer_chemistry.dataset.converters import concat_mols\nfrom chainer_chemistry.dataset.preprocessors import preprocess_method_dict\nfrom chainer_chemistry import datasets as D\nfrom chainer_chemistry.datasets import NumpyTupleDataset\nfrom chainer_chemistry.links.scaler.standard_scaler import StandardScaler\nfrom chainer_chemistry.models.prediction import Regressor\nfrom chainer_chemistry.dataset.splitters.random_splitter import RandomSplitter","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_max_atoms = 29\n\ngraph_conv = WeaveNet(\n weave_channels=[32, 32], hidden_dim=32,\n n_atom=num_max_atoms,\n n_sub_layer=2, n_atom_types=MAX_ATOMIC_NUM)\npredictor = EdgePredictor(graph_conv, coupling_type_scaler)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nbatchsize = 128\n#device = -1 # Use CPU\ndevice = 0 # Use GPU, much faster!\nepoch = 3\nout = 'results'\nseed = 0\n\n\n#splitter = RandomSplitter()\n#train_data, valid_data = splitter.train_valid_split(train, return_index=False, frac_train=0.9, frac_val=0.1)\ntrain_data, valid_data = split_dataset_random(train, int(len(train) * 0.9), seed)\ntrain_iter = iterators.SerialIterator(train_data, batchsize)\nvalid_iter = iterators.SerialIterator(valid_data, batchsize, repeat=False,\n shuffle=False)\n\nregressor = Regressor(predictor, lossfun=F.mean_squared_error,\n #metrics_fun=metrics_fun,\n device=device)\n# Set up the optimizer.\noptimizer = optimizers.Adam()\noptimizer.setup(regressor)\nconverter.return_label = True\n\n# Set up the updater.\nupdater = training.StandardUpdater(train_iter, optimizer, device=device,\n converter=converter)\n\n# Set up the trainer.\ntrainer = training.Trainer(updater, (epoch, 'epoch'), out=out)\n# trainer.extend(E.Evaluator(valid_iter, regressor, device=device,\n# converter=converter))\ntrainer.extend(LogMAEEvaluator(valid_iter, regressor, device=device,\n converter=converter, name='val', eval_func=predictor))\ntrainer.extend(E.snapshot(), trigger=(epoch, 'epoch'))\ntrainer.extend(E.LogReport(trigger=(1, 'epoch'))) # (100, 'iteration') can be used for frequent logging.\ntry:\n # https://github.com/grafi-tt/chaineripy\n # Override PrintReport & ProgressBar to use chaineripy on jupyter notebook\n from chaineripy.extensions import PrintReport, ProgressBar\n trainer.extend(PrintReport([\n 'epoch', 'main/loss', 'validation/main/loss',\n 'val/main/logmae_mean', 'val/main/logmae_0', 'val/main/logmae_1', 'val/main/logmae_2', 'val/main/logmae_3', 'val/main/logmae_4',\n 'val/main/logmae_5', 'val/main/logmae_6', 'val/main/logmae_7',\n 'elapsed_time']))\n trainer.extend(ProgressBar())\nexcept ImportError:\n trainer.extend(E.PrintReport([\n 'epoch', 'main/loss', 'validation/main/loss',\n 'val/main/logmae_mean', 'val/main/logmae_0', 'val/main/logmae_1', 'val/main/logmae_2', 'val/main/logmae_3', 'val/main/logmae_4',\n 'val/main/logmae_5', 'val/main/logmae_6', 'val/main/logmae_7',\n 'elapsed_time']))\n trainer.extend(E.ProgressBar())\n\n\ntrainer.run()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Prediction for test label\n\nconverter.return_label = False\ntest_pred = regressor.predict(test, batchsize=128, converter=converter)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"test_pred.shape, test_pred.dtype","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check test prediction distribution\n\nsns.distplot(test_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.alltrue(test_df.id.values == sample_submission_df.id.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Assign test prediction to submit format.\nsample_submission_df['scalar_coupling_constant'] = test_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Save to csv file and we are ready to submit!\n\nsample_submission_df.to_csv('submission.csv', index=False)\nprint('saved to submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n\n# Going further\n\nAbove example shows the baseline modeling using graph convolution network. Below is the topics we can consider further."},{"metadata":{},"cell_type":"markdown","source":"#### CV strategy\n\nI splitted dataset just randomly in above example. But is this enough for molecular data?\n\nWhen molecular structure differs, its prediction may be more difficult. And **scaffold split** is sometimes used to separate data by molecular \"scaffold\". [Scaffold splitting](https://github.com/pfnet-research/chainer-chemistry/blob/master/chainer_chemistry/dataset/splitters/scaffold_splitter.py) is provided in chainer-chemistry, and we may use this for more conservative validation. \n\n#### Data preprocessing\n\nMany kinds of feature extraction is provided in RDKit modules, and using these feature may help improving your model.\nReference\n\n - [RDKit get started](https://www.rdkit.org/docs/GettingStartedInPython.html)\n - [preprocessor class](https://github.com/pfnet-research/chainer-chemistry/tree/master/chainer_chemistry/dataset/preprocessors) in chainer-chemistry for many kinds of atom feature extraction with RDKit.\n\n\n**Deal with each coupling type separately**\n\nIt seems value is very different. you may consider training different models etc.\n\n\n#### Models\n\nYou may consider using other graph convolutional models, as also discussed in this thread.\n\n - [Which graph CNN is the best?](https://www.kaggle.com/c/champs-scalar-coupling/discussion/93972#latest-546673)\n\nMany models are available in chainer-chemistry as well. Especially `SchNet` deals atom position explicitly, which may be useful. \n \n#### Understanding QM9 Dataset\n\nRules are updated that we can use QM9 Dataset as external dataset.\n\n - [IMPORTANT: Rules Update - Please Read](https://www.kaggle.com/c/champs-scalar-coupling/discussion/95496#latest-554642)\n\nSo understanding original QM9 Dataset is quite important.\n\n[Training example of QM9](https://github.com/pfnet-research/chainer-chemistry/tree/master/examples/qm9) is provided in chainer-chemistry.\nQM9 dataset itself is downloaded by the library as well.\nYou may consider using pre-trained model with this example as well.\n\n - [QM9 dataset explanation notebook](https://github.com/pfnet-research/chainer-chemistry/blob/master/examples/qm9/qm9_dataset_exploration.ipynb)"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}