GCN_chemo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://docs.dgl.ai/tutorials/basics/4_batch.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from collections import namedtuple\n",
    "import dgl\n",
    "from dgl import DGLGraph\n",
    "import dgl.function as fn\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch.utils.data import Dataset, TensorDataset\n",
    "from torch.utils.data import DataLoader\n",
    "import torch.optim as optim\n",
    "\n",
    "import networkx as nx\n",
    "import copy\n",
    "import os\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import RDConfig\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ELEM_LIST = ['C', 'N', 'O', 'S', 'F', 'Si', 'P', 'Cl', 'Br', 'Mg', 'Na', 'Ca', 'Fe', 'Al', 'I', 'B', 'K', 'Se', 'Zn', 'H', 'Cu', 'Mn', 'unknown']\n",
    "\n",
    "ATOM_FDIM = len(ELEM_LIST) + 6 + 5 + 1\n",
    "MAX_ATOMNUM =60\n",
    "BOND_FDIM = 5 \n",
    "MAX_NB = 10\n",
    "\n",
    "PAPER = os.getenv('PAPER', False)\n",
    "\n",
    "def onek_encoding_unk(x, allowable_set):\n",
    "    if x not in allowable_set:\n",
    "        x = allowable_set[-1]\n",
    "    return [x == s for s in allowable_set]\n",
    "\n",
    "# Note that during graph decoding they don't predict stereochemistry-related\n",
    "# characteristics (i.e. Chiral Atoms, E-Z, Cis-Trans).  Instead, they decode\n",
    "# the 2-D graph first, then enumerate all possible 3-D forms and find the\n",
    "# one with highest score.\n",
    "'''\n",
    "def atom_features(atom):\n",
    "    return (torch.Tensor(onek_encoding_unk(atom.GetSymbol(), ELEM_LIST)\n",
    "            + onek_encoding_unk(atom.GetDegree(), [0,1,2,3,4,5])\n",
    "            + onek_encoding_unk(atom.GetFormalCharge(), [-1,-2,1,2,0])\n",
    "            + [atom.GetIsAromatic()]))\n",
    "'''\n",
    "def atom_features(atom):\n",
    "    return (onek_encoding_unk(atom.GetSymbol(), ELEM_LIST)\n",
    "            + onek_encoding_unk(atom.GetDegree(), [0,1,2,3,4,5])\n",
    "            + onek_encoding_unk(atom.GetFormalCharge(), [-1,-2,1,2,0])\n",
    "            + [atom.GetIsAromatic()])\n",
    "\n",
    "def bond_features(bond):\n",
    "    bt = bond.GetBondType()\n",
    "    return (torch.Tensor([bt == Chem.rdchem.BondType.SINGLE, bt == Chem.rdchem.BondType.DOUBLE, bt == Chem.rdchem.BondType.TRIPLE, bt == Chem.rdchem.BondType.AROMATIC, bond.IsInRing()]))\n",
    "\n",
    "def mol2dgl_single(mols):\n",
    "    cand_graphs = []\n",
    "    n_nodes = 0\n",
    "    n_edges = 0\n",
    "    bond_x = []\n",
    "\n",
    "    for mol in mols:\n",
    "        n_atoms = mol.GetNumAtoms()\n",
    "        n_bonds = mol.GetNumBonds()\n",
    "        g = DGLGraph()        \n",
    "        nodeF = []\n",
    "        for i, atom in enumerate(mol.GetAtoms()):\n",
    "            assert i == atom.GetIdx()\n",
    "            nodeF.append(atom_features(atom))\n",
    "        g.add_nodes(n_atoms)\n",
    "\n",
    "        bond_src = []\n",
    "        bond_dst = []\n",
    "        for i, bond in enumerate(mol.GetBonds()):\n",
    "            a1 = bond.GetBeginAtom()\n",
    "            a2 = bond.GetEndAtom()\n",
    "            begin_idx = a1.GetIdx()\n",
    "            end_idx = a2.GetIdx()\n",
    "            features = bond_features(bond)\n",
    "\n",
    "            bond_src.append(begin_idx)\n",
    "            bond_dst.append(end_idx)\n",
    "            bond_x.append(features)\n",
    "            bond_src.append(end_idx)\n",
    "            bond_dst.append(begin_idx)\n",
    "            bond_x.append(features)\n",
    "        g.add_edges(bond_src, bond_dst)\n",
    "        g.ndata['h'] = torch.Tensor(nodeF)\n",
    "        cand_graphs.append(g)\n",
    "    return cand_graphs\n",
    "msg = fn.copy_src(src=\"h\", out=\"m\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "solcls = {'(A) low':0, '(B) medium':1, '(C) high':2}\n",
    "train_mols = [m for m in Chem.SDMolSupplier(os.path.join(RDConfig.RDDocsDir,'Book/data/solubility.train.sdf'))]\n",
    "train_y = [solcls[m.GetProp('SOL_classification')] for m in train_mols]\n",
    "test_mols = [m for m in Chem.SDMolSupplier(os.path.join(RDConfig.RDDocsDir,'Book/data/solubility.test.sdf'))]\n",
    "test_y = [solcls[m.GetProp('SOL_classification')] for m in test_mols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_graphs = mol2dgl_single(train_mols)\n",
    "test_graphs = mol2dgl_single(test_mols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def collate(sample):\n",
    "    graphs, labels = map(list,zip(*sample))\n",
    "    batched_graph = dgl.batch(graphs)\n",
    "    return batched_graph, torch.tensor(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reduce(nodes):\n",
    "    # summazation by avarage is different part\n",
    "    accum = torch.mean(nodes.mailbox['m'], 1)\n",
    "    return {'h': accum}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NodeApplyModule(nn.Module):\n",
    "    def __init__(self, in_feats, out_feats, activation):\n",
    "        super(NodeApplyModule, self).__init__()\n",
    "        self.linear = nn.Linear(in_feats, out_feats)\n",
    "        self.activation = activation\n",
    "    \n",
    "    def forward(self, node):\n",
    "        h = self.linear(node.data['h'])\n",
    "        h = self.activation(h)\n",
    "        return {'h': h}\n",
    "    \n",
    "\n",
    "class GCN(nn.Module):\n",
    "    def __init__(self, in_feats, out_feats, activation):\n",
    "        super(GCN, self).__init__()\n",
    "        self.apply_mod = NodeApplyModule(in_feats, out_feats, activation)\n",
    "    \n",
    "    def forward(self, g, feature):\n",
    "        g.ndata['h'] = feature\n",
    "        g.update_all(msg, reduce)\n",
    "        g.apply_nodes(func=self.apply_mod)\n",
    "        h =  g.ndata.pop('h')\n",
    "        #print(h.shape)\n",
    "        return h\n",
    "    \n",
    "class Classifier(nn.Module):\n",
    "    def __init__(self, in_dim, hidden_dim, n_classes):\n",
    "        super(Classifier, self).__init__()\n",
    "        self.layers = nn.ModuleList([GCN(in_dim, hidden_dim, F.relu),\n",
    "                                    GCN(hidden_dim, hidden_dim, F.relu)])\n",
    "        self.classify = nn.Linear(hidden_dim, n_classes)\n",
    "    def forward(self, g):\n",
    "        h = g.ndata['h']\n",
    "        for conv in self.layers:\n",
    "            h = conv(g, h)\n",
    "        g.ndata['h'] = h\n",
    "        hg = dgl.mean_nodes(g, 'h')\n",
    "        return self.classify(hg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = list(zip(train_graphs, train_y))\n",
    "data_loader = DataLoader(dataset, batch_size=32, shuffle=True, collate_fn=collate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Classifier(ATOM_FDIM, 256, len(solcls))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Classifier(\n",
       "  (layers): ModuleList(\n",
       "    (0): GCN(\n",
       "      (apply_mod): NodeApplyModule(\n",
       "        (linear): Linear(in_features=35, out_features=256, bias=True)\n",
       "      )\n",
       "    )\n",
       "    (1): GCN(\n",
       "      (apply_mod): NodeApplyModule(\n",
       "        (linear): Linear(in_features=256, out_features=256, bias=True)\n",
       "      )\n",
       "    )\n",
       "  )\n",
       "  (classify): Linear(in_features=256, out_features=3, bias=True)\n",
       ")"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loss_func = nn.CrossEntropyLoss()\n",
    "optimizer = optim.Adam(model.parameters(), lr=0.001)\n",
    "model.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 20, loss 0.6104\n",
      "Epoch 40, loss 0.5616\n",
      "Epoch 60, loss 0.5348\n",
      "Epoch 80, loss 0.5095\n",
      "Epoch 100, loss 0.4915\n",
      "Epoch 120, loss 0.5163\n",
      "Epoch 140, loss 0.5348\n",
      "Epoch 160, loss 0.4385\n",
      "Epoch 180, loss 0.4421\n",
      "Epoch 200, loss 0.4318\n"
     ]
    }
   ],
   "source": [
    "epoch_losses = []\n",
    "for epoch in range(200):\n",
    "    epoch_loss = 0\n",
    "    for i, (bg, label) in enumerate(data_loader):\n",
    "        bg.set_e_initializer(dgl.init.zero_initializer)\n",
    "        bg.set_n_initializer(dgl.init.zero_initializer)        \n",
    "        pred = model(bg)\n",
    "        loss = loss_func(pred, label)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        epoch_loss += loss.detach().item()\n",
    "    epoch_loss /= (i + 1)\n",
    "    if (epoch+1) % 20 == 0:\n",
    "        print('Epoch {}, loss {:.4f}'.format(epoch+1, epoch_loss))\n",
    "    epoch_losses.append(epoch_loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a194ae940>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXl4VdXV/78rYYYQphDmMCqgImDEWbEiIkVQbCvUSh1aqtZWrbbCz7c4tLbVvp201LdoreirIFoUXgvFEQcqSkRAkXkOYBgSEQETAvv3x7q759yTc4ckdzz5fp4nz77n3J1zVs69+Z511l57bTHGgBBCSLDISbcBhBBCEg/FnRBCAgjFnRBCAgjFnRBCAgjFnRBCAgjFnRBCAgjFnRBCAgjFnRBCAgjFnRBCAkijdJ24Q4cOpmfPnuk6PSGEZCUffvjhPmNMQax+aRP3nj17oqSkJF2nJ4SQrEREtsXTj2EZQggJIBR3QggJIBR3QggJIBR3QggJIBR3QggJIBR3QggJIBR3QggJIFkn7u++C0ydCnB1QEIIiUzWiXtJCfCb3wAVFem2hBBCMpesE/fCQm3LytJrByGEZDIUd0IICSAUd0IICSAUd0IICSBZJ+7t2gG5uRR3QgiJRtaJe04O0LEjxZ0QQqKRdeIOaGiG4k4IIZGhuBNCSAChuBNCSACJKe4i8oSI7BGRTyK8LyLysIhsFJFVIjI08WaGY8WdJQgIIcSfeDz3JwGMivL+pQD6hX4mA3i0/mZFp7AQqKwEDhxI9pkIISQ7iSnuxpi3AZRH6TIOwFNGWQqgjYh0TpSBfjDXnRBCopOImHtXADtc26WhfUmD4k4IIdFJhLiLzz7faLiITBaREhEp2bt3b51PSHEnhJDoJELcSwF0d213A7DLr6MxZoYxptgYU1xQUFDnE1LcCSEkOokQ9/kAJoWyZs4EcMAYszsBx41Ihw46U5XiTggh/jSK1UFEZgEYDqCDiJQCuAdAYwAwxvwPgAUARgPYCOAwgOuSZawlNxdo3x6oR2SHEEICTUxxN8ZMjPG+AfDDhFkUJ/n5TIUkhJBIZOUMVQBo04biTgghkchacc/PBz7/PN1WEEJIZpK14k7PnRBCIpO14k7PnRBCIpO14k7PnRBCIpO14p6fDxw6BFRXp9sSQgjJPLJW3Nu00ZbeOyGE1CRrxT0/X1uKOyGE1CRrxd167hxUJYSQmmStuNNzJ4SQyGStuNNzJ4SQyGStuNNzJ4SQyGStuNNzJ4SQyGStuLdurS09d0IIqUnWintuLpCXR8+dEEL8yFpxB1jTnRBCIpH14k7PnRBCapLV4s7iYYQQ4k9Wizs9d0II8SerxZ2eOyGE+JPV4k7PnRBC/MlqcbeeuzHptoQQQjKLrBb3/HxdrOPIkXRbQgghmUVc4i4io0RknYhsFJEpPu8XicjrIrJKRBaLSLfEm1oTW4KgoiIVZyOEkOwhpriLSC6A6QAuBTAQwEQRGejp9t8AnjLGDAJwP4BfJ9pQPzp00Hb//lScjRBCsod4PPdhADYaYzYbY6oAzAYwztNnIIDXQ6/f9Hk/KVhx37cvFWcjhJDsIR5x7wpgh2u7NLTPzUoAV4ZeXwEgT0Ta19+86FDcCSHEn3jEXXz2efNT7gRwgYh8BOACADsBVNc4kMhkESkRkZK9e/fW2lgvFHdCCPEnHnEvBdDdtd0NwC53B2PMLmPMeGPMEAB3h/bVmF5kjJlhjCk2xhQXFBTUw2ylXTttKe6EEBJOPOK+DEA/EeklIk0ATAAw391BRDqIiD3WVABPJNZMfxo31oyZBDwEEEJIoIgp7saYagC3AFgEYA2AOcaY1SJyv4iMDXUbDmCdiKwHUAjggSTZW4MOHei5E0KIl0bxdDLGLACwwLNvmuv1CwBeSKxp8UFxJ4SQmmT1DFWA4k4IIX5kvbgXFFDcCSHES9aLu/XcWTyMEEIcAiHuX30FHD6cbksIISRzCIS4AwzNEEKIG4o7IYQEkMCIOycyEUKIQ2DEnZ47IYQ4UNwJISSAZL24t2kD5OZS3AkhxE3Wi3tODtC+PbBnT7otIYSQzCHrxR0ABg0CXnuNE5kIIcQSCHG/5hpgyxZgyZJ0W0IIIZlBIMR9/HigRQvg6afTbQkhhGQGgRD3Vq1U4J97TksREEJIQycQ4g4Al18OHDgAfPxxui0hhJD0Exhx791b2+3b02sHIYRkAoER96IibbdtS68dhBCSCQRG3Nu21dg7PXdCCAmQuIsAPXrQcyeEECBA4g5oaIbiTgghARR3hmUIISRg4t6jB7B/P3DoULotIYSQ9BKXuIvIKBFZJyIbRWSKz/s9RORNEflIRFaJyOjEmxobZswQQogSU9xFJBfAdACXAhgIYKKIDPR0+y8Ac4wxQwBMAPCXRBsaD1bcGZohhDR04vHchwHYaIzZbIypAjAbwDhPHwOgdeh1PoBdiTMxfnr00JaeOyGkodMojj5dAexwbZcCOMPT514Ar4jIjwC0BDAiIdbVki5ddOEOeu6EkIZOPJ67+OzzVk6fCOBJY0w3AKMBPC0iNY4tIpNFpERESvYmYUXr3FygWzd67oQQEo+4lwLo7truhpphlxsAzAEAY8x7AJoB6OA9kDFmhjGm2BhTXFBQUDeLY9C7N7B5c1IOTQghWUM84r4MQD8R6SUiTaADpvM9fbYDuAgARGQAVNwT75rHQd++wMaN6TgzIYRkDjHF3RhTDeAWAIsArIFmxawWkftFZGyo2x0Avi8iKwHMAnCtMelZ9K5vX2DvXi3/SwghDZV4BlRhjFkAYIFn3zTX608BnJNY0+pGv37abtwInHZaem0hhJB0EagZqoB67gBDM4SQhk3gxN0u2kFxJ4Q0ZAIn7i1bar47xZ0Q0pAJnLgDGnffsCHdVhBCSPoIpLgzHZIQ0tAJrLiXlQEHD6bbEkIISQ+BFHebDvn66+m1gySfv/8d+OijdFtBSOYRSHEfNQo45RRg0iRg1ap0W0OSya23Ao8/nm4rCMk8AinuLVsCCxZoe/vt6baGJItjxzT0VlWVbksIyTwCKe6AVoe87DJg5cp0W0KShR1TOXo0vXYQkokEVtwBYOBAXVM1CdWFSQbwxRfaUtwJqUngxR0APv00vXaQ5GCLw1HcCakJxZ1kLfTcCYlMoMW9a1cgL4/iHlQo7oREJtDiLqLeO8U9mDAsQ0hkAi3uAMU9yFjPnamQhNSkQYj7Z58B5eXptoQkGnruhESmQYg7wCnqQYQxd0IiE3hxP+88oF074Pe/T7clJNHQcyckMoEX97w84M47tRzBe++l2xqSSOi5ExKZwIs7APzoR0BBAfCrX6XbEpJIKO6ERKZBiHurVsAVVwBLlgDGpNsakigYliEkMg1C3AHgpJOAigpdxIMEA3ruhEQmLnEXkVEisk5ENorIFJ/3/yAiK0I/60Xk88SbWj9OOknb1avTawdJHPTcCYlMTHEXkVwA0wFcCmAggIkiMtDdxxhzuzFmsDFmMIBHAMxNhrH1geIePOi5ExKZeDz3YQA2GmM2G2OqAMwGMC5K/4kAZiXCuERSWAi0bcvZqkGCnjshkYlH3LsC2OHaLg3tq4GIFAHoBeCN+puWWETUe6fnHgyqqoCvvtLXFHdCahKPuIvPvkg5JxMAvGCMOeZ7IJHJIlIiIiV707CChhV3ZsxkPzYk07IlxZ0QP+IR91IA3V3b3QDsitB3AqKEZIwxM4wxxcaY4oKCgvitTBADB6YuY+aDD3SNT5IcrLh36ABUV/OGTYiXeMR9GYB+ItJLRJpABXy+t5OInAigLYCMnQdqB1VXrUruedauBc44A5g3L7nnaci4xR2g906Il5jiboypBnALgEUA1gCYY4xZLSL3i8hYV9eJAGYbk7k+VHEx0Lo18Ne/Jvc8a9Zou3lzcs/TkLGDqRR3QvyJK8/dGLPAGHOCMaaPMeaB0L5pxpj5rj73GmNq5MBnEvn5wO23A3PnJr5KZFkZcOaZwJYtjqjvihS8IvXGeu7t22tLcScknAYzQ9Vy221AmzbAvfcm9rjLlwPvvw+8+irFPRXQcyckOg1O3Nu0Ae64A5g/H1i2LHHH3bNH208+obinAnruhESnwYk7ANx6q9Z4v+eexB3TZuBQ3FPDl19q27atthR3QsJpkOKelwf87GfAwoUaSkkEVtw//hjYulVf79zJFL1kYScw5eVpS3EnJJwGKe4A8MMfAi1aADNnJuZ4Niyzb5/OnuzfXwXo84wroRYMKiuBnBygWTPdprgTEk6DFfdWrYDRo4EXXwSOH6//8crKtMSB5dxztWVoJjlUVQFNmwKNG+s2xZ2QcBqsuAPA+PHAZ5/p8nvbt4fPKD1yBFixIv5j7dkDDBnibJ93nrYU9+RQWUlxJyQaDVrcv/51oEkT4Ac/AHr2DF9E++GHgWHDnKyMWJSVqbi3aQM0aqQzVAGKe7KguBMSnQYt7q1bAxdfrMXEGjcG/vEP571ly1Qwtm+PfZzjx4G9e7Ws8MknA0VFQPdQNZ50ifuBA8EezKW4ExKdBi3uAPDHPwKzZgF3363FvmzWiw3JxCPu5eUa0iksBH75S30CaNFCvfh0iPsbb6gtf/tb6s+dKqy4N2mi2xT32mOdEhJMGry49+0LTJgAXHaZeroLF2ooZtMmfX/Hjui/DziZMh07AhdcAIwNVdzp0kXTIVPJqlXA5Zer+G3bltpzpxKv515VlV57spH584EePdQ5IcGjwYu7ZfBgoGtX4OWXw6tGRvLcV6/WOjKA4+0XFob36dIl9Z77o4/qTap58/jHC7KRykr12hmWqTs7dmi67r596baEJAOKewgRYMwY4F//At58U/e1bOnvuRujfW++WbetuHfsGN6vqAjYsCExqZbxsns30KuX1lyx9VeCCGPu9efIEW0PHUqvHSQ5UNxd3HwzcPgw8MADQEGBevN+nvuGDToL1VaWtGEZr+c+fLg+8i5fnkyrw9mzR28y+fnBFnfmudcfK+6HD6fXDpIcKO4uBg0CJk5Ur3DwYI1H+nnur7yibVmZimlZGZCbq/Vq3Iwcqe2iRcm1201Zmd5kWrcOfliG4l4/bAkHeu7BhOLu4b77NE+9uNgRd29Y5ZVXdOo7oLVk9uxRTz/HczU7dtTc91SKe0Px3Cnu9Yeee7ChuHvo21fDLVOmaK760aNO2AXQcMCbb+rsVkAHX6237Mcll+gM2FR40YcPa7XEwkIVd3ruJBqMuQcbirsPJ5+sYY0ePXTbHXd/+20V0KuvVhEtKdHKkv36+R/rkkt0AWc7SJtM3CmZrVvTcyfRsWEZeu7BhOIeBTvL1Mbdq6uBn/4U6NwZGDFCY/Rz5qio3nCD/zGGDdP2k0+Sb687JZNhGRILeu7BhuIeBa/n/sc/6szVRx7RqpKnnKKC36uXM3jqpUULvRnYSVHJxO255+erAFZWJv+86YB57vWHMfdgQ3GPQtu2KuKrVgH79+tg65gxTrx90CBtf/CDmoOpbnr31tWZjh0DJk0Cli5Njr1uz711a30d1Lg7Pff6w2yZYNMo3QZkMiIaW//b33Qg9csvgV//2qnbPmYMcNNNKu7R6N1bY+7r1wNPP63Hmj078fZ6PXdAQzMFBYk/V7phnnv9YVgmOmVlmg03YkS6Lakb9NxjcO+9KiLPPgtceaUOtlratwf+8hctEBaNPn20xswHH+j2woXJqYVSVqYee7NmwfbcjXHEPSdHfyjutYdhmej85S8abs3W2jtxibuIjBKRdSKyUUSmROjzLRH5VERWi8iziTUzfXTqBNx1lwrIz39et2P07q2CNH++bn/xBfDWW4mz0WJz3IFwzz1o2Btj06baNmnCwmF1gWGZ6Ozbp/+3776bbkvqRkxxF5FcANMBXApgIICJIjLQ06cfgKkAzjHGnATgtiTYmjbuvhvYuBE49dS6/X6fPtouXKgpk82bA/PmJc4+y549Tr69Ffcgeu52kNiKe+PG9NzrAj336FjH6O2302tHXYnHcx8GYKMxZrMxpgrAbADjPH2+D2C6MaYCAIwxexAgcnI0I6au9O6t7ZEjukLTyJEq7oleTKOszPHcbVgmiJ47xT0xMOYenYYg7l0BuCuslIb2uTkBwAkiskRElorIqEQZGAQKCzUlElDv/8orgdLSxGfN+HnuFHcSCU5iis7nn2u7fDlw8GB6bakL8Yi7+Ozz+pyNAPQDMBzARACPi0iNYUYRmSwiJSJSsrcBLQEj4njvp54KjBung56zZjl9KivrVxq4ulpjhF7PnWEZEgl67tE5cECdpGPHgH//O93W1J54xL0UQHfXdjcA3iUoSgHMM8YcNcZsAbAOKvZhGGNmGGOKjTHFBUHMz4uCW9xbt9bFuefM0S9OZSXQvz9w//11P74d/LHi3rSp/gTZc7dL7FHca48xznWk5+7PgQOaBimSvLkpySQecV8GoJ+I9BKRJgAmAJjv6fMSgAsBQEQ6QMM0mxNpaLZz/vnA0KGO+E6cqDHyxYuBuXO1Pvw779T9+HY5v66ugFlQSxB4s2Uo7rXHhmQAeu6ROHBA/59atszO/6OY4m6MqQZwC4BFANYAmGOMWS0i94tIaLVQLAKwX0Q+BfAmgJ8aY/Yny+hs5I47gA8/dLZHj9b673fdBTz8sO5bubLug6xW3Lt1c/ZZcf/Od4B//rNux60rH3+cvOXbGJapPzYkA9Bz9+P4cQ1p5udrCNV9M8wW4pqhaoxZAGCBZ98012sD4CehHxIHzZsDf/+7xt8BTZHcsAH47DOtRVNbSku1dYt769YaK9yxA1izRkNBqeJrXwO+9S1g+vTEHztTxd3emMVvlCpOFi/WEJ6ta5QsrFi1bUvP3Y+DB/XzzM/X/1X3zTBb4AzVNDJ2LHDPPbqQ9m9+o/tWrdJl/p57rnbHKi3VRUbc67jm5zsVLZcvD39ySCYHDqjXnqxKmJkq7tdcoz/14corgQcfTIw90bBi1b69hrmqq5N/zmzChmHatMlez53inmbuvVerTg4frtsLF+pM2Ouv12Jj8VJaqvFBdwEzmw7Zp496H489liiro2NvKGvXJuf4mSru771Xv7/56FGd6p6scJYbt7gDDM14seJuPXeKO6kTdv3Vrl2BRx91Hu8nT44/Bl9aGh6SAZx0yPHjNUTy7LOpmaZvSyTv2ZOcuhyZKO7HjunfXZ986IqK8DaZWLGy4s7QTDg2x93G3BmWIfVi0CAV37PPBn77W+D114FXX43vd/3E3Xrul1yiTwYHD/ov+J1o3OdYty7xx89Ecd+5U0MbX35Z92PYG2EqxN2KVYcO2tJzD4eeO0kotj78pEm6slNBgVamO3JEs10iTXIyxl/cTzhBB2fPPRfo2VP3bd1aN9uqq3Wxkj//OXZft7gnIzTjl+ee7sJhW7ZoWx/PPR3iTs/dH2/MnZ47qRdjxgDFxcBVV6lX+r3vAf/3fzrwOmaM5sMfOwb87nfA4MFaax5QMThypKa433yzik7TpkBRke6ri7h/9ZXeIG6/HfjFL2L3375dQ0xNmiRH3P2qQqbbc7fi/uWXdU9n3R9KHk5HWIaeezj03ElCOfdcYNkypz68XQTktddUwGbPBp56CrjzTp0ANWsWsHu3fxokoCl5VgC7ddPB1m3bam/XihW6CHjv3ipAscok7NihTwr9+iXXc091WGbtWp1F7Pf3W3E3pu5CaT33zz9PfFE5L96wDD33cIIQc+dKTBlMUREwbZp+uXbuBGbM0PTCU05RoT/pJPXmbcVKr7i7adxY36+L524nSF1wgebmV1Q4Hp8f27dr9ctOnXQyU6JJh7h//DFw0UXA3r06KcyWk7BYcQc0NNOyZe3PYT3348f1GHZAPBkwWyY6Bw6oQ9WsGT13kiTuuUdnsU6YoKK2bh3wk58AAwfqz/PPR/bcvRQV+Xvux45F/z17/CFDtI1W8+34ce3fvbvWy9m0KfHx8HSI+/jxTrhkv8/ca6+41wV3ZlGyQzPMlonOgQPOE7TXcy8v17GwZD9d1ReKe5Zw1lkqzoWFWpcGAL75Ta01/eabGnLp1Cn6MXr2rOm5P/64/t769ZF/b+dOFdL+/XV7T5Rq/Xv3qph37643n2PHEh+aqazUkFOj0HNnssX9yBFdrGXkSN32S+/cssURg2wQd3ru0fn8cyfbzOu5z5oF/PCHTspvpkJxzxJE9Es1d67jsU6cqK9nz9ZZro1iBNl69nRS9gCgpES/pPv2AX/9a+Tf27lTB0jt7Ndonrv9wvfoAZx2mr5O9MzYykr9u+00/2SLu83+GTxYW6+4V1YCu3ZpuAyoezqk+4kgVeLOmLs/ttwvUNNz/+wzbTP9mlHcs4izztIceMuJJ2qY5tln9ScWRUXqSZeWqgiPH69PAiNHAjNnRo4r2tmvtkpzNHG3Qti9uw6o5uUlT9wtscT9llv0JlZXvOLuDcts366P6DaVtT6eu13UJRVhGRHnaSPThSrVuMW9eXP9v7FOUVmZtpl+zSjuWU6PHurBn3de7L421/3993Vw8IsvdNHun/5UBevFF/1/z3ru1svzE3djNI3zG99w7MrJ0Th9usX9lVeAf/0r3NZZs+LPgLBPI3YNXa/nbscxrOdeH3G36+3abI1kceSIM1gowrCMF6/nDjjfF+u5Z/o1o7g3IGyu+w03aBhh4UL1Rr/2NRWVRx5R4Vu2zFl5xk6Qsnnr+fn+4r59u4r4uHHA//6vE8s97TQtZZzIwlRVVc4EJiC6uB87puK7bZvTZ9Uq4Nvfjr842/btKoBFRZrB4hV368nbrKX6hGWsuKciLGOFvUWLzPdCU83nnztPNc2ba2ufbK3nTnEnGUP37vrPfOiQpjSedZbuz8nRDJz33lPBu/hifRKYMUOFrLLSycQpKHDE/fhxrVO/apWT8njnnc7kKkDF/cgRLTmcKGrjue/apTcDW/sFcAqybdgQ3/m2b9dB56ZN9ablDctYsbdleuvjudsnnlSEZaxotWyZ+UKVauLx3DP9hsg89wZE06bAZZdpDvq3vhX+3vXX6+zTq69WsbzgAp1EZb/QdoUnt7j/+9/A73+vfbqHFmI8+eTw49pB1ZISJ2zhx7Fj+pQQa1AY8Bd3Y/QYubnhfd0pips2qWdsM4Y2bYp9LkBj7vbva9eupuduhdj2qYu4V1Xp73XooB5jKjx3K1otW9avJk7QqK5W4XbH3AG9IRpDz51kKPPmAf/v/9Xc36yZeuHHj2vJ4QULgFatVLyBcHG3qZAvvKDtO++o515U5PxDWE44QY/z7rvR7br2Wi1uFk/usJ+4A/7eu1fcASdGHq+4b9/ueOWRxL15cxXJ5s3jE8rXXw8XcPu6fXtdQCNVYRkg/DMlzueXl6et23P/4ovsWXuW4k7+w623ai2bu+7SL/SYMU4ow4ZlOnZUz/34cRX3nBydNfvuu/6eeU6Oxrdnzoy8eMfevZrOuWRJfIOvtRV3W4ahLuJuTLi4+4VlKipUkAG9kcXy3Ldu1YWX3fX17Q2jXbvUiLs7LNOpkxNqIE64pVUrbd2eu/s6ZXpYhuJO/kPjxiroNjRiM19EnKX/Cgo0L37pUs2iueEG3b9jR+Swy69+paGGm27Sm8Lu3cDddzuezzPP6KNw48bxLSjiFXc7uOo3E3bLFr0x9e7tiLkNy1RU6L7u3TWjxo/9+9Vji+a5l5c74p6XF1vcFy7Udteu8PMAjueeqmwZgOLuxXrutoSE23O3IRmAnjvJYkaNUq+lY0fHOy4oUCGeMUNF9Re/cN6LJO7t2wMPPaTe/cSJ6rX+6leOoD75pKZRTpyo+fqxwhq18dw3b9Yslj59nIHUbdt00hegk7dKSyM/MbgnZdm/paIivHiY23PPy4ttvxV3dygk1Z67OyzTubM+PXGpPcV65Fbc3Z47xZ0EgpYtVXCLi519diLTM89o2mNhIXD66bov2oDpddfpOrFz5qi3nJurKZeffKKpktdeC3z/+yqMTz8d3S6vuNv0wXnzavbdssUR902bNGZaUaF5/gDwxBPauv9p3bgnZQEqvsePOyVhAT1eu3b6OlZYprISeOMNfe0Wd+u5t2uXmgFVb1jGGMbdLV5xd3vu9gnHZp1lMhR3EpXHHtM4vMWKe3W1ijGgM1zz83XwNBIiGstfsEBr4ZxyimbQLFqk719+OXDOOcCwYerlR/Miq6rCxd2mbk6bpuJtsWUBrLgfOqQ3FAC48EJtraj6CZsxai8QHpYBataBiTcs8847akdeXvg57bqp7gHVZBam8oZlgIYZmnnzzZrftWgx97IydUwKC+m5kywnJ8ep4QI44t6zp+P9Tp2qeezuiUWRuPRSza8//XQV99de0zIKXbvqee6+W2Pis2dHPkZlZfi5RHQBkz17tLVs26YC2bu3493bmaoDB4YXWvPz3G+6ScNPN97o1NWxk7Pc4l6bmPtrr2kYaezY8HNu3643yNat9VhHjyZXPLxhGUDHQpLBwYNavTTT2LhRJ/DNmRO+P5bnXlCgn3MgxF1ERonIOhHZKCJTfN6/VkT2isiK0M/3Em8qyQSsCH//+yr8gAqtFYh4Of10FcVXX3VuEoAO6J5yCnDffZEfe71hGXu8Sy/VyVk2Hm6rUfbpo7n9eXnA9Om6r6jIEfyTT67puVdUaDx+8mQt72qxnrv1+I8e1VCSO1smWsx9yxY9d69e6q3bcstbtzrlIewNxHrzyeDw4fCwDJA8z/3vf9d5FdEqj6YD+/d6J9jFirkXFuqs3qwXdxHJBTAdwKUABgKYKCIDfbo+Z4wZHPp5PMF2kgyhsFBDGz/7Wf2OY+P4x46Fi3tODvDwwxofv+WW8N954QUNwRw8WFPcAU253LFDZ9oCwFtvab+hQ1Uwp0xxwhGFhWrDSSdpOMjrua9cqe0VV4Q/uXjDMjarxe6P5bnv3q2DuR076k3IHsct7nYhkHhn0NaW8nK9OdlyCYWFjm3JYPXq8DZTsOMa3pTYWNkynTplR8mGeDz3YQA2GmM2G2OqAMwGMC65ZpFM5rTT4ptJGo2TT9Z/GhGdvORm+HCdSPXkkyrY1rN68EENaxw44C/u48bpMWfN0u033lDhtv+ct91awXDJAAAUl0lEQVSmTx49e+p5f/tbLaJWWKhesjv2asXdFguzeMMyViC8YZlI8fLdu/Upx4Z59uzRvm5xt3Xz163zP4abL76ovYdvS0XYKpbNmqn9yfLc7eeXyBIUicB+dhs3hu+P5rl/9pl+X7KhZEM84t4VgGs9e5SG9nm5UkRWicgLItI9IdaRwNK4sYZShg1zvF4306bpAOz8+RpSef11jdHbUIqfuOflaXmFOXPUw1qxQmOqlhYtgJdf1gVKrA0tW6rQGhM+OWnlSo2tehdAscWkbF8r8u6wzPHjkcsn+4l7ebl6i7awW+fOepx4xP3GGzWT56GHYq9ta7Hi7s5uSmauuxX1Tz9NzvHrSiTP3Yq7Lb/cqJE+UXo99yCIu/js8/ol/wegpzFmEIDXAMz0PZDIZBEpEZGSvdGKgpMGwaxZwD/+4f9ebq6mTq5Yodkx48bpvjfe0BWo3KLt5rvf1ZxtWzvH22/wYPXm3diwhDs0s3Kleu3i+fY3aqQDn9E8d8A/NHPwoAqHV9ztpCrruYuo9x6PuC9bpsJz113ASy/F7g9oobf27cPHSTp1Sk5YZt8+58ki0Z77hx/WL2XU/m55efgA+aFD6q3bOkUiuv3ZZ/pdtDH3IIRlSgG4PfFuAHa5Oxhj9htjQhUX8BiA0/wOZIyZYYwpNsYUF9i0C9Jg6drVqVkTib59NVZ+6JCmXPbooZ75JZf49x892ll+sFWr8Bz9SLiFdt069chWr64ZkrEUFjolDKxAuGPugL+4W/H0irs9lhV3QDOIYi1PWFmpE7N+9CMVouXLo/e3fPyxeu3uG1fnzvX33NeuDR98BhxBt39PvE8XsTh6FDj3XOAPf6j7Mdw3Brf3fuhQzQXOmzVzbsKBGVAFsAxAPxHpJSJNAEwAMN/dQUTcuRJjAWRYdI1kM3fdpYJ9112x+4polkuPHjr4amevRsN67mvWaBx6+HAVTrvykpcLLtAniKNHI3vufhkzttxA587qOefk6NOC13MHVAy3b48uIBs2qFieeqqGq+JZq/b4cRV3G2+32LBMfXLrH3xQV7xye7RW3MeP178lUeuObtumoS+7eHtdKC93bnBucf/yy5ri3ry58zl16hSQmLsxphrALQAWQUV7jjFmtYjcLyJjQ91+LCKrRWQlgB8DuDZZBpOGR/Pm6q1fcEF8/du21XDOU0/F19+K+4sv6mO3negUyXMfPVo98yVL/GPugJZW8AqP9dy7dFFht9UYt27V/HYbzwecQdVoGTNWzPv31x+/sMemTeE1d7ZsUfH1zibu1EnF6uKLNRwWjXXrdAzEy+LF2roFfM0a9XJHjXK2E4EdBK1PdLeiQm+i7uMBen3s52hp1sz5u4IUloExZoEx5gRjTB9jzAOhfdOMMfNDr6caY04yxpxqjLnQGJPg9e4JqR1t29b8B41Efr7m6r/1lnr648bpPiuwXi66SPstWKAC0bKl84Rg49g/+5neHNwC7w7LABqaseJeVBQeJrGiEy3u7g55DBigNwJ3xs9jj+k6tu5QyapV2no9d1v18/XXtYJnNO67TydhucVt61bHs/WK+4kn6qQxt831xXra9RX3Ll00NOgVdz/P3dYusuJeVZXZ9Xg4Q5U0eEScvPPiYs2nX7s2ckgnLw84/3xH3N3ZPiefrAK6cKGGdr79bUcAdu/WLB/robvF3R2SAVSURaKHWtas0ZtCixYq7kePOqI3e7ZOwDLGyS9fskTDJvn5mt/v5oortM7O7bfrDSXaRCwbLnr5ZWffW2+Fv++2ccAAXYSkoCDzxL1tW6fukCVSzB3QsY327Z33412HNx1Q3AmBM8B57rmaEeNNgfQyerSK5ttvOyEZyymnaBji0Ue1lszdd+t+mwZpPfSOHdVj3LjRSYO0NG+ugm89bTeTJ2tBNyucgPOUYW8Gc+fquENxsQpXdTXw9a/rjeDtt/3F67rrtOaOMU6evx/2acRdImLxYr3J5eY6A8R79qjQDxmi2wMGJC4dMlFhmbZtddKYrRgKRPbcAf3McnKcNMlMjrtT3AmBE3c/99z4+l97rXrpmzfXFHfLNddoHvpDD2kVTSvuFrvwSatWwPd8CnYMH64Dt7ZEAaAx88cfV2FdscIRddtaz3jtWg0L9e+v4r5pk07+mjatZkjGjV0WMVLmzfHjWse/USN9OrGF2hYvVnu7dnU89/ff1/bMM7UdMEDtS0RBNOtpHzzorIxUW6y4d+2qA9v2OvsNqFrP3X5PrLhnctyd4k4InH/as8+Or3+7dlpRcNgwRxD9+NOftGLlTTdpuMMt7t/4BjBpEvDRR/6DtyNHqgC5a80/+qh6jlagrajn52v8eM0aFan16zXW3aePlmT46CPt5w3HeOncWa9FJHHfs0efAiZOVFGdP1899a1bVdyLihzPfelSvQkMHarbAwbo31Pf0sLHj+tNtXVr3d67V58IapPzfuSI2t+2rf7Nx445+fh+A6reOjz03AnJEq66CrjzTo0Nx0uHDipg7kqUXpo00cqShw6px2sXCQH0KWHmzMi5/rbmjl3U5PBh9dqvuEIXNRk0yCldDDgZM9u2qXD176/ibgzwz386faIhomIcafESG5IZP15j6IsWOfH24cM1FGQ996VL9aZlhTBRg6q7dmka5Bln6PbevZpJdfXV8R/DncJqPxObqhot5m6dAPs+xZ2QDGfUKK01U1u8M1j96N/fWY6wNtUzCwpUaK24//d/qyj9+Mfqga9cGV5Df8gQ3We9dCvugIp7r141RcuPoUPVE/YbLLSLl/TooStqvfqqPsG0b6829eihN4CqKuCDD5yQDOCMD9Q37m5DMvbYn36qXvfChU7RuFi4xd1b8jhazJ1hGUJIGPfeq0LsFrt4GDlSBWvmTOCXv9RwyHnn+fcdMUJF1a5D6xb3iorYIRnLaadpmMIuB+jGeu7dumlOfFmZMwchJ0fDMtXVmlL55Zfhf2/XrppptGaNetu2ouZLL+lyje+8oyWfr702un1ecXcL+rRp2i5bpjZGKqngnllsxX3XLie9MZLnHi0ss3q1Hi9ShtP8+cAjj0T/2xIJxZ2QFNCli8bcI9XEicSkSZo6ee21Kox//GPkvuefr2GgRYs0ZNS+vQ7aWqGyYZFYjB6tIZ8bb9QB3YcfdjzU0lI9R4cOKu6ACpyt7GlXrLLF2WzoBNCnnAED9Mli2DBHxKdMUVG26aUzZ0au/b5sGfDAAxpvt2MdVtyvv16rhm7dqqGinTs1/dMPt+duBXv3bicFNFLMPVpY5t139bh+NX4WL9YxlqlTk7vClhuKOyEZzIABGsN+6inNK7cpm360aOEURbOToESc+vDxeu5Nm2pM/+BBjfvfequWFgBU3Lt2VS+9Wzcn1OIV97lzdTygb9+af8+SJY4AV1ToTe/mm/WJ46OPdBDWPn24WbtWbwDV1fpUUVCgqZcrV2p7zTXa79NPnZuDDVF5cYt706Z6I9y9u2a5X0skz90dlrHntGE0S3m5jlEYo/137EBKoLgTkuE0b67CddZZsftab9o9cGpDM/GKu+07b55m54wZo08M5eUq7nY2K6CDu0VFzrGtuOfnqwfuHZOwN4P8fA3LPPusbo8bp+mgp56qr598UlMp7azXY8d03KJFC91/9tl6g2nfXt8rKnJKKqxZUztxBzQ0s2tXZHGPFHN3e+52NvGSJeH7ly7V89n5Dqmqa09xJyRAjByprVvcTzhBPdtYmTJ+x7rxRg2DHDyoWUFecb/vPvWU7ZKLeXk6C/aZZ7TOvJcLL9Tff/JJ3X74YW3d1TsnT9YB0jPPVBGvrlZP/t//1puMe4KZzW7q10+FvqBAPfx4xT0/X9suXaJ77jZMY+PzfuK+fr0+WVVV6UQxy4oV2n7nO9pS3AkhtWboUF2z9LrrnH0/+YkuDB5PpowfgwYBEyZoNtH27eHi3qiRI3SWP/9ZZ8P6MWyYhiXGjtWxhPXrVZjdJRwuvljr/P/85yq4ixfrQOQZZzgCabGVw/v103bAAM3S2b3bqXTpV8q4okKF3dZs93ru3pj7pEn6JGPttH/z8uUactm3T3Pvr7lGwzzu0MyKFRoa69tXb0CpWrSE4k5IgBDRgUq7HCCgoYQRI+p33OnTVdSrq8PFva7k5DjZLsOGhb8nooI5daqK7NSpKog33FAzzOMV9/79nZIN3/ymts8+q+LsXmmrvDx8ZnGXLnoTsHX4vTfC9u31hmTJzVURnzVLq4n+7ncaHjr1VM1mcmcarVjhlI8eOJCeOyEkg2jXTr3poqKaYlxX7BjC6af7v9+8uQpqSYm+tqtrubHibgdubUwfcMT9jjuAp58Gnn9et48c0cFcOxYBqOdeXe3Mro3nKcf9xPLoo9qeeKKOGaxdqzH4gwe1Do4Vd1tfJxUZMxR3QkhcDBmiA5y1zdWPxIgR6sF7F0h3c9VV2o4f78TH3fh57pbiYg2HtGih/aw3/cgjGhr6r/9y+tpYui1IFo+4t2ypPxMnat0eQMc3xo3T1/Pm6cIoxoSLe3l5/QqexUs917AnhJC6cfbZGipxL1Li5ZJLNBxz++3+748cqWGYXr1023ruPXqot//44xo+eeYZzd7ZtQv49a91TMB9U7ElCOziKPGI+/jx6v2fdJKGZzp21L+lTRvNwX/pJSd2bytjukswREtrTQQUd0JI2ogm7IAKs50Q5cc554QveN69u3rqtiyDrb1TXq6LlowYoYOmNm/fYsXdZrbEs9DLn/6k7dGjGr93l4K4/HKdmHXggMbrbf0gdwmGeFcWqysUd0JIYMjJAW67reZs3Asv1Jm1a9ZoGQdvzn9RkXryixfroK3Na4+Hxo11oRN32Oiqq3SQtapKs5XsQHC3bjrByz1zN1mISdVcWA/FxcWmpKQkLecmhDQ8rrxSa+EsXqwpnF527NBsl6qq6CtRpRsR+dAYUxyrHz13QkiD4PnntRa8n7ADGtKZOzdyueNsg+JOCGkQ5OQ4M2kjMXx49OydbIKpkIQQEkAo7oQQEkDiEncRGSUi60Rko4hMidLvGyJiRCRmsJ8QQkjyiCnuIpILYDqASwEMBDBRRGqU/ReRPAA/BvB+oo0khBBSO+Lx3IcB2GiM2WyMqQIwG8A4n36/APAQgK8SaB8hhJA6EI+4dwXgXjukNLTvP4jIEADdjTEvJ9A2QgghdSQecfdb3/0/M59EJAfAHwDcEfNAIpNFpERESvamonIOIYQ0UOIR91IA7jVVugHY5drOA3AygMUishXAmQDm+w2qGmNmGGOKjTHFBbacGyGEkIQTs/yAiDQCsB7ARQB2AlgG4NvGmNUR+i8GcKcxJmptARHZC2BbHWwGgA4A9tXxd5NNptpGu2oH7ao9mWpb0OwqMsbE9I5jzlA1xlSLyC0AFgHIBfCEMWa1iNwPoMQYM78OxiEe4yIhIiXx1FZIB5lqG+2qHbSr9mSqbQ3VrrjKDxhjFgBY4Nk3LULf4fU3ixBCSH3gDFVCCAkg2SruM9JtQBQy1TbaVTtoV+3JVNsapF1pq+dOCCEkeWSr504IISQKWSfu8RYxS4Ed3UXkTRFZIyKrReTW0P57RWSniKwI/YxOg21bReTj0PlLQvvaicirIrIh1LZNsU0nuq7JChH5QkRuS9f1EpEnRGSPiHzi2ud7jUR5OPSdWyUiQ1Ns129FZG3o3C+KSJvQ/p4icsR17f4nxXZF/OxEZGroeq0TkUuSZVcU255z2bVVRFaE9qfkmkXRh9R9x4wxWfMDTcXcBKA3gCYAVgIYmCZbOgMYGnqdB50LMBDAvdA8/3Rep60AOnj2PQRgSuj1FAAPpvlz/AxAUbquF4DzAQwF8EmsawRgNICF0NnaZwJ4P8V2jQTQKPT6QZddPd390nC9fD+70P/BSgBNAfQK/c/mptI2z/u/AzAtldcsij6k7DuWbZ57vEXMko4xZrcxZnno9UEAa+CpuZNhjAMwM/R6JoDL02jLRQA2GWPqOomt3hhj3gZQ7tkd6RqNA/CUUZYCaCMinVNllzHmFWNMdWhzKXSWeEqJcL0iMQ7AbGNMpTFmC4CN0P/dlNsmIgLgWwBmJev8EWyKpA8p+45lm7jHLGKWDkSkJ4AhcMod3xJ6tHoi1eGPEAbAKyLyoYhMDu0rNMbsBvSLB6BjGuyyTED4P1u6r5cl0jXKpO/d9VAPz9JLRD4SkbdE5Lw02OP32WXS9ToPQJkxZoNrX0qvmUcfUvYdyzZxj1rELB2ISCsA/wBwmzHmCwCPAugDYDCA3dBHwlRzjjFmKLQG/w9F5Pw02OCLiDQBMBbA86FdmXC9YpER3zsRuRtANYBnQrt2A+hhjBkC4CcAnhWR1ik0KdJnlxHXK8REhDsSKb1mPvoQsavPvnpds2wT91hFzFKKiDSGfnDPGGPmAoAxpswYc8wYcxzAY0ji42gkjDG7Qu0eAC+GbCizj3mhdk+q7QpxKYDlxpiykI1pv14uIl2jtH/vROS7AMYAuNqEgrShsMf+0OsPobHtE1JlU5TPLu3XC/hPXazxAJ6z+1J5zfz0ASn8jmWbuC8D0E9EeoU8wAkA6lTbpr6EYnl/A7DGGPN71353nOwKAJ94fzfJdrUUXRULItISOhj3CfQ6fTfU7bsA5qXSLhdhnlS6r5eHSNdoPoBJoYyGMwEcsI/WqUBERgG4C8BYY8xh1/4C0ZXSICK9AfQDsDmFdkX67OYDmCAiTUWkV8iuD1Jll4sRANYaY0rtjlRds0j6gFR+x5I9apzoH+io8nroHffuNNpxLvSxaRWAFaGf0QCeBvBxaP98AJ1TbFdvaKbCSgCr7TUC0B7A6wA2hNp2abhmLQDsB5Dv2peW6wW9wewGcBTqNd0Q6RpBH5mnh75zHwMoTrFdG6HxWPs9+59Q3ytDn/FKAMsBXJZiuyJ+dgDuDl2vdQAuTfVnGdr/JIAbPX1Tcs2i6EPKvmOcoUoIIQEk28IyhBBC4oDiTgghAYTiTgghAYTiTgghAYTiTgghAYTiTgghAYTiTgghAYTiTgghAeT/Awcq91ZSPWcuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(epoch_losses, c='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "test_bg = dgl.batch(test_graphs)\n",
    "test_y_tensor = torch.tensor(test_y).float().view(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_bg.set_e_initializer(dgl.init.zero_initializer)\n",
    "test_bg.set_n_initializer(dgl.init.zero_initializer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit = model(test_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "probs = torch.softmax(logit, 1).detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_y = np.argmax(probs,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.metrics import confusion_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7587548638132295"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(test_y, pred_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.70      0.86      0.78       102\n",
      "           1       0.79      0.64      0.71       115\n",
      "           2       0.87      0.82      0.85        40\n",
      "\n",
      "   micro avg       0.76      0.76      0.76       257\n",
      "   macro avg       0.79      0.78      0.78       257\n",
      "weighted avg       0.77      0.76      0.76       257\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(classification_report(test_y, pred_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[88 14  0]\n",
      " [36 74  5]\n",
      " [ 1  6 33]]\n"
     ]
    }
   ],
   "source": [
    "print(confusion_matrix(test_y, pred_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import AllChem\n",
    "from rdkit.Chem.Descriptors import rdMolDescriptors\n",
    "from sklearn.preprocessing import normalize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate 3D conf\n",
    "train_mols2 = copy.deepcopy(train_mols)\n",
    "test_mols2 = copy.deepcopy(test_mols)\n",
    "\n",
    "ps = AllChem.ETKDGv2()\n",
    "for m in train_mols2:\n",
    "    m = Chem.AddHs(m)\n",
    "    AllChem.EmbedMolecule(m,ps)\n",
    "for m in test_mols2:\n",
    "    m = Chem.AddHs(m)\n",
    "    AllChem.EmbedMolecule(m,ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_dragon_type_desc(mol):\n",
    "    return rdMolDescriptors.CalcAUTOCORR3D(mol) + rdMolDescriptors.CalcMORSE(mol) + \\\n",
    "        rdMolDescriptors.CalcRDF(mol) + rdMolDescriptors.CalcWHIM(mol)\n",
    "train_X = normalize([calc_dragon_type_desc(m) for m in train_mols2])\n",
    "test_X = normalize([calc_dragon_type_desc(m) for m in test_mols2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfc = RandomForestClassifier(n_estimators=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',\n",
       "            max_depth=None, max_features='auto', max_leaf_nodes=None,\n",
       "            min_impurity_decrease=0.0, min_impurity_split=None,\n",
       "            min_samples_leaf=1, min_samples_split=2,\n",
       "            min_weight_fraction_leaf=0.0, n_estimators=100, n_jobs=None,\n",
       "            oob_score=False, random_state=None, verbose=0,\n",
       "            warm_start=False)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rfc.fit(train_X, train_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_pred_y = rfc.predict(test_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7587548638132295"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(test_y, rf_pred_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.77      0.87      0.82       102\n",
      "           1       0.79      0.66      0.72       115\n",
      "           2       0.67      0.75      0.71        40\n",
      "\n",
      "   micro avg       0.76      0.76      0.76       257\n",
      "   macro avg       0.74      0.76      0.75       257\n",
      "weighted avg       0.76      0.76      0.76       257\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(classification_report(test_y, rf_pred_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[89 11  2]\n",
      " [26 76 13]\n",
      " [ 1  9 30]]\n"
     ]
    }
   ],
   "source": [
    "print(confusion_matrix(test_y, rf_pred_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}