{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Intro_before_Intro\n\nThis kernel helped me win Silver and 57th place! It is a slight variation of [SchNet Starter Kit](https://www.kaggle.com/toshik/schnet-starter-kit) by @toshik.\n\n[What else wento into #57 blend is summarized here](https://www.kaggle.com/c/champs-scalar-coupling/discussion/106359#)\n\nOne must run 8 kernels such as this one separately for each coupling type and then combine everything into a single submission.\nI expect LB ~ -1.97 with the parameters that I used (see below). Since this variation is actually ~20% faster version of what I actually used (due to testing only on a single type), one could increase the #epochs listed below by 20-30% and possibly get -2.0 to -2.1 score (just guessing here, I don't want to waste ~50hrs of GPU time to find out) while stil being within GPU kernel limit.\n\nbatch_size and # epochs used (**by no means these are optimal**):\n\n* 1JHN 9 100\n* 2JHN 13 80\n* 3JHN 14 80\n* 1JHC 16 65\n* 2JHC 14 60\n* 3JHC 16 65\n* 2JHH 16 65\n* 3JHH 14 60"},{"metadata":{},"cell_type":"markdown","source":"# Introduction\nI share a simple example which employ SchNet for predicting coupling constants.\nI hope this kernel helps beginners of DNN and will be used as a starter kit.\n\nThe core idea is same as Heng's, which employ GNN as a feature exstractor.\nTwo feature vectors are concatenated and thrown into regression header.\nMore details of his idea and discussions can be read in below pages.\n\n* Which graph CNN is the best (with starter kit at LB -1.469)?<br>https://www.kaggle.com/c/champs-scalar-coupling/discussion/93972#latest-591759\n\nDue to the limitation of Kaggle kernel, the model in not trained completely in this Kernel.\nYou would achieve better score by continuing training procedure longer."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"},"cell_type":"markdown","source":"# Install packages\nIn this example, I use chainer chemistry which offer an implementation of SchNet.\nThis library can be install by PIP.\n* Chainer Chemistry: A Library for Deep Learning in Biology and Chemistry<br>https://github.com/pfnet-research/chainer-chemistry"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"!pip uninstall -y tensorflow\n!pip install chainer-chemistry==0.5.0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Import packages\nNext, I import main packages. Other sub-modules are imported later."},{"metadata":{"trusted":false},"cell_type":"code","source":"import random\nimport numpy as np\nimport pandas as pd\nimport chainer\nimport chainer_chemistry\nfrom IPython.display import display","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load dataset\nIn this example, 90% of training data is used actual training data, and the other 10% is used for validation.\nEach dataset is grouped by molecule_name name for following procedures."},{"metadata":{"trusted":false},"cell_type":"code","source":"TYPE_VAL = '1JHN'  # change to the required type, must run 8 kernels to cover the train set.\ndef load_dataset():\n\n    train = pd.merge(pd.read_csv('../input/champs-scalar-coupling/train.csv'),\n                     pd.read_csv('../input/champs-scalar-coupling/scalar_coupling_contributions.csv'))\n\n    test = pd.read_csv('../input/champs-scalar-coupling/test.csv')\n    train = train[train['type'] == TYPE_VAL].copy()\n    train = train.reset_index(drop=True).copy()\n    test = test[test['type'] == TYPE_VAL].copy()\n    test = test.reset_index(drop=True).copy()\n    print(train.shape)\n    print(test.shape)\n\n    counts = train['molecule_name'].value_counts()\n    moles = list(counts.index)\n    random.seed(999)\n    random.shuffle(moles)\n\n    num_train = int(len(moles) * 0.999)\n    train_moles = sorted(moles[:num_train])\n    valid_moles = sorted(moles[num_train:])\n    test_moles = sorted(list(set(test['molecule_name'])))\n    \n    valid = train.query('molecule_name not in @train_moles')\n    train = train.query('molecule_name in @train_moles')\n\n    train.sort_values('molecule_name', inplace=True)\n    valid.sort_values('molecule_name', inplace=True)\n    test.sort_values('molecule_name', inplace=True)\n\n    return train, valid, test, train_moles, valid_moles, test_moles\n\ntrain, valid, test, train_moles, valid_moles, test_moles = load_dataset()\n\ntrain_gp = train.groupby('molecule_name')\nvalid_gp = valid.groupby('molecule_name')\ntest_gp = test.groupby('molecule_name')\n\nstructures = pd.read_csv('../input/champs-scalar-coupling/structures.csv')\nstructures_groups = structures.groupby('molecule_name')\nprint(structures.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## train data"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(train.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## validation data"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(valid.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## test data"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(test.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## structures"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(structures.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing\nI implemented a class named `Graph` whose instances contain molecules.\nThe distances between atoms are calculated in the initializer of this class.\n## Define Graph class"},{"metadata":{"trusted":false},"cell_type":"code","source":"from scipy.spatial import distance\n\n\nclass Graph:\n\n    def __init__(self, points_df, list_atoms):\n\n        self.points = points_df[['x', 'y', 'z']].values\n\n        self._dists = distance.cdist(self.points, self.points)\n\n        self.adj = self._dists < 1.5\n        self.num_nodes = len(points_df)\n\n        self.atoms = points_df['atom']\n        dict_atoms = {at: i for i, at in enumerate(list_atoms)}\n\n        atom_index = [dict_atoms[atom] for atom in self.atoms]\n        one_hot = np.identity(len(dict_atoms))[atom_index]\n\n        bond = np.sum(self.adj, 1) - 1\n        bonds = np.identity(len(dict_atoms))[bond - 1]\n\n        self._array = np.concatenate([one_hot, bonds], axis=1).astype(np.float32)\n\n    @property\n    def input_array(self):\n        return self._array\n\n    @property\n    def dists(self):\n        return self._dists.astype(np.float32)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Convert into graph object\nEach dataset is represented as a list of Graphs and prediction targets."},{"metadata":{"trusted":false},"cell_type":"code","source":"list_atoms = list(set(structures['atom']))\nprint('list of atoms')\nprint(list_atoms)\n    \ntrain_graphs = list()\ntrain_targets = list()\nprint('preprocess training molecules ...')\nfor mole in train_moles:\n    train_graphs.append(Graph(structures_groups.get_group(mole), list_atoms))\n    train_targets.append(train_gp.get_group(mole))\n\nvalid_graphs = list()\nvalid_targets = list()\nprint('preprocess validation molecules ...')\nfor mole in valid_moles:\n    valid_graphs.append(Graph(structures_groups.get_group(mole), list_atoms))\n    valid_targets.append(valid_gp.get_group(mole))\n\ntest_graphs = list()\ntest_targets = list()\nprint('preprocess test molecules ...')\nfor mole in test_moles:\n    test_graphs.append(Graph(structures_groups.get_group(mole), list_atoms))\n    test_targets.append(test_gp.get_group(mole))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Convert into chainer's dataset\nThis type of dataset can be handled by `DictDataset`.\nGraph objects and prediction targets are merged as a `DictDataset`."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer.datasets.dict_dataset import DictDataset\n\ntrain_dataset = DictDataset(graphs=train_graphs, targets=train_targets)\nvalid_dataset = DictDataset(graphs=valid_graphs, targets=valid_targets)\ntest_dataset = DictDataset(graphs=test_graphs, targets=test_targets)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model\n## Build SchNet model\nThe prediction model is implemented as follows.\nFirst, fully connected layer is applied to input arrays to align dimensions.\nNext, SchNet layer is applied for feature extraction.\nFinally, features vectors are concatenated and thrown into three layers MLP.\nI add batch-normalization layers like ResNet."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer import reporter\nfrom chainer import functions as F\nfrom chainer import links as L\nfrom chainer_chemistry.links import SchNetUpdate\nfrom chainer_chemistry.links import GraphLinear, GraphBatchNormalization\n\nclass SchNetUpdateBN(SchNetUpdate):\n\n    def __init__(self, *args, **kwargs):\n        super(SchNetUpdateBN, self).__init__(*args, **kwargs)\n        with self.init_scope():\n            self.bn = GraphBatchNormalization(args[0])\n\n    def __call__(self, h, adj, **kwargs):\n        v = self.linear[0](h)\n        v = self.cfconv(v, adj)\n        v = self.linear[1](v)\n        v = F.softplus(v)\n        v = self.linear[2](v)\n        return h + self.bn(v)\n\nclass SchNet(chainer.Chain):\n\n    def __init__(self, num_layer=3):\n        super(SchNet, self).__init__()\n\n        self.num_layer = num_layer\n\n        with self.init_scope():\n            self.gn = GraphLinear(512)\n            for l in range(self.num_layer):\n                self.add_link('sch{}'.format(l), SchNetUpdateBN(512))\n\n            self.interaction1 = L.Linear(128)\n            self.interaction2 = L.Linear(128)\n            self.interaction3 = L.Linear(4)\n\n    def __call__(self, input_array, dists, pairs_index, targets):\n\n        out = self.predict(input_array, dists, pairs_index)\n        loss = F.mean_absolute_error(out, targets)\n        reporter.report({'loss': loss}, self)\n        return loss\n\n    def predict(self, input_array, dists, pairs_index, **kwargs):\n\n        h = self.gn(input_array)\n\n        for l in range(self.num_layer):\n            h = self['sch{}'.format(l)](h, dists)\n\n        h = F.concat((h, input_array), axis=2)\n\n        concat = F.concat([\n            h[pairs_index[:, 0], pairs_index[:, 1], :],\n            h[pairs_index[:, 0], pairs_index[:, 2], :],\n            F.expand_dims(dists[pairs_index[:, 0],\n                                pairs_index[:, 1],\n                                pairs_index[:, 2]], 1)\n        ], axis=1)\n\n        h1 = F.leaky_relu(self.interaction1(concat))\n        h2 = F.leaky_relu(self.interaction2(h1))\n        out = self.interaction3(h2)\n\n        return out\n\nmodel = SchNet(num_layer=3)\nmodel.to_gpu(device=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training preparation\n## Make samplers\nFor mini-batch training, I implement a sampler named `SameSizeSampler`.\nThe molecules which have same number of atoms are selected simultaneously."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer.iterators import OrderSampler\n\nclass SameSizeSampler(OrderSampler):\n\n    def __init__(self, structures_groups, moles, batch_size,\n                 random_state=None, use_remainder=False):\n\n        self.structures_groups = structures_groups\n        self.moles = moles\n        self.batch_size = batch_size\n        if random_state is None:\n            random_state = np.random.random.__self__\n        self._random = random_state\n        self.use_remainder = use_remainder\n\n    def __call__(self, current_order, current_position):\n\n        batches = list()\n\n        atom_counts = pd.DataFrame()\n        atom_counts['mol_index'] = np.arange(len(self.moles))\n        atom_counts['molecular_name'] = self.moles\n        atom_counts['num_atom'] = [len(self.structures_groups.get_group(mol))\n                                   for mol in self.moles]\n\n        num_atom_counts = atom_counts['num_atom'].value_counts()\n\n        for count, num_mol in num_atom_counts.to_dict().items():\n            if self.use_remainder:\n                num_batch_for_this = -(-num_mol // self.batch_size)\n            else:\n                num_batch_for_this = num_mol // self.batch_size\n\n            target_mols = atom_counts.query('num_atom==@count')['mol_index'].values\n            random.shuffle(target_mols)\n\n            devider = np.arange(0, len(target_mols), self.batch_size)\n            devider = np.append(devider, 99999)\n\n            if self.use_remainder:\n                target_mols = np.append(\n                    target_mols,\n                    np.repeat(target_mols[-1], -len(target_mols) % self.batch_size))\n\n            for b in range(num_batch_for_this):\n                batches.append(target_mols[devider[b]:devider[b + 1]])\n\n        random.shuffle(batches)\n        batches = np.concatenate(batches).astype(np.int32)\n\n        return batches\n\nbatch_size = 9  #change batch size according to coupling type per the table at the top or to your own value\ntrain_sampler = SameSizeSampler(structures_groups, train_moles, batch_size)\nvalid_sampler = SameSizeSampler(structures_groups, valid_moles, batch_size,\n                                use_remainder=True)\ntest_sampler = SameSizeSampler(structures_groups, test_moles, batch_size,\n                               use_remainder=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Make iterators, oprimizer\nIterators for data feeding is made as below."},{"metadata":{"trusted":false},"cell_type":"code","source":"train_iter = chainer.iterators.SerialIterator(\n    train_dataset, batch_size, order_sampler=train_sampler)\n\nvalid_iter = chainer.iterators.SerialIterator(\n    valid_dataset, batch_size, repeat=False, order_sampler=valid_sampler)\n\ntest_iter = chainer.iterators.SerialIterator(\n    test_dataset, batch_size, repeat=False, order_sampler=test_sampler)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Make optimizer\nAdam is used as an optimizer."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer import optimizers\noptimizer = optimizers.Adam(alpha=1e-3)\noptimizer.setup(model)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Make updator\nSince the model receives input arrays separately, I implement an original converter.\n`input_array` and `dists` are exstracted from `Graph` object and `pair_index` and `targets` are exstracted from `targets` object.\n`targets` is added only for training.\nWhen this converter is used for evaluation, `targets` is not added."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer import training\nfrom chainer.dataset import to_device\n\ndef coupling_converter(batch, device):\n\n    list_array = list()\n    list_dists = list()\n    list_targets = list()\n    list_pairs_index = list()\n\n    with_target = 'fc' in batch[0]['targets'].columns\n\n    for i, d in enumerate(batch):\n        list_array.append(d['graphs'].input_array)\n        list_dists.append(d['graphs'].dists)\n        if with_target:\n            list_targets.append(\n                d['targets'][['fc', 'sd', 'pso', 'dso']].values.astype(np.float32))\n\n        sample_index = np.full((len(d['targets']), 1), i)\n        atom_index = d['targets'][['atom_index_0', 'atom_index_1']].values\n\n        list_pairs_index.append(np.concatenate([sample_index, atom_index], axis=1))\n\n    input_array = to_device(device, np.stack(list_array))\n    dists = to_device(device, np.stack(list_dists))\n    pairs_index = np.concatenate(list_pairs_index)\n\n    array = {'input_array': input_array, 'dists': dists, 'pairs_index': pairs_index}\n\n    if with_target:\n        array['targets'] = to_device(device, np.concatenate(list_targets))\n\n    return array\n\nupdater = training.StandardUpdater(train_iter, optimizer,\n                                   converter=coupling_converter, device=0)\ntrainer = training.Trainer(updater, (1, 'epoch'), out=\"result\")  #change #epochs according to coupling type per the table at the top or to your own va","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training extensions\n## Evaluator\nI implemented an Evaluator which measure validation score during training.\nThe prediction for test data is also calculated in this evaluator and the submision file is generated."},{"metadata":{"trusted":false},"cell_type":"code","source":"from chainer.training.extensions import Evaluator\nfrom chainer import cuda\n\nclass TypeWiseEvaluator(Evaluator):\n\n    def __init__(self, iterator, target, converter, device, name,\n                 is_validate=False, is_submit=False):\n\n        super(TypeWiseEvaluator, self).__init__(\n            iterator, target, converter=converter, device=device)\n\n        self.is_validate = is_validate\n        self.is_submit = is_submit\n        self.name = name\n\n    def calc_score(self, df_truth, pred):\n\n        target_types = list(set(df_truth['type']))\n\n        diff = df_truth['scalar_coupling_constant'] - pred\n\n        scores = 0\n        metrics = {}\n\n        for target_type in target_types:\n\n            target_pair = df_truth['type'] == target_type\n            score_exp = np.mean(np.abs(diff[target_pair]))\n            scores += np.log(score_exp)\n\n            metrics[target_type] = scores\n\n        metrics['ALL_LogMAE'] = scores / len(target_types)\n\n        observation = {}\n        with reporter.report_scope(observation):\n            reporter.report(metrics, self._targets['main'])\n\n        return observation\n\n    def evaluate(self):\n        iterator = self._iterators['main']\n        eval_func = self._targets['main']\n\n        iterator.reset()\n        it = iterator\n\n        y_total = []\n        t_total = []\n\n        for batch in it:\n            in_arrays = self.converter(batch, self.device)\n            with chainer.no_backprop_mode(), chainer.using_config('train', False):\n                y = eval_func.predict(**in_arrays)\n\n            y_data = cuda.to_cpu(y.data)\n            y_total.append(y_data)\n            t_total.extend([d['targets'] for d in batch])\n\n        df_truth = pd.concat(t_total, axis=0)\n        y_pred = np.sum(np.concatenate(y_total), axis=1)\n\n        if self.is_submit:\n            submit = pd.DataFrame()\n            submit['id'] = df_truth['id']\n            submit['scalar_coupling_constant'] = y_pred\n            submit.drop_duplicates(subset='id', inplace=True)\n            submit.sort_values('id', inplace=True)\n            submit.to_csv(f'kernel_schnet_{TYPE_VAL}.csv', index=False)\n\n        if self.is_validate:\n            return self.calc_score(df_truth, y_pred)\n\n        return {}\n\ntrainer.extend(\n    TypeWiseEvaluator(iterator=valid_iter, target=model, converter=coupling_converter, \n                      name='valid', device=0, is_validate=False))  ##Not checking validation speeds things up. Change to True if you need validation score\ntrainer.extend(\n    TypeWiseEvaluator(iterator=test_iter, target=model, converter=coupling_converter,\n                      name='test', device=0, is_submit=True))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Other extensions\nExponentialShift is set as a learning rate scheduler.\nAn extension which turn off training mode is also set to deactivate normalizatoin from second epoch.\n\nLog options are set to report the metrics.\nThis helps us to analyze the result of training."},{"metadata":{"trusted":false},"cell_type":"code","source":"trainer.extend(training.extensions.ExponentialShift('alpha', 0.99999))\n\nfrom chainer.training import make_extension\n\ndef stop_train_mode(trigger):\n    @make_extension(trigger=trigger)\n    def _stop_train_mode(_):\n        chainer.config.train = False\n    return _stop_train_mode\n\ntrainer.extend(stop_train_mode(trigger=(1, 'epoch')))\n\ntrainer.extend(\n    training.extensions.observe_value(\n        'alpha', lambda tr: tr.updater.get_optimizer('main').alpha))\n\ntrainer.extend(training.extensions.LogReport())\ntrainer.extend(training.extensions.PrintReport(\n    ['epoch', 'elapsed_time', 'main/loss', 'valid/main/ALL_LogMAE', 'alpha']))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training\n## Run\nI tuned number of epochs to prevent timeout.\nSchNet tends to be underfitting, longer training makes the model better basically."},{"metadata":{"trusted":false},"cell_type":"code","source":"chainer.config.train = True\ntrainer.run()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Check output"},{"metadata":{"trusted":false},"cell_type":"code","source":"submit = pd.read_csv(f'kernel_schnet_{TYPE_VAL}.csv')\ndisplay(submit.head())\nprint('shape: {}'.format(submit.shape))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# For more improvement\nThis example can be improved by below ways.\n* Train more\n* Tune hyperparameters\n* Add original feature\n* Try different GNNs\n* Blend with Gradient Boosting Machines\n\nYou can start with this kernel and don't forget upvote :)\n\n# References\n\n* SchNet: A continuous-filter convolutional neural network for modeling quantum interactions<br>https://arxiv.org/abs/1706.08566\n* SchNet - a deep learning architecture for molecules and materials<br>https://arxiv.org/abs/1712.06113\n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.6"}},"nbformat":4,"nbformat_minor":1}