tutorial-7-dynamical-graph-convolutional-neural-nets.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Jet tagging with Dynamical Graph Convolutional Neural Nets (DGCNNs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is a graph?\n",
    "\n",
    "Essentially, just a structure that has:\n",
    "- nodes: entities with a set of properties or features\n",
    "- edges: these connect two nodes, they can be uni- or bi-directional.\n",
    "\n",
    "Graph networks model graph structures in a neural network architecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is this useful?\n",
    "\n",
    "Jets are graphs (or, point clouds).  Each constituent is a node, and edges connect neighbouring constituents.  Jet observables are generally constructed by computing quantities between pairs of particles, for example the angular separation between two particles or the invariant mass of two particles.  Since our networks will serve as jet observables, the graph construction is very natural.  \n",
    "\n",
    "The constituents in a jet also have very trivial symmetry, permutation invariance.  This means that any observable we construct should not depend on some arbitrary ordering of the jet constituents.  Graph networks are inherently permutation invariant.  If you remember the regression networks we built in tutorials 2 and 4, there we ordered particles by type and $p_T$ to input them to the network.  This was not a permutation invariant construction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we will learn how to construct graph neural networks using pytorch, and apply them to jet tagging.\n",
    "\n",
    "For this tutorial I am using the following sources:\n",
    "- https://pytorch-geometric.readthedocs.io/en/latest/notes/introduction.html#\n",
    "- https://towardsdatascience.com/program-a-simple-graph-net-in-pytorch-e00b500a642d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outline:\n",
    "- Imports and plotting set-up\n",
    "- Building the dataset\n",
    "    - Preprocessing\n",
    "    - Dataloaders\n",
    "- Data handling in graph networks\n",
    "- Jets, nodes, and edges\n",
    "- Message passing $\\rightarrow$ generalising convolutions\n",
    "- Dynamic edge convolutions\n",
    "- Building a Dynamic Graph Convolutional Neural Network (DGCNN)\n",
    "- Optimising the DGCNN\n",
    "- Study the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some interesting papers:\n",
    "- Dynamic Graph CNN for Learning on Point Clouds\n",
    "    - https://arxiv.org/abs/1801.07829\n",
    "    - This is the paper that introduces the DGCNN as we know it\n",
    "- ParticleNet: Jet Tagging via Particle Clouds\n",
    "    - https://arxiv.org/abs/1902.08570\n",
    "    - In this paper by Huilin Qu and Loukas Gouskos they apply the DGCNN idea to tagging jets\n",
    "- An Efficient Lorentz Equivariant Graph Neural Network for Jet Tagging\n",
    "    - https://arxiv.org/abs/2201.08187v2\n",
    "    - In this paper by Huilin Qu and others they construct a Lorentz equivariant DGCNN for tagging jets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import random\n",
    "import time\n",
    "import numpy as np\n",
    "import torch\n",
    "from torch.utils.data import Dataset\n",
    "from torch.utils.data import DataLoader\n",
    "from torch import nn\n",
    "from torch.nn import functional as F\n",
    "from torch.nn import Sequential, Linear, ReLU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import matplotlib\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "from matplotlib.lines import Line2D\n",
    "from matplotlib.font_manager import FontProperties\n",
    "import matplotlib.colors as mcolors\n",
    "import colorsys\n",
    "\n",
    "labelfont = FontProperties()\n",
    "labelfont.set_family('serif')\n",
    "labelfont.set_name('Times New Roman')\n",
    "labelfont.set_size(14)\n",
    "\n",
    "axislabelfont = FontProperties()\n",
    "axislabelfont.set_family('serif')\n",
    "axislabelfont.set_name('Times New Roman')\n",
    "axislabelfont.set_size(22)\n",
    "\n",
    "tickfont = FontProperties()\n",
    "tickfont.set_family('serif')\n",
    "tickfont.set_name('Times New Roman')\n",
    "tickfont.set_size(16)\n",
    "\n",
    "axisfontsize = 16\n",
    "labelfontsize = 16\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "plt.rcParams[\"mathtext.default\"] = \"rm\"\n",
    "plt.rcParams['text.usetex'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using the top-tagging data again, so let's just load what we saved out before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = np.load( \"tutorial-3-data/TopTagging/x_train_short.npy\")\n",
    "y_train = np.load( \"tutorial-3-data/TopTagging/y_train_short.npy\")\n",
    "X_test = np.load( \"tutorial-3-data/TopTagging/x_test_short.npy\")\n",
    "y_test = np.load( \"tutorial-3-data/TopTagging/y_test_short.npy\")\n",
    "X_val = np.load( \"tutorial-3-data/TopTagging/x_val_short.npy\")\n",
    "y_val = np.load( \"tutorial-3-data/TopTagging/y_val_short.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each jet we'll keep the $50$ constituents with the highest $p_T$, with each constituent consisting of $(p_T,\\eta,\\phi)$.\n",
    "So each constituent in the jet has three features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still need to preprocess the jets in much the same way as we did for the CNN study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some initial settings\n",
    "__n_warning__ = 0.7\n",
    "n_shift_phi, n_shift_eta = 0, 0\n",
    "\n",
    "# Grid settings\n",
    "xpixels = np.arange(-2.6, 2.6, 0.029)\n",
    "ypixels = np.arange(-np.pi, np.pi, 0.035)\n",
    "\n",
    "# Calculate the pseudorapidity of pixel entries\n",
    "def eta (pT, pz):\n",
    "    small = 1e-10\n",
    "    small_pT = (np.abs(pT) < small)\n",
    "    small_pz = (np.abs(pz) < small)\n",
    "    not_small = ~(small_pT | small_pz)\n",
    "    theta = np.arctan(pT[not_small]/pz[not_small])\n",
    "    theta[theta < 0] += np.pi\n",
    "    etas = np.zeros_like(pT)\n",
    "    etas[small_pz] = 0\n",
    "    etas[small_pT] = 1e-10\n",
    "    etas[not_small] = np.log(np.tan(theta/2))\n",
    "    return etas\n",
    "\n",
    "# Calculate the azimuthal angle of pixel entries\n",
    "def phi (px, py):\n",
    "    \"\"\"\n",
    "    phis are returned in rad., np.arctan(0,0)=0 -> zero constituents set to -np.pi\n",
    "    \"\"\"\n",
    "    phis = np.arctan2(py,px)\n",
    "    phis[phis < 0] += 2*np.pi\n",
    "    phis[phis > 2*np.pi] -= 2*np.pi\n",
    "    phis = phis - np.pi \n",
    "    return phis\n",
    "\n",
    "# function to calculate masses\n",
    "def mass (E,px,py,pz):\n",
    "    mass = np.sqrt(np.maximum(0.,E**2-px**2-py**2-pz**2))\n",
    "    return mass\n",
    "\n",
    "# function to return the image momenta for centroid and principal axis\n",
    "def img_mom (x, y, weights, x_power, y_power):\n",
    "    return ((x**x_power)*(y**y_power)*weights).sum()\n",
    "\n",
    "# returns the jet image\n",
    "def orig_image (etas, phis, es):\n",
    "    \"\"\"\n",
    "    Gives the value on grid with minimal distance,\n",
    "    eg. for xpixel = (0,1,2,3,..) eta=1.3 -> xpixel=1, eta=1.6 ->xpixel=2\n",
    "    \"\"\"\n",
    "    z = np.zeros((etas.shape[0],len(xpixels),len(ypixels)))\n",
    "    in_grid = ~((etas < xpixels[0]) | (etas > xpixels[-1]) | (phis < ypixels[0]) | (phis > ypixels[-1]))\n",
    "    xcoords = np.argmin(np.abs(etas[:,None,:] - xpixels[None,:,None]),axis=1)\n",
    "    ycoords = np.argmin(np.abs(phis[:,None,:] - ypixels[None,:,None]),axis=1)\n",
    "    ncoords = np.repeat(np.arange(etas.shape[0])[:,None],etas.shape[1],axis=1)\n",
    "    z[ncoords[in_grid],ycoords[in_grid],xcoords[in_grid]] = es[in_grid]\n",
    "    return z\n",
    "\n",
    "# preprocess the jet\n",
    "def preprocessing( x ,y, weights, rotate=True, flip=True ):\n",
    "    \"\"\"\n",
    "    (x,y) are the coordinates and weights the corresponding values, shifts\n",
    "    centroid to origin, rotates image, so that principal axis is vertical,\n",
    "    flips image, so that most weights lay in (x<0, y>0)-plane.\n",
    "    Method for calculating principal axis (similar to tensor of inertia):\n",
    "    https://en.wikipedia.org/wiki/Image_moment\n",
    "    here: y=phi, phi has modulo 2*np.pi but it's not been taken care of hear,\n",
    "    so possible issues with calculating the centroid\n",
    "    -> pre-shifting of events outside of this function solves the problem\n",
    "    for iamge-data with Delta_phi < 2*np.pi\n",
    "    \"\"\"\n",
    "\n",
    "    # Shift\n",
    "    x_centroid = img_mom(x, y, weights, 1, 0) / weights.sum()\n",
    "    y_centroid = img_mom(x, y, weights, 0, 1)/ weights.sum()\n",
    "    x = x - x_centroid\n",
    "    y = y - y_centroid\n",
    "\n",
    "    # Check if shifting worked, there can be problems with modulo variables like phi (y)\n",
    "    # x and y are sorted after highest weight, 0-comp. gives hottest event\n",
    "    # for Jet-like Images Centroid should be close to hottest constituen (pT-sorted arrays)  \n",
    "    global n_shift_phi\n",
    "    global n_shift_eta\n",
    "    if np.abs(x[0]) > __n_warning__:\n",
    "        n_shift_eta += 1  \n",
    "    if np.abs(y[0]) > __n_warning__:\n",
    "        n_shift_phi += 1       \n",
    "\n",
    "    if rotate:\n",
    "        #Ccovariant matrix, eigenvectors corr. to principal axis\n",
    "        u11 = img_mom(x, y, weights, 1, 1) / weights.sum()\n",
    "        u20 = img_mom(x, y, weights, 2, 0) / weights.sum()\n",
    "        u02 = img_mom(x, y, weights, 0, 2) / weights.sum()\n",
    "        cov = np.array([[u20, u11], [u11, u02]])\n",
    "\n",
    "        # Eigenvalues and eigenvectors of covariant matrix\n",
    "        evals, evecs = np.linalg.eig(cov)\n",
    "\n",
    "        # Sorts the eigenvalues, v1, [::-1] turns array around, \n",
    "        sort_indices = np.argsort(evals)[::-1]\n",
    "        e_1 = evecs[:, sort_indices[0]]  # Eigenvector with largest eigenvalue\n",
    "        e_2 = evecs[:, sort_indices[1]]\n",
    "\n",
    "        # Theta to x_asix, arctan2 gives correct angle\n",
    "        theta = np.arctan2(e_1[0], e_1[1])\n",
    "  \n",
    "        # Rotation, so that princple axis is vertical\n",
    "        # anti-clockwise rotation matrix\n",
    "        rotation = np.matrix([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])\n",
    "        transformed_mat = rotation * np.stack([x,y])\n",
    "        x_rot, y_rot = transformed_mat.A\n",
    "    else: \n",
    "        x_rot, y_rot = x, y\n",
    "  \n",
    "    # Flipping\n",
    "    n_flips = 0\n",
    "    if flip:\n",
    "        if weights[x_rot<0.].sum() < weights[x_rot>0.].sum():\n",
    "            x_rot = -x_rot\n",
    "            n_flips  += 1\n",
    "        if weights[y_rot<0.].sum() > weights[y_rot>0.].sum():\n",
    "            y_rot = -y_rot\n",
    "            n_flips += 1\n",
    "            \n",
    "    return x_rot, y_rot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to convert the jet to a particle cloud\n",
    "def constit_to_cloud( jets, n_constit, norm, rotate, flip ):\n",
    "    \n",
    "    print( \"Crop constituents\" )\n",
    "    jets = jets[:,0:n_constit,:]\n",
    "    \n",
    "    print( \"Calculating pT\" )\n",
    "    E     = jets[:,:,0]\n",
    "    pxs   = jets[:,:,1]\n",
    "    pys   = jets[:,:,2]\n",
    "    pzs   = jets[:,:,3]\n",
    "    pT    = np.sqrt( pxs**2+pys**2 )\n",
    "    \n",
    "    print( \"Calculating eta\" )\n",
    "    etas  = eta( pT, pzs )\n",
    "    \n",
    "    print( \"Calculating phi\" )\n",
    "    phis  = phi( pxs, pys )\n",
    "    \n",
    "    print( \"Calculating the mass\" )\n",
    "    E_tot = E.sum( axis=1 )\n",
    "    px_tot = pxs.sum( axis=1 )\n",
    "    py_tot = pys.sum( axis=1 )\n",
    "    pz_tot = pzs.sum( axis=1 )\n",
    "    j_mass = mass( E_tot, px_tot, py_tot, pz_tot )\n",
    "    \n",
    "    print( \"Pre-shifting the phis\" )\n",
    "    phis = ( phis.T - phis[:,0] ).T\n",
    "    phis[phis < -np.pi] += 2*np.pi\n",
    "    phis[phis > np.pi] -= 2*np.pi\n",
    "    \n",
    "    print( \"Using pT as weight\" )\n",
    "    weights = pT\n",
    "    \n",
    "    print( \"Preprocessing\" )\n",
    "    \n",
    "    for i in range( np.shape(etas)[0] ):\n",
    "        etas[i,:], phis[i,:] = preprocessing( etas[i,:], phis[i,:], weights[i,:], rotate, flip )\n",
    "    \n",
    "    cloud = np.stack( [ pT/np.mean(pT), etas, phis ], axis=2 )\n",
    "    \n",
    "    return cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Crop constituents\n",
      "Calculating pT\n",
      "Calculating eta\n",
      "Calculating phi\n",
      "Calculating the mass\n",
      "Pre-shifting the phis\n",
      "Using pT as weight\n",
      "Preprocessing\n",
      "Crop constituents\n",
      "Calculating pT\n",
      "Calculating eta\n",
      "Calculating phi\n",
      "Calculating the mass\n",
      "Pre-shifting the phis\n",
      "Using pT as weight\n",
      "Preprocessing\n",
      "Crop constituents\n",
      "Calculating pT\n",
      "Calculating eta\n",
      "Calculating phi\n",
      "Calculating the mass\n",
      "Pre-shifting the phis\n",
      "Using pT as weight\n",
      "Preprocessing\n"
     ]
    }
   ],
   "source": [
    "z_train = constit_to_cloud( X_train, 50, True, True, True )\n",
    "z_val = constit_to_cloud( X_val, 50, True, True, True )\n",
    "z_test = constit_to_cloud( X_test, 50, True, True, True )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train data shape: (30000, 50, 3)\n",
      "Val data shape: (30000, 50, 3)\n",
      "Test data shape: (30000, 50, 3)\n"
     ]
    }
   ],
   "source": [
    "print( \"Train data shape: \" + str(z_train.shape) )\n",
    "print( \"Val data shape: \" + str(z_val.shape) )\n",
    "print( \"Test data shape: \" + str(z_test.shape) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_train_p = torch.Tensor( z_train )\n",
    "z_val_p = torch.Tensor( z_val )\n",
    "z_test_p = torch.Tensor( z_test )\n",
    "y_train_p = torch.Tensor( y_train )\n",
    "y_val_p = torch.Tensor( y_val )\n",
    "y_test_p = torch.Tensor( y_test )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataloaders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "class topcloud_dataset( Dataset ):\n",
    "    \n",
    "    def __init__( self, data, labels ):\n",
    "        self.data = data\n",
    "        self.labels = labels\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.labels)\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        return self.data[idx], self.labels[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "trn_dataset = topcloud_dataset( z_train_p, y_train_p.unsqueeze(-1) )\n",
    "val_dataset = topcloud_dataset( z_val_p, y_val_p.unsqueeze(-1) )\n",
    "tst_dataset = topcloud_dataset( z_test_p, y_test_p.unsqueeze(-1) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "trn_dataloader = DataLoader( trn_dataset, batch_size=64, shuffle=True )\n",
    "val_dataloader = DataLoader( val_dataset, batch_size=64, shuffle=True )\n",
    "tst_dataloader = DataLoader( tst_dataset, batch_size=64, shuffle=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data handling in graph networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is handled slightly differently in graph networks than in the networks we studied previously, this is because the data now has two different attributes, the nodes and edges.  The `torch-geometric` package has it's own Data class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch_geometric.data import Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, let's create a graph with 3 nodes and 4 edges.  Remember that edges can generally go in either direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"simple-graph-image.png\" width=\"500\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_index = torch.tensor( [[0, 1, 1, 2], \n",
    "                            [1, 0, 2, 1]], \n",
    "                          dtype=torch.long )\n",
    "\n",
    "x = torch.tensor( [[-1], [0], [1]], dtype=torch.float )\n",
    "\n",
    "data = Data( x=x, edge_index=edge_index )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the data we have defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of nodes: 3\n",
      "number of edges: 4\n",
      "number of node features: 1\n",
      "number of edge features: 0\n",
      "keys: ['x', 'edge_index']\n"
     ]
    }
   ],
   "source": [
    "print( \"number of nodes: \" + str(data.num_nodes) )\n",
    "print( \"number of edges: \" + str(data.num_edges) )\n",
    "print( \"number of node features: \" + str(data.num_node_features) )\n",
    "print( \"number of edge features: \" + str(data.num_edge_features) )\n",
    "print( \"keys: \" + str(data.keys) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is awkward.  Ideally the `edge_index` would be a list of tuples, with each tuple specifying the start and end point of the edge.\n",
    "\n",
    "So we instead define it as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_index = torch.tensor( [[0, 1],\n",
    "                            [1, 0],\n",
    "                            [1, 2],\n",
    "                            [2, 1]], \n",
    "                          dtype=torch.long )\n",
    "\n",
    "x = torch.tensor( [[-1], [0], [1]], dtype=torch.float )\n",
    "\n",
    "data = Data( x=x, edge_index=edge_index.t().contiguous() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `.t()` here just applies a transpose, and `.contiguous()` ensures that the data has the correct memory format.\n",
    "\n",
    "We can again check the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of nodes: 3\n",
      "number of edges: 4\n",
      "number of node features: 1\n",
      "number of edge features: 0\n",
      "keys: ['x', 'edge_index']\n"
     ]
    }
   ],
   "source": [
    "print( \"number of nodes: \" + str(data.num_nodes) )\n",
    "print( \"number of edges: \" + str(data.num_edges) )\n",
    "print( \"number of node features: \" + str(data.num_node_features) )\n",
    "print( \"number of edge features: \" + str(data.num_edge_features) )\n",
    "print( \"keys: \" + str(data.keys) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jets, nodes, and edges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The jets we consider usually have $\\mathcal{O}(10^2)$ constituents, which would mean $\\mathcal{O}(10^4)$ edges.\n",
    "This is a lot, and it would mean we need a graph neural network with many parameters.\n",
    "The way around this is to only consider sub-graphs of neighbouring nodes within the complete graph, i.e. we only connect neighbouring nodes by edges.\n",
    "The nodes which are considered as neughbouring can be determined by any clustering algorithm, in practice the k-means algorithm is usually used.\n",
    "This funcitonality is defined in the `knn` class within the `torch_cluster` package, or in the `knn_graph` class within the `torch_geometric.nn` package.\n",
    "The latter returns a graph in the correct format, so we'll use that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we've discussed the structure of graphs and graph networks, but we haven't explained what they do with the data.\n",
    "The basic function of a graph network is the Edge Convolution, this funcitonality is coded in the `EdgeConv` class in `torch_geometric.nn.conv`.\n",
    "However this requires that we give the network the edges in the format discussed above.\n",
    "We will instead use the Dynamic Edge Convolution function, which automatically performs the k-means clustering to determine which are the neighbouring nodes.\n",
    "This functionality is coded in the `DynamicEdgeConv` class in the same package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Message passing $\\rightarrow$ generalising convolutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But let's take a step back and think about convolutions in more detail.\n",
    "\n",
    "We have studied convolutional networks in the context of jet images, where the convolutional operator scans over teh image pixels in a discrete manner.\n",
    "Now we work with point clouds, or graphs, where the data is less structured and so the idea of a convolution becomes more complicated.\n",
    "\n",
    "The generalisation of a convolution to irregular domains is usually encoded in 'message passing' (or 'neighbourhood aggregation') schemes.  \n",
    "\n",
    "The graph network is constructed of layers labelled by $k$.\n",
    "At each layer the data is again treated as a graph, so that as we move through the network the graph is being altered by each layer.\n",
    "We denote the features of node $i$ in layer $(k-1)$ of the graph network as ${\\bf x}_i^{(k-1)}$, and the edge features from node $j$ to $i$ as ${\\bf e}_{j,i}$.\n",
    "The message passing in a graph network is then expressed as:\n",
    "\n",
    "$\n",
    "{\\bf x}_i^{(k)} = \\gamma^{(k)} \\left( {\\bf x}_i^{(k-1)}, \\square_{j\\in\\mathcal{N}(i)} \\phi^{(k)}\\left( {\\bf x}_i^{(k)}, {\\bf x}_j^{(k)}, {\\bf e}_{j,i} \\right) \\right)\n",
    "$\n",
    "\n",
    "where $\\gamma$ and $\\phi$ are some differentiable functions, i.e. neural networks, and $\\square$ is a permutation invariant function like a sum, mean, or a max.\n",
    "$\\mathcal{N}(i)$ is the function which defines the nearest neighbours, i.e. k-means or something similar.\n",
    "\n",
    "The `MessagePassing` class in `torch_geometric.nn` is the base class for the graph convolution functions.\n",
    "To define a message passing layer the user needs to define the networks $\\phi$ (`.message()`) and $\\gamma$ (`.update()`), and the aggregation function $\\square$ (e.g. `aggr=\"mean\"`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic edge convolutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Edge convolutions are just a specific message passing choice.\n",
    "\n",
    "The choice is:\n",
    "\n",
    "$\n",
    "{\\bf x}_i^{(k)} = \\text{mean}_{j\\in \\mathcal{N}(i)} h_{\\Theta}\\left( {\\bf x}_i^{(k-1)}, {\\bf x}_j^{(k-1)} - {\\bf x}_i^{(k-1)} \\right)\n",
    "$\n",
    "\n",
    "where $h_{\\Theta}$ is just a neural network.\n",
    "\n",
    "We could simply load the classes for `EdgeConv` and `DynamicEdgeConv`, but let's have a look at them here.\n",
    "The edge convolutions in pytorch are implemented using the `MessagePassing` class as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Callable, Optional, Union\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "\n",
    "from torch_geometric.nn.inits import reset\n",
    "from torch_geometric.nn.conv import MessagePassing\n",
    "from torch_geometric.typing import Adj, OptTensor, PairOptTensor, PairTensor\n",
    "    \n",
    "class EdgeConv( MessagePassing ):\n",
    "    r\"\"\"The edge convolutional operator from the `\"Dynamic Graph CNN for\n",
    "    Learning on Point Clouds\" <https://arxiv.org/abs/1801.07829>`_ paper\n",
    "\n",
    "    .. math::\n",
    "        \\mathbf{x}^{\\prime}_i = \\sum_{j \\in \\mathcal{N}(i)}\n",
    "        h_{\\mathbf{\\Theta}}(\\mathbf{x}_i \\, \\Vert \\,\n",
    "        \\mathbf{x}_j - \\mathbf{x}_i),\n",
    "\n",
    "    where :math:`h_{\\mathbf{\\Theta}}` denotes a neural network, *.i.e.* a MLP.\n",
    "\n",
    "    Args:\n",
    "        nn (torch.nn.Module): A neural network :math:`h_{\\mathbf{\\Theta}}` that\n",
    "            maps pair-wise concatenated node features :obj:`x` of shape\n",
    "            :obj:`[-1, 2 * in_channels]` to shape :obj:`[-1, out_channels]`,\n",
    "            *e.g.*, defined by :class:`torch.nn.Sequential`.\n",
    "        aggr (string, optional): The aggregation scheme to use\n",
    "            (:obj:`\"add\"`, :obj:`\"mean\"`, :obj:`\"max\"`).\n",
    "            (default: :obj:`\"max\"`)\n",
    "        **kwargs (optional): Additional arguments of\n",
    "            :class:`torch_geometric.nn.conv.MessagePassing`.\n",
    "\n",
    "    Shapes:\n",
    "        - **input:**\n",
    "          node features :math:`(|\\mathcal{V}|, F_{in})` or\n",
    "          :math:`((|\\mathcal{V}|, F_{in}), (|\\mathcal{V}|, F_{in}))`\n",
    "          if bipartite,\n",
    "          edge indices :math:`(2, |\\mathcal{E}|)`\n",
    "        - **output:** node features :math:`(|\\mathcal{V}|, F_{out})` or\n",
    "          :math:`(|\\mathcal{V}_t|, F_{out})` if bipartite\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, nn: Callable, aggr: str = 'max', **kwargs):\n",
    "        super().__init__(aggr=aggr, **kwargs)\n",
    "        self.nn = nn\n",
    "        self.reset_parameters()\n",
    "\n",
    "    def reset_parameters(self):\n",
    "        reset(self.nn)\n",
    "\n",
    "    def forward(self, x: Union[Tensor, PairTensor], edge_index: Adj) -> Tensor:\n",
    "        \"\"\"\"\"\"\n",
    "        if isinstance(x, Tensor):\n",
    "            x: PairTensor = (x, x)\n",
    "        # propagate_type: (x: PairTensor)\n",
    "        return self.propagate(edge_index, x=x, size=None)\n",
    "\n",
    "\n",
    "    def message(self, x_i: Tensor, x_j: Tensor) -> Tensor:\n",
    "        return self.nn(torch.cat([x_i, x_j - x_i], dim=-1))\n",
    "\n",
    "    def __repr__(self) -> str:\n",
    "        return f'{self.__class__.__name__}(nn={self.nn})'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This requires us to provide the `edge_index` information in the forward pass.\n",
    "However as we discussed before we want to determine the nearest neighbours on the fly using dynamics edge convolutions.\n",
    "Pytorch does this as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch_cluster import knn\n",
    "\n",
    "class DynamicEdgeConv( MessagePassing ):\n",
    "    r\"\"\"The dynamic edge convolutional operator from the `\"Dynamic Graph CNN\n",
    "    for Learning on Point Clouds\" <https://arxiv.org/abs/1801.07829>`_ paper\n",
    "    (see :class:`torch_geometric.nn.conv.EdgeConv`), where the graph is\n",
    "    dynamically constructed using nearest neighbors in the feature space.\n",
    "\n",
    "    Args:\n",
    "        nn (torch.nn.Module): A neural network :math:`h_{\\mathbf{\\Theta}}` that\n",
    "            maps pair-wise concatenated node features :obj:`x` of shape\n",
    "            `:obj:`[-1, 2 * in_channels]` to shape :obj:`[-1, out_channels]`,\n",
    "            *e.g.* defined by :class:`torch.nn.Sequential`.\n",
    "        k (int): Number of nearest neighbors.\n",
    "        aggr (string): The aggregation operator to use (:obj:`\"add\"`,\n",
    "            :obj:`\"mean\"`, :obj:`\"max\"`). (default: :obj:`\"max\"`)\n",
    "        num_workers (int): Number of workers to use for k-NN computation.\n",
    "            Has no effect in case :obj:`batch` is not :obj:`None`, or the input\n",
    "            lies on the GPU. (default: :obj:`1`)\n",
    "        **kwargs (optional): Additional arguments of\n",
    "            :class:`torch_geometric.nn.conv.MessagePassing`.\n",
    "\n",
    "    Shapes:\n",
    "        - **input:**\n",
    "          node features :math:`(|\\mathcal{V}|, F_{in})` or\n",
    "          :math:`((|\\mathcal{V}|, F_{in}), (|\\mathcal{V}|, F_{in}))`\n",
    "          if bipartite,\n",
    "          batch vector :math:`(|\\mathcal{V}|)` or\n",
    "          :math:`((|\\mathcal{V}|), (|\\mathcal{V}|))`\n",
    "          if bipartite *(optional)*\n",
    "        - **output:** node features :math:`(|\\mathcal{V}|, F_{out})` or\n",
    "          :math:`(|\\mathcal{V}_t|, F_{out})` if bipartite\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, nn: Callable, k: int, aggr: str = 'max',\n",
    "                 num_workers: int = 1, **kwargs):\n",
    "        super().__init__(aggr=aggr, flow='source_to_target', **kwargs)\n",
    "\n",
    "        if knn is None:\n",
    "            raise ImportError('`DynamicEdgeConv` requires `torch-cluster`.')\n",
    "\n",
    "        self.nn = nn\n",
    "        self.k = k\n",
    "        self.num_workers = num_workers\n",
    "        self.reset_parameters()\n",
    "\n",
    "    def reset_parameters(self):\n",
    "        reset(self.nn)\n",
    "\n",
    "    def forward(\n",
    "            self, x: Union[Tensor, PairTensor],\n",
    "            batch: Union[OptTensor, Optional[PairTensor]] = None) -> Tensor:\n",
    "        # type: (Tensor, OptTensor) -> Tensor  # noqa\n",
    "        # type: (PairTensor, Optional[PairTensor]) -> Tensor  # noqa\n",
    "        \"\"\"\"\"\"\n",
    "        if isinstance(x, Tensor):\n",
    "            x: PairTensor = (x, x)\n",
    "\n",
    "        if x[0].dim() != 2:\n",
    "            raise ValueError(\"Static graphs not supported in DynamicEdgeConv\")\n",
    "\n",
    "        b: PairOptTensor = (None, None)\n",
    "        if isinstance(batch, Tensor):\n",
    "            b = (batch, batch)\n",
    "        elif isinstance(batch, tuple):\n",
    "            assert batch is not None\n",
    "            b = (batch[0], batch[1])\n",
    "\n",
    "        edge_index = knn(x[0], x[1], self.k, b[0], b[1]).flip([0])\n",
    "\n",
    "        # propagate_type: (x: PairTensor)\n",
    "        return self.propagate(edge_index, x=x, size=None)\n",
    "\n",
    "\n",
    "    def message(self, x_i: Tensor, x_j: Tensor) -> Tensor:\n",
    "        return self.nn(torch.cat([x_i, x_j - x_i], dim=-1))\n",
    "\n",
    "    def __repr__(self) -> str:\n",
    "        return f'{self.__class__.__name__}(nn={self.nn}, k={self.k})'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test it.\n",
    "\n",
    "To define the dynamic edge convolutional layer we need to supply it with a neural network ($h_{\\Theta}$) and a number of clusters ($k$) for the k-means clustering.  The dimensions of the neural network must be ( 2 x input_features, output_features ), because we concatenate $x_i$ and $x_j - x_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnt = nn.Linear( 2*3, 12 )\n",
    "dec = DynamicEdgeConv( nnt, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DynamicEdgeConv(nn=Linear(in_features=6, out_features=12, bias=True), k=5)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's pass through a random 'jet' with 40 constituents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt = torch.rand( 50, 3 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([50, 12])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dec( xt ).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the output has the dimension specified by the NN we defined."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do the same for EdgeConv.\n",
    "\n",
    "But here we need to supply edges, for this we'll just generate a random tensor of integers between 0 and 50 with shape (2,50)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnt = nn.Linear( 2*3, 12 )\n",
    "ec = EdgeConv( nnt, aggr='mean' )\n",
    "edges = torch.randint( 50, [2,10] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[18, 29, 36, 42, 25, 38, 36, 21, 37, 33],\n",
       "        [ 4,  2, 31, 47, 12,  2, 26, 29, 21,  3]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([50, 12])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec( xt, edges ).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the EdgeConv and DynamicEdgeConv return a point cloud.\n",
    "In a graph neural network each EdgeConv layer will return another point cloud, with which we will determine where to place the edges, and then pass through another EndgeConv layer to get another point cloud.\n",
    "This whole procedure is completely permutation equivariant, i.e. a change in the order of the inputs will only result is a change in the order of the outputs.\n",
    "In the end we want to aggregate over the point cloud to get a permutation invariant result.\n",
    "In this example we would simply aggregate over the dimension 50 row to get a 12 dimensional output.\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([ 0.0028,  0.0043, -0.0206, -0.0314, -0.0160, -0.0774,  0.0048,  0.0691,\n",
       "        -0.0618,  0.0452,  0.0293,  0.0545], grad_fn=<MeanBackward1>)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec( xt, edges ).mean( dim=0 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see that this is invariant we take the simpler example from before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_index = torch.tensor( [[0, 1],\n",
    "                            [1, 0],\n",
    "                            [1, 2],\n",
    "                            [2, 1]], \n",
    "                          dtype=torch.long )\n",
    "\n",
    "x = torch.tensor( [[-1], [0], [1]], dtype=torch.float )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[ 1.0610,  0.6449,  0.2553, -0.6052],\n",
       "        [ 0.2668, -0.2579,  0.5094,  0.0332],\n",
       "        [-0.5274, -1.1607,  0.7636,  0.6716]], grad_fn=<DivBackward0>)"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec2( x, edge_index.t().contiguous() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now permute the first two nodes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_p = torch.tensor( [[0], [-1], [1]], dtype=torch.float )\n",
    "edge_index_p = torch.tensor( [[1, 0],\n",
    "                            [0, 1],\n",
    "                            [0, 2],\n",
    "                            [2, 0]], \n",
    "                          dtype=torch.long )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[ 0.2668, -0.2579,  0.5094,  0.0332],\n",
       "        [ 1.0610,  0.6449,  0.2553, -0.6052],\n",
       "        [-0.5274, -1.1607,  0.7636,  0.6716]], grad_fn=<DivBackward0>)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec2( x_p, edge_index_p.t().contiguous() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we see that the first two rows have been permuted, so the result is permutation equivariant, as expected.\n",
    "\n",
    "And if we aggregate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([ 0.2668, -0.2579,  0.5094,  0.0332], grad_fn=<MeanBackward1>)"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec2( x, edge_index.t().contiguous() ).mean( dim=0 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([ 0.2668, -0.2579,  0.5094,  0.0332], grad_fn=<MeanBackward1>)"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ec2( x_p, edge_index_p.t().contiguous() ).mean( dim=0 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the result is permutation invariant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the graph neural network this would then be passed to a fully connected network and propagated to the final layer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a Dynamic Graph Convolutional Neural Network (DGCNN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "class JetNet( torch.nn.Module ):\n",
    "    \n",
    "    def __init__( self ):\n",
    "        super( JetNet, self ).__init__()\n",
    "        \n",
    "        self.net1 = nn.Sequential(\n",
    "            nn.Linear( 2*3, 32 ),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear( 32, 32 ),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear( 32, 16 )\n",
    "        )\n",
    "\n",
    "        self.net2 = nn.Sequential(\n",
    "            nn.Linear( 2*16, 32 ),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear( 32, 32 ),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear( 32, 32 )\n",
    "        )\n",
    "\n",
    "        self.conv1 = DynamicEdgeConv( self.net1, 10, aggr='mean' )\n",
    "        \n",
    "        self.conv2 = DynamicEdgeConv( self.net2, 10, aggr='mean' )\n",
    "        \n",
    "        self.net3 = nn.Sequential(\n",
    "            nn.Linear( 32, 10 ),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear( 10, 1 )\n",
    "        )\n",
    "\n",
    "    def forward( self, data ):\n",
    "\n",
    "        x = self.conv1( data )\n",
    "        x = F.relu( x )\n",
    "        x = F.dropout( x, p=0.1, training=self.training )\n",
    "        x = self.conv2( x )\n",
    "        x = F.relu( x )\n",
    "        x = x.mean(0)\n",
    "        x = self.net3( x )\n",
    "\n",
    "        return F.sigmoid( x )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's just pass some data through it to make sure this works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "testnet = JetNet()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt = torch.rand( 50, 3 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([1])"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testnet( xt ).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the dimension of the output that we expect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can view the network elements by just calling it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JetNet(\n",
       "  (net1): Sequential(\n",
       "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
       "  )\n",
       "  (net2): Sequential(\n",
       "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
       "  )\n",
       "  (conv1): DynamicEdgeConv(nn=Sequential(\n",
       "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
       "  ), k=10)\n",
       "  (conv2): DynamicEdgeConv(nn=Sequential(\n",
       "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
       "  ), k=10)\n",
       "  (net3): Sequential(\n",
       "    (0): Linear(in_features=32, out_features=10, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=10, out_features=1, bias=True)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testnet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimising the DGCNN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial I have used a standard pytorch Dataloader to store the data, just as we've done in previous tutorials.\n",
    "But to batch the data for DGCNNs constructed using `torch_geometry` we would need to use the `Data` package within `torch_geometry`.\n",
    "I have resolved this below (in an unstatisfactory way) by manually looping over each jet each batch.\n",
    "This fixes the problem, but would certainly slow the code down if we were to use a GPU.\n",
    "In practice we should use the correct Dataloaders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss_fn = nn.BCELoss()\n",
    "batch_size = 64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_epoch( dataloader, model, loss_fn, optimizer ):\n",
    "    \n",
    "    size = len( dataloader.dataset )\n",
    "    \n",
    "    for batch, (X_batch, y_batch) in enumerate( dataloader ):\n",
    "        \n",
    "        loss = 0.0\n",
    "        \n",
    "        for (X, y) in zip( X_batch, y_batch ):\n",
    "\n",
    "            # pass data through network\n",
    "            pred = model(X)\n",
    "\n",
    "            # compute loss\n",
    "            loss += loss_fn( pred, y ) / batch_size\n",
    "\n",
    "        # reset gradients in optimizer\n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        # compute gradients\n",
    "        loss.backward()\n",
    "        \n",
    "        # update weights with optimizer\n",
    "        optimizer.step()\n",
    "        \n",
    "        # print the training loss every 100 updates\n",
    "        if batch % 100 == 0:\n",
    "            #loss, current = loss.item(), batch * len(X)\n",
    "            current = batch * len( X_batch )\n",
    "            print( f\"current batch loss: {loss:>7f}  [{current:>5d}/{size:>5d}]\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def val_pass( dataloader, model, loss_fn ):\n",
    "    \n",
    "    size = len( dataloader.dataset )\n",
    "    num_batches = len( dataloader )\n",
    "    vl = 0.0\n",
    "\n",
    "    with torch.no_grad():\n",
    "        for X_batch, y_batch in dataloader:\n",
    "            for (X, y) in zip(X_batch,y_batch):\n",
    "                pred = model( X )\n",
    "                vl += loss_fn( pred, y ).item() / batch_size\n",
    "\n",
    "    vl /= num_batches\n",
    "    print( f\"avg val loss per batch: {vl:>8f}\" )\n",
    "    \n",
    "    return vl\n",
    "\n",
    "def trn_pass( dataloader, model, loss_fn ):\n",
    "    \n",
    "    size = len( dataloader.dataset )\n",
    "    num_batches = len( dataloader )\n",
    "    tl = 0.0\n",
    "\n",
    "    with torch.no_grad():\n",
    "        for X_batch, y_batch in dataloader:\n",
    "            for (X, y) in zip(X_batch,y_batch):\n",
    "                pred = model( X )\n",
    "                tl += loss_fn( pred, y ).item() / batch_size\n",
    "\n",
    "    tl /= num_batches\n",
    "    print( f\"avg trn loss per batch: {tl:>8f}\" )\n",
    "    \n",
    "    return tl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can train the network.\n",
    "Graph networks can be slow, especially on a laptop CPU.\n",
    "So I only train it here for 10 epochs.\n",
    "But you'll see that even with 10 epochs we get quite a good performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------\n",
      "model architecture \n",
      "-----------------------------------------------\n",
      "JetNet(\n",
      "  (net1): Sequential(\n",
      "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
      "    (1): ReLU()\n",
      "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (3): ReLU()\n",
      "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
      "  )\n",
      "  (net2): Sequential(\n",
      "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (1): ReLU()\n",
      "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (3): ReLU()\n",
      "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
      "  )\n",
      "  (conv1): DynamicEdgeConv(nn=Sequential(\n",
      "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
      "    (1): ReLU()\n",
      "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (3): ReLU()\n",
      "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
      "  ), k=10)\n",
      "  (conv2): DynamicEdgeConv(nn=Sequential(\n",
      "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (1): ReLU()\n",
      "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
      "    (3): ReLU()\n",
      "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
      "  ), k=10)\n",
      "  (net3): Sequential(\n",
      "    (0): Linear(in_features=32, out_features=10, bias=True)\n",
      "    (1): ReLU()\n",
      "    (2): Linear(in_features=10, out_features=1, bias=True)\n",
      "  )\n",
      ")\n",
      "-----------------------------------------------\n",
      "Epoch 1\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.688897  [    0/30000]\n",
      "current batch loss: 0.643746  [ 6400/30000]\n",
      "current batch loss: 0.562593  [12800/30000]\n",
      "current batch loss: 0.575729  [19200/30000]\n",
      "current batch loss: 0.502195  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.557880\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.547772\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 2\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.461260  [    0/30000]\n",
      "current batch loss: 0.501755  [ 6400/30000]\n",
      "current batch loss: 0.514454  [12800/30000]\n",
      "current batch loss: 0.577440  [19200/30000]\n",
      "current batch loss: 0.595107  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.537895\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.525343\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 3\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.556473  [    0/30000]\n",
      "current batch loss: 0.635104  [ 6400/30000]\n",
      "current batch loss: 0.375902  [12800/30000]\n",
      "current batch loss: 0.259500  [19200/30000]\n",
      "current batch loss: 0.253679  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.306987\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.302586\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 4\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.375568  [    0/30000]\n",
      "current batch loss: 0.383267  [ 6400/30000]\n",
      "current batch loss: 0.441775  [12800/30000]\n",
      "current batch loss: 0.237086  [19200/30000]\n",
      "current batch loss: 0.342605  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.278998\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.276466\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 5\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.271998  [    0/30000]\n",
      "current batch loss: 0.382415  [ 6400/30000]\n",
      "current batch loss: 0.288617  [12800/30000]\n",
      "current batch loss: 0.299773  [19200/30000]\n",
      "current batch loss: 0.276738  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.269817\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.266514\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 6\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.292407  [    0/30000]\n",
      "current batch loss: 0.279011  [ 6400/30000]\n",
      "current batch loss: 0.345525  [12800/30000]\n",
      "current batch loss: 0.269051  [19200/30000]\n",
      "current batch loss: 0.278432  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.264490\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.263246\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 7\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.319298  [    0/30000]\n",
      "current batch loss: 0.333138  [ 6400/30000]\n",
      "current batch loss: 0.232179  [12800/30000]\n",
      "current batch loss: 0.285207  [19200/30000]\n",
      "current batch loss: 0.220867  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.255705\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.254744\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 8\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.227725  [    0/30000]\n",
      "current batch loss: 0.140448  [ 6400/30000]\n",
      "current batch loss: 0.386471  [12800/30000]\n",
      "current batch loss: 0.353087  [19200/30000]\n",
      "current batch loss: 0.258750  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.252674\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.251380\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 9\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.352660  [    0/30000]\n",
      "current batch loss: 0.259513  [ 6400/30000]\n",
      "current batch loss: 0.252328  [12800/30000]\n",
      "current batch loss: 0.140094  [19200/30000]\n",
      "current batch loss: 0.196019  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.302173\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.301284\n",
      "-----------------------------------------------\n",
      "|\n",
      "-----------------------------------------------\n",
      "Epoch 10\n",
      "-----------------------------------------------\n",
      "current batch loss: 0.326917  [    0/30000]\n",
      "current batch loss: 0.210762  [ 6400/30000]\n",
      "current batch loss: 0.271036  [12800/30000]\n",
      "current batch loss: 0.060773  [19200/30000]\n",
      "current batch loss: 0.139073  [25600/30000]\n",
      "-----------------------------------------------\n",
      "avg trn loss per batch: 0.249912\n",
      "-----------------------------------------------\n",
      "avg val loss per batch: 0.247864\n",
      "-----------------------------------------------\n",
      "|\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# a useful function to present things clearer\n",
    "def separator():\n",
    "    print( \"-----------------------------------------------\" )\n",
    "\n",
    "epochs = 10\n",
    "\n",
    "# re-initialise the model and the optimizer\n",
    "model = JetNet()\n",
    "learning_rate = 0.001\n",
    "optimizer = torch.optim.Adam( model.parameters(), lr=learning_rate )\n",
    "separator()\n",
    "print( \"model architecture \")\n",
    "separator()\n",
    "print( model )\n",
    "\n",
    "# track train and val losses\n",
    "trn_losses = []\n",
    "val_losses = []\n",
    "\n",
    "for t in range(epochs):\n",
    "    separator()\n",
    "    print( f\"Epoch {t+1}\" )\n",
    "    separator()\n",
    "    train_epoch( trn_dataloader, model, loss_fn, optimizer )\n",
    "    separator()\n",
    "    model.eval()\n",
    "    trn_loss = trn_pass( trn_dataloader, model, loss_fn )\n",
    "    trn_losses.append( trn_loss )\n",
    "    separator()\n",
    "    val_loss = val_pass( val_dataloader, model, loss_fn )\n",
    "    val_losses.append( val_loss )\n",
    "    model.train()\n",
    "    separator()\n",
    "    print( \"|\" )\n",
    "    \n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Study the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JetNet(\n",
       "  (net1): Sequential(\n",
       "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
       "  )\n",
       "  (net2): Sequential(\n",
       "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
       "  )\n",
       "  (conv1): DynamicEdgeConv(nn=Sequential(\n",
       "    (0): Linear(in_features=6, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=16, bias=True)\n",
       "  ), k=10)\n",
       "  (conv2): DynamicEdgeConv(nn=Sequential(\n",
       "    (0): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=32, out_features=32, bias=True)\n",
       "    (3): ReLU()\n",
       "    (4): Linear(in_features=32, out_features=32, bias=True)\n",
       "  ), k=10)\n",
       "  (net3): Sequential(\n",
       "    (0): Linear(in_features=32, out_features=10, bias=True)\n",
       "    (1): ReLU()\n",
       "    (2): Linear(in_features=10, out_features=1, bias=True)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avg trn loss per batch: 0.249912\n",
      "avg val loss per batch: 0.247864\n",
      "avg val loss per batch: 0.249365\n",
      "-----------------------------------------------\n",
      "Final evaluation of performance:\n",
      "BCE-loss on train dataset: 0.24991207184104408\n",
      "BCE-loss on validation dataset: 0.24786406089310592\n",
      "BCE-loss on test dataset: 0.2493650476925356\n",
      "-----------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "trn_loss = trn_pass( trn_dataloader, model, loss_fn )\n",
    "val_loss = val_pass( val_dataloader, model, loss_fn )\n",
    "test_loss = val_pass( tst_dataloader, model, loss_fn )\n",
    "\n",
    "separator()\n",
    "print(\"Final evaluation of performance:\")\n",
    "print(\"BCE-loss on train dataset: {}\".format(trn_loss))\n",
    "print(\"BCE-loss on validation dataset: {}\".format(val_loss))\n",
    "print(\"BCE-loss on test dataset: {}\".format(test_loss))\n",
    "separator()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "samp_size = 10000\n",
    "\n",
    "train_pred = []\n",
    "for z in z_train[:samp_size]:\n",
    "    train_pred.append( model( torch.Tensor(z) ).detach() )\n",
    "\n",
    "val_pred = []\n",
    "for z in z_val[:samp_size]:\n",
    "    val_pred.append( model( torch.Tensor(z) ).detach() )\n",
    "\n",
    "test_pred = []\n",
    "for z in z_test[:samp_size]:\n",
    "    test_pred.append( model( torch.Tensor(z) ).detach() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pred_lab = torch.round( torch.Tensor( train_pred[:samp_size] ) )\n",
    "val_pred_lab = torch.round( torch.Tensor( val_pred[:samp_size] ) )\n",
    "test_pred_lab = torch.round( torch.Tensor( test_pred[:samp_size] ) )\n",
    "\n",
    "train_correct = ( train_pred_lab == torch.Tensor( y_train[:samp_size] ) ).sum().item()\n",
    "val_correct = ( val_pred_lab == torch.Tensor( y_val[:samp_size] ) ).sum().item()\n",
    "test_correct = ( test_pred_lab == torch.Tensor( y_test[:samp_size] ) ).sum().item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------\n",
      "Evaluation of accuracy: \n",
      "Accuracy on training dataset: 89.9 %\n",
      "Accuracy on validation dataset: 90.3 %\n",
      "Accuracy on test dataset: 89.9 %\n",
      "-----------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "separator()\n",
    "print('Evaluation of accuracy: ')\n",
    "print('Accuracy on training dataset: {:.1f} %'.format(train_correct*100/len(train_pred)))\n",
    "print('Accuracy on validation dataset: {:.1f} %'.format(val_correct*100/len(val_pred)))\n",
    "print('Accuracy on test dataset: {:.1f} %'.format(test_correct*100/len(test_pred)))\n",
    "separator()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pred = np.array( train_pred )\n",
    "val_pred = np.array( val_pred )\n",
    "test_pred = np.array( test_pred )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve, roc_auc_score\n",
    "\n",
    "def closest_point(array, tpr_p=0.3):\n",
    "    dist = ((array-tpr_p)**2)\n",
    "    return np.argmin(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAFMCAYAAAD89+yjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAACVrklEQVR4nOzdd1hUV/rA8e+hKooMYI19MGpMTCJgetE4pG6yKaC/7G56hGzfbIGY7ZVgdrMtuxswyW52s0WZVNMZE9KLMKbaIqMxJlZwFFSkzPn9cWeQMsAAM9wZ5v08D4/Mre9ckLnvPee8R2mtEUIIIYQQQgjRvRizAxBCCCGEEEKIcCeJkxBCCCGEEEL0QhInIYQQQgghhOiFJE5CCCGEEEII0QtJnIQQQgghhBCiF5I4CSGEEEIIIUQv4swOYLCMHj1aT5s2rd/7Hzp0iBEjRgQvoAgT7e8f5BpE+/sHuQYw8GtQXV29T2s9JoghDRn9/ZyS30u5BiDXAOQagFwDCO3nVNQkTtOmTaOqqqrf+1dWVrJgwYLgBRRhov39g1yDaH//INcABn4NlFKfBC+aoaW/n1PyeynXAOQagFwDkGsAof2ckq56QgghhBBCCNELSZyEEEIIEymlLldKlR04cMDsUIQQQvRAEichhBDCRFrr1Vrr/JSUFLNDEUII0QNJnIQQQgghhBCiF5I4CSGEEEIIIUQvTK+qp5SyAPlAuta6KIDtCwEXkAagtS4LaYBCCCFECCmlLgcunzFjhtmhCCGE6IGpLU5KKRtgAzIASwDblwAurbXdmzBlKKVyQxulEEIIEToyxkkIISKDqYmT1tqhtbYD7gB3yfdu77MSKAh6YEIIIYQQQgjRTsSMcVJKZfpZ7MZosRJCCCGEEEKIkDF9jFMfpAF1nZZ1fi1EwJpaPBxuaqHVo40vrWlpNb5vONpCc6sHjwatNR4NG+taSajZh9bg0brLv3vrjxIfG0OLx0Nzq6al1cPBxhb2H25ieHxsh/P4vv+k9jDJw+KIUapDbBrtN2btf3H3y/0cp/ttez52bW0jD39S1eMefT92N++zhzgGemz/23Z37I4r9u8/wootb/Xp59Dnn2W3QXZdtHlPPRNShqO6rgq6mS2biNfNTJs2kwWDcD7RNzf9/R32NTTxo8tO4HRrutnhiF4cPXqUuro66uvraW1tDcoxU1JS2LBhQ1COFankGsg1AOMabNq0iWHDhjFy5EhSU1OJiQlOW1EkJU6W7lYopSxaa7ef5fkYhScYN24clZWV/T55Q0PDgPaPdGa/f601LRoOHtUcOKo53KJp9sC+IxoFtGqMhKcZDhzV7D1iJD2tGlo8UNfoodVjHMuDscwT+H31Me+83e/3EB8DMarjlwIONsFxI/zc+nZzN9zdTXJfbp6V6tuttgI8nlZqj+ztdJy+HaMvy/2t6GuC0O05+xl3q6eVPfv2D+gYAW0f4M9+XCLU19czJin0nQcmt3zCF5qe4bmjP4vqv4WhEIziEOkjE3lp017e2VoniVOYO3r0KNu3byc1NZVp06YRHx/f57/J/tTX15OcnByECCOXXAO5BgAHDx5kxIgRHD58GLfbzcGDB5k8eTJxcQNPeyIpcXLjraTXTufXHXgLSJQBZGdn6wULFvT75JWVlQxk/0gXzPd/pKmVHfsPc7TFw9GWVnbsP8Ln7kY8WrN13yE+qT2E1nC0xcOm3fU0t3q6fTrfneRhcXg8muxpaSTExZAQG0PdoSaOHzeSYfGxxMcqtIYYpRg7KpEYpYiLUcTEGP/GxihaWjWjkxO9SY7i/ffeY968U1FATIwyEh+liFGqbZuEuBhGJsYRF6uIi4khLlaRGBdDYlxsUK6dmaL9/wBE4TU4Wg+fVYN1AXAxNC1j5xvvRNc1GARa69XA6uzs7KX9PcZdV8/FXr0jiFGJUKmrqyM1NZXRo0ebHYoQQ5JSitjYWJKTkxk5ciQ7d+6krq6OsWPHDvjYkZQ41dG11ckC4K+1SZhHa83RFg879h+msdnDvoajOLe7cazfTe2ho+w+eLTXY4wemcDciSlMtIzlYGMz2VNTSYyP5WiLh1njkomLVUxIGcaw+FhGJMYxMjGOhNgY4mONpCcYT+86a94RyxnyJFdEiz0bYNX1cGAHfOcDGDEaEpLMjkqIiFdfX8+0adPMDkOIqKCUIj09ne3bt0dX4qS1diql3J0WpwEOE8KJeoebWqhYv5uDjS3sqz/KMx/spKnVw+GmVvbWd58YzZ2YwtWZkxiREMvMcckkxMUQG6OYnJrEmOREkhJiQ5L0CCH64P1VsPrbkDACvrTSSJqEEEHR2tpKfHy82WEIETUSEhJoaWkJyrHCOnFSSlmBzHYlyFcppXLbvc4BSs2JLro0tWqeev9zXtq4l8/ch3nL1bUux9T0JC6YNZa4WEV8bAxpIxKYPT6ZxPhYUobHc8qkFEmKhAhnWsPT34OqB2DKWZD7IIyaYHZUQgw58lkoxOAJ5v83UxMnb4lxG5DrfV0IOLTWTu8muRjJkR1Aa12glCr0TpxrBWo6zeskgsR9uIl3P3WzfudB1n9+kKfePwysa1tvO2Ecs8cnc1XmRMaNGsbw+FhiY+SDQIiIphQMGwVnfQsW/RRiw/rZmhBCCDGoTP1U9CZITmB5N+uXd17nXSZC5HP3EX742Ae8tKlj9bTkBLjp3OO59rTJjB81TJ6WCTGUbHoWhqfClDOMhEn+fwshhBBdyONEwduuWv77znbWbNhD/dFjfUAXZ0/ilnOsTE1P4q3XX2XBgpkmRimECLrWFnjxl/D6H2DmxcZ4JkmaBl0wypELMZQ5nU4yMzPNDkMISZyi1Qc7DnD3C5t4ZfOxlqWp6UlcdNJ4lsyfzPxpPVZ6F0JEuvpdYL8FPnkNsm6Ci+8yO6KoFYxy5EIMVW63m6ysLPbv34/FYul2u6KiItauXctLL71EZmYm2dnZlJZ2HAafl5eHw2HUFFu8eHGX9T5lZWWUl5djsVhIS0tr29dms+F2u1m1ahX5+fnBeYP94Ha7KS4uJiMjA4CamhpKSkoC3r+oqKjDa3/7ZmRkUFRUxOLFiwFwOBxUVFRQUlLS489hqJPEKYporXFu38+3/vsun7mPtC0/Z8ZofvHFE7GOGWlidEKIQbN/GzxwITQehKtK4ZT/MzsiIYTwa9WqVYCRzBQWFna7XUlJCe+//z6nnHIKJSUl2Gy2LtuUl5dTVlaG1Wr1u97lcrUlSBUVFR3WORwO7HY7FRUVbQmLWfLy8igtLcVqtQJG3Dk5OV1i9icjI4Py8vK2Fjy73e53X5fLRUFBAQUFBQBYrVYqKiqiOmkCCP108yIsbNx1kOnLnuGav73JZ+4jfOn0KTz/nfPYdtdlPHzr6ZI0CRFNUqbACVfA0hclaRJChDW3201ubi4rV64MyvGsVmtbwtGey+UiKyuLkpISvy0wNpsNq9VKWVlZUOLoL7vd3uU9+L6323uul1ZUVERmZmaHbo+5ubm4XK62ljifkpISqqurqaiooLq6mpqaGr/XLdpI4hQFlj+3kYv/8CoAcyaM4omvn81vrprLrPHJJkcmhBg0R/bDY7cZE9rGxMBlv4Vxc8yOSgghuuVyubBarSxZsgSn04nL5QrZuQoKCli8eLHfliifzMxMU7voAaxcuZKsrKwuy3NycrrteuhTVlZGTk5Ol+U2m43y8vIuyzMzM7HZbDK+rB3pqjeEvbJ5Lzf+/R082nhdel0WF5043tyghBCD7/N1sOp6OLjTKAKRMsnsiIQQolcOh6MtUbFYLJSWlvZpLE+g7HY7DoeDmpqaXrfNy8vD6XT2ul2oOByOtu5z7VmtVqqqqrrdz+1243a728ZstZeRkdFr0iUM0uI0RJVXfcr1DxpJ05wJo3jxe+dL0iREtNEaqh40xjNpDTc/BydeaXZUQggRELfb3fb94sWLe+2K1l++8UKBdEXzddkzQ0/Jj8Vi6XC9uuNvjJLFYum2Nc/pdHbpxhfNJHEagu545H1+YH8fgGWXzOaZb58rY5iEiEZr74enbofp50PBKzAp2+yIhBAiIC6Xq0MXsYKCAlwuV0hae6qqqvqUDOXm5gY9hkDU1dX1uk13yZMvYfK33res/bq1a9dit9uxWCzYbDaKiopMH98VDqSr3hDzy6fW87+1nwLw8C2nc87xo02OSAgx6LQ25mM6eQloD8xfaoxrEmFJ5nESoqv23fTAGG9jtVpZuXJl0MfcuN3uoFeLy8vLC6gFqL2CgoKQJmW5ublUVFR0OcfatWv9xtJ+vFdJSQkZGRndViSMFpI4DSFNLR4eeG0rAFU/sjF6ZKLJEQkhBt0HdqOl6brHYNgoOL1rX3gRXmQeJyECk5ubS1lZWdDHOQXaza0v/BVbMNuKFSu6FJZwOp1trW3tk0d/yVFubi5FRUVUV1eHNM5wJonTENHY3ErO718G4ILZYyVpEiLatByF538Ia1fA5DPgaAPEDzc7KiGE6DOn00l1dbXfIghutxuHwxHUVo/s7Ow+VezzVfsbbP7GNnXWU8uZxWKhurqa5cuXY7VaqaurIzs7u60lqTcZGRksX768LyEPOZI4DQEtrR7OW/4Se+qPctr0NO6/XsYxCBFV3Nth1Q3wuRPO/AbYfgax8WZHJYQQ/eJwOLqt8uZwOCgvL/ebOE2fPh3ofpwPGElPdnbH+6S8vDwKCgoC6rLndrs7tNIMJovF0lbIoXN3RZfLFVB3Q4vF0mUi4ZUrV3Z4PxkZGZSUlJg2liucSeI0BDg27GZP/VGso0ewquBMs8MRQgy2J74OtVtg8b9gzhVmRyNM9LuKzXx94QxiYpTZoQgREgUFBRQVFXWbWFmt1h5bj/wlR/n5+ZSWllJcXNxrN8BAW7tCNcYpOzvbb5GImpqafrfCORyODu/bYrH4HUdWU1MT9XM6SeIU4T6pPcRtDxsVZp785jkmRyOEGDSeVqN7XkISXP5HoyBEeobZUQmTxMYoRg2L42BjC/sPN5Eu3bVFhHI4HD0mD75xNna73e92JSUlFBUVdWlVgZ6LQJSXl5OVlUVOTk63CYjb7aauri6glp1QjXHKy8ujoqKiy0S8DoeDZcuW9bivr+tj+6TTV2q8/XtesmSJ3xY1u91OUVFRv2MfCqTMUgTbeeAI599dCcAPLprFyETJg4WICg174F9XwhNfMxKmNKskTVFOKcX3L5pldhhCDFhJSUmP3eB88y111+KUm5uLzWbrMj7K5XJRXFzcJeFof9zq6mqKior8Jgcul4uysrJu9x8s+fn5uFyuDq1qTqeTtLS0LolkXl5elzmY8vLy2r53u92UlJR0SfJsNluX0uO+cVFmv3+zyZ12BPvfO0bZ8R9eegJLzzNnMjYhxCD75A0ovwka3XDZ74yy40IIEeEcDkfbXE1ZWVmsWbPGb8uObxuXy9U2NqlzC1FpaSl2u528vDzS0tKwWCykp6f32g3PlzyVlZWRk5ODxWLBarWSnp6O1Wr124plhjVr1lBcXExGhvHArKamhoqKii7bORwOcnJy2l6XlJRQXFzcYdvy8vIu19nXHc+XQLrdbjIyMvyeI9pI4hShmls9/PnFj4mNUVx35lSzwxFChJrW8MafwfEzSJ0KX7HD+LlmRyWEEEFhs9moqanpdbvS0tJuW5vay83N7Xdxg/z8/LBuWbFYLAGVZN+/f3+/9gMjeYr28Uz+SOIUoX799AY8Gm5fdDzD4mPNDkcIEWqH9sJrv4fZl8EX74VhKWZHJIQQQkQVSZwi0N76o/zjjW0AfOWMKeYGI4QIrX1bjDFMI8dCfiVYpkj3PCGEEMIEUhwiAn131bsA/DbvFKmcJMRQpTVUPwR/Owvevs9YljpVkiYhhBDCJNLiFGFcext49eN9AORmTTI5GiFESDQdhqe/B+/9B6wL4eTFZkckQkgpdTlw+YwZM8wORQghRA+kxSnC3PJQFQD/WXq6yZEIIUJi3xa4fxG8919YsAy+8giMGG12VCKEtNartdb5KSkybk0IIcKZtDhFkOseeJut+w4xd2IKZ2XIjZQQQ9KhvXC41kiYZiwyOxohxCD7+eqPWP/5wT7v19raSmysOcWi5hw3ip9efqIp527P7Xazd+/eHueBMpPL5RrU2Ab7fNFAWpwixKd1h9u66P03/wyToxFCBFVLE2x6zvh+6pnw7fckaRJCiD665557cDqdZofh1/Lly0lLSxvUczqdzi4T4IqBkRanCPGt/60DYMX12YxMlB+bEEPGgR1QfiPsqIKvvQVjZ0P8cLOjEhEoxls4JOtXHW+URibG8ey3z2VyWpIZYYk+6m/LTX19PcnJyUGOxhx2u52VK1dSXl7ep/0ef/xxfv/734coqv4rKyvDZrN1mWjW7XZ3mcg20HmWXC5Xh/msXC4XJSUljBkzpm1Zbm4uRUVFWK3WsGt58k2u6+PvffuuT3p6OrW1tbjd7rb3Yxa5A48Ab9TsY912NxMtw8mZM87scIQQwbLFAY8shdZmyPuHkTQJ0U8XnzSe/YeaaPbotmWf1h3msXWfsetgoyROIuwVFBQARhJQV1fXp32dTiennnpqCKIaGLfbTXV1td8JdfPy8igtLW1LBFwuFzk5OVRUVPR4TJfLhd1u75Bs2O12srKyePnllzn55JPbli9btoylS5f2OQkNpYyMDMrLy9sm2LXb7V3ety9p6pxQ5eTkdLhmg0266kWAB1/bCsDKAumiJ8SQ8crd8HAuJE8w5mc68UqzIxIRbvTIRL656Hi+mzOz7euaTKm+KiJHaWkppaWl5OXl9Wvfm266KQRRDUxRUVGX1hUwkoXOLUG+7+12e4/HbN/S5JObm4vb7ebvf/97h+UWiwWr1Ro2XfaKiorIzMxsS5rAiN3lcnWIsbi4mGXLlnXZv6SkxO/7HyySOIW5hqMtODbsAWBSqjwtFGLISEqHU78EtzpgtJShFkKIgXA4HCxcuNDsMLrorkDDypUrycrK6rLc16LSm5UrV3ZZZrFYOHDgQJflBQUFpiYb7ZWVlZGTk9Nluc1m69Aq5nK5cLlcXbazWCy43e5QhtgjSZzC3KLfVQJw7WlTzA1ECDFw29+GDU8Z32fdBF/8CyTIA5FQU0pNV0rdqpSaZnYsQojgczgc2Gw2s8PowtcFzR+Hw+E3obJarVRVVfV43JKSEqqrqzssc7vduN1urrzySr/HDIeiGb4Y/RXJyMjI6NDiZLVaWbp0aZckyW6396tFMlgkcQpjJc9tZPfBowD85qqTTI5GCNFvWsObf4F/XAqVd4HHA0oZXyLktNZbtdb3A+F3ZyWEGLDS0tK28VHdWb58OXl5eeTl5fW6bbBUVFR06JLm01MC0d8WlaVLl5Kfn99tq5vNZguL5AnoUiTDt6x9C1NJSQkul4vp06e3JVS++M1MkiVxCmN/q6wB4N2f5KDkBkuIyNR4AFZ+BZ6/E2ZeDDc+BTHypzdYlFKn9qElKSOUsQghzOF0Ov0mKGB0+crIyMBisVBeXk55efmgdVurqqoiOzu7y/JACl8Ekjw5nc62hHDJkiU9vq/OLTpm8CVM/t6bb1n7dVu3bsVqtZKTk0NOTg51dXUUFhaGPtAeSFW9MPXch7sAGDcqEUtSgsnRCCH6pfEAlC2A/Z/Ahb+CM78hrUxB4k2WqgGL93Wp1vprndZnAvMBq/f78HjcKoQIGrvdTm5ubrfrc3JyyM3N9VvVLtTcbrff1pVg8RVZ8JUir6ur49prr/W7rdVq7bVaX3t5eXl9bvkqKCjo8WcBRiGIioqKLtutXbu2y7YWi4UlS5aQnZ1NWVkZLpeLiooKKUcuOmpu9XDbw0bf1X/cdJrJ0Qgh+m1YCpyUCxkXGBPbimByAg6gAkgF8pVSt2itH1BK3Qd0vkty+FkWUkopG1CH0UXQobU2LXF771M386cN7uSbQgwGXyU+f4qKinC5XH6rs/k4HA4KCgooKCjo0ppRVlZGUVERa9as6bZFqyd9LaneX1arldLSUlJTU9m9ezc//vGPu2xjsVj6FE+oypevWLGiS1EMp9PZlgy1b5UqKipq+9kWFBSQl5dHVlZWv38ewSCJUxj605qPAVg4awwnTBhlcjRCiD5pPgLPLYOsG+G4U+GCH5od0ZCjlPo+sFRr/Ui7xcuVUquUUvuBNCALcAForbuWmQp9jFagSGudo5RKA5YBgz6iecbYkQDs2H9ksE8tRMi53W7q6uq6bYGw2+1YLBaWLl3atiwtLa1DomWz2bBarX7HzeTn57eVzw4mf2ObOutPS1V+fj4/+clP/CZOEFj3v1CzWCxUV1ezfPlyrFYrdXV1ZGdnk5GR0eHnmJeX1yF5y8zMpKampm2cWk1NjRnhS+IUbj5zH+HPL24BYHnuKSZHI4Tok9oaWHUD7P4Axsw2EicRCvO11r/1s7wUKNRaXzTYAXWmtXYBvnJaVqBrP5RBMD5lGKOGyUe9GJpWrVrVY6EHl8tFfn5+r2Oaqqqq/CZH7VtC+qO7RMVisbQVQ+h8XpfL1WPS5Ha7ycvLo6SkpMu+6enpbcformJfOLBYLF1a91auXNkWn6+Lo7/rUF5eTkZGRrfvMdTkr2mYebR6BwCFF89iTHKiydEIIQK2/kl44usQEwtftsPx/kvQiqDY383yKowueUGjlLJgdPFL11p3mcVSKVWI0bKVBqC1Luu0PhfI8LevEGJgSktLWbNmTbfru7v5bq+n5GigZc57Ond2drbfrnM1NTU9ntM3Uay/pKu2thbw36LV1/FWoRrj1B2Hw0FJSQlgvMeeWuXM6qYHkjiFnd9VbAYg/9zweCoghAjApudg1XUwMQvyHgLLZLMjGuq034VaH1BKdZ0xsZ+8Y5QsdFONTylVAqzVWtt9r5VSub7X3pjsSimrUsqmtTa3pJUQQ4jv5rqnZMBms/V68++bT8lut7e1Yvhu/isqKigqMp552O32tm6B7RMbX5czMBK59gUY0tLSuk1Y8vLyqKio6FK0wuFw9DgmKzMzk8LCQr8Jim8SYH/nc7lcbS1SgQjVGCdfC2H7VkBftT/fdc3MzOxxLiu3221a65nUxA0jVduMJw9T05OIi5UfjRBhT3vv32fY4OK74KbnJGkyn9+kCkApVdynA2nt8CZB7m42yW+fJAErgQLvuSze1iowWsEGp/6xEEOAb56jntjt9l7nYyopKWHVqlVdjtV+vqCKigrmz59Pbm4uhYWFbYkSHCsn7mt5qqmp6TAXUk5ODjabjdzcXL83+zabrdsEID8/H5fL1SEWp9NJWlpal6QoLy+vQynx+fPnU1bWoXG7rRXqD3/4g9/z1dTUmNpS0177CWzdbjclJSVdErWCgoIOPwuf5cuX+10+WMKixam3rg7dbO/2vrRorZeHNMBBUvaK8Z/nD0tONTcQIUTval6Eip/CVx6FkWPgjK+aHVE0yVZKLcR/l7353bQ6WfGWLg8GpZS/OxA3xybZzQfSgSLvctNK2rV4NA1HW8w6vRABKyoqwu12tyU7eXl5pKWlUVBQ0OWmf+XKlVRXV/d4PKvVSnV1NUVFRW1zObnd7g6tPFVVVR1aierq6nC73W0tWr7kyWKxtHUlAyNxg2PdxlwuV5cudjk5OTidzm673q1Zs4bi4mIyMoxG7ZqaGr8lwx0OBzk5x7p/5+bm4nQ6KSoqIj09ndraWlwuF1u3biU2NtbvuaqqqgZt/qqelJSUUFxc3OF9lpeXd2kly8/Pb6t4aLFY2t5nQUFBdJcjD6SrQ6ftC9snSkqpzM7LIpHWmhfW7wbglEkWc4MRQnTP44FX7obKYhgzC44eNBInMZiyMFpxupsUy9/jSA30+FCuj9IwSo231/51GWDzdvfLw4SKej6Hm1qp/qS7YWFChA9fYtLbDX5fijb4SnX707nAgK/1x2KxUFVVhcViwel0thWZaG/lypUdWk6cTmeH5AaMBCcnJ6fbSVs7J2Pd2b+/6/9f3xxOndXX13dZZmbXts4Cfc9gtNgNZIxZKJieOGF0dWj/IbcSKAH8Jk7AEqAtSdJaO5VS3XcGjRB/ecmopHfRieOIiZEJMoUIR/FNB+HfuVCzBk5eAl/4PSSMMDusaOSk74lIKrC0160CZ+luhVLKorV2c+xzrNuxTUqpfLzzS40bN47Kyso+B9LQ0NDjfuOSFLEtR/p17EjR2zUIJykpKX5vbgeqtbU1JMcNR/feey/XXXddl/fb12vw4Ycfcsopp7Tt8+9//5urrrqK+vp6nn32Wb7zne9w5ZVXct555zF+/HgWLlzYNmZpy5YtfOtb32rb97nnnuPnP/95l/NPnjyZ999/n+nTpw/wXQfG3zX485//7Pd6DVX+rkFjY2NQ/kaYmjgF0NXBnzqlVLnWOs97jHyMZCtiaa3599vbAbg7T0qQCxGurK5/wN5XjYQp6yZQ8pDDJA6t9dY+7rM1mIUj8N/9rs/d8bxd08sAsrOz9YIFC/ocSGVlJT3tN2PLWxxt9rBgwVl9Pnak6O0ahJMNGzaQnJwc9OPW19eH5Ljh6OWXX+bBBx/ssryv1+Ckk07i6aefbtvn8ccfZ82aNSQnJ1NZWcnf//53kpOTGT16NJMmTeLtt98mOzub5ORkZsyYwYgRI0hOTsblcrFu3TrOPffcLue45557OkzkGmqdr4Hb7eb999/vdm6nocjf78GwYcOYN2/egI9tdotTb10d/CkAqr2THBYDrh669Q34SZ5PKJ9mvfxpMzsPNHG5NR7nW6+H5BwDFUlP80Il2q9B1L5/rYltPUJrXBJHx+fx2cTLaGiwwssvmx2ZKcLh90BrfUc/97s7iGHU0bXVyeI9jzuI5xFCtDPQEuHtWa1WMjIy2irq+cbauN3utnFNYBQqqKqqwmq1ti0rKSmhtLSUuro6KioqyM7O9nsOi8VCVlYWTqfTlOIMxcXFAXeNE70zO3GydLeiXVeHDrTWLm9lpByMLn3L6aZbXzCe5PmE8mnWt3/+AgA/uOYcpqQnheQcAxVJT/NCJdqvQVS+/8aD8OQ3oX4n3Pg0la++Tna0XYNOovL3wA9vN3F3p8Vp9GMeKaXU5cDlM2bMCEZoQgxppaWlPZbr7it/448sFkuH4gX+Sn9brda2hKS3pCg/P7+tbHlf5lIaKF9RiXAZ3zQUmJ04ueljVwelVClQqrVe7h10W66Usvq67kWarfsOceBIM7PGJYdt0iREVNr9Eay6Huq2wqKfgPJfqUiYRyl1KkbX7nTgf1rr9wY5hFWdihnl0I+y41rr1cDq7OzsYI7BaqNQrPvUzam/eKHD8m8snMGtMmegiDDLli0zvax253mHVq5c2eu8R4WFhbhcrkFNnKxWqyRNQWZ24tSnrg7eMVFurbXTu41DKTUd6Gtf97Dxzze3AfDDy04wNxAhxDHv/gee+i4MGwU3PAnTzjE7ItGJUur7tCsUBBR6K6z+LojnyMRIzHK9rwsxxlf5PoMKlFKF3od4VqCmu67jZrrt/AwyxnQsYvLous94b8cBkyISov/MTprAqL7nmxuqoqKCFStWBJSgDHYSI0lT8JmaOPWjq0MaUNvpGG6lVETOxu7xaP7++jYAzpsp5YyFCAvNjfDKb2FSNlzzACSPMzsi0YlSah5wG0ZlPd/f//nAfUopR7BanrwJkpOOCVrnbQY8FUaou+qdc/xozjl+dIdlr368LyTnEiIatC9NHm7lskVoxZgdAN6uDu1ed+jqoJSy+tZrrR3e9bRbb8GYPDfibK87DEDmFIu5gQghYP82aD4C8cPghtVw3eOSNIWvO4AcrfUjWusD3i8HkA3caXJsfaa1Xq21zk9JSTE7FCGEED0wPXHSWhcAVqWUzVsFr3NXh1yMSno+Bd5JcvO92y/uNA9UxPjzi8bcTdeeNsXkSISIchuegvvOA8fPjNcpEyHW7J7MogfKXzlybxfviO26Pdg+3l3PQ29so7G51exQhBAiIoTFnUFPXR2869pPeOvC/6zwEWfDzoMAXDVvosmRCBGlWpthzc/hjT/DhFPhjK+aHZEITE/TVkgftABMSU+ictNefvrkR0xJS2Lh7LFmhySEEGHP9BanaHXoaAvrdx7k/JljiIuVH4MQg+7gTnjoCiNpyr4ZbnkBUqeZHZUIjDY7gGBSSl2ulCo7cGDwijU8eMN87LedCUCLZ0hdTiGECBm5YzfJb1/YBMClc8ebHIkQUarpENS54Or74Qu/h7hEsyMSIeadAzDsmDHGKSZGMSxeSuwLIURfhEVXvWjkq6aXmzXZ3ECEiCYeD2x8Ck64HEbPgG+/ZxSDEJEmWym1ENjvZ12GUuqCbvazAcGbOXOI+Ndbn/DK5r3cfM50po8e0fsOQggRpSRxMoH7cBNglCCPjVEmRyNElDhcB48VwMcvwJcfgeNtkjRFriyMMuTd/QEt6Ga59Elr5zjLcKxjRrD+8wO8snkvY5MT+eai480OSwghwpYkTiZ4efNeAC6cI6WOhRgUn1XDqhuhfidc+luYscjsiMTAODHmcOoLBdwXglgiVtqIBF783gJaPZqMO5+RrFIIIXohiZMJ3t5qFIS6OlOq6QkRcs5/wVO3Q/IEuOV5mJhldkRi4Bz+ypH3RilVEYpgBirUE+AGam/9URqbW2XskxBCdEOKQ5hgp/sIAEkJkrcKEXLJE2CGDQpelqRpiNBa39HP/e4OdizBYPYEuApIiIvhX299wk+f+MiUGIQQIhLInbsJtu47xDkzRpsdhhBD156NsOMdyLzeGMt0vM3siEQIKKVGAWmAVWv9otnxRKqYGMWjXz2L/H9WceBIs9nhCCFCzO12U1xcTEZGRtvrwsJCk6OKDNLiNMiONLWyrfYwo4ZLzipESLy/ClYshJd+A0cbzI5GhIBS6mOlVC1QAliBPnfbEx2dNDGF5GHxHGpqYXvtYT5zH0FrGfUkBk9RURFZWVkopcjKyqKgoKDtKy8vj+XLl4fs3GVlZWRlZZGamhq0Y7pcLnJyckhNTcXhcATtuH05f15eHm63u8Nyt9vNokWLWLZsGfn5+eTn52OxWEJ6fYcSuXsfZLevfBeASalJ5gYixFDTchSeuwOqHoQpZ0Hug5A40uyoRGhkADla6zVmBzKUJMTF8OrH+zjv7pcA+F3eKVyTNcnkqES0KCkpweVykZGRQUlJCTZbx54Cdrud1NRUqqursVqtQT13fn4+2dnZLFoUvMJBVquVioqKtladweZ2u3E4HNTV1WGxWNqWFxUVsWTJkg7LysvLycnJGfwgI5AkToPsxY17ALjdNtPkSIQYQlpb4B+XwY61cNa3YNFPIVb+vA1hjqGUNIVLcYjf5p3Ch58d4GiLhzsf+4DaQ0dNjUdEr7S0tC7LcnNzWblyJVlZWezf728Kt4Fpn0hEwnF7k5mZ6fc6VVVVUVDQccaGioqwrJsTlqSr3iBqafXQ1OphQsowhidI1SIhgiY2DuYuhv/7D1z4S0mahj5X+xdKqRSl1FKl1CpvN761SqnvK6VONSm+PjG7OITPrPHJXJM1iS+eehwAjvV7pLueCCs5OTm43W6cTqfZoUSszl33RN/I3cUg2lZ7GIBFJ4w1ORIhhoDWFnjpVzDlTJh5EZyeb3ZEYvC427/QWh8AVgArlFJVQL7Wep0ZgQ0F8bHGM9V3ttXh2neIjDHS5VWEB99Nv1mtOEJI4jSI3v3UDcDZGVJRT4gBqd8N9pvhk9fgrFYjcRLRpKdmkKrukial1AVSfa93CXExLM89mUL7+zS1eMwOR4g2paWl5Ofndxjj5HK5KCgooKqqioceeoikpCTcbjd1dXVUVFSwYsUKv4lWUVERGRkZpKWlUVdXR3Z2dp/j6VydDmDx4sW9JnYOhwOXy0VaWhpr164lJyeny5guh8PRlijW1dW1bVtSUtLj+h/96Ee43W7y8vKoqqqivLwcm82G3W5n7dq1uFwuiouLsVqtbdfJ4XBgs9koLy/vEIPT6cThcGC1Wqmrq+tQfc/pdLJ06VJcLhdr1qzB5XJRV1dHeXn5kO76J4nTIHrPmzidMtliahxCRLRtrxlJU+NBuPI+OPVasyMSg8/Sw7qeBj/kAZI4BWDUMOP24JZ/rCUhrvte/XGxMZRcM5esqV3HpIgB+PtlXZedeCWcthSaDsO/87qsjpt9NZx5CxyqhVXXd91//s1w0jVwYAc8WtB1/VnfgFmXwL6PYfV3uq4/7/uQsRB2vg/PLeu47qanA3pb/eF2u6mqqqK0tJSioiLy8zv2LvAVYUhNTeW9997jy1/+cofEKi8vr8ONvK+qXHl5eYftOo/7CSSurKwsKioqOhxn+fLlPZb2drmMnsa+95Gbm0tWVhbLli0jNze3bRun09nhOC6Xi9LS0oDWWyyWLoUpcnNzyc3NxW63s2zZMjIzMzvEvHbt2g5xOhwOSkpKOly75cuXU1BQQGlpKZmZmVRXV7dVDfS9n6KiItxu95BtFZTEaRC9vmUfSsGElGFmhyJEZNr5Pjx0OaRZ4brHYNyJZkckzJGtlFqIMXdrZ9Zu1lkAmdArQNnT0liSPZnGltZut2lq8fDsh7t4f8cBSZxEUK1cubItwXC5XFRUVFBQUNCWWPiTlpbG1q1bOyQx2dnZXRKipUuXYrPZulTmy8vLY9WqVQHHuHTpUnJzczscx+FwUFRU1GPiZLfbKS0tpaampm2ZLxnxvT+Hw9El8bBarW3JTm/rfQJNXvxt54upvcLCQpRSlJSUtO2TlpZGbW1t2+tQFO4IJ5I4DZKWVg+ufYc4eVIKSvn7rBdCdMvjgZgYGD8XLlkOp/wfJCabHZUwTxbgwH/iBEbLkj9S6SBAo0cmUpJ7co/buA838eyHu6hYv5sJKcO4+KQJgxRdFOipBSchye/6lvp645sR6T3vnzKp5/Wjj+95/YSTQ9rCBLBkyZIOSUBhYWFby1Hnm3kfi8XCqaee2mVZZ3a7nerq6i7L/VXy64ndbu/SJS07O7vb+Hxyc3O7xOXrLtj+OFlZWbjd7rZ5luBYq1hv6wfK6XTicrn8dl+0Wq1UVVV16Fo4f/78oJw3EkjiNEg27DT+oJ17vIxvEqJPPl8HT3wD8v5hfKCfttTsiIT5nHSfHHVHAfeFIJYBC5dy5H2VlBCHdfQI3nLV8kntYUmcREgtW7aMrKwsioqKup3HqbcJbH3V+Abajcx3nM5xWCyWLl0JO7NarW3buN1uXC5Xl25ymZmZlJeXs3Tp0rb3W1BQ0NaS1dP6el8CPQC+1j5/E/eWlJR0SaiGarc8fyRxGiRvb60FjO4PQogAaA3Vf4dni2DEWGhqMDsiET4cWuutfd1JKRWWNYy11quB1dnZ2RH1VCAhLoYXv7+AH5S/R3n1DjweTUyM9KgQoeFLUpxO54AnwO1r61Kw+VqrsrKysNlszJ8/v0uS4huT5CvQUFpaytq1a9sKOHS3/sEHHxxwfL5EqKeuke2ZfT0Hk8zjNEj+/fZ2YmMU5x8/xuxQhAh/TYfgsQJ46naYfh7c9iocN8/sqESY0FrfMZj7iZ4dbjLGQa3fedDkSEQ06Nw60xe+7n++FhUzjlNUVERpaanfCoE+y5cv73CuwsJCampqcDqduN3uXtcPlK9FaaDXaSiSxGmQfO4+Qqs8jRMiMG/cC++vgoU/gi+VQ1L0PM0SItLkZk0C4LkPd1HbcNTkaMRQ5WsFaX8z358b+9zcXL9d0Pp6rNzc3C7lu33H6WmC3uXLl1NUVNRhma/UN0BZWRlgtEp1ZrPZ2sZCdbc+GMUZLBZLWwW+zpxOZ1RPQCyJ0yA43NTC0RYPthPGmR2KEOGt8YDx79nfhpuehfN/YBSFEEKErdEjE1EK7n1pC/94Y5vZ4YgI5kseuktibDZbh5v29glQoC0tK1asYOXKlV22r6io6FNrzYoVK3A4HF2SCLvd3qGwhdvt7nJcf+fxJUS+bm/FxcV+t/G1UHW3fvr06T2eu66urkMhiu7iKSkpobS0tMvPwuFwdHh//o43lMkYp0Hw8FufAHDZyeNNjkSIMNVyFJ7/IWxxQMHLMCwFpp5pdlRCiADMnZTCuh/nMO+XFbzw0W6+d+Ess0MSEWj58uXU1NSQn59PRUVFh8lefdoXRMjIyMBms+FyuSgpKcHlcnHPPfewc+dOCgsL28p+g1FqvKCgAJvNhsViYc2aNRQXF7dVg6urq6OgoICysjJycnIoKSnpUtq7M4vFQnV1dVssmZmZHSaIbR9XcXExLpeLwsJCqqur2xIS3zny8/OpqamhqKiInJwcLBYL5eXllJWVtSVSvuP5zt3T+s7nBqNLX0lJCW63m6KiIgoKCsjPz6eoqAi73d52DXylxq1WK9XV1RQXF5Oent6WsHV+f77j2Wy2Lj+voUhpHR3VWbOzs3VVVVW/96+srGTBggX92ves4jV8fqCRj35+ESMSIzNXHcj7Hyqi/RqE7P27t8OqG+BzJ5z5DbD9DGLjg3+eIIj23wEY+DVQSlVrrbvWuBX9/pwKl9/LaXc8TUJcDJt/dcmgnztcrkEgNmzYwAknnBD049bX15OcHN3TNMg1kGsA/q9BX/7f9fQ5FZl38RFm/+FmgIhNmoQImc0vwKNLQXtg8b9gzhVmRySE6Kcbz5rGP9/cxjV/e6PD8qXnTpdS5UKIIUEGD4TYp3WHOdLc2jZ4VgjhpTW8/gdImQz5lZI0iaillLpcKVV24MABs0MZkAvnjOPsGaMZHh/b9vXBZwd48r3PaWn1mB2eEEIMmDSBhNiX7n8LgAWzpAy5EAA07AEVa8xuv/ifkDAC4oebHZUYIpRSo4BsoEprfdC77FSt9bumBtaDSJ3HqbOzZozmrBkdJ3lf+NtKnvlgF7Ex7/Hna2VKASFEZJPEKcTivRXBLpsr3RSE4JM3oPwmmJgF1/4HRozufR8hAqSUug/IB2qAIuDRY6vU97XWvzUtuCj1x/87lW//712qttXx66fXd1inlGJx9iRmjI3u8RhCiMgR9MRJKTUdWIQxs/u2YB8/krR6NK59h7CdMA6lZP4mEcW0hjf+BI6fQ+o0WHin2RGJIUYp9QOgRmsd4319jW+d1nodsE4pdbXW+tHujiGC7+RJFi6YPZb/vrOdf7+9vcO6w02tHGlq5fsXzSI5MU7mORRChL2gJ05a663A/UqpW4H7g338SPJI9Q4AZo4baXIkQpjoiBse/xpsehrmfBGuuBeGjTI7KjH0uLXWK9q9jo6SsRHgx1+Yw4+/MKfL8sxfVvCvtz7hX299Qm7WJH6bd4oJ0QkhROD6VBxCKXWqUmpagJtn9D2coeWtrbUAXHfmVJMjEcJEnlbY8xFcVAx5D0nSJEKlNoBtrCGPQgTsj/93Kj/+whwmpw3HXr2DxuZWs0MSQogeBZQ4KaWmKaVqgWqgRin1Vz/rr1ZKFSulViqlPkY+oHjmg51kTU1lQooMfBdRRmvYsBpam40iEF97G878GkiXVRE6nR/Wdfhl8z70k0F1YeTc48dwyznTmekd4+Tcvt/kiAZPtMyhKUQ4COb/t0BbnJzAGuA2YBlwoVLqFmgbjOsC7BiDcfOArRgDdKNaY7OHpIRYs8MQYnA1HTa65q38Cqz7l7Esfpi5MYlo4FBKPa+UWuitrKeh7cHeD4AK4DemRij8+uoCI+d9/sNdtLR6hnxSERsbS3Nzs9lhCBE1mpqaiIsLzuikXo+ilPo+sFRr/Ui7xcuVUquUUvuBNCALI3lCax3ZE1EEydEWo8vB2GS5YRRRZN8WWHU97FkP598BmTeYHZGIElrrdUqpu4EVwHSgfVEeO3Chrzy5CC9ZU1MBeOjNT3jozU+YPy2V8tvOMjmq0ElOTubgwYOMHi0NoEKEmtaa2tpaUlJSgnK8QNKv+d2UcC0FCrXWFwUlkiHmuQ93AXC8FIYQ0WLz82C/BWLj4SuPwIxFZkckoozW2gHMUEplYiRPboz5nOSBXhhTSnHfVzL5eHcDjg27cW53c+VfXm9bf9W8idxw1jTzAgyytLQ0tm83KgyOGjWK+Ph4qbwrRBBprWltbeXw4cO43W5aWloYO3ZsUI4dSOLUXafjKsARlCiGoBfW7wbgohPHmxyJEINk1HEwcR5c+TdImWR2NCKKaa2dGF3MUUqNUkpNi/bpMcLdxSdN4OKTIHNqKmWvuNpKIr67fT+/enr9kEqcEhMTmTJlCnV1dWzbto3W1uAUxWhsbGTYsOju5SLXQK4BGNcgKSmJ4cOHM2LECFJTU4mJ6VM9vG4Fkjj57WystT6glHIFJYohRmvNmg27SU2KZ/roEWaHI0TouD+F9Y/DWd+E8XPhhtVmRySilFLqLmAeRsJUqrXeppRaCWQCa5RSqUCRJFDh7ewZozl7xrEubDn3vMzHexr4W2UNANYxI4bEA8nExEQmTJjAhAkTgnbMyspK5s2bF7TjRSK5BnINILTXYKAjpbodwamUKtZaLwvkIEqpQowxUmkAWuuyXra3YBSpqPEuqvI+YQwLH31+kMZmD/83f4rZoQgROlsc8MhSo3LenCvBMtnsiER0W4uRMG2FtglxM7XWx/s28I7Z9df13FRKqcuBy2fMmGF2KGHnujOn8pMnPqLkuY0AJMTFsPlXl5gclRAiWgWSOGUrpRbiv8ve/G5anayAJZAAlFIlwFqttd33WimV63vtZ3sLUK61zvG+zsdIovICOd9geG+HG4CLT4r8p2JCdOFphZdL4OXlMHYOLP6nJE0iHKT6kiavAuC+TttsJQxprVcDq7Ozs5eaHUu4uf7MaSzONv6+/N6xmbJXXDz34U4S4mI4Z8YYEuKC0/1GCCECEUjilIUxlqm7kYtFfpZpoMdWo3bytdbtj7ESKMGoguTPCozCFD6rCLOxVp+7jwAw5ziZ6FMMQeU3woYn4ZQvwWW/g4QksyMSAto93FNKpWA8wOv82TC061wPUcPijWk9Ro9IRGu47WGjg8nV8yZyydwJzBgrRZiEEIMjkMTJSd9bc1KBXp+ceSsfdeYGbD3slgssVUpZAYu3i567j/GF1LufulEKRiQEp2a8EGHlpGtghg0yr5cJbUU4aZ8U5QNurfW7nbZJH7xwRLDdfM50zps5hr31R/nKA2/z6LrPeHTdZwB849RERn92gJMmBqfksBBC+BPInb2jU/eHQGwNsHBEGlDXaVnn123aJVrZeOeNUkqVY8wz5e5jjCHz7nY3s8YlExsjN5ViCNAa3vorx322HVgAJ15pckBC+HXAO4bpAEavhVzfCqXUNcAdhFGXbtF3sTGKWeOTmTU+mVcLF3LgSDO/fGo9b2+t4953j3Lvu6+xtfhSKe0thAiZXhMnrfUd/Tmw1vruADazdLdCKWXxkwxZ2x3flzitxOi+1+UD0Tv+KR9g3LhxVFZWBhCSfw0NDQHvf6iplZbGQwM6X7jpy/sfqqLxGsS2HGL2xj8zZt+bjEg9g8qXXorqVqZo/B3oLFyvgdZ6jfeBnQ3I0lqvg7YiEWB067YB95sUogiiyWlJTAYevvV0tuxp4PZ/vc7GOg9Fj7zPglljuXRu8KrVCSGEj9l9ydx4K+m10/l15+3BmEPKx0W7J4vteavzlQFkZ2frBQsW9CdGwChtGMj+dYea4LkKzpw9mQULTur3+cJNoO9/KIu6a7DrA1h1O+z/BC78NR8fPZEFCxeaHZWpou53wI8wvwa1WusV7RcE+BBPRKj42BhOmDCKJbMS+Pmbjayq2sGqqh18fWEG382ZJT0/hBBBFXDipJQ6FeNpXTrwP631e0E4fx1dW50sAN10vXP5Wef2xuevhWrQvbx5DwAzx8lgVRHB6nfDAxfCsBS48WmYeiaEYSuDED5KqRcwihnJOKYoND0llo2/vJg1G/bw9f84+ctLNTz7wS6sY3xzKSpuPmcaZ2WM7vE4QgjRk4ASJ2+/8eXtFhUqpQq11r8byMm11k6llLvT4jS6qZKntXYppdydkiQLxiDgzscxxYOvbQMgZ46UIhcRyOOBmBhIHmdUzJuRAyPHmB2VEIEoR8YwRbVh8bFcdvIETjxuAd8rf48jTa3sPNAIGPMrOjbsZmp6Eh6tiYuJ4ZdfPInMqRYUiuEJsSZHL4SIBL0mTkqpecBtGB9IvoRmPnCfUsoRhJanVZ3mbcqhXblxb/W8zHbri4HFHCt3vsS7LCx8vKcegPEpw0yORIg+qq0B+01w0W9g2jlw6pfMjkiIvqijl3LjfZmYXUSuaaNH8MhXz+qw7B+vb+XdT90AHDjSzEub9vKVB95uW5+bNYmSa06Wrn1CiB4F0uJ0B5DTqbKeQymVjZHgLBlIAFrrAqVUoVLKhlH8oabT5Le5GMmU3bv9cu/2hd71tVrr5YSB5lYPjc0ezj1eugKICLP+SXji6xATCy1HzY5GiP6oAfKVUunAWoxu3O2rtKZhdDeXxCkK3Xj29A6vn3j3M3Z5W6OKn92IvXoHG3Ye5BdfPImsqalmhCiEiACBJE7KXzlyrbVbKRWUWdh7Sny865b7WRZ2Nu0yWptOn95TfQshwkhrMzh+Bm/eCxOzIO8fYJlidlRC9MeL3n/r8N9lLw2QSX4EAF88dWLb9+fNHMMlf3yVjz4/yDV/e4NL547nrmtOZtSweBMjFEKEo5gAtul2XiVgX7ACGQrWbDAKQ5wpg09FpPig3EiaTsuHm56VpElEMpfWOk1rPaObrzSMqSsGjVIq3/tV6u12LsLQCRNGUfObS/n5FSeiFDzzwS5O/tkLOLfvNzs0IUSYCSRx6rHPuDimYsMuAE6YkGxyJEL04ojb+Pfk/4MbVsOld0NcoqkhCTFASwPYpiTkUXh5J2yv8k6LUU67sbsi/MTGKG44axobfnExeVmTALj6r28YU4wIIYRXIIlTvyilwqZgw2BRKJIT40hKMHt6LCG64fFAZQn8ORPcnxoV9KafZ3ZUQgxYuwlvRymlLlBKjfKt806ngb9u5yFkBQq831cB2YN4btFPw+JjuTvvFCZahgOQ+csKXt68l6ptdXz42QG0lmfJQkSzQO7ws5VSCwF/bdYZSqkLutkv6gbhfvDZAb5wssxWLsLUoVp4LB+2OODkJZAkY/HE0KKUug/IxygUUQQ8emyV+r7W+rf9OKbFe8x0rXWRn/WFGHMMpkHbxOt4ixz5Ch1l03HidhHmXilcyGV/epWNu+q54cF32pYXXz2Xa0+TLs1CRKtAEqcsjDLk3dXoLOhmeVQ9lmk42gLA+FFShlyEoU/XQvmNcGgPfOH3kHUTKCm7K4YOpdQPMKqyxnhfX+Nb522NWqeUulpr/Wh3x/BzTBvGXIEZ3awvAdb6KsEqpUo6Ta/hU4DMMRVRYmMUT33zHKo/2U9Tq4e6Q018+3/v8tvnN0niJEQUCyRxctL3P/gKuK/v4USu973zQ0xKHW5uIEL4s/Z+o9T4LS/AcfPMjkaIUHBrrdsXfxjwwzuttQNAKTUfI4HqLL9TK9RKjHFUbYmTUiofKAqXSdpF4OJiYzjdmt72+s5HP6D2UBP3v+riujOnkhgnk+YKEW0CSZwc/ekXrpSq6Ec8EevpD3YCMHeSxdxAhPBpPAiNbqNS3mW/A08zDJf5ScSQVRvANkGrbOct/tCZG6Obum8bG8ZnqEspZfMlYiIy/eqqk7h95Xv86ukN/OrpDdx2fgYXzB7LaTIFiRBRo9fiEFrrO/pzYK313f3ZL1J9vKcBgMwpFnMDEQJg90ewYiH870tGQYjEkZI0iaGuc3e6Dn1RlVLTgGDOFZFG1+k62l57E6tyoEIpVUP33dqDx9N6rGKmCLqr5k3iue+cy/FjRwJw38s1LC59k5ZWj8mRCSEGi5R/C6LjUoahZNyIMNu7/4Wnbodho+DyPxqV84QY+hxKqeeBu4BqvF31vAlTHkaBh6wgns/S3QqllEVr7QR6fVrh7cqXDzBu3DgqKyv7HEhDQwOVlZVkVX2Xo4lpfDj3R30+RqTzXYPB8MNMgBH8/M0jbD3gYcYPn2X8CIXWMHFkDAUnJ5IYN/j3AoN5DcKVXAO5BhDaa9Br4uQtK27xvqwG6voyuDZavLO1juyp8kRfmKi5EZ4tBOdDMPUcyH0QkseZHZUQg0JrvU4pdTfGJLfTgfYPsuzAhVrrg0E8pRtvJb12+txny1uFrwwgOztbL1iwoM+BVFZWsmDBAtgzl+Sd79GfY0S6tmswiOad3syPHv+wrUT5U+/vZPfhVgoch8mamorWmkvnTuDWcwdn7mMzrkG4kWsg1wBCew0CaXEqAAq11veHJIIhxJKUYHYIItrtfBfO+S4s/CHESoOyiB5KqWneMUQzlFLzMMYzuTEmoT0QglPW0bXVyQJgWiGIMbNhw5PQfATipVBRqKUMj+fP1x4rtvObq5u5/X/vcrTF6Lr32pZ9OLe7ufns6cTESG8UIYaCQO6sXJI09exoSysAlqR4kyMRUenjCph8GgxLgVsqIC7R7IiEMEM5MB+OlR8P5cm01k6llLvT4jSM6Tv6RCl1OXD5jBkzBhbUmFmgPVC7BcbPHdixRJ+NGhbPAzfOb3t9/YPv8Mrmvbj2HWKGd1yUECKyBTL4oV+T9vUwMe6Qc+BwMwBzJozqZUshgqi1BV74Mfw7F177g7FMkiYRvbKUUs8rpa4exHOuUkrltnudA5T29SBa69Va6/yUlJSBRTNmtvHv3k0DO44IiiXZkwGw3fMyazbsNjkaIUQwBNLi1N+5MPKAF/u5b0TZeaARAGmJF4OmfhfYb4ZPXof5t8KCfhW/FGIoKdJa362Umq6UWgqkAHat9bb+HtBbGc8G5HpfF2KUF3cCaK0LlFKF3rLjVowJeDtPfjt40jNgwZ0wdo5pIYhjLj5pPF85YwoPv7WdWx6q6tDqpLWmxaP5+oIZjBoez0UnjpPiUkJEgEASp/5OUDA4oyHDgK8U+dTRI0yORESFHVXw32uhqQGuvh9O7uv81EIMPb4pMLzzDq4AUEotUkrlYDwAXNXX4hDeBMkJLO9hm27XBSpoXfXiEmFBUe/biUERG6P41ZVzmT1+FG/WdJxmzDf3Y+Ej7wMwLT2JF7+3QMZCCRHmAkmccpRSrSGPJIK97TL+IM4cl2xyJCIqjDoORh8Pl90DY2ebHY0QYUtrvQZYo5S6BtiqlHJorZeYHVdnWuvVwOrs7OylAz7Y4Tqoc8Gk7IEHJoLiK2dM5StnTO2w7F6t2XmgkQNHmrnkj6+yrfYwz320i0mpx4p6JCXEkTFmhLRECRFGAioOAazs43FHA7f2PZzIVLl5LwDjkmV8iQiRw3XwThmc9wMjcbrpGbMjEiKseedvKsA7RxJGK1Sfxx9FnDfvhdf/CHfuhDip9BqulFIcZxnOcZbh/P3G+dz0j7V87d9Ov9s6vnseM8bKg1khwkEgiVOVrwtEXyilBjjKNTJ4PJq99UeZkDKMuFiZaFSEwI5qKL8BGnbD8TkwMZhzeAoxNCilVmI8sFsC3AbMw5i/abG35Sk6jJkNnhaj1UlapCPCuceP5t+3nk5j87HOPQcbm7nz0Q850tyK7Z5XSIiLoanFw/D4WB756llMHz2C4QmxJkYtRHQKZXEIdz/3iyj1jS0ALJg11uRIxJCjNay9H55bBskT4ObnJGkSont5GEUcnBgtS6tCNH9T0AVtjBMYJckB9m6UxClCxMXGcPaM0V2WXzVvEo9U72gbR33fyzUcaW7l0j+9CsAz3zqXEyYkS1c+IQZRyIpDaK2josxX7aGjAMwaJ3M0iCB7/ofw1l/g+AvhqlJI6m+dFiGiggvI887hFFGCOsYp/XhASUnyIeKarElt33/vwpm8+vFevvO/dznY2NKWQJVel0VCXAxb9rdyvtaSSAkRQoEkTplKqWStdX3Io4lArr2HABg1XCa/FUE254tGsnTOdyFGuoEK0YvS3pImpdQFWuuhPU1GQhKkTjVanMSQEh8bwwWzx/H+zy7ixY27ue1hJ00tHgr+Vd22zb3vV/CFkydwxSnHMW9KKglx8tkhRDAFkjiNBl5UShUDzoHMiTEUtXg8AEyXUuQiGN5fBfs+hgt+CFNON76EEL0KcCxuAdEwv+AXfg8jpPv4UHbB7HFs+uXFbNxVT2NzKw1HW7jugXc4cKSZf7+9nX+/vZ3UpHiqf5QjJc6FCKJAEqdpGN31LECGUipTa/1oKIOKJC9u3ANA8jBpcRID0NwIzy+Dqgdh6tnQ0iQVsYToA6XUKIxiEIu624T+j9kNqaCOcQLIuCA4xxFhTSnFCRNGtb3+x8UjOP/883lnax1Lyt5i/+FmrHc+wz9vPo3zZo4xMVIhho5e23C11ge01lu11uu01mskaeqotqEJMCavE6Jf9m+DBy8ykqazvw3XPylJkxB9dz9QDmQDM7r5Csvqelrr1Vrr/JSUIBWjPVQL75dDw57gHE9EDKUUp1vT+eBnF5LiHUJw/YPvkP0rBx5PWD43ECKiSOfXAdq4q54paUlSilz0T/MReOAiqNsK//cfyPkFxAbSECyE6KRCa73C+5Bvq58vF9EwjxMYD2MevRU+fcfsSIRJkofF895PL+S/S88AYF/DUax3PsMrm/fyad1hk6MTInLJ3f4AfeY+QmqSdNMTfeQdG0f8cLjsd1DwMsy+zNyYhIhsdb1toLV+ZDACMd2Ymca/UiAi6p2Zkc7Hv76k7fX1D77Ductf4ksr3uK9T93mBSZEhJLEaQCONBmT1VnHSCly0Qf1u+Ghy+G9lcbrE74AadPNjUmIyOdWSk3raQOl1PcHKRZzJSbDqElSklwARjW+rcWX8sTXzyb/PCsAb9TU8sW/vM7i+97kLy9tYV/DUZOjFCIySJ+gAdi46yAAs8cnmxyJiBjbXgP7zdB4ELJuNDsaIYYSDeQqpTKAavy3QC0BfjuoUZllzCzYu8HsKESYUEpxymQLp0y2UHTxbF75eC83/X0t72yr451tddz9/CaWnjudK06ZyInHjZJKfEJ0QxKnAfh4tzGb98xxkjiJXng88MYfYc0vIM0K1z0O4+aYHZUQQ4nd+28dMN/PegsQlk27Qa+qB3DcqfD6H+FoAyRKrwhxTGyMYuGssWz59SU0tXq44cF3WLttPyte3cqKV7cCcP/12SycPZZYSaCE6EASpwFoajXGqcwYKx9KohefvA6On8GcK+GL9xpdaYQQwVSltb6wpw2UUvcNVjB9obVeDazOzs5eGrSDnpZvfEnSJLoRFxtDXGwM5bedRcPRFl77eB+F9vc42NjCrf+sAuBrCzL4bs5MKYAlhFfQEielVLHWelmwjhcJNuw0uuqljpDS0aIbR/bD8FSYfi7c+LQxR5OSJ3hChEBBANuUhDyKcJE83uwIRAQZmRjHxSeN5+KTxrP+84Nc+dfXaWrx8NfKGv5aWcPw+FheuP08JqfJ1CsiugXzEUJhEI8VEfbUHyUpIZaRidJwJzrRGtY+AL+fC585jWXTzpGkSYgg8E5224HWemtv+wWyzZDysQMe/7rx90iIAM05bhSbf3UJ7/30Qr50+hQAjjS3cu7yl/j6v51s3HWQZm+PGyGiTTATp6i7I9y0q57UJGltEp00HYJH8+Hp78KU0yF1mtkRCTHU9KvlSCn1t2AHEtYObId3H4Z9H5sdiYhAKcPj+c1Vc9lafCk/uGgWAE9/sJOL//Aqx//wWU766fPc8OA7fPT5AbQk5yJKBLOpJOr+12yvO8xp09LMDkOEk72bYdV1RhnghT+Cc78HMdI3XIggy1ZKLaTvD+yyQxFM2MpYZPy7xXFsbich+kgpxdcXziAvaxLO7ft5s6aWivW7+fxAIy9v3svLm/cCcNncCfz52nlSkU8MadLHrJ+aWoxm6gmWYSZHIsLKh3Y4tA+uewwyFpodjRBDVRbgoO+JU1g+4AtJVT2A1KmQfryROJ35teAeW0SdsaOGcfFJE7j4pAn8/IsncaSplbe21vLM+zspr97B0x/sNL6+dQ4nHpdidrhChEQwE6cDQTxW2Nt1oBGAuRPlj0PUazkK+z8xnuieVwjZt0DyOLOjEmIocwJ5fdxHAdFTVc9nhg2q/w7NRyB+eNAPL6LX8IRYFs4ay8JZY7njktks/G0lBxtbuOxPr5E5xUJiXCy/W3wKx1nk904MHUFLnLTWUdVnzbl9PwDJw6TRLqq5t8OqG+Dg5/AtJySMkKRJiNCr6k+hB6WUMxTBhLXjbfDp21C/05hDTogQSB+ZyPs/u4i/vLSFN2tqeW3LPgDOuutFAG4+ezrxsYrsaWmMTU7klMkWE6MVov/C4q5fKVUIuIA0AK11WR/2LdVaB1KGNqjWeROn82eOHexTi3Cx+QV4dCloD3zxL0bSJIQIOa31bf3c745gxxL2ZtiMLyEGwdcXzuDrC2fgPtzES5v2cPvK9wB48HXjOUfpKy4AZo4byezxo8ielsr1Z04zK1wh+sz0xEkpVQKs1Vrbfa+VUrm+1wHsa8ojtLXbjMRpfIqMcYo6nlaoLIZX7oZxc2HxQ5CeYXZUQgjRveZGiJfPKzE4LEkJXDVvElfNmwRAq0ezaVc96z7dz70vbmHz7gY2727gyfc+5ydPfMSfrp3HFaccZ3LUQvTO9MQJyNdaF7V7vRKj1GyPiZNSKjOkUfUiPlZhSYo3MwRhGgW7PoB5X4FLfyvjBoQQ4e2jx+GxAvjaW5A23exoRBSKjVHMOW4Uc44bxZdPnwoYY8XPKF4DwLf+u45v/XcdV8+byGnT07js5AkkD5N7LBF+TK2T3E3y4wYC6VeQDVQENaA+2F53mPlSijyqpLjXg/tTo7z44n8a3fMkaRJChLuJWUYRm3f/Y3YkQrQZnzKMbXddxn9uPb1t2aPrPuOORz9g7s9e4NaH1raVOhciXJg9wUwaUNdpWefXXSilcoFVIYkoQC0ezf5DTWaGIAaL1vD6nzj13R+C42fGsrhEU0MSQoiAWSbDjEWw7mGjq7EQYeSsGaPZdtdlbC2+lMrvL2grHOHYsIcbHnyHs+96kZc27qHVE5azCYgoY3ZXPUt3K5RSFq21299ywK21divV8xQeSql8IB9g3LhxVFZW9jvQhoaGDvu3tLSQousHdMxI0vn9R4u45gZmbfoTY/a9za7U+WyxXE1rFF4HiN7fgfbkGsg1iFiZNxiTc29xwMyLzI5GiC6UUkwbPYInvn42Wmte+Xgft/2rms/cR7jpH2sBmDUumYdvPZ0xyfLwUpgj6ImTUmqa1npbgJu78VbSa6e3/m+LA626592uDCA7O1svWLAgwLC6qqysxLe/x6M58twzTJ8ymQUL5vT7mJGk/fuPGrU18PDVcGAHXFTMpsYTWLAweie1jcrfgU7kGsg1CIWQTYDb3qxLYMRYqH5IEicR9pRSnD9zDBt+eTEffX6Arz7sZHvdYTbtrmf+rx0AvPKDhUxJTzI5UhFtgtpVTyl1DVDk/f4HSqlapdTflFKjutmljq6tThaAblqbMjFmizfVvkNHgTCdgl4Ez8ixkDoNbnwazvwa9NLCKYQQ/aG1Xq21zk9JCeGE6rHxcOVfIecXoTuHECFw4nEpvFK4kJrfXMqPLjuB5ETjmf95d7/EtDue5sn3PqexWbqgisER7Ban/VrrryqlpgN3ATla6xeVUrcC93feWGvtVEq5Oy1Oo/vkKA2wteuiNx+weueBsmutXcF4E73ZsLMegKnypGPoaToMr/0ezrkdEpPh+ifMjkgIIYLj+ByzIxCi32JjFLeea+XWc63c88Im/vTiFsCoyOezYFIcnyRs4/ozp9LbcA4h+iPYiZMVeBEoALZqrV/0Lu9phvdVneZtygFKfSuVUlYgU2tt11o7aJdUeccwWbXWy4P5Jnqz+0AjACdNDOHTQTH49m2BVdfDnvVw3Kkw+zKzIxJCDIBSqlhrvczsOMLKjmp492G49HdGhVAhItB3L5zFdy+cxcZdB3ny3c8pe8VFi0dTuaOFyh0f8dMnP2LmuJGMHpnIvV/KJG1EgtkhiyEi2H811yilXsAoyJDfbnm3vdq01gUYrUY2byJU02ny21yMRKwD77Z53n0LvUUjBsWhphYARo+QwYlDxkePQdkCqN8JX7FL0iTE0FBodgBhp84FVQ/C1kqzIxFiwGaPH0XhxbPZ8ptLqfnNpfzu/OGMTU5k1LA4Nu9u4I2aWjJ/WcEFv6vk07rDZocrhoCgtjhprbcCF/peK6VSgGqMbnsv9rBfty1G3nVd1rcv/DDYfBUxU2QC3KHhzb/C88tg0nzI+wekTDI7IiFEcEhfnc5OuNwoEvHa7yHjArOjESJoYmMU6cNjeOeHxu91c6uHHz72AauqduDae4hzl7/ESRNHYb/tLIbFx5ocrYhUIS1HrrU+oJTK8SZUQ8aWPQ0AxMfKZ/KQMOtiaNgFC38EcdKcL8QQIjV8OosfZozhfH4ZbH0Vpp9rdkRChER8bAzLc0/h11fN5bZ/VbNm4x4+/Owgs3/8HGkjEjjDmkbGmJGcPMnCglljiI+VrquidyGfx2moJU0AIxONJxXD5YlF5NqyBjY8CV/4A6RZpdKUECJ6ZN8Er/8RKoth2jlSMVQMafGxMTxw43y01txTsZk/v7iFukNNPPPBrg7bfXvR8XzHdrwUlRA96lfipJT6m9b6qz2snw7Uaq0P9juyMNbqgeTEOPnPFYk8rfDycni5BMaeAEf2Q1JvU4cJISLUAbMDCEvxw2HRj6F+F2gPKHkIKIY+pRTfu3AW37twFlprDje18v6OA1y74i0A/rjmY/645mPmTbHwHdtMzp0xmpgYuc8THfW3XTK9/Qul1N/av/a2MuUzRLV6PPKfKRId2gcPXwMv3wWnXAu3rpGkSYghTGst/8G7M+8rcN73IUaSJhF9lFKMSIzjzIx0tt11Ge/+JIeTJxmVktdtd3PDg+9gvfMZHl/3mcmRinATrK566X6WDdknfXvqj0rPhkijNfzrKti7Ca74M8y7TrqnCCGim8cD6x+D5Akw9SyzoxHCNJakBJ78xjm0tHp4Z2sdX7r/bQC+s/JdvrPyXWaMHckflpwq09CIfrc4rfVOatuTIfvb9dxHu2hu8ZgdhgiE1sbNgVJw0W/g1grIvF6SJiFE2FBKXa6UKjtwYJCfN3paoOJn8GwhtLYM7rmFCENxsTGcNWM02+66jL9+ObNt+ZY9DXzhz69xxm/WcNezG2n1SN2ZaNWvxElrfTdwh1LqN0qpUXSqXORdNiMI8YWlkYlxjBs1zOwwRG8aDxgT2r76O+P19HNhwinmxiSEEJ1orVdrrfNTUgb5eWNcAlz4S9j1AbxT2vv2QkSRS+dOYNtdl7Htrsv4yRfmYB0zgl0HG7nv5Roy7nyGRb+r5GdPfsSrH++luVUepkeLgXTVuxB4ASgC3N5xTm7AAiwGsgYaXLhSwHkzx5gdhujJrg+MpGn/JzDlDLOjEUKI8DTni3D8RfDir43vZR47Ibq4+Zzp3HzOdJpbPfz1pRoqNuziw88OUrP3EP94YxsA86el8n/zp3D2jNGMT5GH60NVvxMnrbULmKGUKgGuAQowEqcKIFtrvS0YAYYjj4YY6eoVvpz/gme+D8NT4canYeqZZkckhBDhSSm49G74y+nwTCFc+x+zIxIibMXHxvBt2/F823Y8La0ePqk7jL16B3+rrGHttv2s3bYfgFHD4rh07gTOnzmGzKmpjBmZKEXFhogBF4fQWhdhtDpFjVaPRuZJC1O1NfDUd2Dq2XDNAzBSWgaFEKJHqVPh4mIYbjHGhcqDQSF6FRcbQ8aYkRRdPJuii2fjPtzE61tqecS5gzdravnf2k/539pP27afaBnONVmTyMuaxOS0JBMjFwMR8glwh6JWreXJQbg5st9oYUrPMFqZJs2XMrtCCBGo7JuOfS/JkxB9ZklK4LKTJ3DZyRMA2LH/MG+76vjgswO8vHkvW/cd4k9rPuZPaz5mWHwMv8s7lYtOHEecPImPKJI49VFzq4emFg+x8qESPtY/CU98A666D2ZfKmOahBAdKKVewBh368DoTl6ltX7X1KDC1Vt/gx1rjRZ7+ZwTot8mpSYxKSuJa7KMcYOHm1q4r7KGP724hcZmD1//jxOAuBhFjFJccepxTElL4oYzp5GSFG9m6KIHkjj1UZW3/+qIRLl0pmttBsfP4M17YWIWjJ9rdkRCiPBUCjgBK5AJLFdKZQFV3uWlQ3lcbp+0NsGHjxjzOs3vbdYRIUSgkhLi+O6Fs/juhbPYdaCRJ9/7jENHW3l03Q4UCnv1DgDuqdgMwEUnjqP46pNJG5FgZtiiE7n776M3XbUAnHe8jJ0x1cHPofxG+PRtOK0ALvyVUVpXCCG60lrrrcBWYA1wt1LKCuQDdUC5UmqptEIBZ34Ttr4Cz90Jk8+A8SeZHZEQQ874lGHkn5cBwO05MwE42tLK4+s+46n3d/Lqx/t4/qPdPP9RBbPHJ3PVvIksOmEcM8aONDNsgSROfaa1MWXV3ElDdn7fyFDzIuz+CHIfhJOuMTsaIUR4y1BKndo+MdJau5RSL2itX8RIpIqBd7s7QNSIiYEr74P7zgH7TbD0JUiUmzUhQi0xLpYl86ewZP4UtNbYq3fw4Ovb2LDzIMXPbqT42Y2kJsVzVsZoLpk7notOHE+8jI8adJI49ZH7cDOxUhjCHB4P7N0A406EU78MGYtg1ASzoxJChDmt9d1KqReUUhoox+iiB5ADvOj93mFKcOFo5Bi4ugwevhq2vQqzLjE7IiGiilKKvOzJ5GVPpqnFwyPOHTz0xjY27qrn6Q928vQHOwG4dO54lsyfwvkyt+igkcSpj3YfbKTVo80OI/ocqoVHl8L2N+EbVZAyUZImIUTAtNYXKqWuAZYAtwE1eKfSUEotAuZhdOMTANbz4VvrwDLF7EiEiGoJcTFce9oUrj1tCi2tHlz7DvGL1et5vWYfz3ywi2c+2AXAhXPGcdc1J5sc7dAniVMfjUiMI1WqnQyuT9dC+Q1waC9cUgKjjjM7IiFEBNJaPwI80s3qA4MZS0TwJU1bHNCwB079krnxCBHl4mJjmDkumYdvPR2ALXvq+fnq9bz68T5eWL+bF9ZXMCwWvt76Md+4YAZKKmMGnXSO7KN9DUdJlQong+et++DvF0NMHNzyAmTfLCVyhRB9opRapJT6frvXo5RSo3yvtdZrtNYrQnBei1KqRCmVGexjDxqt4Z0V8OQ3oeYls6MRQrQzY2wy/7rldLYWX8pvrpqLdcwIGlvhdxWbmb7sGe56diPb9h0yO8whRRKnPnLtPcShoy1mhxE9arfAjBwoeBmOm2d2NEKIyGQFRvteaK0PAvOVUheE+LzZ3nNHLqXg6hUwehb878uw/S2zIxJCdKKU4kunT+HF7y2gNCeJszLSAbjv5RoW/LaSaXc8zdl3vUj1J3V4ZLjJgEhXvT7aU9/ISROlol5I7V4PnhaYcDJcXAwq1qj0JIQQ/VMHrPS98CZMtRiJzYvd7TRQWmuHUiovVMcfNMNGwXWPwT8uhX/nwfVPwMTIbUQTYihLjFX8Z+kZaK15Yf1unv9oFy9t3MNn7iNc87c3AUgfkUDB+da2kugicJI49ZFHI+UfQ+nd/8JTtxtJ083PQ6yMJxNC9I2f0uOPeLvrrcUoDlGAkTiVBnAsC8Z8T+la6yI/6wsBF5DmPVdZMN5D2EkeB9c/CX+/BD56VBInIcKcUoqLTjTKlgO8s7WOZz/cyZs1tWzcVc9vntnIPRWbuXDOeH63+BS5tw2QJE59sP9QE60eTbqMcQq+5kZ4thCcD8G0c+GaB2QskxCiv170lh53YZQZrwDWYrQwpWitswM5iFLKBlgAv49llVIlwFqttd33WimV63s95KRMhKUvQpLRDQhPK8TEmhuTECIgp01P47TpaQC8sWUfRY++z6d1R3jyvc958r3Pufa0yVx56kROm54mRSV6IOllH7iPNANwymSLuYEMNQ174YEcI2k657tw3ePG000hhOifIq11OkZLUR1wB7ANKAOWKKWual8cojtaa4c3CXJ3s0l+pyRpJUZr1tA1YrTxUGv/J/DXM2Hrq2ZHJIToo7NmjObVwgv4+NeXkDPHuN/67zufsqTsLaYve4Y/rfmYllaPyVGGJ2lx6gNfUYjJqUkmRzLEDE81yt4u/CHMutjsaIQQEc5XIU9rvQ5YB9wNoJSah9HqdC3wgFJqrdb6ov6co5tKeW7A1p/jRZzYeCOBevgayPsHzL7U7IiEEH0UHxvDiuuz8Xg0r27Zx+8rNvPup27uqdjMPRWbyc2axIVzxnFmRjrJw2ToBEji1Cefu48A0oMsKFqb4dV7IPsmGDkW/u/fZkfUJx6Ph/3799PQ0EBjYyMez9B/MpOSksKGDRvMDsNUcg2Ma7Bp0yaGDRvGyJEjSU1NJSZyirfUeJOpFQBKqYFU+knDaM1qr8Nrb1e/bMCtlEJr7ex8EKVUPkbLGOPGjaOysrLPgTQ0NPRrv4GKm/lDTn7/FyT/78t8fPxSPp9oXvJk1jUIJ3IN5BrAwK7Bd04E94zh/PXdo2ze78FevQN79Q4AxgxXXDcngbmjY8O+K18ofw8kceqDbbVGLfxp6SNMjiTCHdwJ9pth+xtGa9Pp+WZH1CctLS18+umnxMXFkZaWRlJSEjExMWH/h2Sg6uvrSU5ONjsMU8k1gIMHDzJixAgOHz6M2+3m4MGDTJ48mbi48P04UUpdA5QAad6xTysxuvMNZNJbSw/ns2it3VprB5DV00G8xSTKALKzs/WCBQv6HEhlZSX92S8ozrsAHrmFmZtLmTl7DmTdaEoYpl6DMCHXQK4BBOcaXOlth/+07jCPr/uMB1/fyt7DzdxTfRSAhbPG8JcvZ5KUEJ5/90P5exAxjwnDwc4DjQBMTZeuev3mehlKz4Wd7xpzg0RY0gRQV1dHYmIikyZNIjk5mdjY8H/6IkSwKKWIjY0lOTmZSZMmkZiYSF1d54aXsGPVWs/QWqdhtABtwyggMW0Ax3TjraTXTufXQ1/iSPi//0DOL+Gka8yORggRRJPTkvjmouNZ95MLeeOOC7j57OkAvLRpL3N+8jxX/uV1KjftQevomRtKEqc+SPCWahyRGJ4Zdtj76DH415VGK9PSl+DkxWZH1C8HDhwgPT1dkiUR9ZRSpKenc+DAQBpugk8p9X2l1KntFtX4vtFab9VaL9dazwdyB3CaOrq2Olm853D35UBKqcuVUmXhdh0DFhMLZ38LEpOh+Qis/ArsfN/sqIQQQXScZTg/uXwONb+5lDsvnU3K8Hje/dTNjX9fy/Rlz3DDg+9w/6sus8MMOckA+mD9zoPEx8rNcr9NPx9OK4ALfmQ8pYxQLS0tJCRISXohABISEmhpaTE7jM5uAy5USmVjJE0upZRVa/3bTtv1+1Nea+1USrk7LU7DKH/e12OtBlZnZ2cv7W88YcO9HXZUGZVSL/0tZF5ndkRCiCCKjVHkn5dB/nkZbNx1kL+/to23ttby8ua9vLx5L796egPWMSPIzZrE5Scfx6TU4UPqQbMkTn1gSUogLnIGQYeHz6rhjXvhqlJISoNL7jI7oqAYSn8EhBiIMP2/UKC1XgNtlfRsGInUnRiJVJV3u5pu9g/Uqk7zNuUQwKS6Q9qYWVDwKjxyCzz5Dfjkdbj0bqM1SggxpMweP4qS3JMBOHCkmV8+tZ7KTXtw7T3E8uc2sfy5TcTFKB792lmcPMlibrBBIolTHxxtbmVy2nCzw4gMWsPa++G5ZZA8Hg7ugDSr2VEJIaKAL2nyft+5JPl0wAq4tNZbezqOt+S4DW+XPqVUIeDwVcfTWhcopQq91fOsGFX7hubkt30xcgxc9xi8XAKv3A1Nh2DJv8yOSggRQinD4/lt3ikANDa38vbWOu589AM+cx/hintfZ9yoRM6ZMYaii2cxdtQwk6PtP0mc+mDjrnrpqheIow2w+lvw4SNw/IXHWpuEEMJk3mSpx4Sp3bZOwAks72GbbtcFSil1OXD5jBkzBnqo8BETCwvvhIwLICndWNZ4AOKGQ5x0dRZiKBsWH8v5M8fw+h0X8NLGPfzpxY/ZvKueR5w7eMS5g2tPm8K3Fs1gQkrkNUZI4tQHY5IT2X+oyewwwt+j+bD5WVj0Ezj7dpDujUKIEFNK/U1r/dUe1k8HarXWBwcxrIAMqTFOnU0549j3T90OezbCF++Fif7mDxZCDDULZ49l4eyxaK2pWL+bHz3+If99Zzv/fWc7s8cnc2ZGOndeegLxsZFxrxgZUYYJj9ZMTI287HjQeFqNfy/4EVz3OJz7PUmaRBuns8vcn0IEU3r7F0qpv7V/7W1pirz5D4aSuXlwpA7uXwQv/BiaDpsdkRBikCiluPDE8bzzQxvLc0/m2tOmsHFXPX9/fRvH//BZ7NU7IqKsudzV9oHHo4kJz4HQ5mpuNJ4krv6W8XrcHLCeb25MIqy43W6ysrJwu909bldUVEROTg5KKbKysigoKOiyTV5eHqmpqaSmpvpd71NWVkZOTg55eXkUFBRQUFCAw+Foi6esrGxA72mg3G43RUVFlJWVUVZWRlFRUZ/2Lyoq6vDlT0ZGBmVlZbjdbtxuN3a7nYKCgl5/DkNEup9lEVrve4iYdQl87S2Y9xV440/w19PhM3mgIkS0WZw9meKr57K1+FJuOHMqAN8vf4+5P3uBHfvD+4GKdNXrg1atiY2RxKmD/dtg1Q3GhLZnfQs8HmllEl2sWrUKMJKZwsLCbrcrKSnB5XKRkZFBSUkJNputyzbl5eWUlZVhtVr9rne5XOTl5WGz2aioqOiwzuFwYLfbqaioICMjY4DvamDy8vIoLS3FajWKprhcLnJycrrE7E9GRgbl5eVkZhrdnex2u999XS5XW9IIYLVaqaiowGKxBPfNhIe1Sqlbtdb397BNyqBF0wdDcoxTd4Zb4Io/w8lLYM0vIGWSsVw+O4SIOkopfv7Fk/huzixuX/UuL27cwzklL3G7bSbfth1vdnh+yV+pPvB4kBan9jY9C6XnQd1WWPJvuPCX8sEn/HK73eTm5rJy5cqgHM9qtbYlHO25XC6ysrIoKSmhpKSky3qbzYbVajW9tclut3d5D77v7faei7IVFRWRmZnZljQB5Obm4nK52lrUfEpKSqiurqaiooLq6mpqamr8XrehQGt9N3CHUuo3SqlRQIc+H95lYZmZaK1Xa63zU1LCMq8LjWnnwC0vwMixRhXWf14Bjp/D0XqzIxNCDLKUpHgevHE+37rA+BP9e8dmFv2ukoajYTdHYHgkTt5yrrlKqXylVI990JVSFu/2hUqp8t62D6aNuw4iDU5eR/bDowVgmQoFL8MJXzA7IhGmXC4XVquVJUuW4HQ6cblCN7N4QUEBixcv9tsS5ZOZmUl+vrlDXVauXElWVlaX5Tk5OZSW9jwNkK8LYmc2m43y8vIuyzMzM7HZbB0SrSHsQmAxsB+wKaX+ppQq9o532goMjYnkhprmwzBqIrx2D/wpE6ofgtbwu2ESQoTWdy+cxZvLLmDOhFHU7D3EST99np0HjpgdVgemJ05KqRKM+TTsWusyIEMpldvDLsu01su9X3lA0WAlT6lJCRxsbB6MU4Wvw3XG08HhqXD943BLBaRNNzsqEcYcDge5ubnk5uZisVh6TQz6y26343A4AhorlJeXF5IYAuVwOPy2/FitVqqqqvzsYfCNVUpL61rePyMjo0uLU7TRWru01jOA32IkTwXerzQgW2u9zcTwRHcSRsDVpXDri5A6zRgv+7ezoHag8xMLISLNhJThPPXNc1g4awwAX/jTa2HV8hQOY5zytdbt73RWAiVAl/4qSikLxiSD7ZUCRUDI+97UHmri8vGjQn2a8LXtNbDfbFTLO71Aysl24+erP2L952FX8bhHc44bxU8vPzEkx25fiGDx4sXY7Xa/3egGyjdeKJCuaDabzbQCCT0lPxaLJaC4/I1Rslgs3bbmOZ1O6urqemyJG0q8nyl9q7Zhoqga49STSVlG972NT0HVg0YrFED9Lhg5DqSrvBBRISZG8febTuPB17byi6fWs+DuSt6+c1FY1BkwtcXJOyt7Z26Mmdq7Y1NKtb8zctM1mQq6Fo/RXT6cst5B4/Ewefsj8NDlkJhs9E0XIgAul6tDF7GCggJcLldISpNXVVX1afxObm5PDduhU1dX1+s23SVPvoTJ33rfsvbr1q5di91ux2KxYLPZ2qr4ifASlWOcuqMUnHA5XPcYxA+DliZYsQgevBhqXjJ6PAghosLN50znlMkW9jUc5dqyt8wOBzC/xSkN6HwX0e1dhdbaDaR2WpwDhLx/ymFvvjR7fHKoTxVejuyHx75KhutZmHOlUQ1pWBS3ugUgVC03kcjhcHQYT5SZmYnVamXlypVBH3PjdruDXi0uLy+vQyLS2tpKbGxsj/sUFBSENCnLzc2loqKiyznWrl3rN5b2rUwlJSVkZGR0W5Ew3HknsS0H/geUheNktiLIlIJzvwuv/g7+dSVMOs3o9TDzImmBEiIKPP61s7jw96/wzrY6vvnfdfz52nmmxmN24mTpboVSyuJNlLrl7bpnAxZ1sz4f74SH48aNo7Kysp9hwl73IUCxecsWKlu39/s4kcay/z1O/tjB+snXs2/M1fBW9M650dDQQGVlJSkpKdTXR1/lp9bW1j6/78bGxi77XHHFFZSVlfGjH/3I7z4NDQ0AHD58uNvzHT58mIaGhg7rLRYLtbW1Qf3ZPPjggx1eB5I4AT3G4Ht/hw4d6rLd4cOH2/bv7jz33HMP5513Xod93333XSZNMso6x8bGtq07/fTT/V7/H/zgB7zyyiu9vg9//P0eNDY2Dujva6C8k9hmK6WuAexKKY2RSK2SJGqIio2H+bcYcz+t+xe89kf47xK4/gmwLjA7OiFEiCml+M/SM5j/awer3/ucLXsaePbb55oWj9mJkxuj1am9rh3/u7cCyNNa+72b9xabKAPIzs7WCxYs6EeIhs9Wvwgc4fysuSw4eUK/jxMRtIZd78OEU4AFYFvCvuqNDOT6DQWVlZUsWLCADRs2kJwcZS2PGDfzfXnfTqeTjz76iO9///td1rndbt5++22/rR4jR44EICkpqdvzJSUlMXLkyA7rs7OzcblcAcfoq/bXF329Bv5MmTIFgBEjRnQ5VlJSEgCTJ0/udv/k5GTWrVvH3/72N6xWK3V1dWRnZ3PCCSdgtVp7je+EE07gD3/4Q7/fh79rMGzYMObNG7yngFrrR4BHlFIpGFX0JIka6uISYf6tkHmDMRXGdO8k62/8GVDEtgzNMvtCCBiTnMhbyxZxRvEaNuw8yNZ9h5g+eoQpsZidONXRtdXJAm3d8rqllCoESrXWg1JGqsXbrzocBqaFVNMheOp2+KAclr4Ix82D5PHARrMjExHG4XB0W0HP4XBQXl7uN3HyJTM9FUlwuVxkZ2d3WJaXl0dBQUFAXfbcbjdOp9OUOY0sFktbIYfO3RVdLldA3Q0tFkuXiYRXrlzZ4f34JhE2ayzXYNBaH8B4gLYikpMoKQ7RB7HxMOeKY68/eRM2Pc2ZsUnQfAOcthTSzZ3cWggRfONThnH/9dnc+s8qFv62klcLFzI5LWnQ4zA1cdJaO5VS7k6L0+hlzJK3XLnTlzQppWyhTqB2H/INSB3CA1P3boZV18PejbDwThh/itkRiSGqoKCAoqKibhMrq9Xa43xP/pKj/Px8SktLKS4u7rVqn8PhCGiMT6jGOGVnZ/stElFTU9PvsUcOh6PD+7ZYLH7HkdXU1AzJOZ26SaJeVErVYjxke9TUAHugtV4NrM7Ozl5qdiwR59r/wGfV1D7xM8atXQFv32dMxn7WN82OTAgRZLY547jxrGn8441tXPfA26z53oJBb9Awu8UJYJVSKldr7Ss/noNRYhwAbwW9TN96pZQNb3LlHeOUBmQS4gIRvjGoEy2Dn90Oig8fgSe+aVQxuu4xyFhodkQigvnmbupObm4uRUVF2O12v9uVlJRQVFTUpVUFei4CUV5eTlZWFjk5Od0mIG63m7q6uoBadjpPKBuMrnpgJGQVFRVdJuJ1OBwsW7asx30LCgoAOiSdvvmb2r/nJUuW+G1Rs9vtAc11Fcn8JFH5SqkqoAZYGc5JlOiHiVlsmPM9xmWWQfXfj1V+3f0RfPwCnPoVGDnG3BiFEEHxsytO5JPaQ7y0aS8/KH+Pe5acOqjnN30CXK11AWBVStm8xRxq2iVRALkYExj6ikFUYCRW+71fNcD8UMe5v9FoaUqIM/2ShcbBz2H8XLjtNUmaxICVlJT02A3ON99Sdy1Oubm52Gy2tiTBx+VyUVxc3CXhaH/c6upqioqK/CYHLpeLsrKybvcfLPn5+bhcrg6tak6nk7S0tC6JZF5eXpeJbdtP4Ot2uykpKemS5Nlsti6lx5cvX47VajX9/Q8mrfUBrfXdWuts4A6MSdarlFIrlVJXmx2fCKJRE4zeEsd5x9ttWQOOn8E9Jxi9KbY4wOMxNUQhxMA9cINx2//ous9oahnc/9Ph0OKE1np5L+uWe793A6YMMvLlS0kJvVfUihju7bD/E5h+Lpz5DTj9NqP/uBD95HA42uZqysrKYs2aNX5bdnzbuFyutrFJnVuISktLsdvt5OXlkZaWhsViIT09vddueL7kqaysjJycHCwWC1arlfT0dKxWq99WLDOsWbOG4uJiMjKM8Rg1NTVUVFR02c7hcJCTk9P2uqSkhOLi4g7blpeXd7nOvu54vgTS7XaTkZHh9xzRwluV727gbm9p81yl1Fpgqdb6XVODE8F39rfg+AuNanzv/gfWPwFjT4Svvi6lzIWIYDExiq+cMYWH39pOxfrdXDaIRdvCInGKBK3ehHZk4hC5ZB9XwKNLISEZvuU0EiZJmsQA2Ww2ampqet2utLS029am9nJzc/td3CA/Pz+sW1YsFkuvSSDA/v37+7UfGMnTUBzPFAztkyizYxEhNHY2XPRrWPQT2PgUHKo1kiat4bECmHImnHglDO88RaQQIpwVnJfBw29t55kPdg5q4jRE+50F367DRuYUFxvhT6k8rfDir+DfuTBqElz/uCRMQgghhra4RDjpGjjd+zDl0D74fB089R347Uz435fho8eh+YiZUQohAjQ5LYlh8TE8/cFO9hxsHLTzSuIUoARvwjQiIYJbnJqPGDOvv3K3MZngrRVStlUIETWUUqeaHYM/SqnLlVJlBw4cMDuU6DFyDHz9HcivhPlL4dN3oPwG2Pycsf5oPbQcNTVEIUTPfnjZHACue+CdQTunJE4B+rzBw8jEOGIieR6nuGGQZoUv/sX4ih9udkRCCDGYCnrfZPBprVdrrfNTUlLMDiW6KGUUkrj4N/C9jXD9E3D8Rca6t0vh7uPhsdtg03PQPHhPtIUQgbnujKlkTU1l0+56mlsHp0iEJE4Bio+BhqMtZofRd1rDG/fC3k3Gh8TlfzRam4QQIvqkmR2ACFMxsWBdAAneKUemng0nfAE2PQP/XQJ3Z8CjBcZnqhAibFx04jgA1m7rOjdiKERwv7PB9VmDhxMmjDI7jL454oYnvm4MiG3YbUwKKIQQEU4ptQqY3sfdLED3NfKFaG/qmcZXyx9g6yuw4Umj656vGt/zPzS6us+8xCiDLoQwxZWnTuQ3z2zkT2s+5qyM0SE/nyROARoep6g7FEH9nXe+Z8xbcWAHXFQMZ3zV7IiEECJYVgH5QHlvG7aTCgztmX9F8MUlwPE248un+QhsehbqaoDbYcIpRgJ10jUwZqZpoQoRjcaOGkbysDjectXR0uohLja0nekkcQqQ1jAnUlqcPnkT/vlFSEqHG5+BKaebHZEQQgSN1tqulFqktV7Rl/2UUtLiJAYufjh8sxr2bDCKSWx+Dl4ugRGjjcTpUC1sfdno+pckvUOFCLWzM0bz3Ee7eL2mlvNnjgnpuSRxClCrhthIKQwxMRPm3wLnfs/4Qy6EEENPf/4gu4MdhIhSSsG4OcbXud81ypv7pvb4+AV4/DZQMTAxC2bYIOMCOC4TYuW2S4hgu/Xc6Tz30S5e2rgn5ImTFIcIkEdrYsJ5pvF9W+C/XzLGNcUlwsXFkjQJIYYsrfVt/djnjlDEIgQjRsMwb1XEkxfDzS/AeT8wuqtU3gUP5BhjjcFoqaqtkUITQgTJKZMtALy0aU/IzyWPPgK0/6gO3xanjx6HJ75hPO2qqzGecAkhhIgISqnLgctnzJhhdigiGGJijS7yU06HhXfC4TpjnqiUicb6l34NG1Ybk9Bbz4dp58L0cyFlkrlxCxGh4mNjGD9qGJ/UHqa+sTmk55IWpwC1eKC2ocnsMDpqaYJn7zAm7Rs7G257VZImIYSIMDKP0xCXlAazLj72OucXcNk9MCnLKDLx+G3wvy8dW1/zEhz4bPDjFCKC3XzONAA27aoP6XmkxSlAsQpmjBtpdhgdVfwY3r4Pzvga2H5uVP8RQgghRPhKsxpf828Bjwf2fASNB411LU1GEtV8GFKnGfNJTTnTKDRhmWxm1EKEtTkTjAdP63ceZEoIzyOJU4A0kBgXJg10nlajK8A5txt/VOdcYXZEQgghhOirmBgYP7fd6zi4+XnY9ip88obRIvXuv2HBnbCgyEiw1j1sdAMcf/KxghRCRLlTJhuJ08e7G5hiCd15JHEK0JEWiDN7jJOn1Sh5umMtfNkOyeMlaRJCCCGGipgYmHCy8XXm140WqX2bINE7Hcrn6+D5Zcb3ccON7vmTT4OsG8yLWYgwMDLRSGlaPJ6QnkcSpwAcaWoFoOFoi3lBHNoHj9wCrko45UvQ2my0OgkhopLb7aa4uJiMjIy214WFhSZHJYQIqpgYGHvCsdfW8+H29fDp28e+3vgTnHiVsX7Tc/DRozAx2xhDNW6udOMXUUEpxXEpw1i/s56LQjh9miROAdh9sBGASalJ5gSw/W0ovxEO18IVf4Z51xlzSAgRZoqKinA4HDidTjIzM8nOzm5bV1dXx/z580N2c19WVkZpaSkul4v9+/cH5Zgul4uCggKqqqooLy/n9NMHdzJpl8tFUVERK1aswGKxtC13u90sWrSINWvWtC0vKytj+fLlkjwJMdSlTISUq+Gkq43XTYeNaUg2vQoNu4wHrO+vNNbFJhqtV9c9Dokjja5+CSONhEyIIeZIcytJ8bFA6CrrSeIUgOZWo9lvomX44J+8tRkeKzCeGN1aARNOGfwYhAhQSUkJLpeLjIwMSkpKsNlsHdbb7XZSU1Oprq7GarUG9dz5+flkZ2ezaNGioB3TarVSUVHR1qoz2NxuNw6Hg7q6ug6JU1FREUuWLOmwrLy8nJycnMEPUghhroR2D3WzboTMG+DADvisCnZUwf5tRtIEsPpbsOVFOO4UOG7esa/UaSYELkRwzZuSyrufumFm6Mb+SeIUAI93jrqEwSwO0XgQ4oYZCdO1/4XkCTDcMnjnF2KA0tK6tpXn5uaycuVKsrKygtYq1F77RCISjtubzMxMv9epqqqKgoKCDssqKioGKywhRDhTyqjAZ5l8rAufz4lXw/BUY6zUm38FT7NRZOK2V431ax8wJvKdcKpR+U9apkQEmWgZzosb99DsCV16I4lTADze2b0HrTbErg9g1fUw61K46Ncd+zcLEeFycnKw2+1t3flE37ndbrNDEEJEojlXHCsq1XIUdn9klD4H0NqYnPdwrfE6YSSMOwlOWQLZN3v3aZIxUyJspSYZLU21R3TIziGJUwB8iZMajHFF6x6Gp78Hwyww+7LQn0+IQea76TerFUcIIQTGuKiJ7R5eKQXf3Qh7N8LO94yvXR/A4TpjfeNBWG6F0ccbCdW4E2H8SXBcpjHJrxAms44xuqQ2tUriZCpv3kRMKBOn5iPwzPeNxGn6eXDNAzBybOjOJ4RJSktLyc/P7zDGqXMRBjASrLq6OioqKroUR/ApKioiIyODtLQ06urqOhSjCFTn6nQAixcv7jWxczgcuFwu0tLSWLt2LTk5OV3GdDkcjrZEsa6urm3bkpKSXte73W7y8vLaronNZsNut7N27VpcLhfFxcVYrda26+RwOLDZbG3Xz8fpdOJwOLBardTV1XWovud0Olm6dCkul4s1a9bgcrmoq6ujvLxcuv4JEY3iEo6VQ+e6jus8LXD2t41kattr8MEqY/mlv4XTloJ7u9H9b9wcGHsijJ0NCSMG/S2I6JWUYFSb9oQub5LEKRCD0lVv/yfw4WNw7vdh4Z1SajzS/d1Pa+GJVxofLk2H4d95Xdef+iWY92U4VGt01exs/s1w0jXGoN9HC7quP+sbMOsS2PcxrP5O1/XnfR8yFsLO9+G5ZR3X3fR0IO+q39xuN1VVVZSWllJUVER+fn6H9b4iDKmpqTidTnJzczskVnl5eR1u5H1V5crLyzts13ncTyBxZWVlUVFR0eE4vVWnc7lcAG3vIzc3l6ysLJYtW0Zubm7bNk6ns8NxXC4XpaWlAa23WCxdClPk5uaSm5uL3W5n2bJlHbo6Ll++nLVr13aI0+FwUFJS0uHaLV++nIKCAkpLS8nMzKS6uprU1FQcDkfb+ykqKsLtdkur4CBRSl0OXD5jxgyzQxGie0lpsOjHx14frjO6+qV5/3bW1oDzoWNd/1CQOtV4EDwpG+p3w6E9MHqm0dolRJDFem/UQzmTkyROAfCEssXp83fhuFONJzPfWgfJ44J/DiFMsHLlyrYEw+VyUVFRQUFBQVti4U9aWho1NTUdkpjs7OwuCdHSpUux2WxdKvPl5eWxatWqgGNcunRplyTN4XBQVFTUY+Jkt9spLS2lpqambZkvGfG9P4fD0SXxsFqtbclOb+t9Ak1e/G3ni6m9wsJClFKUlJS07ZOWlkZtbW3b61AU7hDd01qvBlZnZ2cvNTsWIQKWlAbTzz32OmMhLPsM3Ntg93rY4/3y9Z756DF4rghUrJFsjZlljOE+8+tGwQqtZaoVMSBtiZO0OJnLN4+Tryx5ULQ2g+Nn8Oa9sOTfcMIXJGkaSnpqwUlI6nn9iPSe16dM6nn96ON7Xj/h5JC3MAEsWbKkQxJQWFjY1nLU+Wbex2KxkJWV1WVZZ3a7nerq6i7L/VXy64ndbu/SJS07O7vb+Hxyc3O7xOXrLtj+OFlZWbjdbvLz89u29yWBva0fKKfTicvl8tt90Wq1UlVV1aFr4fz584NyXiFEFIuJMZKiNKtxX9PeCV+AEaONMVR7NsDeTbD5eTjndmP9mp8bydXoWTBmpvHv6Jkw+TRJqERAfInTgaMyxslUvpamMclBalo++Lkxoe2nb8Np+XC8zL0iosOyZcvIysqiqKio23mcekt+nE4nMPDiEr7jdI7DYrF06UrYmdVqbdvG7Xbjcrm6dJPLzMykvLycpUuXtr3fgoKCtpas3tYPlK+1z+FwdFlXUlLSJaGSbnlCiJBKmQRzO/U4aF+lb9xJxpxTezcbk/i2HjVaogq3GutfvhsOfGo8HBw9E9JngGUqxMqtrDDEecvnhzLNlt+2ANQ2HAUgKSEIl8tVCfZboKURch80xqwIESV8SYrT6RzwBLh9bV0KNl9rVVZWFjabjfnz53dJUnxjknwFGkpLS1m7dm1bAYfe1g+ELxHqqWtke2ZfTyFEFGpf2nxu7rHEytMK7k+McVG+1ib3Ntj07LFy6QATToGCV4zv31nBhM+3wrY4I6kaOU5aqqJM8jDjPj2EPfUkcQpEw9EW4Fi1jgFpPAAjxsDifxpN0UJEobVr1wZ8Q9+Zr/ufy+Ua0DxQ7Y/T1ySuqKgIp9PZoZufrwXLp32BiczMTDIzMyksLCQjIwO3201ZWVmP6wfaAuRrUerP+xNCCFPFxB7r8ufzxb8Y/x6ug9otRiGk2HaJ12t/YNbBHbDZu13CSDj1y3DpcuP1h49A8nHGMUeOlaRqCBqMH6lMCd0HKd6JtfrsUK3xlARgzheNGbolaRJRyJcM+LqRdf4+ULm5uX67oPX1WLm5uX5bd3wV77qzfPlyioqKOizzlfoGKCsrA4xWqc5sNlvbWKje1g+ExWJpq8DXmdPp7PH9CSFE2EpKM8Y9zfsynNyuQu133ufNM1bAVx41SqTP+4ox1xRAc6PR2+fvF8PvZsJvJsLfzoHqfxjrW1vA9bJRUt3TOuhvSQSH8nbS0yFscpLEKQBb9jQAkBDbj8v16VooPQ8eWXpsErnYfiZgQkQAX/LQXRJjs9k63LS3T4B8+/ZmxYoVrFy5ssv2FRUVAR/DdxyHw9ElibDb7R1as9xud5fj+juPL+HxdXsrLi72u42vBai39d2du66urkty5S+ekpISSktLu/wsHA5Hh/fn73hCCBFRYmI5OmwszFhkTP1xSQlk3WCsi02Ab1bDl+1wyd2QeR2MmnCsxerAdvjnFfCHufDr8fDnLHg4F7asMdY3HYJdH8LRenPemwiIr8VJuuqZbGSicZmGxfehq57W8E4ZPP9DGHUc3LhaZtYWQ97y5cupqakhPz+fioqKDpO9+rQviJCRkYHNZsPlclFSUtI2savL5aKwsLCt7DcYpcYLCgqw2WxYLBbWrFlDcXFxWzW4uro6CgoKKCsrIycnh5KSkl678lksFqqrq9tiyczM7DBBbOe4NmzYwI9//GOqq6vbEhLfOfLz86mpqaGoqIicnBwsFgvl5eWUlZW1JVK+4/nO3dP6zucGo0ufb3LcoqIiCgoKyM/Pp6ioCLvd3nYNfKXGrVYr1dXVFBcXk56e3paQdX5/vuPZbLYuPy8hhIh4MTGQnmF8+TNyPFz/JOzfCnVbj/3bYlRV5vN18A/v/IxJ6ZA6zShMcc7tRqXaI25j7FXK5I7jtsSgCsm0QZ0oHcr2rDCSnZ2tq6qq+rXvz1d/xP/e3saGX/mZ1NQfjwceuQU+ehRmXQpX/tWoDBPBKisrWbBggdlhmMp3DTZs2MAJJ5xgdjiDrr6+nuTkZLPDMJVcA//XoC//J5RS1VrrrjXSRb8/p+Tvs1wDkGsAIbwGh/bB1leMqn/uT4x/92+Dq++HyfPhvZXwWD6gjIfllqlgmQIX/Agsk6FhDzQ1wKhJIU+sovn3YPPuei78/St87ZRECq+19b5DN3r6nJIWpwB4PJqYQJNY3wRuY2aD7edw9rdlAKIQQgghRKQaMRpOurr79VPOgC/+1Rgj5f4E9n8Cn7wOyjvEY92/YM0vAAXJE4xkKmUyfOH3MGyUkYQ1Nxol2xNHDsY7GpJ8d9vSVc9kLYEmTu/9z3jSMP08WFDU+/ZCCCGEECKypU41vroz+3KjO6D7E3B/asxH9bkT4pOM9a//CaoeML4fnmokUJapsORh4+H75+uMOa9SJhrHkbmr/JIxTmHipY17ep5Mq7kRni0E50Nw4lVG4iSEECKqKaXyARdgBRxa676XkBRCRL4xM3uupnxaPkw9y0iofIlVU8OxTKDyLtj8nPG9ijVarSbOMxIrMCo3tzZDykQSju43hozERGP9t9BnTpI4BSB9ZCKfH2j0v7LOBatugF3vG4MEF/5ocIMTQggRdpRSViBDa13mfV0O5PW8lxAiKo2dbXx156LfwPxbjYTqwGdw8DNIGHFs/Uu/hl0fAHAWwNvxMOviY4nVG/cac2ONOs6Yy2rUhCHZciUtTmGixaOZN9ZPRb3aGihbaCS41640fkmFEEIIyAVq2r3u/2zNQojo1lNFQIDrnjCSqoOfs7m6kpnjkoxWKZ83/wL1n3fcZ86VsPgh4/tHlhpjrZIneJOr8TB6ptFlMIL4quoN+cRJKVWI0Z0hDcD3hC5Y2w9Uq8dDrL++emlWmH8zZN1olKYUQggxpCilLEA+kK617jJ4tYfPo3TvciGECK0R6cbXcafy+a4kZnauqvfd9XBkPxz83GitOvi5UZwCjMl/d39kLG90H9vnjK/Dxb+BpsNw3zlGQjVynJFUJU+A6efChFOMCYNbGju2gJmkrThECCuGm544KaVKgLVaa7vvtVIq1/d6oNsHQ2OzB4uveuTBnfD0d+HiYiNZsv0sVKcVQghhIqWUDbAAfh/1mvF5JIQQfaaUMZdoUhqMP6njutg4+NobxvdNh6Fhl3GvO2KMsayl0Zirqn4XfFZl/NvSaHQfnHCKMd/VvVmQOKpdYjUesm+BqWcac1ztet9YN3IcDEsJWbVp32GbWkNyeCAMEicgv9NTvJVACdDdB09ftx+w7XWHGT021qjhb7/ZmEF672ZpZYpiWmuUlJkXIqRP9symtXYAKKXmYyRQnfX0eVSLtxVKCCEiQkKS0ZsqzXpsWVIa5P3j2GutjZapGG8KMWwULPopNOyG+p1Qvxt2rDW6AgLsfBf++cVj+8cNMxKoL/7FaLXauwnWP3EssRo59tj3fRyDFR9rFMQ42DREW5yUUv76fLsBv7NW9XX7YEmKVyxuegz+uQrSZ8ANT/U8iE8MaXFxcTQ1NZGYmGh2KEKYrqmpibi4cHgGN7gC+DyyAwXt1jlDHZMQQoScUkbJdJ+RY+Hc73a//YRT4fonjcSqYbfRYtWw25gbC2Dn+0Zxi86WvggTs2DDaninzEikRoyFkWOMf0+43Ejamo9ATDzExjEi0fgsGhYXugfbZn/apQF1nZZ1fj2Q7YPiOp7h/xr/ByddA5f/SSYni3IpKSnU1tYyYcIEaXUSUU1rTW1tLSkpKWaHYoYeP4+01i6l1Fpvdz8rIJP7CSGiz3ALWM/vfv3JeTDnCmjY4/3aZfzra/XSHmg5arRiNeyB5sPG8owLjMTp9T9BZTEkpTFy+GiWxc1gGzeF7O2YnThZuluhlLJord0D2d47h0Y+wLhx46isrOxXkK8nLWJV0gjGpl8Mb1b16xiRrqGhod/Xb6hofw1GjBjBoUOHGDVqFImJicTExAz5JKq1tZX6+nqzwzCVXANoaWnB7XZz9OhRDh48iNvt5tChQ6xfv97s0AabpbsVvs+jQMY6BeNzSv4+yzUAuQYg1wCGwjUYAUyHt9/zvk6BjDvb1sa2HCG+2U1j9XpQm0hxjyR16hLimw8Qe3Q/CYlJJOqjIbsGZidObrr2Ae+pT3iftvdWNyoDyM7O1gs6VxkJ0IIFUFmZSn/3HwoqKyuj+v1Dx2vg8XjYv38/hw4dwu124/F4zA1uEDQ2NjJs2DCzwzCVXAPjGiQlJTF8+HDS0tLIyMggJionWuzz55dfwfickr/Pcg1ArgHINYBovAYLOry6idBeA7MTpzq6PrWzAPhpberP9kKERExMDOnp6aSnp5sdyqCprKxk3rx5ZodhKrkGcg3akc8jIYSIMqY+JtT6/9u7n+S2kSuO479XNQfAyMk6NZw5gSyfIORydpJ9gkg3kGpOoKJuIPkEtrzLUrxBJN1AnMo6iYw5wcuiGzIGJoU/JAiA+H6qWDb+iegW2U/deGj4o8KoXd6BpMU29gcAoA3bjEdm9quZ3fzxxx/bODUAQEv6kF/x2cyOc8szSdfZgplNCttf3R8AgB3ZSjxy93+6++lIJ9kAgMHovOPk7meSJmY2jTfJPhVuqD1WbkrXCvsDALAxMzs0s3OFODQ1s/P8NOTEIwAYl67vcZIkuftVybarFesAAGhNTMd7VCEGFfbZOB6Z2a+Sfv3ll182/VEAgBZ1fsUJAIAxI1UPAIaBjhMAAAAAlKDjBAAAAAAl6DgBANAhpiMHgGGg4wQAQIe4xwkAhsHcvetz2Akz+4+kf2/wI/4i6b9bOp0hGnv5Jepg7OWXqANp8zr4m7v/dVsns082iFN8LqkDiTqQqAOJOpBajFOj6Thtyszu3f2o6/PoytjLL1EHYy+/RB1I1EEf8TuhDiTqQKIOJOpAarcOSNUDAAAAgBJ0nAAAAACgBB2n6m66PoGOjb38EnUw9vJL1IFEHfQRvxPqQKIOJOpAog6kFuuAe5wAAAAAoARXnAAAAACgBB0nAAAAACjxQ9cn0Bdmdi5pKelAktz91fzIuvv3XZ3ymFki6TQuvpN0N/TyS5v9Ts3s2t3P2jq3XWjwHUgk/SbpKa66d/fHNs+xbQ3bgTQuJu5+1eoJtiz33X7j7hcV9t+rdrCvxh6fJGKURIySiFMScarzOOXuo39Jmks6Xre86f59fzUpf2H5SdJp1+XYZR2sOPau6zLs+DOQ5Mus0Ijddl2OHdfBeWH5sLhuSC9JU0nHkq4lXW+7vng1/r2MOj41rYPCMjFq4DGq4eeAOEWc2np72Hkl9OEl6Wth+fC1Rqbu/n1/1SlPbIhuC+vOJT11XY5dfgYK+w0+KDX4DtwWGqNE0qTrcuy4Dh5W1UvX5dhCPcwrBqRa9cWr8e+jVj3v4++FGEWMalIHxCniVBvt4ejvcTKzwxWrU4Ve7cb7913D8kzNbFLYf7Jm397b8Hd6JOluqye0Yw3LfyxpYWYTMzt099Tdl62c4A40rINnM7vN/YxTSZ+2fGq9tG/tYF+NPT5JxCiJGCURpyTiVF1ttYej7zgp5Dw+F9YVlzfZv+9qlSc2PD8WGp+ZpEUbJ7cjjX6nZnYs6XMrZ7Rbtcqfa4yOcutuY97xUDX5DJwp/IH2NeZQP7v7l1bOrn/2rR3sq7HHJ4kYJRGjJOKURJyqq5X2kI5TuHS70povWN39+y5Zt6FKeeI+U0mlN+j1WLJuw7o6iOtTd09bOaPdStZtWFP+l5Fbd196uNH2k6SPWz+z3UnWbVj3GYh/mF1KuldIG3jXxon1VLJuw0Dbwb5K1m0YSXySiFESMUoiTknEqbqSdRs2aQ/pOIXLdgeFdcXlTfbvu1SbleejpBMf9iw1qerXwXt3H/IIZl6q+t8BKTTEmaVCWsRQpar5GTCza0kLd58pjGif5lMi9lyq/WoH+yrVuOOTRIySiFEScUoiTtWVqoX2kOnIw2W7pLAukcIl/y3s33eNyxMv+17vQeNcqw5iCsDQy5xX9zOwXLEtlcIozhi+B/EzkGZ/jLn7wsx+kvR7q2fZH/vWDvbV2OOTRIySiFEScUoiTtXVSns4+o6Tuz+aWVpYfaA1jU7d/fuuaXli7vRjFpDMbDrU4NSgDg4Ucoaz5XeSJjFIfxnazacNvgNLM0sLwSfRgNNCGn4G/lf4GamZDfI7UNe+tYN9Nfb4JBGjJGKURJySiFN1tdUekqoXfI6NbGamMEe8JCnOyHJcdf8BqlV+M5sqfPjuzSyJsxetmr1kSCrXgbsv3P0qeynMWJTG5cEFpKjud+BS0vvc8oe4bshqfQbiduW2J4qjnPtoBO1gX409PknEKIkYJRGnJOLUq3bRHlqc13z04kjMo+INhZ57snDcNos5oqX7D1HV8scv3dcVP+KLu5/s4lzbUvczENefSjpRmLnnUtLNUEezGn4HXvjAn0Yu1auD+MfYmb49kX7Q7UBM65gqlEkKwWWRpXmMoR3sq7HHJ4kYJRGjJOKURJxSx3GKjhMAAAAAlCBVDwAAAABK0HECAAAAgBJ0nAAAAACgBB0nAAAAAChBxwkAAAAAStBxAgAAAIASP3R9AgCaMbO5wvMMDhWeUXAvKY2bE4VnFkwlLd3951eOyRzEfy+zZyKUvFf+uH/tw/MxAAAA1uE5TsCAxYfbPSk88G2xYvtU0q27/1jlmLjtTtJ1sSNUctyxpI+S3g78yfQAgIriA4d/l/SsMKj2HDcdKQzeLVasu3H3i3g8g3kYFK44AXvM3Rdm9tnMkhVPi39esf/SzC4k3ZrZohisXjnui5l9kPQg6cfvDwEA7KEDSUtJf8/HGDM7lzR395P8zmZ2KulttuzuF7lBuYt1g3lm9qfBvArHHZvZVzGYhy3jHidg/93q28hdFVkAmtZ8nztJiZkd1jwOADBMicIVobTKzu5+88rmlYN5ki4kzV+JLSsH8xRi2UOV8wKqouME7KGYPpFZKqRH1PWm5v7Ze6YN3gsAMDwHCqlyddTtzDCYh94gVQ+jFVMJ3sXFZ3c/6/J8tuy9pBvpZcSuTqrCUfz3U833PFPIXSctAgDGIWnQ5n93hagiBvPQOTpOGJ3cBAjf5V/vkZlix6mOeKXqWiFnvHQUMe5/pNBpmpekYQAAtqjrAcCYEtf2MQzmoTfoOGGM7iR9WfdHfpyJ7lqrZ5Y7lTRXuBG2bnpCm87MbKYwwnakbyNtZcc85ZbfKKT0nZSU7UPsfCruP1Ooq9oBFABQ30gGABnMQ+/QccKoxClMJ5Iu1+0TZ6Jb6ltedX7bjZnNe9ZpkkLHZSFJMZ/7Y4VjbldNYV7Bp0L5r8zs1sxme5buCAB99eoA4EAxmIfeo+OEsTlWyHf+aGbZulXpDUerGurYKen1ZX93fzSzuikNm7qU9BA7lb2uHwAYsioDgAPFYB56j44TxiZ7+N7axrSkczTViitRm4jPmkgKq6+yBwQ2tNVzrCCrr953LAFg4KoOAI4Zg3loBR0njE2q8hl2ppKWZnas0NFa5i75zxTucVLcfhC3v3RU4s26WUN95u6z197M3bf+wNgOUwnfSSI9AgDaUzoA2ERLg3hdYTAPraDjhLFZqHzihJmku6yzFHOus87AkaT7OIHEQtJv8edl9xfdKd7EGnOujzQC7p7Gkc+X50WZ2YSRPgDYulQVptiOMWiq0HFIyu75aWMQrwcYzMNW8QBcjM2FpPeFB8Qqd2OpFO5vys+md2Bm2UP0nhU7T+6euvtFtm+8ApW/2jPR7lLmkh2+x7qH6S4URvcydR9WCAAoV2UAUAoTJNzEjIhZIc7tNXdP43//NJjXzdlgn9BxwqjEKyBvJc3N7NzMTmNq3bP00rC+XCXJGtrYCB8pjPIdKjxgtuiDpNvc8qHCzEetiDcIX8fFuZldxythrx1zXjhmXvG9zhWfiaEQgFcddyLp0czmcdr2Xd9nBQBjUGUAUJKOYns8cfezHWcA1H1YbRNJ/JfBPOyMuXvX5wD0Rux4nGS547HD8LO7n5nZrcLsPV/M7EEhJW9hZklMVXuQ9I/silM+ba+r8gAA9k/sJF1IelIY0EsU7ntKc/tkzyScKMSiq+9+0PbPax7P5b2+pbEvtWbGvBhjZwqdmkdJiyr3VGWxObcqLR4XO5Yf4/s/xZ9N+jg2QscJyMmCUa7j9KDwsNs03jj7U/z/nULQOlBM24sdq8vc/U0Pe5ozDgDosWxAL/5/qtBxIR4BG2JyCCDH3Zdm9pSbUe8kdooSxQ5S3PVaIXVvmVt3ofAAvwOFEbT7nZ48AGD04sDdneIVmZgZ8bnbswL2A1ecgBbENII3A53GFQAwYPE+0ywtbSLpcz6ND0AzdJyALcgmZchyuGOK3wn51AAAAPuBVD1gOyaSkvgso5nCJBF0mgAAAPYEV5wAAAAAoATPcQIAAACAEnScAAAAAKAEHScAAAAAKEHHCQAAAABK0HECAAAAgBJ0nAAAAACgxP8BmXRA+QkiYTkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14,5))\n",
    "\n",
    "fpr, tpr, th = roc_curve(y_test[:samp_size], test_pred)\n",
    "auc_score = roc_auc_score(y_test[:samp_size], test_pred)\n",
    "rnd_class = np.linspace(0, 1, 100)\n",
    "\n",
    "ax[0].plot(fpr, tpr, label='AUC = {:.2f}'.format(auc_score) )\n",
    "ax[0].plot(rnd_class, rnd_class, '--', label='Rnd classifier')\n",
    "ax[1].plot(tpr, 1/fpr, label='AUC = {:.2f}\\n $1/\\epsilon_{{bkg}}$(0.3) = {:.0f}'.format(auc_score, 1/fpr[closest_point(tpr, tpr_p=0.3)]))\n",
    "ax[1].plot(rnd_class, 1/rnd_class, '--', label='Rnd classifier')\n",
    "ax[1].set_yscale('log')\n",
    "\n",
    "ax[0].set_xlabel('$\\epsilon_{bkg}$ - FPR', fontproperties=axislabelfont)\n",
    "ax[0].set_ylabel('$\\epsilon_{s}$ - TPR', fontproperties=axislabelfont)\n",
    "\n",
    "ax[1].set_xlabel('$\\epsilon_{s}$ - TPR', fontproperties=axislabelfont)\n",
    "ax[1].set_ylabel('1/$\\epsilon_{bkg}$ - Inverse FPR', fontproperties=axislabelfont)\n",
    "for i in range(len(ax)):\n",
    "    ax[i].legend(prop=axislabelfont)\n",
    "    ax[i].tick_params(labelsize=axisfontsize)\n",
    "    ax[i].grid('on')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFgCAYAAABEyiulAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABEZElEQVR4nO3deXycZb3//9eVvU3STtJM2tI9oXubqW3g6PfgAiYoIKiQluMGiJAI7h5txHPOT0WlJK5HFEwKBxCO2iYgbigkiqhHRdpA9zVTWtrSkjaZdMue6/fH3BOmaZY2meSeJO/n45FHey9zz6dpm/dcy33dxlqLiIiIjCwxbhcgIiIiF04BLiIiMgIpwEVEREYgBbiIiMgIpAAXEREZgeLcLiAaZGRk2NmzZ7tdhoiIyDk2btx4zFrr7b5fAQ7Mnj2bDRs2uF2GiIjIOYwx+3vary50ERGREUgBLiIiMgIpwEVEREYgBbiIiMgIpAAXEREZgRTgIiIiI5BuIxMRiXItLS3U19dz8uRJOjo63C5HIiA2NpbU1FTS09NJTEwc0DUU4CIiUaylpYUDBw6QlpbG7NmziY+PxxjjdlkyCNZa2traOHHiBAcOHGDmzJkDCnF1oYuIRLH6+nrS0tLIyMggISFB4T0KGGNISEggIyODtLQ06uvrB3QdBbiISBQ7efIkEyZMcLsMGSITJkzg5MmTA3qtAjzC2hsa3C5BREaRjo4O4uPj3S5Dhkh8fPyA5zUowCOo49RpavPy2X/zLZx4+mlsa6vbJYnIKKBu89FrMH+3CvCIsky6/XbaDh7k0Of/nT3vuJyj3/oWrft7XIdeRERkwBTgERSbkkLGx4vIrnqWGWvXMn7FcuofeZTad72b/R/9KCd+/3u1ykVEJCKi4jYyY8xqwA+kA1hry/s4txBYAVQ4u1YCJdZa//kcHw4mJoaUt15Gylsvo+3o6zQ++QQNFRUc+uzniJ00Cc/11+NZtZKEGTOGqyQRERllXG+BG2NKAL+1ttIJ7mxjTEE/L1sFVAElQFkP4dzf8WETPzmTjDvu4OKqKmaU/ZhxPh/HH3qI2vwrOfCx2zjxzLPYtja3yhMRGfGKioooKioa1DVqamrIz88nPz8/QlUNPdcDHCi01laGba8D+vybsNamWWuNtXaFtbbmQo+7wcTGkvL2tzPj/h9x8R//QMYnP0lLbS2HPvMZ9lxxBa9/7/u0HjzkdpkiIsOitLQ0YtdauXIlK1euHNQ1li9fTnFxcYQqGh6udqEbY5b3sDsA5A1zKcMqfsoUvJ/8BBl3fJxTf/4zgXXrOb52LcfLy0m+7DI8q1aSevnlmLioGOEQEYm4F198MWLXyssb1ZHRK7cTIh3ovgRNv0vSOOPc9fQyZt7f8WhhYmNJvfxyUi+/nLbDhwlUPkGgspJDn/o0cV4vEwtuIK2ggPhp09wuVUSizJF77qFlx05Xa0hcuIApX/7yBb9upLV0o5XbAe7p7YAxxmOtDfRwaAMQCJu0VmGMqQ/rhu/veOj6hUAhwMyZMwf9Bwmp3l/NpVMvZULCha2cFH/RRXg//Sky7ryDU88/T8O6dRz/cRnHf1xG8tveStqNN5LytrepVS4iI1p1dTV+v5+amhpKS0vxeDwUFhZSXV1NcXExubm55Ofns27dOm688UYKCgqoqamhvr6eQCDAiy++SH5+fleru6ampusDQVVVVde2x+PpGhevqqoiOzubwsLCC6o1EAhQXl5OVlYWAH6/n9WrV3cdLy0tZfny5QQCAaqqqiguLiYrK6vX/RFnrXXti2BXeUO3fVmABTzneY3VwMaBHrfWsmLFChsJr516zS59ZKm99PFL7b0v3GsPnjw4qOu1Hjxoj37/+3b3ZW+12+cvsLvf9nb7+g/us62HD0ekXhGJftu3b3e7hIirqKiwBQUF5+wvKyuzWVlZtqGhwW7cuNFu3LjRWmvt8uXLbUVFRdd5Ho/nrNdVVVXZvLy8s66flZVla2trrbXWNjQ02GDc9a37dZYvX24bGhrOum5hYWFXreE1VVRU2I0bN/a6vy/9/R0DG2wP2eV2c66ec1vhHgDbc+sbY0yetbY6bJcfWH6+x4fSlOQprHvPOh7d/ig/3/lzfrbzZ+TPyufmxTezJGPJBV8vfto0Mj/zGbx33snJP/2JwLr1HLv/fo498AApb3sbnhtXBVvlsbFD8KcRERle6enpAHg8HpYvf+PHdkVFxTkt2EAggMfj6fE6of2h14S2+3pNd5WVlWe9FqCgoICVK1dSUlJCVlYWRUVF1NfXk5eXR0FB8Oap+vr6HvcPBVdnodvgDPFAt93pQPW5Z4MxJguoMsZ4uh3yn8/x4bBw0kLufeu9/O6G3/GRRR/hr4f+ygd++wFu+f0tPHfgOTpt5wVf08THMyE/n5kPriW76lkm3XYbTVu3cvCOO9mbl0/dj35E29GjQ/CnEREZXuHBHZKenk5paSnl5eVUVwfjob8neA22y9rv93d9oAjn8Xjw+/3k5eVRVlZGVVUVK1asYMWKFQQCgV73D4VouI1sfbf7vvOBstCGMSYrdNwGx7WLu7XObyR4v3e/x4fTlOQp/Hvuv1NVUMUXc7/I4VOH+fRzn+a9T72X9bvW09zePKDrJkyfTubnP8fc5/7ItP/+bxLnzOHYfT9k7+VX8Oqdn+DU889jB7gwvojIcPP7/fj9fbexVqxYQV5eHoWFhWfNOB+qYITgB4CePiQEAgGysrKorq4mLy+PiooKGhoayM3N7fqA0dP+oeB6gFtri4AsY0yeM7Gs1p494ayAs+8LrzTGrHa+SoAqe/Ys8/6OD6uUhBRuWnwTT1//NKVvK2V8/Hi+/o+vc2Xlldz/8v0cbzo+oOua+HgmvOtKZv7PQ2Q/+wyTPnYrTZs28WrRx6nNv5JjDzxA29HXI/ynEREZvKysrK7w9fv9fbaWQxPYQi3z8NAOtcZ7MtBnbIcUFBQQCATO+nBRWVlJQUEBHo+Hmpqas94/NGGut/1DwQTHx8e23Nxcu2HDhmF5L2stG45u4CfbfsKfDv6JhJgErrv4Om5adBNzJs4Z3LVbWzn5hz/QsH49Z/7+D4iNJfWKy/GsupHkf/1/mBjXP6+JyAXasWMHCxcudLuMiCsqKmLFihUAXbPQS0pK8Pv9FBUVUVhY2DX+XFRURHZ29lkhXlVV1RWOa9as6Xp9bm4uxcXFbNiwgZKSElatWsWaNWsoLS2lsLCw1xnhNTU1Z12nsLCQQCDAmjVruOSSS7pmwYdmoYda1aFu9tAM9d7296W/v2NjzEZrbe45+xXgwxvg4fyNfh7b/hi/2vsrWjtbecf0d3DT4pvInZw76McHtr7yCg0VFTQ++Qs6GhqInzYNz8qVeG64njivN0J/AhEZaqM1wOUNCvBBcCvAQ443HWfdrnX8fOfPaWhpYNGkRdyy+BbyZ+UTFzO4GwU6W1s5WVVFYN16zvzznxAXR+oVV+C5cRXJb3mLWuUiUU4BPvopwAfB7QAPaW5v5le1v+Kx7Y/xyolXmJo8lQ8v/DDXz72elISUQV+/xb+PwPr1NP7iF3Q0NpIwaxYzH32E+ClTIlC9iAwFBfjoN9AAV/MriiTFJbFq/ip++b5f8oPLf8BFKRfxrQ3fIr8yn+9s+A5HTh8Z1PUTs+Yw+UvFXPzn55m6Zg2t+/dz8plnIlS9iIgMJwV4FIoxMVw+83Ieefcj/Oyan3HZtMv4yfafcNUTV/Glv3yJHcd3DO76iYl43v8+4qZMoWnT5ghVLSIiw0kBHuWWZCzhW2//Fk9f/zT/tuDfeO7Ac6z6zSpue+Y2/nLwLwxmCGScz0fTZgW4iMhIpAAfIaalTKP40mKqVlbxuRWfY9+Jfdz5hzt5/y/fz5N7nqSlo+WCrzkuJ4e2gwdpPz6we9FFRMQ9CvARZkLCBG5dciu/v/733HPZPcTFxPGVv32Fd1W+i7JNZQSaA+d9rXG+HAB1o4uIjEAK8BEqPjaea7OvpeLaCtZeuZYFkxbww5d/SH5lPt/4xzc4cOJAv9dIWrwYYmNp2rxpGCoWEZFIcvtpZDJIxhjePPXNvHnqm9nTsIefbP8JT+55kvW71nPFzCu4efHNLPMu63FhmJhx40icN49mjYOLiIw4aoGPInPT5vL1f/06z9zwDLctvY0Xj7zITb+7iQ8//WGefeVZOjrPfcjJuJwcmjZvwXZe+FPSRETEPQrwUcg73sunl3+aqoIq7rr0Luqb6/n35/+da35xDf+7438503am69xxOTl0njpF6759LlYsIiIXSgE+io2PH88HF36Q37z/N3zvHd8jY1wG9/7zXvIr8/nvmv+m7kwd45b5AE1kExEZaRTgY0BsTCx5s/J4/OrHeeyqx7h0yqU8tOUh3v3Eu3nF00ZMaipNmzSRTUSiW01NDfn5+eTn50fkvJFOAT7GLMtcxvcu/x5Pve8p2jrb+NPB5xm3dIkWdBGRqLd8+XKKi4sjdt5IpwAfo7ImZpE1MYuX614mKSeHlt276WxqcrssERE5T7qNbAzzZfqo3l9N4tIC6Oigeds2xuee88AbEYlCJf8sYWf9TldrWJC+gOJLR39LN1qpBT6GLfMu40TrCV6/OA3QRDYRGXqVlZVkZ2eTnZ1NTU1N1760tDSKioqoqamhurqayspKiouLqa6ujsj7BgIBSktLqayspLKyktLS0rOOl5aWdr1vUVERfr+/z/3RQC3wMcznDc5A39q6n6XTpmkcXGQEGakt34KCAjweD8XFxSxfvrxrX319PYWFhaxYsYK77rqLgoICCgoKSEtLo6GhYdDv+853vpM//OEPeDwegK5ALisro7y8nKysLPLy8rrODwQCve6PFgrwMWz2xNmkJqSyqW4Tl/p8nHnpJbdLEpExIC8vD7/fT01NDcuXL6e6uppVq1YBUFFRQVZW1lnnBwKBruAdiMrKSoCzrlFQUMDKlSspKSkhKyuLoqIi6uvrycvLo6CgAID6+voe90cLdaGPYTEmhhxvDpvqNjHOl0P7a6/RdvR1t8sSkTHgrrvuoqysDAC/398Vrunp6ZSWllJeXt7VfV5fXz+o9/L7/aSnp5+z3+Px4Pf7ycvLo6ysjKqqKlasWMGKFSsIBAK97o8WCvAxzuf1URuopX3xxQB6sImIDIvCwkLWr1+P3+8/q8W9YsUK8vLyKCwsjFjXdVZWVo8fAgKBAFlZWVRXV5OXl0dFRQUNDQ3k5uZ2fYDoaX+0UICPcT6vD4tld0Y7xMfrwSYiMiw8Hg+5ubkUFRV1BXVNTQ319fVdY+PhoT2YyWwFBQUEAoGzJqBVVlZ2jceHJs6FFBUVddXT0/5ooTHwMS4nIweDYUtgO++eP18z0UVk2BQXF3fNRIfgAiyrVq2itLS0K8TXrl1LSUlJ1wz1srIyNmzYQHl5OYWFhT1et6fzNm7cyJo1a7jkkkuor68nEAhQUVEBvNGVHhor9/v9rF69mvLy8h73RwtjrXW7Btfl5ubaDRs2uF2Ga67/1fVkjsvkq3+fRuNTTzHvxX9iYmPdLktEgB07drBw4UK3y5Ah1N/fsTFmo7X2nEU61IUu+Lw+NtdtJjFnCZ1nztCyt9btkkREpB8KcMHn9XGy7SRH504CNJFNRGQkUIDLGwu6xL9O7MSJejKZiMgIoAAXZk+YzcTEiWyu20xSTg7NmsgmIhL1FOCCMQaf1xdc0CUnh5a9e+k4ddrtskTEocnGo9dg/m4V4AIEu9H9jX7allwM1tK8davbJYkIkJCQQJMe9TtqNTU1kZiYOKDXRn2AG2NWG2MKjDGFxpieb/p749xCY0yZMSbP+SozxmT19RoJCo2D77kouK0Hm4hEh4yMDA4ePEh9fT1tbW1qjY8C1lra2tqor6/n4MGDTJo0aUDXieqFXIwxJcCL1trK0LYxpiC03YtVQCFQA9xurY2eZ79FsaUZS4kxMWxtquVds2ZqJrpIlJg4cSKJiYnU1dVx/Phx2tvb3S5JIiAuLo6kpCRmzpxJUlLSwK4R4ZoirdBaG/7MvHVACdBrgFtr04a8qlFofPx45nrmsun1TbzP5+P03/+OtRZjjNuliYx5SUlJzJgxw+0yJMpEbRe6MWZ5D7sDQF4P+yUCfF4fm49tJmHpUjrqjtH+2mtulyQiIr2I2gAH0oHuj4/p95lyzjj4eY2Zy9mWZS7jdNtpjswPPnZP4+AiItErmgPc09sBY0xvxzYA1dbaSmttOZBvjOnxCexOwG8wxmyoq6sbdLGjQWgi286Uk5iEBD3YREQkikVzgAcItsLDnftE9jDW2ppuk9ZeBO7q5dxya22utTbX6/UOqtDRYkbqDNIS09hUv5WkhQvVAhcRiWLRHOD1nNsK9wBYawM9vcAY03183A/0NJYuPQhf0CXJl0Pztm3Ytja3yxIRkR5EbYBba2sItsLDpQM9PtXdud+7qofudd1GdgF8mT5eOfEKrUsuxjY307Jnj9sliYhID6I2wB3ru41h5wNloQ1jTFbouNN1XtytdX4jwdvO5Dx1LegyM3iHoR5sIiISnaI6wK21RUCWs6paIVDbbRGXAqAobLvSWblttbMITJUzmU3O0+JJi4k1sWzrPEhseromsomIRKloX8gFa21pP8dKw7b94dty4cbHj2de2jw2123m/Tk5msgmIhKloroFLu7weX1sObaFBN8SWv1+Ok6ccLskERHpRgEu5/Bl+jjTfobD8zMAaNqyxeWKRESkOwW4nCM0kW3XpGYAmtWNLiISdRTgco7pKdNJT0pn84ldJGRnayKbiEgUUoDLOYwxLPMuY1PdJsbl5NC0aZOeQSwiEmUU4NIjX6aPAycP0JKTTUdDA20HD7pdkoiIhFGAS4+6FnSZlQCgbnQRkSijAJceLZ60mDgTx/aEY5ikJJo2a0U2EZFoogCXHiXFJTE/fT6bj28hafFimtUCFxGJKgpw6ZXP62Prsa0k5CyheccOOltb3S5JREQcCnDp1bLMZTS1N3FokRfb2krLzp1ulyQiIg4FuPSqa0GXzHZAE9lERKKJAlx6NTV5Kt5xXra27CPO69WDTUREoogCXHpljMHn9bGpbhNJvhzNRBcRiSIKcOmTz+vj4KmDNOXMpW3/AdobGtwuSUREUIBLP3yZwXHwvXMSAWjWk8lERKKCAlz6tGjSIuJi4tiR0gjGaCKbiEiUUIBLnxJjE1mUvojNge0kzp1L0yaNg4uIRAMFuPQrx5vD9uPbifctoWnLFj2ZTEQkCijApV++TB/NHc0cXJxJZ2Mjra+84nZJIiJjngJc+rXMuwyA3VODLe9m3Q8uIuI6Bbj0a0ryFDLHZ7LVHiJm/HhNZBMRiQIKcDkvy7zL2HxsM0lLlmhFNhGRKKAAl/Pi8/o4fPowZ5ZdTPPOnXQ2N7tdkojImKYAl/PStaDL3GRob6d5+w6XKxIRGdsU4HJeFqYvJD4mnh0TTwNoXXQREZcpwOW8JMQmsGjSIrae3kPc1KmaiS4i4jIFuJw3n9fHtmPbiPMt0Ux0ERGXKcDlvPm8Plo7Wzm4dDJthw7Rfvy42yWJiIxZCnA5b8sylwGwe7oBUCtcRMRFCnA5b5njM5maPJXt8XUQG6sHm4iIuCjO7QIAjDGrAT+QDmCtLb+A15ZZa4vCtguBFUCFs2slUGKt9Ueu4rHL5/Xxct3LJM6fp5noIiIucr0FbowpAfzW2konuLONMQUX8NqsHg6tAqqAEqBM4R05Pq+PI6ePcOpNF9O8ZSu2s9PtkkRExiTXAxwotNZWhm2vA4p6OznEGLO8t2PW2jRrrbHWrrDW1kSiSAnyeYMLutTOS6Xz1Cla/fpsJCLiBlcDvJcQDgB55/HyXIKtbBlGC9IXkBibyM704FKqmsgmIuIOt1vg6UB9t33dt8/hdLGv7+N4oTGmwPm1cJA1Spj42HgWT1rM1tZXiElN1YNNRERc4naAe3o7YIzp8ZizP2CtDfTy0g1AddiYen5PY+pOuG8wxmyoq6u70LrHNJ/Xx476HcTlLFKAi4i4xO0AD+DMPA/Tfbu7Vdba6t4OWmtruk1aexG4q4fzyq21udbaXK/Xe771CsEAb+ts48CyqbTs2kXnmTNulyQiMua4HeD1nNsK9wD01MJ2xsx7DW/nnO7j536g1wlvcuFCTybbMzMOOjtp3rbN5YpERMYeV+8Dt9bWGGMC3Xan03tIpwN5xpjQ9iVAlnMfeWgme5UxJq3bBwBNlY6gjHEZTEuZxo5xDbwVaNq8mfGXXOJ2WSIiY0o0LOSy3hhTEHYrWT5QFjpojMkCljtj2tWEhbszQS3LWlsatq+4W3jfSPB+cImgHG8OG49sJG76NM1EFxFxgdtd6DirqGUZY/KcQK7tdl94AT3cF+6cu9J57eqwSW+VzvZqZ6GXqgtZ2U3OzzLvMl5vep1TK+ZpIpuIiAuioQVOeAu6l2PnHHdC+Zxgdiaw9Xo9iYzQOHjtooks/uUR2o6+TvzkTJerEhEZO1xvgcvINC9tHkmxSezMaAXQuugiIsNMAS4DEh8Tz+KMxWzrPATx8TTryWQiIsNKAS4D5vP62BnYhVk0TxPZRESGmQJcBszn9dHe2c7BFdNo2rYN29HhdkkiImOGAlwGLPRksj2zE7BnztCyd6/LFYmIjB0KcBmwSeMmMSN1BjuSGwFo0ji4iMiwUYDLoPi8Prac3kPMxAm6H1xEZBgpwGVQfF4fx5qOceKS+TRrIpuIyLBRgMughMbBa5ek07J3Lx2nTrlckYjI2KAAl0GZmzaXcXHj2DW5HayleetWt0sSERkTFOAyKHExcSzNWMp2cwRA94OLiAwTBbgMms/rY/eJWmzWDE1kExEZJgpwGTSf10e7befApTNp2rwJa63bJYmIjHoKcBm0HG8OAHuzxtFRd4z2115zuSIRkdFPAS6DlpaUxqwJs9gx4SSgBV1ERIZDxALcGDPBGDM7UteTkcXn9bGl2Q8J8ZrIJiIyDC44wI0x9xpjnjHGrAkFtjFmHbAR+JIxZp2CfOzxeX00tDQQyL1YE9lERIZB3ABe8yJQZq3dB2CM+SKw3Fo7N3SCMeYLwLcjU6KMBKEFXfw5GaQ9/E9sWxsmPt7lqkRERq+BdKGnhcLbUQSUdTtnHzKmXOy5mOT4ZHZNtdiWFpp373a7JBGRUW0gAd4Q+o0xZiKQBVR3O0f3EY0xsTGxLMlYwvb4OgCa1Y0uIjKkBhLg4eFcCASstS93O2fSgCuSEcvn9bH31Cu0ZqZpIpuIyBAbyBh4ozPG3QiUAAWhA8aYG4AvASsjU56MJD6vjw7bwatvnkWKbiUTERlSF9wCt9b+AXjC2VxhrX0SuiazZQHrgbyIVSgjRmgi2955ybTu20dHY6PLFYmIjF4DaYHjTGJb2213GZBurX1lsEXJyDQxcSJzJs5hZ/sZ3gU0bdlKymX/6nZZIiKjUqTvAy/WfeBjm8/rY3vbq1hjaNqsbnQRkaGi+8AlonxeH0/tfYrjOTNI1UQ2EZEho/vAJaK6FnRZlknT5s16MpmIyBDRfeASUdmebFLiU9g93dDR0EDbwYNulyQiMirpPnCJqBgTQ443hx2JxwFoelnj4CIiQ2EgAd5ojPmCMeZ2gveB3xY6YIy5wRjzIlAVqQJl5PF5fdQ2H6Q5NVEPNhERGSK6D1wizuf10Wk7OfCWWZqJLiIyRAZ1H7jzDPArnN1l1toTkSstyBizGvAD6c57l1/Aa8ustUWRrkn6ttS7FIC98ycw/7lNdLa2EpOQ4HJVIiKjy0C60HGCez0QIDiBrRpocO4PnxWp4owxJYDfWlvpBHe2Maagv9eFvTYrUrXI+ZuQMIHsidnsSm/GtrXRsnOn2yWJiIw6A1nIZSJQSXCcO9taGwOkARcDfwAqjTETIlRfobW2Mmx7HcHb1vqrcXmE3l8GyJfpY1vnITpBDzYRERkCA2mB3w6stNauDd0Pbq1ttNbus9aWAvnAXYMtrJcQDnB+4+u5aCKdq5Z5l3Gi/SSvX5xOkx5sIiIScQOahW6t7fUpFdbaAMEx68FKB+q77eu+fQ6ni319BN5fBqFrQZfcqZqJLiIyBAZ7H/hgzumPp7cDxpgejzn7A86HiD4ZYwqNMRuMMRvq6uoGWKL0ZvbE2aQmpLJnZhxtBw7Q3tDQ/4tEROS8DWgp1b4OOuPfFw+snLMEcGaeh+m+3d0qa233VeF6ZK0tt9bmWmtzvV7vQOqTPoQWdNk5PgBAs1rhIiIRNZAALzfGPGuMeX/4ZDVnZvptBCey3ROB2uo5txXuga5u+rM4Y+bnFd4yPHxeH/6Ww5xJitFENhGRCLvg+8CttY3GmI8DPwaeMMaEd5fXEGwFD/p+cGttjTEm0G13Or2HdDqQZ4wJbV8CZDn3kVdaayMxLi8XwOf1YbG8knsRXrXARUQiaqALufiBK40xc4DQbPGabk8pi4T1xpiCsFvJ8gl78pkxJovgo0wrna7z6rBjhUCWMzNeXJCTkYPBULvIw5J1W7CdnZiYAS09ICIi3Qzqp6lz69gTzlfEHyHqrKKWZYzJcwK5ttt94QX0cF+4c+5K57Wre5v0JkMrJSGFi9MuZpe3jc7GRlpf2e92SSIio8aAWuDhjDFvInjfdTZwjODkM7+19o+DvTZAXy1o59g5x51V2857yVUZOj6vj9+feDq4oMvmTSRmzXG7JBGRUWHQ/ZnW2pecRV2+BLxEMFCfHXRlMir4vD5OdZzhtenjNBNdRCSCIjog6TypLC/S15WRK7Sgy77cizQTXUQkgiIetNbaGnQ7lzhmT5jNxMSJ7JmTQPOuXXQ2N7tdkojIqDBULWXdsiUAGGPweX3sSDkB7e00b9/hdkkiIqNCnwEe9qzvC6V1M6WLz+vjlfajnEoKTmQTEZHB668FvnKA1/UM8HUyCnU92GTJJD2ZTEQkQvq7jSzPGPN+oNenj/Uid4D1yCi0NGMpMSYG/5J0VjyriWwiIpHQX4BnA08M4LqReBqZjBLj48cz1zOX3Z3NtB0+TPuxY8RlZLhdlojIiNZfgFfTw0pn/TDoedzSjc/r47eNv6bTQNPmzaReMdDpFSIiAv0H+IDWNzfG6DYyOcuyzGWs372eg5mxeDcpwEVEBqvPSWzO6moXbKCvk9GrayLb8smaiS4iEgFaMU2GxYzUGaQlprE3exzNm7dgOzrcLklEZERTgMuwCC3osnPiaTpPn6bVr7V+REQGQwEuw8aX6eNA5zFnQRfdTiYiMhgKcBk2oXHwvReP14NNREQGSQEuw2bxpMXEmlhqczLUAhcRGSQFuAyb8fHjmZc2j90XQcvu3XSeOeN2SSIiI5YCXIaVz+tjV8JxOm0nzdu2uV2OiMiIpQCXYeXL9HHGtnDAix5sIiIyCApwGVahiWy1i9M1kU1EZBAU4DKspqdMJz0pndp5KZrIJiIyCApwGVbGGJZ5l7Ej/QztR4/SdvSo2yWJiIxICnAZdr5MH4cIcGKcxsFFRAZKAS7DLjQOvntmHM3qRhcRGZD+HicqEnGLJy0mzsThXzpRE9lERAZILXAZdklxScxPn8+eGbE0bd2KbW93uyQRkRFHAS6u8Hl97EpqoL25iZa9e90uR0RkxFGAiyuWZS6jmTb2Z6JudBGRAVCAiyvOejLZZs1EFxG5UApwccXU5Kl4x3nZu2CCZqKLiAyAAlxcYYwJjoNntNKyt5aOU6fcLklEZESJitvIjDGrAT+QDmCtLe/jXA+wytnMds4vDjteCKwAKpxdK4ESa60/4oXLoPi8PqoPVNM4ztK8dSvJb36z2yWJiIwYrrfAjTElgN9aW+kEd7YxpqCPl5QA66215U5w5zkfAMKtAqqcc8sU3tHJl+ks6DLN0PSyxsFFRC6E6wEOFFprK8O21wFFfZyfC+SFbfuBS8JPsNamWWuNtXaFtbYmcqVKJC2atIi4mDj2LpqoB5uIiFwgVwPcGLO8h90Bzg7oszihHB74ywm2tmWESYxNZFH6IvbMiqdp82astW6XJCIyYrg9Bp4O1Hfb1327V07XeXX3MXNnHLye8xhTF3fleHOoOLadlvpm2g8fJn7aNLdLEhEZEdzuQvf0dsCZrNbrMSekAWq7Hd5AMNRDY+r5PY2pG2MKjTEbjDEb6urqLrxyiQhfpo8W2oMLuqgbXUTkvLkd4AGcVnKY7tvnsNYGnElspQQDuiLsWE23SWsvAnf1cI1ya22utTbX6/UOrHoZtGXeZQDBbnStyCYict7cDvB6zm2FeyAY0t1Pdlre3WecVwEFYed0Hz/3Exwnlyg0JXkKmeMz2btggp4NLiJyAVwNcGeGeKDb7nSgupeX5AIlvXWvG2OygKoejus2sii2zLuMXZntNG/fjm1rc7scEZERwe0WOMD6bmPU+UBZaMMYkxU6bq2tBoq7tc7zgVLnuL+H4zcSvB9copTP6+No3Gnq41po3rXb7XJEREYEt2ehY60tMsasdrq+s4DabreJFRAM6dC+yrBu9ElAlTMWTh/HNQs9ioUv6LJg8ybGLVnsckUiItHP9QAH6BbAPR0rDdv2h2/3cH6fxyX6LExfSHxMPHsuTqB502b44AfdLklEJOpFQxe6jHEJsQksmrSIvXOSdCuZiMh5UoBLVPB5fexJPc2Z/fvoaGx0uxwRkainAJeosCxzGW2mg31ToGnzFrfLERGJegpwiQo+b2giWwxNm3U/uIhIfxTgEhUyx2cyNXkqe+enaBxcROQ8KMAlavi8PnZPsTRv0pPJRET6owCXqOHz+qiLb6KuPUDbq6+6XY6ISFRTgEvUeGMc3OjBJiIi/VCAS9RYkL6AxNhEds+K1zi4iEg/FOASNeJj41k8aTF7spP0ZDIRkX4owCWq+Lw+/BObObV7O52trW6XIyIStRTgElV8Xh9tphP/pHZaduxwuxwRkailAJeo0vVksos0kU1EpC8KcIkqGeMymJYyLTgOrolsIiK9UoBL1PF5feyeZjijJVVFRHqlAJeo4/P6OJ7QypHAq7TX17tdjohIVFKAS9TpGgefZtSNLiLSCwW4RJ15afNIik1k9/QYmhXgIiI9UoBL1ImPiWdxxhL2ZCVpJrqISC8U4BKVfF4f+9Jaady+GdvZ6XY5IiJRRwEuUWmZdxntxrJ3/ElaX9nvdjkiIlFHAS5RKcebA8Du6YYm3U4mInIOBbhEpUnjJjEjZQa7Z8bpwSYiIj1QgEvU8mX62DMjVgEuItIDBbhELZ/XR0NiO68e2UNnc7Pb5YiIRBUFuEQtn9dZ0GVKB83bt7tcjYhIdFGAS9SamzaXcbFJwRXZdD+4iMhZFOASteJi4ljqzWHP7ATNRBcR6UYBLlHN5/XxyqQOGrcqwEVEwinAJar5vD46jGUnr9FQUYHt6HC7JBGRqKAAl6gWWtBlX+5FHPmv/499BSs5/Y8XXK5KRMR9cW4X0B9jzGrAD6QDWGvL+zjXA6xyNrOd84uHuEQZQmlJacyaMIsD+VlcdNnnef073+HALbeQcsUVZH7xCyTOmeN2iSIirojqFrgxpgTwW2srneDONsYU9PGSEmC9tbbcCe485wOAjGA+r49NdZuYcPXVZD/9NN7Pf54zL7yA/9rrOPLNe2hvaHC7RBGRYRfVAQ4UWmsrw7bXAUV9nJ8L5IVt+4FLhqIwGT5vynwT9c31/Nf//RdH2xvIKLyd7Gd+j+eGG2j43/+l9t1XUf/oo9jWVrdLFZEh0LJ3L4eLi2k9eNDtUqJK1Aa4MWZ5D7sDnB3QZ7HWrugW+MuBqgiXJsPsuuzr+PDCD/P0vqe55hfXUPLPEhqTDVO/9lXmPPULxi1ZwtE191J77bWcrK7GWut2ySISIS179vDzr3yQm6b+ll9/6YO0HT3qdklRI2oDnOCYd323fd23e+V0nVf3NWYuI0NCbALFlxbz2/f/luuyr+NnO3/GVU9exX0v3Ufr7CnMfOhBZqwtx8THc/CTn+LATTfTtG2b22WLyCA1797Nz7/6Ib6b18SxtFjufXsDVV/4EO3Hj7tdWlSI5gD39HbAmazW6zFjTKGzWdvHeYXGmA3GmA11dXUDLlKGz9SUqXz1/32VX7z3F7xt+tso31zOVU9cxUNbHiL2LblkPfUUU776FVpqa3mlYCWHv3SXPq2LjFDNu3azzgnvpZ6F/PL9vyIjJZO7LzvKXz77EToaG90u0XUmWrsbjTF5QIW1Ni1sXxbBUE6z1gbO4xpVQMBau7Kv83Jzc+2GDRsGWbEMtx3Hd3DfS/fxl0N/IWNcBkU5Rdww9wZizjRzvLyc+kcehdhYJn3sY0z62K3EjB/vdskich6ad+1i/dc+zHfzm1niWUjZex4mOT6ZQ6cO8eGnbqSjsZHSF+dyyQ8fIzYl2e1yh5wxZqO1Nrf7/mhugddzbivcA9BTeDst7+4zzquAvmatywi2cNJC7s+7n0fe/QgzU2fyzRe+ybVPXcvTr/+JSZ/7LFm/e5qUy9/BsR/9iNp3X0XgyV9gOzvdLltE+tC8axcVXw2G92LPAn58zf+QHB8M6Wkp03jwPY/SkTqe/3jTbjZ99rYx/aTCqA1wa20NwUlr4dKB6l5ekguU9NW9LqPTiskreOTdj3D/O+9nQsIEvvzXL1Pw6wL+0rmLad/9LrN++lPipkzhtS9/mX0FBZx+4Z9ulywiPWjeuZOKr36Y71wZDO+yax4mJSHlrHOyPdn8+KoHOelJ5D/nb2HHv985Zu9AidoAd6zvdt93PlAW2jDGZIWOW2urgeJurfN8oHQ4ChV3GWN46/S38vP3/Jxvv/3btHe289nnPsuHnv4Qm6e2MPvnP+Oib3+bjkCAAzffzKuf/CStr7zidtki4mjesaMrvBf1Et4hS71Lue/K+zmaGcd/Tfsne4s/j21vH+aK3Re1Y+AhTrd4DZAFZ6/E5hzLt9bmO9tZvNFlPgk4bq3tN8A1Bj76tHe28+vaX3P/pvs5cvoI/zLlX/j08k+zJHUe9Y/+hONlZXS2tpL+oQ+ScccdxHo8bpcsMmY1b99Oxd038Z13tbDQM5/yax4mNSG139f98cAf+fxzn2XRvg7uab6aWfeUYGKivV164XobA4/6AB8OCvDRq6WjhYpdFazdspb65noun3E5n3rTp5jTkUbdD+4jUFlJTGoq3k/cSdq//RsmIcHtkkXGlObt26n42k18590XFt4hv9z7S/7z//6TS3d2cve4lVz0X1/BGDOEFQ8/BXgfFOCj3+m20zy+/XEe2fYIp9tOc03WNdy57E68h5t4vaSE03/7GwmzZpG5+oukXHHFqPsBIBKNmrZt44m7b+Y772phftp81l5geIc8tv0xSl8s5R2bOvmPqTcz+QtfHFX/hxXgfVCAjx2NLY08tPUhfrrjp3R0dnDDvBu4fentJG/cxdGSUlpraxl/6aVkFq9m3OLFbpcrMmo1bd3Gk3ffxLff3Tqo8A750Us/4sebf8x7Xujkc0vuJPPOT0SwWncpwPugAB97Xj/zOuWby3li9xPExcTxgYUf4KMLboJfVVH3g/voCASY+L734f3sZ4ifPNntckVGlaYtW3ny6zc74T2P8mseZkLChEFd01rLvS+s4ae7fsa/Pd/Bx9+2mkm33BKZgl2mAO+DAnzsevXkqzzw8gP8xv8bkuOTuWXxLXxgxntp/p/HqX/0JxAXF1wI5taPaiEYkQho2rKVJ+++hW9f1RKx8A7ptJ38x1++zG/2/ZaPPdPBze/9Kmk3rur/hVFOAd4HBbjsadjDD1/6IX989Y+kJ6Vz29LbeF/yv9L4/fs4+bvfE5eZifdzn2Pie68blbNcRYZD05YtTni3Mi9tLmuveSRi4R3S1tnG5//4OZ4/+Dyf+nUnq265l4nXXhvR9xhuCvA+KMAlZHPdZn7w0g944bUXmJI8hTt8d5B3YgbHS75N8+bNJC1aRGZxMcn/cqnbpYqMKE2bN/OLu2/hW1e3Mdczl7XXPMzExIlD8l4tHS18/NkiXjqykS8+abnuM98nNa/XB1lGPQV4HxTg0t0/XvsHP6j5AVuObWH2hNl8wncn/7K1jWPf+z7tr71GSt47mfyFL5Awe7bbpYpEvaZNm3jq7o9Sek0bc9PmsvbqoQvvkNNtp/nY7z7KnmM7+Y8Ky5VffoCUt142pO85VBTgfVCAS0+stTz36nPc99J97A3sZWH6Qj65pIiFz+yhvnytFoIROQ9NL7/MU1+/dVjDO6ShuYGbf/sRjja8ylfWG97xzbWMv+SSYXnvSFKA90EBLn3p6Ozg6X1P86OXf8ShU4dYnrmcT2TdxIzHnyfwxBNaCEakF2deeolffv1WSt/TzsVpc3lwGMM75MjpI9z0249wuuF1vl4Ry1u+9yjjli4d1hoGSwHeBwW4nI+2jjae3PMkZZvLqGuq47Jpl/Hx9OtI+2EFp//29+BCMMWrSbn88lG1iITIQJypeYlffeNjlL6njey0i3nw6keGPbxD9p/Yz02//TAxDSf4+hOJ5N7/OEnz57lSy0AowPugAJcL0dTexM92/oyHtjzEidYTvGvWu/ho+6UkffdRWv1+kpYuZfzyN5E4fwFJCxeQkJ1NjFrmMoacqXmJX999KyXXtbse3iE763fy0d/dwoTjTXzjVynkrH2cxDlzXK3pfCnA+6AAl4E40XqCR7c9ymPbH6O1o5Xrsq7lA/svIumXf6Jl925s6DnFcXEkZmWRuGA+SQsWkrRgPokLFhCXnu7uH0BkCJypqeHXd3+MkuvayfJk89A1j7oe3iEvvf4Stz9zGxcdbefuZzwsevinJEyf5nZZ/VKA90EBLoNxvOk4D255kHW71gFwXfZ1LPe+iXktHibvP0nbrj0079pJy46dtL/+etfr4rxeEhcuIGn+AifcF5AwezYmNtatP4rIoJzZuJFff/02J7wv5sGrH8aT5HG7rLP89dBf+dQfPsncQ5av/Hky8x59nPjMTLfL6pMCvA8KcImE1069xgObHuB3+35Hc0ew9T0ubhwL0xeyaNIiFk1axPz46Uw5eIa2nbtp2bWT5p27aKmtBedZxiYpicS5c4OtdKcLPnH+fGJTen4uski0OLNhA7/+xu2UXNfOHE82D139SNSFd8jv9/2e1X9ezZv2wZc3ziD7J48Tl5bmdlm9UoD3QQEukdTe2Y6/0c/249u7vnbV7+o11BdOmMfUunbadu2hZedOmnftomXHDjoaG7uuGT99erCVHgr1BQuInzZNk+UkKpx58UV+/c1CSq/tYHZaVlSHd0jF7gru/vvd/OtOwxd2ZTPnkUeJnRDZVeEiRQHeBwW4DLULDfVF6YuY1pzkhPoumnfupGXnTlr37wfn/2xMSgqJ84Nd76Eu+MS5c4lJSnLzjypjzJkXX+Q33yik5LoOZqXN4aGrHyEtKXpbs+Ee2vIQ36/5Ple+BJ88uoRZDz0Ulc88UID3QQEubrjQUF88aTEz4r207/XTvGNncFx95y5adu2i88yZ4EVjYkiYM4ek+cGJcl0T5rxetdYl4k7/85/89huFlLy3c8SFd8h3N36Xh7c+zPV/s9zaegkzfvxjYhIT3S7rLArwPijAJVr0F+rj48azIH3BWaE+M2UGnYcOB8fTnXH15p07aD/8Wtd1Y9PTu42rLyAxaw4mPt6tP6qMcKdf+CdPf7OQe9/XySzPyAxvCK64+LW/f40n9jzBR/7QwQfGv53p9/0gqv5vKMD7oACXaNZTqO+s30lLRwvQc6jPmjALTp4KjqeHdcG37N2LbW0FwMTHEz9rJvGTpxA3eTLxUyYTlzmZuCmTiZ88mbjJk4lNS1PLXc5x+h8vOOFtmZU2hwevfpj0pJF7W2RHZwer/7yaZ/c/y8d/28H7pr+bi771rai5I0QB3gcFuIw03UN92/Ft7Krf1XuoZyxm9oTZmPYOWvbto2VXMNRb9++n/ejrtB85QvuxY13j6yEmIYG4zMxgqGdOJm7KFOInZxLnBHz85MnB7vkoaq3I0Dr9j3/w9DeLuPd9lplps3no6kdGdHiHtHW08ck/fpJ/HPo7n/1FO+9a/H6mfuPrUfH4YAV4HxTgMhpcSKgvzljMokmLmJU6i9iYYCvDtrXRfuwY7UeP0nbkKO2vH6Xt6FHajxwN7jsa/DXUgu9iDLEZk95oyU/OJG7yFOImZxI/ZQpxmcF9McnJw/0tkQg7/fe/8/Q9H+8K7wevephJ4ya5XVbEnGk7Q1FVEdte30zxz9t4+2UfYvJ/fNn1XigFeB8U4DJa9RfqBkNKQgoTEiYEvxInvPH7btupCamkJqSS3AzJgWaSjp/CHnUC/+iRYEveCfrOEyfOqSUmNTUY6uqyH5FO/+1v/O6eO1jzfsuMtNk8NMrCO6SxpZFbn7mVA8f9/OdjzbzlPYVkfv5zrtakAO+DAlzGkvBQP3jyICdbT3Ki9UTwq+XEWb9v7Wzt81rj48b3GPopZhwprYbkJkg+1cq4xlbGHz9NYl0j4440kvRqHeb14/132YeCfvJkEmbODK4rH2UzhMeCU//3fzxzz53cc/3oDu+QY03HuOl3NxEIHOErDzex4kOfJePjRa7VowDvgwJcpGctHS3nhHrX73sI/BOtJ7o+EDS1N/V57aTYJFLjkkk1SaR0xJPcFktykyX5dDvjG1sY19DEuLqTJJ9qJ7nZMukkZJyKJWnWLBLnzyNp3jwSna/4adOiYqxyNDr11//jmTV3suZ6y7S0WfzPVY+M6vAOOXTqEDc9/RHaTpzg7gdPk3Pnl0m/6SOu1KIA74MCXCTy2jrazgr7k60n+wz88H2n2k71eM1xNp6Zp5OY8Vo7F+0/xcw6mFFnSbfjSJr7RqAHv+ZG9fKYI8Gpv/w1GN43wLS0WTx01cNkjMtwu6xhUxuo5ebf3cy4ky18tfwUi+/6Bp4bbhj2OhTgfVCAi0SX9s52Tred7gr0xpZGDp46yN7A3uBXw14aWhq6zk/tSGDmiXimH2ph+qEWZtRZZhyDtNTMrkBPmh/8VY93PT+n/vIXnr3nE9xTABd5ZvI/Vz8ypsI7ZEvdFm579ja8jZavrD3FvG98i4nXXDOsNSjA+6AAFxl5jjcdfyPQnVDfG9h7Vus9vS2RGQ0xTD/QxIyjHUw/ZpnREMvEi2aTOG9ucMW6UDf8RRepG94RCu81BTB1DId3yAuvvcCd1Xcypz6O/3jkDBd/5wekXnH5sL2/ArwPCnCR0cFay9EzR7sCfU9gD3sDe/EH/F2r2QFktiQyow6mvdrEzGOWGXWW6WfGMSF7/lld8Enz5hHr8bj3B3LBqT//map7PsE9K43CO8wfD/yRz//p8yypS6L48TNk/ejHJP+//zcs760A74MCXGR06+js4PCpw12BHvraF9hHuw0+yjXGGqacSWD60XZmHG5zuuEtM+IySZ47vyvQE+fPJyEra1R2w596/nmq1nxS4d2LX+79Jf/5f//JWw4l85nKFuY8+BDjly8f8vdVgPdBAS4yNrV1tnHgxIGzuuH3BPbw6olX6aQTgLhOw7ST8Uw/3Mr01zuYUQcz62OYljaL5HnzSZwX1g0/7aIRex/7yeeeo7rk09yz0jDFM4P/ueoRvOO9bpcVdR7f/jglL5bwTv94in7bwexHHmXcksVD+p4K8D4owEUkXEtHC/sa97GnYc9Z4X749OGucxI7YpgWiGH64VZm1tlgsJ8ex7Spc4mfPIXYCanEpKSe/WvqBGJTU4iZMIHYlOCvMcnJro+9n3zuOarv/TRrVsUw2TNd4d2P+1++nwc2PcC128dz8/OG2Y89RuLcuUP2fr0FeNyQveMFMMasBvxAOoC1tryPcz1AobN5CVAVfr4xphBYAVQ4u1YCJdZaf+QrF5HRKDE2kQXpC1iQvuCs/afbTlMbqD0r1HdM28Ofm485Z7Qyvn0H6Wd2BRexOd1B8rEOkpshpQlSmm3w982Q3GRJaYaUFkNqbDJJyROISU0lNjU1GPCpKcSkTiAmNYXYsF+DHwSc85xfzSC680/+MdjyXrMqhkzPNIX3ebjDdweNLY38lJ+SeiaJG269ldmPP07CrFnDWofrAW6MKQFetNZWhraNMQWh7R7cZa0tDnt9rTGme+ivIhjyNcDtCm8RiYTk+GRyvDnkeHPO2t/Y0njWTPjjzce77ms/0tLIiZYTnGrv+d72oGYSO1tJaQ+Q0hoTDPkzluRTHYw/0kbKmc6u4E9phuTQB4EmSG6GuMSknoM+rAcgZkJY6E+YQExKCq179/LH/y5WeF8gYwzFlxZzsvUkP+XXjG+yXP3RjzL78ceJv+iiYavD9QAHCsMDGVgHlADnBLjT+s7qtrsMKAa6Atxaq9UbRGTYTEycyIrJK1gxeUWv57R3tnOq9VTXfe2hxWu6ft9ygsbWxq573xtaG9nv/L6/Ve3GdxpSO1pJbqsnpTVASpMl+Uwn40+1k3y41VnNLvQB4I3wfzUD1nwgAa8T3pnjMyP9rRm1YkwMX/vXr3Gy7SQPvvVPjG8KcMVHb2XW448R5x2eD0GuBrgxpqfpewEgr4+X5RljssJa1QHODXURkagSFxOHJ8mDJ8lzwa9t62gLhnvYanXhwd/9w8DRsA8DrZ0W6P251jNTpvGwwntA4mPi+fbbv80d1Xfwoys3Mv5XRzG3foyZP3l0WFYBdLsFng7Ud9vXfbuLtTYAdP+u5APV4TuccfB6zmNMXUQk2sXHxpMxLmNAt3Q1tzf32uJv7WzlvdnvVbf5ICTGJnLfFffxsWc+xvfeu5ukn7+Cub2QmY88TGxKypC+t9sB7untgDHG4wR2r5wu9TzgnWG7NwCBUAvdGFNhjKnvPqbuhHwhwMyZMwdSu4hI1EuKSyIpLonJyZPdLmXUSo5P5oG8B7j59zdT+m+v8c3nzjCjuRmGOMDdXjcwgNNKDtN9uy9rgZXW2prQDmttTbdJay8Cd3V/obW23Fqba63N9Q7TeIWIiIxOaUlplOeX4xmfzt3vOsmRxOb+XzRIbgd4Pee2wj3Q1V3eK+fWszJrbffu8+7j535g6JfKERGRMW1K8hTKryznqjlXDUuPh6sB7rScA912p9NtTLs7Y0wBUBMK71BoG2OygCqnaz2cbiMTEZEhN2vCLO76l7uIj4kf8vdyuwUOsN4J5JB8greGAcFQDj/uhHU6sMEY43FCezmA03Ve3K31fiPB29JERERGDbcnsWGtLTLGrHaCOQuo7TbhrIBgqFc6LesqZ39Z2Dnh51c63esAk+i2UpuIiMhooLXQ0VroIiISvXpbCz0autBFRETkAinARURERiAFuIiIyAikABcRERmBFOAiIiIjkAJcRERkBFKAi4iIjEAKcBERkRFIC7kAxpg6YH8EL5kBHIvg9UTf06Gi7+vQ0Pc18sby93SWtfacx2YqwIeAMWZDT6vmyMDpezo09H0dGvq+Rp6+p+dSF7qIiMgIpAAXEREZgRTgQ0NPP4s8fU+Hhr6vQ0Pf18jT97QbjYGLiIiMQGqBi4iIjEAKcBERkREozu0CRhNjzGrAD6QDWGs1ZjNIxhgPUOhsXgJU6fsaWcaYMmttkdt1jAbOv9e7gFpn1wZrbY17FY0Ozs/WgLPpsdaWulhO1FCAR4gxpgR40VpbGdo2xhSEtmXA7rLWFoc2jDG1xhh9OIoQ599tltt1jAZOeFdYa/Od7UKCYb7SzbpGOmPM6vDANsYs775vrFIXeuQUdgvrdYBaNYPg/EDsHi5lQPG5Z8uFMsYsd7uGUWYtwX+fIevRv9VIuDF8w+nRuMSlWqKKAjwCevlBGADyhrmU0SjPGBMe4gHUYoyUXKDK7SJGkQKg2hiTZYxZbq0NWGv9bhc1CtQbYypCG07PxjoX64kaCvDISAfqu+3rvi0XyPkBmNbth2A+UO1WTaOFMaaAYAtRIiDsQ3xu2L4KpxdJBqeI4Af5BmcsvF5Dk0EK8Mjw9HZA/4Ejx/le5qFuyUFxvo8Ba23A5VJGk65eIWut3+nmXUewW10GwfkAvwbYAJSg7vMuCvDICODMPA/TfVsGby2wUrN6B22VtVa9GJEVcH7dELbPT7BbXQbBGFMGVDuTA/OBwvAu9bFMs9Ajo55zW+EeCHYDD3Mto5LTdVam4Bkcp6tX38PI88M5/98DEOzx0M+BgXH+vQZCH9qttdXGmDnAPncriw4K8Aiw1tYYYwLddqejH5QR4YzX1oTC2xiTpyAfsHSC44mh7UuALOcDUqUmXQ2MtdZvjAl0C2sPGqoYrHTgePgOa23AGKP//6gLPZLWO0ETks/Zt5TIABhj8gj+J95gjPE4M9J1+9MAWWurrbWloS+Cs9ADzrbCe3DWAKvCtm909skAOR/U88P3OXM49G8VPcwkopxWTA3OhBYtNjI4zn/Uhh4OVVprtTjGIDm346wkOHN6DVCu1uLgOD8DumixkcFzPrQX8cbqdvrZ6lCAi4iIjEDqQhcRERmBFOAiIiIjkAJcRERkBFKAi4iIjEAKcBERkRFIAS4iIjICaSU2kTHAGFNC8EEwoaVU/byxfrcHZ4U25/fZbi7q4tz3W0bw/vSVWnVPpGcKcJExwFpb7ARjLVDSUyg6x6sILkTkWoA7Hx7yjTG1/Z4sMoapC11EgK7gLKKPx+MOs4DbBYhEMwW4iHRxWuZZ/Z4oIq5TgIuMcc4jG8MF3KhDRC6MxsBFJI/gQ3iANx4U4QR7CcHJZO/kjZZ5FpBNcCz9rLFy5wE0hZw9Qa7Hh6Q4595F2EMqgPXdz3WeSOchONEuH7hdD10RUYCLjGlOSBcB5zw1y1pbQ3AyWQNQZK0tCntdFrDRGLMiFOKh2ePW2vyw8zxAhTFmZXjoOvs3AvnhHwKcp3mF15LvXDP0HgAVdHvEpMhYpC50kbGn2BhTZowpA9bS/5h3PcHZ6V2cQC3n7Gfel3XbxgntCud9wq0l+FjY8PDOI9jiD+fp1srfQLDHQGTMUwtcZOw56zYyY0zFebwm0MO+dcBqpzUduo+8p+e0VwNlxhhPWCu8gHNb0RsI9gaE23gedYiMSWqBi8hZrWtjjMcJ5f6EWsa5BBeIoaex6bAWdK5z/dCkOX+38wKh8fcw9edRh8iYpAAXGeN6CM1VkZwkdp4fBkTkAinARaQ7z3meFxo734Azi72XsE4POy80OS789SIyAApwEenukh72eXrYVwRUO13ffoIhvqqH8wpC54Xtq6SH8XJjTFYP96WLSA8U4CLSxZnQ1tO4843dzlvOuZPWVgJF4a1w5/c3cu7ktNuBvB7CuiCshe5BRHqlWegiY0DY08gASowxG+j9aWTFPVyizBhT4Pw+HVgBrAhvVVtr/caYdwJ3hT2IJBt4Z/cxdWttwBizwqmllmDr3WOtLXXuJy8m2MV+lzEmy9lfgPNBwPmgUaYnlclYZqy1btcgIlHMCdgihaVIdFEXuoiIyAikABcRERmBFOAi0p90NKFMJOoowEWkR84tXRUEw/su50EjIhIlNIlNRERkBFILXEREZARSgIuIiIxACnAREZERSAEuIiIyAinARURERqD/H2SFlrlVCIxTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots( 1, 1, figsize=(7,5) )\n",
    "\n",
    "c1 = 'tab:red'\n",
    "c2 = 'tab:green'\n",
    "\n",
    "axs.plot( trn_losses, label=\"train loss\", color=c1 )\n",
    "axs.plot( val_losses, label=\"val   loss\", color=c2 )\n",
    "\n",
    "#axs.set_yscale('log')\n",
    "\n",
    "axs.set_xlabel( \"Epoch\", fontproperties=axislabelfont )\n",
    "axs.set_ylabel( \"Loss\", fontproperties=axislabelfont )\n",
    "\n",
    "xticks = [ int(x) for x in axs.get_xticks() ]\n",
    "axs.set_xticklabels( xticks, fontproperties=tickfont )\n",
    "\n",
    "yticks = [ np.round(x,2) for x in axs.get_yticks() ]\n",
    "axs.set_yticklabels( yticks, fontproperties=tickfont )\n",
    "\n",
    "axs.legend( loc='best', prop=tickfont )\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This training could definitely benefit from more epochs, some regularisation, and maybe early stopping.  But here we just want to show how to get the DGCNN up and running."
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}