data/SBMs/generate_SBM_PATTERN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook for generating and saving SBM PATTERN graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import pickle\n",
    "import time\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.sparse\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate SBM PATTERN graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'nb_clusters': 10, 'size_min': 5, 'size_max': 15, 'p': 0.5, 'q': 0.25, 'p_pattern': 0.5, 'q_pattern': 0.25, 'vocab_size': 3, 'size_subgraph': 10, 'W0': array([[0., 1., 1., 1., 0., 1., 0., 0., 0., 0.],\n",
      "       [1., 0., 0., 0., 1., 0., 1., 0., 1., 1.],\n",
      "       [1., 0., 0., 0., 0., 1., 0., 0., 0., 1.],\n",
      "       [1., 0., 0., 0., 0., 0., 1., 0., 0., 1.],\n",
      "       [0., 1., 0., 0., 0., 1., 0., 1., 1., 1.],\n",
      "       [1., 0., 1., 0., 1., 0., 0., 0., 0., 0.],\n",
      "       [0., 1., 0., 1., 0., 0., 0., 0., 1., 1.],\n",
      "       [0., 0., 0., 0., 1., 0., 0., 0., 0., 1.],\n",
      "       [0., 1., 0., 0., 1., 0., 1., 0., 0., 0.],\n",
      "       [0., 1., 1., 1., 1., 0., 1., 1., 0., 0.]]), 'u0': array([0, 2, 2, 2, 0, 0, 2, 2, 0, 0])}\n",
      "<__main__.generate_SBM_graph object at 0x12ce33890>\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def schuffle(W,c):\n",
    "    # relabel the vertices at random\n",
    "    idx=np.random.permutation( W.shape[0] )\n",
    "    #idx2=np.argsort(idx) # for index ordering wrt classes\n",
    "    W_new=W[idx,:]\n",
    "    W_new=W_new[:,idx]\n",
    "    c_new=c[idx]\n",
    "    return W_new , c_new , idx \n",
    "\n",
    "\n",
    "def block_model(c,p,q):\n",
    "    n=len(c)\n",
    "    W=np.zeros((n,n))\n",
    "    for i in range(n):\n",
    "        for j in range(i+1,n):\n",
    "            if c[i]==c[j]:\n",
    "                prob=p\n",
    "            else:\n",
    "                prob=q\n",
    "            if np.random.binomial(1,prob)==1:\n",
    "                W[i,j]=1\n",
    "                W[j,i]=1     \n",
    "    return W\n",
    "\n",
    "\n",
    "def unbalanced_block_model(nb_of_clust, clust_size_min, clust_size_max, p, q):  \n",
    "    c = []\n",
    "    for r in range(nb_of_clust):\n",
    "        if clust_size_max==clust_size_min:\n",
    "            clust_size_r = clust_size_max\n",
    "        else:\n",
    "            clust_size_r = np.random.randint(clust_size_min,clust_size_max,size=1)[0]\n",
    "        val_r = np.repeat(r,clust_size_r,axis=0)\n",
    "        c.append(val_r)\n",
    "    c = np.concatenate(c)  \n",
    "    W = block_model(c,p,q)  \n",
    "    return W,c\n",
    "\n",
    "\n",
    "def random_pattern(n,p):\n",
    "    W=np.zeros((n,n))\n",
    "    for i in range(n):\n",
    "        for j in range(i+1,n):\n",
    "            if np.random.binomial(1,p)==1:\n",
    "                W[i,j]=1\n",
    "                W[j,i]=1     \n",
    "    return W    \n",
    "\n",
    "\n",
    "    \n",
    "def add_pattern(W0,W,c,nb_of_clust,q):\n",
    "    n=W.shape[0]\n",
    "    n0=W0.shape[0]\n",
    "    V=(np.random.rand(n0,n) < q).astype(float)\n",
    "    W_up=np.concatenate(  ( W , V.T ) , axis=1 )\n",
    "    W_low=np.concatenate( ( V , W0  ) , axis=1 )\n",
    "    W_new=np.concatenate( (W_up,W_low)  , axis=0)\n",
    "    c0=np.full(n0,nb_of_clust)\n",
    "    c_new=np.concatenate( (c, c0),axis=0)\n",
    "    return W_new,c_new\n",
    "\n",
    "\n",
    "class generate_SBM_graph():\n",
    "\n",
    "    def __init__(self, SBM_parameters): \n",
    "\n",
    "        # parameters\n",
    "        nb_of_clust = SBM_parameters['nb_clusters']\n",
    "        clust_size_min = SBM_parameters['size_min']\n",
    "        clust_size_max = SBM_parameters['size_max']\n",
    "        p = SBM_parameters['p']\n",
    "        q = SBM_parameters['q']\n",
    "        p_pattern = SBM_parameters['p_pattern']\n",
    "        q_pattern = SBM_parameters['q_pattern']\n",
    "        vocab_size = SBM_parameters['vocab_size']\n",
    "        W0 = SBM_parameters['W0']\n",
    "        u0 = SBM_parameters['u0']\n",
    "\n",
    "        # block model\n",
    "        W, c = unbalanced_block_model(nb_of_clust, clust_size_min, clust_size_max, p, q)\n",
    "        \n",
    "        # signal on block model\n",
    "        u = np.random.randint(vocab_size, size=W.shape[0])\n",
    "        \n",
    "        # add the subgraph to be detected\n",
    "        W, c = add_pattern(W0,W,c,nb_of_clust,q_pattern)\n",
    "        u = np.concatenate((u,u0),axis=0)\n",
    "        \n",
    "        # shuffle\n",
    "        W, c, idx = schuffle(W,c)\n",
    "        u = u[idx]\n",
    "    \n",
    "        # target\n",
    "        target = (c==nb_of_clust).astype(float)\n",
    "        \n",
    "        # convert to pytorch\n",
    "        W = torch.from_numpy(W)\n",
    "        W = W.to(torch.int8)\n",
    "        idx = torch.from_numpy(idx) \n",
    "        idx = idx.to(torch.int16)\n",
    "        u = torch.from_numpy(u) \n",
    "        u = u.to(torch.int16)                      \n",
    "        target = torch.from_numpy(target)\n",
    "        target = target.to(torch.int16)\n",
    "        \n",
    "        # attributes\n",
    "        self.nb_nodes = W.size(0)\n",
    "        self.W = W\n",
    "        self.rand_idx = idx\n",
    "        self.node_feat = u\n",
    "        self.node_label = target\n",
    "        \n",
    "        \n",
    "# configuration\n",
    "SBM_parameters = {}\n",
    "SBM_parameters['nb_clusters'] = 10\n",
    "SBM_parameters['size_min'] = 5\n",
    "SBM_parameters['size_max'] = 15 # 25\n",
    "SBM_parameters['p'] = 0.5 # 0.5\n",
    "SBM_parameters['q'] = 0.25 # 0.1\n",
    "SBM_parameters['p_pattern'] = 0.5 # 0.5\n",
    "SBM_parameters['q_pattern'] = 0.25 # 0.1    \n",
    "SBM_parameters['vocab_size'] = 3\n",
    "SBM_parameters['size_subgraph'] = 10\n",
    "SBM_parameters['W0'] = random_pattern(SBM_parameters['size_subgraph'],SBM_parameters['p_pattern'])\n",
    "SBM_parameters['u0'] = np.random.randint(SBM_parameters['vocab_size'],size=SBM_parameters['size_subgraph'])\n",
    "        \n",
    "print(SBM_parameters)\n",
    "\n",
    "\n",
    "\n",
    "data = generate_SBM_graph(SBM_parameters)\n",
    "\n",
    "print(data)\n",
    "#print(data.nb_nodes)\n",
    "#print(data.W)\n",
    "#print(data.rand_idx)\n",
    "#print(data.node_feat)\n",
    "#print(data.node_label)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO19f9CfVXXn5xISSSKEJEgMhEiEUNdKVUhFSmdHS1m7bKf+A7S7U8dl7cTObFuLO1NwO9aXP9ZR24KMbrWpWN1Wt7LEWR3sKAy7/qHd4hLaoVULuIaJCeGHCYVCqC+Yu3+8732573nPOfece+/z/T6a75nJfPM+z73nnnvu8zz33nM+59wQY8SMZjSjH386adoCzGhGM5oMzV72Gc3oBKHZyz6jGZ0gNHvZZzSjE4RmL/uMZnSC0Oxln9GMThCa2MseQviFEMIDIYTvhBBumFS7VgohnBNC+N8hhG+HEL4ZQnjX4vVNIYS7QggPLf5unLasOYUQVoUQ/iaEcMfi36OVN4Rwegjh9hDCPyzq+dKRy3vd4rPw9yGE/x5COGXM8pZoIi97CGEVgP8K4F8DeDWAfxtCePUk2nbQCwD+U4zxXwB4I4D/uCjjDQDujjHuBHD34t9joncB+Hb295jlvQXAl2OMrwLwWizIPUp5QwhnA/gtALtijK8BsArAr2Ck8pooxjj4PwCXAvhK9vd7ALxnEm03yPwFAFcAeADA1sVrWwE8MG3ZMhm3YeGB+zkAdyxeG6W8AE4DsB9AINfHKu/ZAL4HYBOAkwHcAeBfjVVey79JLeOT4hIdXLw2SgohnAvg9QDuAbAlxngYABZ/z5yeZCvowwB+B8Dx7NpY5X0lgCcA/OnituMTIYT1GKm8McZDAP4AwAEAhwE8FWO8EyOV10KTetkDc22UON0QwksB7AXw2zHGp6ctj0QhhF8E8HiMcd+0ZTHSyQAuAvCxGOPrATyLES+BF/fibwWwA8BZANaHEH51ulK10aRe9oMAzsn+3gbgkQm1baYQwmosvOifiTF+fvHyYyGErYv3twJ4fFryEboMwC+FEB4G8BcAfi6E8OcYr7wHARyMMd6z+PftWHj5xyrvzwPYH2N8Isb4PIDPA/gZjFfeIk3qZf+/AHaGEHaEENZgwdDxxQm1baIQQgBwK4Bvxxhvym59EcDbF///dizs5adOMcb3xBi3xRjPxYI+/1eM8VcxXnkfBfC9EMJPLF66HMC3MFJ5sbB8f2MIYd3is3E5FgyKY5W3TBM0eFwJ4EEA/w/A707bWMHI97NY2FrcD+BvF/9dCWAzFoxgDy3+bpq2rIzsb8KLBrrRygvgdQDuXdTx/wSwceTy3gjgHwD8PYA/A/CSMctb+hcWOzWjGc3ox5xmCLoZzegEodnLPqMZnSA0e9lnNKMThGYv+4xmdILQ7GWf0YxOEGp62Wsi2UIIu1vanDTN5B2WZvJOjqpf9oZIth81Zc3kHZZm8k6IWmb2NwD4TozxuzHGeSxANt/aR6wZzWhGvakaVBNCuArAL8QYf23x77cBuCTG+BtSnTPOOCOuX78eL3vZy1bcm3/hOJ48No+N69Zgzckrv0Gl+1by8nniiSewYeNmsQ7lx/GX2kzXX/qSk/HMD17AxnVrAEAtq8k9/8JxHPr+P+LsM05364iTpSS/pgcrPfHEE0vPQ+sYt9S31s3l7S2DRT46PrS9ffv2fT/GyAvYACW8GsAnsr/fBuAjTLndWIBI3rt9+/Yo0U13PhBfcf0d8aY7H6i6b6UaPlodeo8rK9VP16/5+F8t3S+V1eRu0REny6Ta7sWjR/+n8Xx5+NLxoe0BuDcK7+zJDR8bUyRbjHEPgD0AsGvXLnEZcdn5m/HX3z2CnVvW4+a7HsRVF2/DOZvWLd2/6uJty36t9L2jx3D7voNL/FI7l52/uVg2/b1zy3pcsmMTW4fKxcmZ/n/Z+Ztx810P4rLzN+Pr33lRhvQ3VyfJkMrmZai8UjtUl6W6kiyJL8fPMz607UTc2EhlOdLGlvICsKTXfCyo3vN2LbL00IPGd+eW9Ut9tfQ5p5aXfSmSDcAhLERe/btaZl//zhHcs/8oACz9XnfFBUv3z9m0btnfVrp930HccvdDS/xSO1//zhG8YcdmtWz6+5Idm8Q6VC5OznTt5rsexC13P4S//u6RFX3M+eb1qUyavJZ2SnWpLJz8HD/P+Eh94sZG6z8lbWwpLwCqjrh2LbL00IPG9+a7HlzWx1Kfc6p+2WOML4QQfgPAV7CQn+uTMcZv1vKTvlwW0r6Qlpm3VDbJxK06ar7O3OxZqsOVpV91uhJ526Xb8cZXbmZXA3S1YpnJNFl6zMDcdalNrr3SLEfv//V3j+DKC7cs3cv5cispzyyqySn1TVpZ5isROm6elUTLzI4Y418C+MsWHomkL5eFtC+kZeYtlU0yAStXHTVfZ2Dl7GmpQ4l+1elK5I2v3CzO6HS1YpnJNFl6zMDcdalNrr3SLJffB7BsTDU9lOTWyPNsSivLRNy4eVYSTS/7EFSzN6/dz3v5a/vZodq2yEV/tZWDVEazNbSuQKxlW3mU6nP3LXqoka+mjmUV6lkVUppoPPtPve6ieO3v/4VqOPJQvuRKRpb8N1/+0La0ZavEj5Pbs3yt6VtvvjUySP2vNWJR/lpZaWnb+xnSng8A1W1qS3PJaGppR+ITQtgXY9zF1ZnozP7ksfmi4chDaZmT+NHfRNwySlu2Svw4uT3L15q+9eZbI4PU/1ojFuWvlZWWtr2fIcrHY8zz8LcYXC3t1DwfE33ZN65bg2sv36kavBJZv4CJD4Alg1Ru5Nu6YS0A2ZWVu5OSASbxSV9Yyt+ydGz9gtcsGVtWA5xMdMlPjaeWpTTHF8CyMdDcXaWtiuYOtJCk59wYl56hmiW0tjSXDKAWgzD3/Jb6P9GXfc3JJxUNXok8LqfELxmkciPfdVdcoM7ouTspyZMbtt6wY/MK/jl5DEj0nvYFr3E1tqwGJJk4PWjGIcnoVLva0gysSaaWFZCk59wYR/XQwj//m8peYxD29H/V3NycS/gW2rNnz9zu3buxbeNanHbKavzyT5+Ds05fi6su3oYNa1cDWPhyffJr+3HZ+Ztx1ulrcdn5m3HbvQexbePapTKrTgIOHD2GX/7pc3D2xnVL/BKf/P7xCHzya/uX6tOy+d8b169exjcRrWOhvM7Tzz3PypD6z/UxUdIHd48S1QvHg8pC61554RY8N398BY+SHiQ5ubF+9Vmn4bRTVi/1e9VJwG33LriVjjwzz7ZNiT4nmkyJv0e/XF89Y2EhabxS229+1ctW6EN7fjesXY0bb7zx8Nzc3B6uvalY4zUXlOXLRV0g9OtJXSya+yj/+/Z9B1nXSs0sa/mCp/5rX+ceLi1u/0n5UVCTBUBkkVMaaw78o4GXrO1xZTzgIq2vve0o0nhpbmiPK5nSaF1vOeQzv86VofuVfF90+KnnzECIHm40y77T064HyCHtl3MgRr7/zHVncdtp5IGqUr1IXhStPqczah/QwEUlnU3CjlKyb3DPvof/CpJA80P8u/jii93A/9bgi6ECEybV3jSDW3rKOYkgF9pfrt892vGQh582Xlb+UAJhJrpn/+gffXzun897s2nPY9knpzJ0z5vvrdL+MPGx7OMksuzZeu/1KD+NF72n2UYs+/sSf03PtA63t5ba4q57noe0133bpdux88xTq20tpToWXWl2CPosJh1x42WVV9uzT/Rlf/9NH5n7xikX4bRTVuPS8/Rl6Ya1C2W0zqYyt927sJdKfD/5tf1Lf7/lNS9fxifdO3D0GPbed8gkS6Kcr1SHk9tSr9RHKj/Hi95Ldc/euG6FTLfdexB77zuEs05fu4yPh7+mZ1rnrNPX4rorLjDphbvueR6+dP+j2HvfIew889QVbZbI0k5J/vzekWfmcc/+oyv0nJdJz2LSETdeVnlHY6BLfvaWPTFHFpghLdvDZzp0PS8vTztSWQ9/S3s18rbqq6e+a9uxPGetdhI3Sev7If5Z9uwHjjwbb7rzgXjPd78fb7rzgXjgyLPFOi380n2tHUuZFn4e/q1te+pK/IaWV+Pfu48txD1b0vNm0W9Lv9N9rFp9fxwgecUgNBQc0hO3LPGwyNICAbXwb23bU9cCJR0aqtrT/dXbdaYBhmqgxS39TvdXrdtwhiTvVEA1HFFDjgVwYiHNSJXftxh+LAYTyfDFlXnzq16GL93/6DIjDgUBUdKMQhZjXglEkutbMrpJhlGONHmpHqihitN3DbCHtmcB7WgkPas5YIg+b1zb9Br9O++PBPiiffvmHbc+8r73/u5HOblHM7NzX67ecEgOqGEBJXiAC5aYZwpgycEkAA96SdQSH63Vl7LRcDJYMtdY5KV6sKzmWmLqa+LRre1wmYZKsfr0mpSfIJElVv+kl6w7TZJ7qjM79+VKX3nqPqqd4TV3UU+XWyLuC15ysXAzg7SK0GC4WtmSy02rI1HrqoiuujwuJ45a+uZ5TrRnUuLDueA0mDeFgZeei9HP7NyXK/+lwQE1M3yP/XErZFWzG6QyUg46ShoMVyuryVeqY5Glpgy36mqZcVv65nlOtGdS4sNBgTWYNwcD13Q9+pmdC0KRgBA1wShcPfpV5sA1EjhFa1sDj6SvLu2btqdOcmngn9KqxwNO8YBptP5bZLDW5e55bBaUh0WXUoAJJxM31hIfLegryaPt2a2rrC9/6uaH537vvR9jC0lm+iH+aa63acFae0FJPbnltes1kM+WXPMt8veSoTUvf4lqYKgWfr359OgrxgKX1azxkmW2Za8NyDODBk2ssdpqM4H0xU77uNxOQa25ln0s1Z1nRVIjP6dbCZasWdaltjW+NZ6ApGcLfNay2qJ9qn02S6sJj30myfCVz/7xEWlmH401vmSZrfWRWvZidH9XY7VtyQKb+pxk8u5jpZz7vfbfltBZarn3jFfJi2DxGnBE9cwlH5Hq5Hx7YDUs/fYkupDk1vzso3nZS+GOvaCqrRDPGpLatoR2evgPBbvsAa3t0V5tmdY8/SU+k4DnltpI19/9B099X2Qire+H+Jf27B7oIAdzLcEMW+GWFsijB9ZbA9NsqZNk0mRr0dk04bIWHp5naWg5e0N0S4SBznqrJg900JLRtSfsMC+jbSUAe8bRod19ktyeDC0enU0TLmvh0TsjbYucY8gSnGiqOehqXRalHHQWIIfFYETb9rgKOX5a3H2NnCW5NeNeSYeaLJJBEKgDIFE+Wg7BRB6wCo1zz/l6XJwe46NUxqMfS+5AWnZ0BjpLVlIN5loCH3iOeCqVkXLRpXx1FqOPZFTyQF8t5AGplHSoyaIdwlkDQKJ8tByCiTxgFS6nW+JrWQVZYMiUeq9aJD3QsqMJhEmZamqyfVrcPNyMWfoiWtrWZPDAeOkMb3FLeeTtkYnW0h7ntkx6tvSp5P7kYKJUzx6wSiKOr2UVVAOyksBbnmy4+arlkh2bTTDq0YBqtu38ySK4oQacoYERegFLpPu9ATgeagFhTAo4NIm2e/Afip8H0FMDqKKEsRjoUqYazoVBT7qgLgYte2nuGillpNVIOgL5qovljKb0miXbp8VV4+Gj9VuiGncRbS+XzZMFt6TnnLgMudz9JJv2DGmkyVBzVDOtqx2fndqUTiTq5dKbyokwQDlfPCUN6FICWtSGp1KZPCegaO1a9uMePjVBQi02AUsuf41Kes5J25vn9xMNkdiiJTQ21dWOz05tcmVbAoMojQYuW4JBWqCDrckJPFZ9Sh54qCdLrcWC3xoGXJKtBBcFfLqn8moZWOneXNq7W8JAtf5qJxDVBA9R+S1eDss4lmxamjV+NIEwdN8ydOBHbxpqT+nZl/fovyU4p1fbNQEqQ+/Hp2kbsPAtPfurTj3jkTimQBht9uCCQ6w+aE8oak3eeI003zMlS/gknbk0C75nVSGRJWxTSxU2dIAK1Z0lcURNfzn+0krSU5a2W7OqA8rnw2nW+NEh6KTgkJyk/aZnL9wroWUizffskVOyE/Q6D04iq89fShU2dIAKlcWSOKKmvxx/gLcFeMqW2rX2sXQ+3Lt/+Py81Ofiyx5COAfAfwPwcgDHAeyJMd4SQtgE4HMAzgXwMIBrYoxPlvgBtmCDoYM6evPvxXfSQSc1PFpl6aEri2fEW9/Kv3dZj4xNYy2t79M/AFsBXLT4/1MBPAjg1QA+BOCGxes3APhgiZe2Z+8RMKDl8e7B1xOY4Qm68ORsr2nLEsDTO/jEUke6VyvLNMapB1lkseoILX72GONhAIcX//9PIYRvAzgbwFsBvGmx2KcBfBXA9f7PzQL1XIpyebx78LXwqNkmeLLA1rRlCeDpHXxiqSPdGwpSytVrHaceZJGlRzCSa88eQjgXwOsB3ANgy+KHADHGwyGEM618OABDC3AhEbc8lI5s9hx1W7P0qo2hLoGLatrK+Sd9XHnhlqW66ZfqiR4h7AEXcTyk/tC6mix0vGhZL7CH6sFSZ+eW9exR4RpZjm62HhXlBVAlMr/sIYSXAtgL4LdjjE+HEKz1dgPYDQDbt28HwH8he+T05gJBUsBKKcuKlW+NDJ46HoCMtS0ugAdYnofckxXXAi5KVAMy0mShfGjQjBfY01JH6xMljwGUo9Ysy4DxZQ8hrMbCi/6ZGOPnFy8/FkLYujirbwXwOFc3xrgHwB4A+KnXXRRvvuvBJehj/oX0fDXpV9Jz4H3pOsefuw7ALIOHWmYPiXLZkuw7t6xfag+wG03TrJLP9HkbllUAJ1feR8/KIZ/Zt25YK7Yl6aFmBUZ1p7WV606TTZKT6leCDVvIYo0PAG4F8O0Y403ZrS8CeDuADyz+fqHE68lj86J7Lf+/NxlETahozcki2v6w136uZfaQiMvplod7SiGukstNSyhiWQVIcnHtateA5TN76ptHD6lvFuJCZaW6LS5ei2u6ZgVsmdkvA/A2AH8XQvjbxWv/GQsv+W0hhHcAOADg6hIjGghDv5BSwADAf5WlX6medQam/Lh9Z9r7phmYfrlrVwHSjCPNnBaq2WN795aeWY/jJ60YSuTZ31tsL5agHG2sUx3LyqHGJkJXWR5dWazxXwMgbdAvN7WySDQQhiYT0AAW0lcZ6HfSh8RP+tJqM3DtKkDav7XMFL322Nre0jPrcXVq++bZ37cGIHkSUlj24bU2kdozEEdzIoznjDMPFJTCWGsgpRwUVjqtk9bJEw9IkFdLsg1PsoYWSHBN8I8laESrJ52YUytvzbNigRxLMFYL/JkbCy5dVn6uQF6HPhcSbPbGG288PDc3t4fr41RTSXvOsqoJyQT67IElHtqXm9tLSl9ji69Ymymk1UuvVYBEpRnNWk8LA62Rt+ZZaYGxWlaW3FhIqyELXqTGezWavPFDUY3VtQcPzz5ZK+uRj/NXDwE5ltqtrVfLpyd5ZKgpa/GhW8avSWcStG6Ifxe+9vUi5G9oKGKPnOe1EFZP3vESnNWSR7+GWvmX+tgLwlwDLbXwmzSc2pLnvkYmjCUtVXK9ASsNVUA7rJWjHjBDiZdWxgLZ1VwsEpzVkke/hlr5l4xXvSDMNdBSC78aWXq11/MZ1WgqOei0pUjLUk5zgUi/mouM8vVAWHMXVAJ7UBdTDlVNbhRpCcflMqNy1bhjEmnwUQnoo7lDqT4sEGaOJHAK5+IrucY4116NC1IDzJRcu5yrkF7jZOoB2ppaDjrAZ6CzkMVgpLnVAPsMLBFndKEGOjorUYMPZ4jjjFi0XguUUoOPSgZKzR1K9ZH3R4IwcyTN5JyLz7O6KBkULS44jl9pDDjDmiWHf4/Zfip54z2uIEtONEtObolvXied8uI5+tiSo5xmGaEuJkt+Pc7tIx1nXOM2k9xAVrlrjipOfdFchFLbmjtUco1xfbOMHx0TzQ1a0gnH15KdSTr2m+pOc71N9GV//00fmfvGKRfhtFNW49LzbO6CT35tP265+6FlddK1A0ePYe99h5Z+zzp9La674gLTg5d45HVuu/cg9t53CEeemcfe+w4ttblh7Yu/FvkAsHUS/51nnrpMzrxsid9t9x5cuv+W17x8WRuanCW5U90v3f/oki5T+xa5PUT7ksaPey6ktrmxSvrIdVTqm2X8JF2dvXHdiv6XdMLxPfLMPO7Zf3TFc8fpIclO34FUZzQv+6du/cTcr7/znS7AhQaasJzmIZHGl5sBqGyer7LUJjerlLK01s6mtA9SRldNL3Qm5mZkC8CH6o7O2pbTf9KKzDO71gJwSuOmrT6lLLAcMEt67jgZpNXFaF72Wz/xJ3N/+sH3qA+p9BXlvp7py8p9YUuk8eVmACrbtw4/bf4qS21yswr9glvk9hCdTWg7ml7oTMzNyLRP3AxEdUdn7ZwHXb3QFZlndtX6pumzNG7a6pM+D1zfUh+k546TQVpdjOZlt+zZPfnSE0l7eEu2T00GCwQ2fY1roJ7cLO7JUkv7b8loKsF8PfniU5+vvHALnps/bjpBNW+PykDHmjvFVbLLaFBjC2y4JiMtXRVp9hm6x87r0hNqW87fs+zZp+pn56jmxFNqHW316XrCKZNVuAbqqVlmATu8t2SptYR2enz+qc9JRssJqrQ9LaiDO8VV8qtbAo56++RpmCn3vEin3uR1qbw1EFiP/FP1s2vkgVRy4X8SzLCFPBDYVn6t0FzPfa8sveG4lrGW2rKMiQeqapHT0u9JQWBddSRo3RD/UnbZ3tDYFjhnrzYnlVW2NUsrlVPLMjsUlLSku9osqz2h0K18LRDpmmemRBgLXDZRbzhgC5yzV5s1kWYW6KSnHxbYpbblGQLC7IGHSjJZ+9gib2++nmxHvQ8skWgq8eyaIcJiMCkdNOg5/LDGQMO5hGpisi2uIIsrUuKnGY6ScU0DykhurxZjFicnlSk3/FEDnWYco6AaTzw/fSZzvhLYinsuSi5UzsirHavlpdEY6BJphojaYJMSvLXG8Ce1m8rm5WtismuOgNZip0vQ4FxOgM8MpEGYexizuDY4o6RkoCsZxzRYskYaZNWTlShdl55xKWdebfYZD03F9abBWi3uh1KGD242KR3kVzsTS3w5l5BlprEcjyyVtcwI1BWmzYI1K6jSGHFEV0V5ne2b1i27p8ktgXUsM6bm8pR0xq22qAwW0BK9x8lrda9qRzaPzvVmcT+UMnxws0lLXjKtbO+9mWfvX2Mn6JExttT3vB1LTjq6KirlJpTkBrBs/HN+1ucutZeXLeV/4/L+S1lgLSuzmkMr0/1V6zacwTLHhGf2P9mzZ+7it1xtmrW1r3EJ3MDNJi1BIlqZNNtpX/tSvjqu/zVla+wRGg8P+EUKLLEEEWmrLQ9cNu2t6YxuCZDy6F1bfZZWnZbZ2gPrpfrQjmyeqOtt286f7HLQfbp3zcf/Sj2Y3tKOt22pTGubPaiHTBa9a/w9bU5Dvh7U6/noIR/lAcX1NrqZPZFlbyN9PWsyyALybM3tZ6WstdpekvLxWLc1eaVZ1LKakTK9eoIxPDMQtxf2WPste+s0BtK+OZdHsndY7BHa6rNGV1qAjfQMUR6jwcZ/8MMfnXtk68+KQR45aQEKUsAKF7bqIVpfC8GkIZZanRT4QflIwRJeeaWgIUuQB+WhBeJY+Ettcn0tBZhwJIUh5/zSGEjBPnlbUnitJgttRyvj0ZUWYCM9Q5THaFxvHrishTyQ2hZ+GsTRAyWVyvaStwePHjwt7bSWsZ54Wluml7weaoENm0ha3w/xL8FlOZoUNNPSnifbZ4/Mq95stV7+vTKxevjWkIdvTdba3plde+n3xzK7rEZDZfn0HNkj3as5EshSpgeAyMK/V9SXh28NtbodAR3W2juzay/9/lhml01kyQJbc5Ahl7lTWmrlyziakdWztK1Z7lmy1XKZY+mBfsDCgyId48stVanuJdk0fVv40n7k17nsrNbMsYkPl223lLU276s0Bp6luZaRtyYjsSZDy8GeiUYTCCPBRFsgjxzfRBqAwQI7tdyTytSuBiigA1gOItEAHOl66egii2wWvhq/UuCLJXMsB08uZa3lwC9UXg/ISsvIW5ORWJOhx2pqdHDZRB6QA63jcblxUEdLxhcJUuqBrNYANzSwjgUWWmq7RodaPc3lJOUQtLipPNlhNNeppSylljx4GvWAH2vW+NGCaqZBJZCDBoyQAD5DydKrzo8TSUAWbmw8ZUvtDCV/TRmMBVRjyS4rkQduWMu3dGy0ljPOEnThkdcC3LDIJ7XbQ3ceqs3AWrO6oDBZbuWToLUWqLXUDgeGoeAXD5CqBqzjyUF3klmTHSidCFNzfE3as9y+76B6rYVv2itJ8tHcaPm1hx57dmkvKcnkkTfJovGzyCe120N3HvLIUCsb1VkKI33Djs0rdEl1leo+9NizrA65dtJzksub/v+Hdz647Df1ReubZcy1tks0VWu8x7LIWSg1i3qpbckazZXxyKBZXy2Wb8nqrJ2LluqUrPIUnMHxtZxxVmMJ1nTn8VLUtMXpJ52/5wHT0P5r8qZnm3oWWsE60jlz6fnDqtVrRMbS+n6IfwlU03uf2zv4YlL7ZUtwhEdejz4lvp5Ak9409F649XlrGZvethzal3R91alnPBJbQTUhhFUA7gVwKMb4iyGETQA+B+BcAA8DuCbG+KSFlwXO6CGPb3QoGGSvOh4fP63TIyNtK6aghYbi3+t5axmbXn2S+pJ+3/0HT31frCx9Beg/AO8G8FkAdyz+/SEANyz+/wYAHyzx0OCyPagGbqllXm2FsWptadet/K2wSktG01YIbI9suxb+Wjbckkwt/akt42nL0zeJ0AqXDSFsA/BvAPyXxZceAN4K4E2L//80gK8CuN7CbyiqgVvmZYbKEtMC2e0BdbVkzWkFbUiQz16QWolfK5TZ2p/aMp62PH2rIesy/sMAfgfAqdm1LTHGwwAQYzwcQjizxGT+hePLDGkWKCVXFoBqQMuNIiWjSm7c8yz3aFnNSCgZEjUDo8Q/l81qZJMMdD0PfpCWthw/aUw0o2AyrlFoqmXbwem3ZGzsvUTnnmsA5r61GEcTFV/2EMIvAng8xrgvhPAmbwMhhN0AdgPAprN3FAH/lq8cwAc8cHnPaFnaJoVNUhioRKW8ZKgzqHUAACAASURBVFxZCapaguxqOdc4eCjto5QxtgSt9ZD12CdOPsvsnCDBtM9D5XSrgUFrpAXuWPrWI/OsZWa/DMAvhRCuBHAKgNNCCH8O4LEQwtbFWX0rgMe5yjHGPQD2AMBPve6ieO3lO5e+tJyrqGSA4IxZlPKyNDiixsjCkWX21OSyXOfKaLH1La4rbeZpmU24GSmtSHZuWb/MXWld4ZRWUlyfKP8a16Onj5S4oJnk/qN986xS6Lh1c71hYY+eDHS/j+UGug+V6g/letNomtDGSVOLTNyY9OijB2I8tAs1/7vG9dhDlqHcf4lvF9cbQx8AcFsI4R0ADgC42lqRC0+UyAPy0GYRDZRSQzV8a8JAPTNMzR6Srkxye0c+81hBS5S/ZhuhgBM642v6qFlJSbYLzo7i3YdbZbGEBWtUso10cb31+Edn9tYvuAXkMKaZ3TObTAvA0kuWmrHRZr8xjWMvmYboE8aWqcby9SxBQDk+2r6+Zv+qrRh6zgDSjGPdU3pmCI8s2uxcWmXRPSbXnseC74ENS3poHUdKratGzVtA5QWwAibrTWIxlZfdYsWkVlguIYEl6N9jMaVtA7K1uMYSm9cvySidjKN5MCT+tbJQC7BmEZY8ABYvBfc3HWvLSUGlvrWOIyWLTBqVdKThI2rwC6PJQUep1e/bo+1Wi73Er6VO6+qlpaxWt8XTYKEWb0dvWXrz8+jOYrkXSVrfD/GPwmV7ZW2dZpbTEo9aOK7Gx1umNZMulbtV3yU9TAIC64Ew12Sc7Q0btmagxdj27IlaoY69IZk9+bbCcTU+3jIePVvkbtV3SQ+TgMACPDCr9zap1zPqASJJNNV4dovxjYOJSoALbkkjZTL1ACA4V1BN37RspBLRpRxntJH0aHENeYycniWkZiyTXK818nqMk1rfLGUkQItm8JPcjJrsHr6eZfxUs8tajG/akcIWiGfNrMQd3+udibm+adlIJdKyvwJ6dlnL199r5Ez8S3JrbXOZYWvl9cxwEmxYK5Nfo4Y0i8GPg3CXMtB6+Hpoojno9uzZM7d79+6mY4m1spb6NW3WZLjl6ljy1ZWynHqOguYynHryu7XkqdOyq1oy51r1oR1EqfXHk0k4laUHRVra1o72lnIeWo7Cltr5ymf/+Ih0ZPPUXW/egJMebdW0aW3b4j5KM5oWoCGtQCxBLVLZnD/ly1HL/libeSSXVU2+dI/7y7Jn18rSFVSNO1AKvCrpx2qXWbVuwxlS/6eSN17LtGnJ3+3JSlozk0mZQS11NHnT1z19ybkveM1KRjvqt5SdVKPaHPKabIB+7LIkQ8r6muokfeez7fEIdaxzfV+yY7N6tHTON5WV8tt7ssF6zgzI9SOtAuiz9eVP3fzwKGb2J4/NV4e45uTZ69TMZJ79fYu8pRNsrKsJj6XWs9drAZxIsgG8TcQSpprXSfrOZ1tAn61zfUuzKmdzyctyK5Fa4JBEkn4snpZ3//D5eYnvRGf2P9mzZ+7it1y9Ilc3t6+98sIteG7++IoTO/KvHP2iWk4SocTxlfbq2iyV5OVmK2kG4K7TvaQ2S3OzUL6X5GRp2YdbZJFsDNxe1ZKjnfKhY5OPUSnvv+X0H2pX0WZgacWm6cOif251J8lLn28tb/xEX/YPfvijc49s/Vmcdfpa9iB5YGFptve+Q3hu/jju2X90qWx+MP1bXvNyXHrewsOTrkn8uGs5cXzP3rhu2W+qm5e99LzNrLxHnpnH3vsOLSuTZJDkza9/6/DTy8pwbVLZjzwzv0xXSSZOFo1fiSyy0HtcH5Oev3T/o9h736EluTVKfOjY5GNUGuv8viRv0t3OM0/FdVdcoH4QE4+zTl+7oqzE36J/2ldNXtpn7WWf6DJ+47o1uPbynSbIoyd5RQtcsRVq6vFBt0BgNX6eFFMtOquBy2r3e8NXPdQD3jtJ+HAXXUnQuiH+aXDZFkjiUDTG7KGe9obKwMr1xwPZndT4tcgwSXk9kOWSDPhRgMsCunFFqtcDHjt0O5a8ej37UZPjzSODlk9N4lM71j2p1fUmle0pV8lw3SLDVF1vuXGBHrSnkWR049xgNccwJ/KAJjR3nWT08QBCPCQdcMgZelqOufYYxyxur0Qe8ItUjzNucoAkyfjKldVcbJIMlvt0DDSXrPTMWw52nOjL/v6bPjL3jVMuYg1qyTDiMdZQZSYjxoGjx5YMU9TgRUkzmFhkom3mbVMjHjX6ePrsIWoUyw1ItL8lo5bG32Mcy/uaDHQWo2lp/KR6nHEzl5HKaylLDawWGSz3qSzcc1F65hO/0bzs6chm7gtZAlrUuixKX3BtdvWAJpIbiYPCSjJwEErPS1fSCSd/aQXVG07Luack92JaFeVlpVmVW8VJxzBzKwitzdLzpc2upVWABcJsgctKbr/RWOPTkc0cCKEEtPAGSeTgB42PBnn0gCYS6EOCwmoy1OYEt2adyfkDehBObzgtV0YKKJGASVoQipZHXwOw1ISeSvxoH7UxtUCYLXDZmj38VANhOPCLJaijZh9Hv+Dafpbuuy17a23fRWGcFgilZyVT2m97bCMePVtmTm7GlHSlAVksdTyrgFIQEUfW2TWX17Jq02wsnn09oPvZp+p6q6Wa7LSWHOU98plzMg6Vi7yGeudmn1Y7NVmHx35kc4+yUFxvU5nZOfJYMUv7Iq4shWZ6VheWmVfbd23ftI79CmtBNENb6i0zGZXBs5+3zGyS7YLTh2ePLc2UlllcC6aSApm0gB7POHoChKSxGV2IK0fe1EravkgqW9rPSXt+yz5L23cB+pls3L6xNXOpRJ7gFipDTaIIT0ZaLVGJZ4+ttV2TdKMUyATIAT2ecfQECEljo4W4juZl7wm3nAQ0U4KzWuC99O/esNZeZOmTl4eFv6YPT1bVXvKWnrMaiLTWtqWPku5GcyLMha99vRvGOVQmVi1T6tBwUw8NnV1WKztpWGtveHILn2lAens8d1D27CcVPzcdKcWz377vIIAXl0fp75zSvT+880GxjFRHK0v5cvwpH4ucFvlqyNMnqYxHxrzs0H3T2p42n0n3nWuz93M31ai3/Fc7PqjHso07yohGqUnZTbUDB6U2a45k4spaloE1WyAPr5atBHcUND22Kh8La3uW7LUtxzHXHjfVcuSz9NxZjj5LPLQjm6cCqkmkHTVUk7HFAqJI/HO+6VcyJmmGEyvQQqPWzDKlMi3ZaXplqmk1ukl8uTqeo7Mkqj1uqiWwiNaxAL4oj9HkoNNcbzUBGRJZgg08MpQAP97+aEEXNX2vyYoyVMYaShpgxpNt1wJkSSTpU3OD1YyJ5r6l0GguUCqBliTXowViTOUdTQ46jXrkO0tkybvtkWHoVUZr3z05+Sz3atvjSMoRz93Tsu2ma17YNNceN1PWjInlIFEtryGwPHRWcz2WctandkeTg06b2RO1ZI61fPU95An08MibcsZZwBMWGSSAhRZoI81AnqAWS957zwxsCe30rJgoaXU1fUgwaku7dFWoBWlpK0irHkYTCGOhlsyxvZNB1AZ6lORNmUuBMnjCIkMJ/MLNmDUnq1j41oBfpH5wZbx2Ga097R6XvKLmbLrSGQRS32qOJy/R6F52aoHUzv/qAfrgyHLmmCSvdp16GDhPg+SVoFZXzSOQt5POySv1g9ZJf2vWck7f0jluGmnn43moNP4e67nUt6SH3CujeZVo25azBq2ya9c5Gt3L7oFo9vjacdTDamvNCe/dO1rsERaoKiUpHFab0Th9a3t0ieiKpxYiXBp/j/VcO3lH8xyVINy1p7nW2GUomV72EMLpAD4B4DUAIoD/AOABAJ8DcC6AhwFcE2N80iy9QEPN1i0y1JStlXton3lNnUlBVD3t1JCGIWjVUeleax97PGdWBN0tAL4cY3wVgNcC+DaAGwDcHWPcCeDuxb9N9L2jx3DzXQ/iG/uP4Oa7HsT3jh4Ty6avpmfJU9t2ugeg2GYqe/ip5wAAh596rtiXVMfTX1pH0wftG9cPWob+cnW2bli7bOYaerx66ZKWTWOVXwPKY01J66N0L13fumGtuR0PJZ01gWpCCKcB+JcA/j0AxBjnAcyHEN4K4E2LxT4N4KsArrcIZgFa5Nd6Uq+2NZeKZCzsCbjw9q0kt0UfQ41NSZdaHYssnF6A6We47bmMT31rjXp7JYAnAPxpCOG1APYBeBeALTHGwwAQYzwcQjjTKji3pOHgrNLB96mTFmOHBMnkjC0emKVkbEv8OGNTDbRWMlhq0NqdW9Yv/VJ4LzWgJbm5OqUtlWbw8hiiqC6TLJedv3lFG9pzAkA1tmnRdDVUA4n2GCFz/qVlfNLZgR8ce1riZ3nZTwZwEYDfjDHeE0K4BY4lewhhN4DdALB9+3YAvKGKGjY0EAGgG45y0mY7amwB9LjzvK5kbEv8OGNTDbS2xmBJD0HMy3IGtCR36TBBi7GwxhBFdZnnHgCgri44gyJt05KbsIZqVhceIyTlr7kTk85Oesm60yR+lpf9IICDMcZ7kgxYeNkfCyFsXZzVtwJ4nKscY9wDYA8A7Nq1K0qNcG6jfEbjvmz5DADoM7BUNv/Spv2U5FqxzCKc64lbXdBf2jfKT9NDIinYx2KQshiQLIErkuuNW5lJM6LFkFZy/yWXI3Vx1tgSPKsLyp+Wpc+hRlzfpNVE0nvTzB5jfDSE8L0Qwk/EGB8AcDmAby3+ezuADyz+fqEovUIlgAXnCqHZZ7QZWCrLHc0ruVYsswg3c1IYpGdvWtJDTlqwj6Rn7nop8EOTW3K9cSszqd9aME7J3ZqPk0fPGpVsF95sPNaZ3ZIhKVHSe+vMDgC/CeAzIYQ1AL4L4FosWPJvCyG8A8ABAFeXmMy/cLz4xfXsbTyznjQTeGa/vL20CrC4Qujq4soLtwAACzzxhH1qYcG5nrUZxzPbcft9DoDDyUDHsxQ+WkP02aF2CU3Pmh6k50yzz0jPcasrWbKbJP5aphrTyx5j/FsAu5hbl3sETckrLBZgy97GM+tpM4J19svbK+2huGt0X8wBTzx7XgmAowVUSHUtpAWuAHpQRz6eqf+98+vRZ4faJbiyLasrzT4jPcetwK+S3WQ0gTDpRBguKKCU4ZULTJCCLDxBDBxJWV89oaicDCkA5soLt+C5+eNsIAw9YYaW4TLnUnkoD+5kEU+wUCm7Knd+mxYAknLXWwKBOH16TgSyluWCemhZSc95e3SstVz4UiCXFg4rBX0l+UeTXTZPXiHBRaUZ17Pnq7Xc0/q0juerzMnABcBIyTAk63iNhVY7WcTbF8muolnu6Vjfvu+gORCoJAPX77wtiTgPgBWWrbVHZ3sNNizZAizhsJLd4Eciu2yJtL2wpWwN5LMFvqnJoAXClGTw7PksVu3avtTyyssPDcOtkamVv6dvJTuSxZ5Ef0eTXTadCNNy2LxGnmycvbOHatlwpbY9mXM9Mnj49cheO0m+tG5vHXpkG0MGWkoYS3bZRGkJcvs+OaNrL77SvR7tcfw92Wo9mXM9Mnj4WeqMiS+t21uHHtl6P0O95JJodAc7asa2lrxv9EigZEDRjuXpnWsttWU5yNAjA82AQzPWaMZJy5FD0lFOlgwtFiNZzZh7jHGeXHycMVIy8tYcC2aRS6tTGq/RZaqx5O5K5DEoWVxvmrukpW1LrrXUlifm2yKD5ObxgGAA2VgmxcdboLsatYy5xxin8bME+0iuTAvoxdO2RV7PEVGURpODTvqaWb6eJR459XC5cHJTtxTnIpNWFZybzuPuS1/75NKjh1dyriHKn3PTSbrjXG+lbK2a7ixjbqkrua4sefDS2HBla1x6Gj86bnSW1vpPx5rWHd3MzpH0NfNABmtmhBqXS40s2qpCcytaZKBfew7AIQUcaW46TXd5u5ZsrR7d1epdc11ZV3xc2RqXniVLkDRLa/3n6lhn+KnO7PnXOQEtLMcDl0AO3GxtmUXorGch+nXWZhNpVZGXpYATCxiox2rFcrSwJXOsBbTjAfZ4MghTPUjAJM84cqstSTYrv5oVQyINrDTqmZ3OaLW55+h1S/5xjpcnZ1si+qXVZhPLqoICTrxZWmtXK54TUFphuDWhoams9bnQwnY50sYR4EFcnr01pzPPikHqo6fuVF/2UngivVfDVwpIKCWk4AJ2uKAJCoTgQhiloAsuWQYNmrFkaS1lpNXqSIEqHEnhq5LuJaJ6rglCsfYNsIWVegAtWl+lMGOLDFIorSdYSaOpvuxaeGL6u5WvZZ/MQVWtIala4gUKUaX1tQy0iU+Nxd6z5+UCVUqBR5pMHiu8Z1/rCSSpCSstzZgluKzUduJjkcGziqmh0RjoEvWGRdbATqWvfA0EslWGoaCklpmsRzsWGXq3WYKUDkkt0OKh5Z4ogi7Fs9OMpnmG0FJ21ZwsWWoTvzfs2LyML0dcFlLgxeyqGo+WbKVaNtV0z9JHLQ+fpOfUJy3rqaQXC9Vkga19HmgdAMUxkfhZ2vH0jcpkGUdtXKUswaM5stkSz07JAkrodeRvjeumRU5N/hZZPLK1yN/ads2ytcW45ylTW7dUr9fSXHo+RnNksyWefYWACjzQ4rpocd1oMfWSq4lze0muMS7mmcbqU1k0t5pELdBVTX4N0GIBiNRAmKVDLGtBUZLL0QLmssC+LXUoWWDOFAyVQDbfvOPWR9733t/9KMtYipAZ4l+KevPQTXc+EF9x/R3xmo//VXzF9XfEm+58oKp+qkf/buFVy1/rU6l+i/y9SNNDzTh56vYcT2/9mrZa6nj0kcquOvWMR6Lw/o0GLkuJBndYgkZoXS6biyfbDCXLl9wyI9BAmHympOAi2icNyGIJvqghaXbK5U569gBEEnkCg+hKjwtYsYBfSsAero4FeERJA3qV9EHhz5xcdNX15U/d/PAoMtV4SMsnZq0LrMzm0pIDzBLAY3HHlHK5WXPXS/d6u2wk1x51nXkBIok8gUFaIEjpUMX8es3e2gI8ouRxM0p91DINUZehloNutC+7BcBBiQPQ9AYmlNrm2illJ+UARLS/HABHyl/eAkjiiLqAaLvedqiuuL5JdShQRjvBRsvvXnJraYAnWsfyjGntSRmFOXcolcvzfI/2Zfd87RNxoZ7AZM70skAnpeyktI4VzsqBR4boq+UkGA9JKx0P4MRygk1JP5rsGuCpJBtHlnz/0qlFmlyesZjKnt1jbdQs7LQ+t/eryWRKybIHtngNpMAdixeB2/vR/b0lyYSUpZQbC6nfltBWyk8LeqrxFlg8MVq4aSnoyWPbaQ2m6hlCO7pAGIvftuZLyK0GaGAJ15ZVXq2utpe0ZCcttSmdSmNZMWghtBZfP+23JbRVO9VWklsjT95/WoZb8ZRWQR7bTmsw1VAhtJSm8rKnfUc6GcV7MogHfkoDS6z7QutZbCWZvKTtw2nACj2VxnIKyeGnnlu256O2kXxfK+0PNX1YxobKQPvey77SsmevlaUG4io9d3nbpWCndH80CLpEmrXRQtKX0HIaiyeEs+Uk0lqS9plcwIq0R9VOIUkrHSlFlrZyaEkXpckg9b2VWvbstbLUeHssz10p2CndH13e+F6zYE2bntl5kvJRGeivRSaLXkvBFtrqqFdgRkuwSEs7Hv7TCJ6xhNdK8qTro80bn+f8nlRed8pPy2FvyUlukdvTl5750Celw151PXn/tfpD5eOnfC3PcY92PAQlb/xoDHSAfrxNrzYpP0tAjCaDBdBSE7zResywt92h+fUKxrHU76E7C1/Lc9yjnV40UdfbR//o43P/fN6bVwR5XHXxtiZYq8c1RgMoarOJSvU1SG1yl3AH95XaroHAtkCDLX211PEEgGjQUkmHrTndSnJrfLnrEtzZQpr81vHXXG8TXcZv2/mTgwRxeIJPphFAQeukoAVP0Mg0A2AmFQBiqdsaGNXS9iTqt/KFsoyfSoirZ7by5ILX+EqgDE+AghZgYyEp+6kl77rl5BaLzko597UZ0xIsQvtaM7tqwURpdWHJc+8hT1ZczwqyRgZNn6XnVZvZJ5qpJh3ZnHyHaY9y+76287RoRg+ujpQF5OvfOVLkz8lSyg6jyZmywzz02LPLYL1av5OL7M/+z4EmnUn30nXu7DSpr1o7NfrR6qa2kstR010NUf6c3JZnMQcc1cqg6dPzvFKaKlyWgy96kh8kqtkfJvJAHT3w1lo5S7YAy37OkjCCroK0RBoWfUjj1xpuK+mDm0GlUF8LFLgGsmtZQXr6mMZAW72VVg6jhcsCK+GLlvO0SvyGgjp64K21ckp8NNCOp/8SGIMLJfboQxq/3kEziSyBQR4ocA1kl5OtBVSTxgCQod01YbaJTC97COE6AL8GIAL4OwDXAlgH4HMAzgXwMIBrYoxPWvhxAIGakFaJXw4tBLAMZijlcM9l8MJ387a1XOiJfx6WScvUgDk8daTwUk7vJX3kvCRAjzd8VcvLn+vMEoKqgYws41UiTe8euG2qnyDdCUbOPSctz2jxZQ8hnA3gtwC8Osb4XAjhNgC/AuDVAO6OMX4ghHADgBsAXG9plPv61YS0SvxqAz9avpqWAAXPiak1bVtICi+1nD5b4sWtPFpPmpH0aglB1VZDLQElHH+pTxa+FNINyOe4DT6zL5ZbG0J4Hgsz+iMA3gPgTYv3Pw3gqzC+7Bz1hCdaIJ81kMTatum9acBwqQwWCKkVommBIHtgyjUw31oaChbbskJrgc2qJPnk8n8A3gXgGQBPAPjM4rV/JGWeLPGZFFy2FSZagkd6+Q4FBZ4UNLNF/mn0dVJ902Tr0e8aHmiBy4YQNgJ4K4AdAP4RwP8IIfyq9WMSQtgNYDcAbN++HcDwcNleRqFexiYLpLaGeuioFyR4CBlr+Uyqb5psPfrdG+5sWcb/PID9McYnACCE8HkAPwPgsRDC1hjj4RDCVgCPc5VjjHsA7AGAXbt2RcC3lNaWK9JBeFy8dYvBJBlKcr6lQ/jyvy056KhhqtRXq45Kfaw59snTrsUQ1prDTWtzyL5psg29JZSeNy2e3bKEvwTAN7GwVw9Y2J//JoDfB3DDYpkbAHyoxKsmb7xGHghs7/zdpbbzv6W2W/Kl/6hR77E50Ul63rS88cWZPcZ4TwjhdgD3AXgBwN9gYaZ+KYDbQgjvAHAAwNXGD1YVce40S1bVUnYVbVbxrECk2Zuro/H3ZMvRXFal+9I9yyxr4Uvrc/rQMsmU+HHXa7LLePh7yFO/JVMNffa1eHaTNT7G+D4A7yOXfwDgckv9HqS50zT3SWnfo923uG4SlY4WLoFSEn/puGitb1IZC6ikxi7hAatI+uLKcuMo8eu1T/bw91CLjcEDAqI6+5HMG08pnxloDjM6a3BgD4v7qPQ11mYTKovGy7KaoAALLk+4tPLgZkwNVJJmk1LOeSl3uwZssa6cbr7rwRW58zh+2nWprCenv8arxu6j6UEaJzomdKzz1Q+93rRn7/mv1559qPPQaviW9uO1e1RpP18TzmuhocKEW/rai3rZC3qH7Q7xvI3mrLeUvKL1DLLSSZ614ZWlIANLYgpPMIoWrigFvlhCKGsCdix1LIk6LPnuS33lypZCc7mc+BZ909Dh2rMMSqQ9OxI/bqylOqmsdorrRF/299/0kblvnHIRTjtlNS49rz4D6233HsTe+w7hrNPX4tLzNuOTX9uPW+5+aInvhrUv/rbwpcTxTdduu/egWQapDlfm7I3rlvEpySjJmYjqylPnrNPX4rorLlgqw9WRxsLTV66sJHd+/VuHn1b7xvFN+jzyzDz23ncIB44ew977Dqly1kxU2rMj8ePGWqqTyv7T33zp+d97z+/8Icdvonv2jevW4NrLdzb5HoFxQid7+WeHqDPJNj17bA+Pmr27h+80MwpLVKOz0WWXraEcOtgL0jhkXU7eBN+Ufi1ZVXtnjC21Z73Xs22LPnpnkK2Ra5qwbOm5wKrV98daP/tYSHO9eetPAn7LyUtdKx73WrrXG0Kpye11Y/Zqu8XdOAT1hk/34Cs9F9ohEVPJVJPIYlyhWTx++afPwSU7NhczhtD6ybikZSeVyJIVhVJuSEn56pIBSvrlDF1ahpbjEe4sMFI2lxrDoqZDT4YaLVtOTQ661uw4lErZZWsMdhzfpFf6rFqMhenvL3/q5ofnfu+9H+Pam+rMbpmtaRYP7dgjWpe7J2VU0ciSFYWSBKpJ9aXfUp9ojrPWFU5NlhuLDmtAJVy2HImfdpxX7xWI57ixFr4UIKOBxaTY/dGCanpBSjkwhtSGJyOOBrX1ZAypAevUGKRK8Nlcbi64x0oSmIdrhwMDlUAlHD9pHLlsLhIAZSiqgcZagFkWYI+n7am+7L0gpfQgQ60NT0YcbYbwZAypgexKs4amM8t+1pItp0SWQyypfqzw1tLxy5Z+9MhC46HWkN8WiLGn7anu2XPy7LPoXsdyAkhNXvCW/WvOK50SIgE3ep1+k/bzV164Bc/NH2f3syUAC8e/Zg9M7RycrcGSTVUibh9N7T69TsSxyuJ9/pK80jNpyV5L9Ty67LIceb5Q3MzSOrtyZVr2r9Qekc9CNRl0Lf2wHIVd2if2CjChMztna7BkU5VIW/ElPq17aivVZEAGytl7LSs+zwpzqnnjPWd5cZZ7T47yHjm/LeeAUYvqVRe/eP5XmnGptblmVcPJmL7yib9mza05q84zQ9IZx+Kd0PSh5YC36GZIsqyGOMhu0oM0TjUQ5tHN7NoXyzpTev2TPXJ+l+pwM2Wql89g1F5Qu6qhRO0RNaG/HP+aGdIa8sutzLTVhWaPmNRMTsmyGpLsCNo49cpem2gqL7sGA5QsvVwdS5kaGVrqWCCZLVBSj3wWPQ8FC+0NPbZ4AKZFveC9Q0GYl0iC1g3xrxdcdmiS2vJAYFvhpj2zqU4yK+404byeOiU+rX216G4IXaElu+xYaGiopqUtDwSWk7PVRVMqW8pYY+GlyVBj5OxNvbLRtGQwspSxGS0xMQAABbxJREFU6G6SzzTwIwSX5VwWyVhDfy2HK2plpZhxi5EpwThzdxJ1sdQYzlpcbrWGKwmSyRlRPX0sUauhlRI3ntIR2Bw8W4Ila7kFqK4491qLQVF6jkdnoEvkgcvm16U8bZ484RbgCXVnWIxM3BE+tG81hrMWl1ut4UqCZHJGVE8fS9RqaKXEjacEyuHg2bRvGl8qpwfy6qGagKCpzuz5l40CT9KXSsvWYQkooWRxT1lmD+nLyoE90ipAmiE1F6Ql44sUqKERlV871piW4YJRaHCSJaOORJ7AIy1bDpU3LyOthhK/vG/bN61j+6L10QLiaiFpzEc7s+dfNuqeAha+VKUD/EoBJZQs7inL7CF9WSUoqTZDcm1LAJGSvNbD/jwZTaUyeTBKciNRPdccQOgJPNIgpR7YMHUR5n2jZSx9HHo/XjPmo9mzp9mP7nk98FZLqGu+F6Ngj1IOM25m4MIxraTNSlK/tTrSvtUCBuJCRkuwZE532mxqkSvXrSX3GsejBjas6VladXKrrZLtojb81vqsj3Zmp1+//IsKlCGUHmsxtxejYA8PuKHliGmNrzd01lLHCgbyrKDo/ZrQYUsYr3Qcs8ajBjZMybLq5FZbgG67mKRnhNJUs8ta9uPajC7N0loWUW0mluprswfNcOtJiuHhq/HT+k15UsuylNyD85BYrOSWAJtSAIhnv6/Bey3PUE3fqC0n74dkn9HgyJbno8Q3kTazTzW7rJY9k2ZVzSllFD3yzDzu2X90qX6qq2UR3XnmqbjuigtYvlJ9LrOplOFWy04q9cPCV+On9ZvypBlYaebYnIdUtiZDLJcFlmar5eQt6TC1m/PwPEM1fUv8uKy7VJ7S3xrlMr3lNS838RnNMp5ml62FVNbAT1ugtBZ+NXBOD98WCHBNO5ayHlk0vtZ+9KQWPUs8elN3/hK0boh/F7729UVo6SRhsbQ9CeI4jUymFhkssF6t3FjIAju1wJFLdWvrDA2trSGJH8YCl33y2LzZ3ZNfG5I0YM80M5laZLACcCatUy9ZDKse4JRUt7YOoOf6mwaMePQGuk/d+om5X3/nO1lXTsloMxRxEFgKceSMTiV4bw1MlMojgXQkY57EozUjLe0bB8TxGPNKfdZcspo7TYoh146VkmTRQFGWjEY1ZSxGw1K25K989o+PjCK77JqTTxJdOZPKF0aJg8BystBYdcAO722RR3MRWWG9rRlpE2mzLOXbEqNvcclagCyWXAOaLBIoypLRqKZMTfAMdXmOJm98cr1x7imP22toktwa3CrA4jaqAb1IpLljJH4WuT2597mVWck1VBv0RJ8PLV98Kbil16qx5OrMy9D+588+1ZmW3Ud6Tyh0eXQHO3LuKY/ba2iS3Br5dakM5zYq9cHTR80dI/GzyK3JQF1MyQWVu6JKLibO9Ubb4lxZVPdfuv9R0S1HD2mU3Iutk0bJ1an1P3/26fjRZ0dzK6a2U5+SXrSDHae6Z9eALRp4oBTE4cmQOtTpIZZQxpoZxwLv9YBSLDJ4oKqWgBUp227NSSg5lYJbPBBjj844W44lQMrTN0p0FZNkGc2JMNyenZK2Z01UstB6EgUMeXqIdq22bQu8t2dOe6/8loAVCj/lAmxqwkJLwS3efksk7ZsTeQKkrH2jJIXoaifChAXX3GQohPAEgGcByMfKWmjV6jWr1m044/gPjj190kvWnZZ+f3jsqe/jh8/Pp/tLf2d1ll3Trr9IZzTLW+iH0nZN/TOwavXTZr4tMnj0zNQ7/oNjT5+0Zu3Lj88/92hx/BzyrHge+vZ7pX5zXgDMMrSQ3NdXxBhfxlWZ6MsOACGEe2OMuybaaAPN5B2WZvJOjk6atgAzmtGMJkOzl31GMzpBaBovOxuRM2KayTsszeSdEE18zz6jGc1oOjRbxs9oRicIzV72Gc3oBKHZyz6jGZ0gNHvZZzSjE4RmL/uMZnSC0P8H6RjKkx6drbIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO19f7BfxXXfZxESSIB+8iMSgiCDSOoYbIMTm5Dp2CHUKc3EmQ44bccex3UGZ6ZJXdyODc0YxD8d48YQJm7iyoBNE6c1BU+dwR0HD63/MElIJeIhsR1+1GKEhAy2hCGyFD+wtn+8t8/7zjvn7Dm7e7/fq+p7Zt689+7dPXv27L13d8/5nLMhxogZzWhG///TSdMWYEYzmtFkaPayz2hGJwjNXvYZzegEodnLPqMZnSA0e9lnNKMThGYv+4xmdILQxF72EMIvhhCeCCE8HUK4cVLtWimEcF4I4X+HEL4ZQvh6COEDC9c3hhC+HEJ4auH3hmnLmlMIYUUI4a9CCA8u/D9aeUMI60MI94cQ/nZBz1eMXN4bFp6Fvwkh/NcQwqljlrdEE3nZQwgrAPwnAP8YwGsB/PMQwmsn0baDXgXwb2OM/wDAWwD8qwUZbwTwcIxxO4CHF/4fE30AwDez/8cs750AvhRj/EkAr8e83KOUN4RwLoB/DeBNMcbXAVgB4J9hpPKaKMY4+A+AKwD8afb/TQBumkTbDTJ/AcDVAJ4AsHnh2mYAT0xbtkzGrZh/4H4ewIML10YpL4C1APYACOT6WOU9F8CzADYCOBnAgwD+0VjltfxMahmfFJdo38K1UVII4QIAbwTwKIBzYowHAGDh99nTk2wZ/S6ADwE4ll0bq7yvAfAdAJ9e2HbcFUI4DSOVN8a4H8DvANgL4ACAl2KMD2Gk8lpoUi97YK6NEqcbQjgdwAMA/k2M8eVpyyNRCOGXALwQY9w9bVmMdDKAywD8QYzxjQC+jxEvgRf24u8AsA3AFgCnhRDeNV2p2mhSL/s+AOdl/28F8NyE2jZTCGEl5l/0z8YYP79w+fkQwuaF+5sBvDAt+QhdCeCXQwjPAPhvAH4+hPBHGK+8+wDsizE+uvD//Zh/+ccq7y8A2BNj/E6M8RUAnwfwsxivvEWa1Mv+fwBsDyFsCyGswryh408m1LaJQggBwN0AvhljvD279ScA3rPw93swv5efOsUYb4oxbo0xXoB5ff6vGOO7MF55vw3g2RDCTyxcugrANzBSeTG/fH9LCGHNwrNxFeYNimOVt0wTNHhcA+BJAP8XwG9P21jByPdzmN9aPA7gaws/1wDYhHkj2FMLvzdOW1ZG9rfiRwa60coL4A0Adi3o+H8A2DByeW8F8LcA/gbAHwI4Zczyln7CQqdmNKMZ/X9OMwTdjGZ0gtDsZZ/RjE4Qmr3sM5rRCUKzl31GMzpBaPayz2hGJwg1vew1kWwhhOtb2pw0zeQdlmbyTo6qX/aGSLbjTVkzeYelmbwTopaZ/WcAPB1j/FaMcQ7zkM139BFrRjOaUW+qBtWEEK4F8Isxxl9f+P/dAN4cY/xNqc6ZZ54ZTzvtNKzbsAkvHpnDhjWrsOrkpd+buVePifdqKPE7/ZSTcfgHry7+tvCfe/UY9n/3ezj3zPVdZJFk6tlXTl5NpzX3Wsco7/8LL76sygtgSVstbXN1LfzyMi+9eBBnnXXWsn54nqkaGay0e/fu78YYz2JvNkAJrwNwV/b/uwH8HlPuesxDJHedf/75McYYb3/oifjjH34w3v7QE5GSdq+GEr93fvLPlvy28O8tiyRT775SfrX6ruHnkVPqf86fttXSNlfXwq+kh9ZnqudzBmBXFN7Zkxs+IqZIthjjTgA7AeDSN1wW7/jyk7jyok0AgGsv37qM6ZUXbcJffOvgYplnDx3B/bv34drLt+K8jWvYa+n/Ky/ahEeePrikbOJ3zSXnAADefcX5eMtrNrFtU76pDFdWIou8id/2c05b/H3Hl59cvE7rW9tKfU38kj62n3Ma3rxt46JOc6J9zHlK97Txs1CSUxoLTu9XXrQJpWeHPgf0eeD4WsZYKkPHkdOv1Pe8rEUG7Rm3UsvLvhjJBmA/5iOv/oVW4cUjc7jz4acAADdcfTFb5pGnD+LRPYfwyNMH8TPbNuH+3fuW1aHX0v9/8a35unnZxA8AHt1zCG95zSaxbcr3vI1rxLISWeRNfO/48pNLZEtU0pHUFu1r0sebt21cotOcaB8pT+1eLSU5pbGgMiVdldqmzwF9HrjxtIyxVIaOI6dfqe95WYsM2jNupeqXPcb4agjhNwH8Kebzc90TY/y6VmfDmlV471XbXV/R9Dt92bkZh84UXNn8i5hTPpPRmbHm6+mZPbQZLu+DJAPVTZot6MyW9700Q3Cy1szo3AqnpI+aOlyZNNOm1Zw242pt0jLSioGbrUsy5mNhec5Kz7GJpPX9ED+XX3559V7Es7f07IG4/WHvvbSlbc+9HmWHtl1Mqo7Gx9JHz55dsvvUyj2ETQgD7dknStqXXVoNWL5+0v6w+uvpIE+favlIZT19rLFdTKqOxsfSR8+KgZvZW+Tu1V8zSV+BIX4sM/veg9+Ptz/0RNx78Pvs/5Y6Fv6Pfuu7y+pIbXNlPXytMmmyeKimbkt7Gh9OL6W2tPue56OlrKePHj75/V7PbU44nmZ2yfiW/rfUsfDnDB0ew5+Hr7fP3j5Z+A1Rx8KH0wugGyE1WTzPR0tZTx89fPL7gN8YOxUDXQ3NvXpsiUuIc0/lRiaPqyUvYzH0cO4SamjRloOS28tiFKLEGW2kPmmuF02Hkr4tbkbJ1amNH2dEpTqUjJCakTC5EXP3olRWkqV126S5hzVXZsn9R4njm28lvEbkib7syfVmmVU9rhZaRvu6au4S6hbJXSLUpVJye1ncMFQmAGK/LV92SYcWfWvuH8sMKelc0qE2xpwsiX9yIwIw60Hja7lHqeQe1p5Fzq0okcT3Z7ZtMr0flFbs2LHDVLAHfebuu3b8xvvfj1/96fOwZf1qXHv5VqxbvRIAsHXDaqw9dSWuvGgT7ts1/1XOyzx76Aju+eoebN2werHOipOAvYeO4G0/eRa++Pi3F++l67/60+fh3A0/mnny+rS9rRtW47Vb1mLtqSvVNhNJ8r7tJ8/CwcNzS9pOJPHLrycZcrlyeTndJaL9Tv+/+4rzsf3sM1R9p3ZyWV4++gru+eqexb5JY8O1reku3aN8ubJU30m/11xyDo7OHVs29nlZOo4rTsKydjQ5JaLPXa6PpDPtOZP0QeWj/eB0R5+DW2+99cCOHTt2cnJPdGZfdfJJ4kyZz7jW2ZrOpukeB1zQgC01IBKpvgawsK5EOL7aKiMR7bcGXpHk1/aU2mxiAUNJ/bXMUhoQic7wElDIssKxkPTcActXZqXnWpNPW+F45E000Zl9586dO66//nr2a1r6YuUzx7GIJWXTbJe+ntvPOW3Z7Mp9Jel16avMzaK0vjbr0T7Q2YhrxzMjSHy4FU5J/lx3529cs6S+pg96z1KWtsnN0pK8acy5lQ7VS9JDWg3k+pDGxLKq41Zx0rhp42lZtUl6pqTN7FN52e/56h7c+fBTWHvqSlxx4fwsla5tWb8aN1x98bKO3LdrHx54bD+2rF+Nbxx4eUnZczeswRUXzivxzoefwsHDc3h0zyFsWb96kf+61fNtUb75dSqXVIern9qW5M/7cPDwHB54bL/aDuWb6ypRSd5cZ7Su1E6uu2cPHV1SX9MHvWcpS9uketHkTWOefuftUL0kPRydO7bsuZDGhHtGqQxffPzby/QrjZs2nrRPludNKjO6l92zX077mHzG3HDaSnZ/yO2b0yrAsiezzKaJSrMqV5auRLQvOJXJskenpN2n+k2za76/l+wHFqL60WwBnlla4q/pTJs5qV2DW6EleemeWloVcDKUVq5WKtkYRrNnT+TZL0v7LG1/mO+bAbsv07K/ovKlMh5rNiDvuyWZOOJsE9b7VL/Jyk3391Z9UNIs+MDyMUny1do5cqI60+wd1K7B9VXyamh791KAUS218JkqqMbi9+T8syV4rMeX6ZXLU6amrIdKfD361aClPaCvtWMitd0bWmuRU4LJ9oLl1sjrIglaN8RPgsv2hremexYIrIdqILAe6Gsv/tOiWllK+hiqbx6o6pj0rBGVE2ODy/aGt0rxy962PG17+lZairbynxbVylKzNO9BHqjqmPSskUfOqbzsFmhmvrwEli6VKBwyQV+5uPB8qcXBCy3Q2tYIMQlCyfWtRDlU0xMPnagUm63VoWU8cdw5SePP8bNAda3yWrYSrdl4SvqtGTONPMv6qRroONJmu2RckYAsHvAIbY9ep3K2GNQsxhor/9zoBtiNj7RtbTUk1aFlSgZCiaTxt4ChWgKjLFDV1hm9pN/eKwYPzHcqLzsX1AGAna25YAMu1xqgB8J4DD29v76UtBmxlF8vzye3ed3qRX7WoIjUT201RGWhwSepHYvuPDMxlyuPrsy0fHpSHU/ettrVCtc2zRZE5ZdWAR55PTTVPbsW9khhp3S/pQVDcF/PUh4xTj7KtxdZXGKpbclF9sjTB0VorUYUdmrJyScFn1h055mJuVx5dGWm5dPj+ui1idSuVmjbwFJ3Iic/AHYV0BLGqtFUXnYt9xrn1qCz94GXjrqzk3pI2mPnX9qava9FPqnfmouspr+WGYzOUpYsqtz4eduxugo9OfqGzMrDkcX2lIjqtyZMOrWHFStXiYUkM/0QP8n11ppfbah87jUy1OSjr22zJx1veeUmxbcXteT6q3mWUt0VZ5z5XBTev9HAZSk8ksI4c0iiJxTVE8IoldUCVRI0lwsjtRLXrge6a+1HThrUU+KnhR3TIKJElmAiD6U2uWCnUh0uxNXTpgTZ5fh6glsS1dShY/Olz9zxzI6bP/IHXNlRwWUBGcYJLN0vWi3qPdINaRZ2y963RJqNoQdUlSMN6mnlZ/Fb1+Te10jb35fq1O6Fa3ASrZ4cax0q2wd/+MqcVHZ0cFnNmlkD5xwa1tpjj9dLfk+dHlBPi9+6N7Xswz11uPpS/2v59iDX8yGt74f4ueT1b+ySgbUG6liTKdYD5/RkwR0K3uvtz1D1Wuta+LVk8+0lZy8ZPG2VMuViLHBZLQddiVqhjp6lXA2c0+NiGgreK1FvWOvQdS38WpfmPeTsJYOnrdRGjdwTfdlPP+VkXLptY9WxPNxyxZKdtAaaWgPn1IA9AFgZLG40D8BHA+Rw/SnxoG5AD8hD60fSRws/bgldc/ihZRlc0mvtMl6SVwOFca5oK030ZT/8g1fNGVhL8eKa8aqUrbVk/GiBc3LyA/b47ZY46B5ff1qn1kio9QPww3wtMeo1IBrPoYqJXwvcmeNrgdS2GGwTTdT19qmdO3dc/vbrWDeVJ/NLIksZKVuLxR2jyUSzrVB3VH6dugol/rV9lMpactBZ2uuRrTXnQTPg1LrEKGnZdL2kZfxNY2txHVvcn1RebdxactBNZWbn3FSlmZwjSxkJ/mjZb2kyeQNraoMuPK4rWrYG+ukJ5PHYHLSVWa89b+koaA9ZVjgW13F+zyqvNm4trsyJvuzakc00uKVXEIC0d5QCbri63J6V7t+1oBEA7J7PEtShkbSfl06G4eC+Hj1TuTnbiyeQJOlGGwNPYI00Jlofa8JhKRQ4H2tLnyR5tbZ7BGedVFWrklLeeE7Y9DX7wz/fizsffmoe59uB0pcwtZm+vCmQ5Knnv78kMIGr+8jTB5fJRENNE98//PO9y/qR7qX6tCzXtoUoX3o99ZGTX6praS/J/T//+vll8lO9aJT0q42BpDtObmlMtD5KZehzw13jnllLnyR5tbZrxovSVOCyibh9Ec0B78nwaiHtNI9S3nENskv3X3mG2zdv22TOdurpoyUXfikvvQY5lmwWSW4ub7oHhlvqR37Psl/25K7vke1VG0ep7VynNEsyJe79KMmp7dmnemSzFsxgCXRoCczQgg0o35Y6rYEakwr4qAk4sgQIDRUY1KqXaQXS5O3W6LdEGAuoJhG3Z9VCAqU9vOQjbU01ZfHp0rao37MHjLaWT00CCcselfLnfL0eCGnJ5sD5mXvpV9ova756j7yWdqXkI5x+uyRUkb4CQ/zQENf86y/NCGMKq5x2W1YaapUh8e/NZxI6LbXNrUx6yGtZQbWEdaNlZg8hnAfgvwD4MQDHAOyMMd4ZQtgI4HMALgDwDIB3xhhftHxgtK8/vVbzBe81q46tLSv1ngVL/HvzmYROS217koR45LWsoDyrLBdJX4H0A2AzgMsW/j4DwJMAXgvgYwBuXLh+I4DbSrzonr0mf7w3CGVS+cEtgQpDBdhIdaaRj77XmEplanTVayw88reMm4VfTSBM0fUWYzwQY3xs4e+/A/BNAOcCeAeAexeK3QvgV7wfGo87QStbe88rg0c+jm/JVVYjvybLxx96sot+PdRrTKUyNbrqNRYe+VvGzcKvhr/LQBdCuADAGwE8CuCcGOMBYP6DEEI4u1R/7tVjYnbSklFJK6stcSRQjaUuJS2YgwJwOOOKFLxQ27cSIMQC7EjUmlWV64u1TQ1I5dmSSGOtGRKlsp4cd5wsEsBHC9LRQEFdtmbSlE9/AJwOYDeAf7rw//fI/ReFetcD2AVg18Zzt3VxXbUYcHrV9RhrWg04tWWHarcXDeWmo/xr+t/btTcp9y1aXW8hhJUAHgDw2Rjj5xcuPx9C2BznZ/XNAF4QPiY7AewEgEvfcFkswWUlF5ZWNpHHxZTPwDXhkCnDrQR1zGeGHrMd1zeL+4jKAvDhpZYVhEVXlrI9QkRrxtr7XKTrpT55QlLz1VZpFddrtZXIYo0PAO4G8M0Y4+3ZrT8B8B4AH134/YUSrwSX5YiC/zXAfym4BSgHlrSGQ96/ex8rgxSKaA1ekHLDcX2z6iGXBfDnikv8PIk/PAFGgD9EtGasa8uWnhVOFomfdpx4j0AmjSwz+5UA3g3gr0MIX1u49u8x/5LfF0J4H4C9AK5rEaTVfSFdl77Kln2tFZzjSaBRmhmp60dLVmANpOD04klI4dFVAkppJ814ZtySTNxYU341wTnc/p6uQCzgF23FkMtXeoZKfTSRtL4f4oe63iZBpb1Sq23AU6YlH/gYbBcWWWrypA8N9hnKHuGxz/Qax1JZKHv2qQbCcFTK0e0NjKHBCjRhghYIYQnqsOR3l2SoCRLRgiM8OpKSblgCYTRdeRJHWIJxPCTJxyWDoG3VnD3AtUfbsgTctCQoobKNJnmFhaS9mGWPpu2dgKVnbnmSCwBybnXLvlCSQWubq8v1sYeOpOv0msWO4kkc0ZKCy8KPylZ7Oqx0j2uPtuWxLdT00aOz0b3snv24tS5XxmIBrinbYnOwkMfP7Kmv7VWH7lNPHla+Hr91ix6G6lONbKPds/eGidZAM2tk6g3nrMkxP1Q+80lCi4eWZew6svKi/cCKlY/HFj/7NKhmSae5fVqWuB6ZapaD2hK6Jse8x1XmodZltsSnZaxrZRm7jqy8aD9WrFl3psRvNAY6mpVTy0giEWckKmWBzbOA0iyimkFGOsiQM+qlsvQwQs1448kKQzOd1BjLLFlRa7LVclTKKGMxFlqMblrfaGYkj4482YA92WZrMizTsR7dwY4cSV9ajzGLMxKVssDS9rQ2LdllOaMeMA+eoIcRal9uzqhHjUwS6KLGWGaZ6XqBPLQc8LksuQz0Wk0Of26sgfpDFC19qwEkpTKe7MlprEd3sKMG/UwZSxNclAMhtGQEtYA/LKAPDiyRZ9jJASgpI0m6pvVN0pUlmCNRDczSYowcythExzOX35MJqGR85EBBluwzAJY8Mx69pr5oICNJ/toMQxJN5WXXoJ8Aiu4ujytE+tKmWZab/TzuNMmVl8/E+YrBekxyqa9Sf4H2fPFSnZac5RppszZgh/dK8uVj1XKij+eY6ESeVVaNW80zJlN52S2QxN4zjJQvnguv9IQ70llDm7VLsEuNP1cW8Ae1SHrxwC9r8rRpZHFXSXnaLO1osFY6TtpKqiZgx7JKLNW1rGpNJJnph/iZBlw2kSfksCbUcChIrFa2Bwy0BTbrydPWSi3hn71grT3lH4oPxgaXrYFF5nUkS7hG1HqrWaw9ucopf4sFXPI8WOrkZUveg5IerRbfUv54TneUn2fMNdgzteDnfafPBdVZXlbK2d7qcZBguNQT49VJouMOLtvqVwX8p4DSfay2l9LCE6V2PXuzGh+vpDOvx6LG4qvVsZ5F5hnzEuw5/18L3y2F+mqejFqPg4Qh4Pb7Ne/BcQeXrdlza/v8mvo10NLekErPHnAoOKeFeljhh5LRYpW3WPBr2rbIo431pPU62j17LzhjTyhmLazV006PzLAeqK2lnRp4aw3/GkiwFxJtlaFV7qHaLPE6ruGyrXDGnlDMWlirp50eS3wP1LYG7FED46wt2xKl56GacWx9RlueHYmXBped6MtOs8tqpAEhEtVk7KR1Le4u7rqUR80D3LC4ID2Am+TeoYAhTiaLGxBYCibxAJ08/ahZ6nLjKpXVxtozjh65LVRyN2tuVknfH/ydl74rNihN+UP8bN3+UxPPjjKUi8XDd9LuqJYjtDjXXuJXo+ehqFd2l0m74DykuVkl+TAW19undu7ccfnbr1sMErEEB3CuHOpS0YIZpMwhnCvEQ5ILhHOjeY4xLgVOaJlqJPef5k7SAkzO37hmCT9Nz6WMPTWZYLQy2nNBdaf13zOOHnktVBprbixSHyR9j8b1dvgHry6Bi1r2KJwrh9ubSW4SKXNIDfRRk0vbL1oy3lA+1kymXBAOdf9p7qQSVJXyK0FpS4FHXN9rQoe154LqTut/jSurd6ivNNbcWNDsyx6360Rn9s/cfdeO33j/+xfDNblQVEtOsBLwQZv9NEAIV98KANFyunnCKVty5klfe8vqiOonB+1o8kry0zoeGSx8uTGSVjaeHG+evIMe4BQ3e5eeRctYUFm0mX0qrjfLvrZmn2W9X6Ie+8He8MjeWVvHtiftUbcnfLj1lBqPvaMnYSx79gSX5WZmCaKqzWDSV96SZEIjKotmW6iZySyrl0RJV9dccg6Ozh1j4ZZUPsvMUwOTpdfzdlv0XAOt1Z4hiZ82jnSVkScLORZhTkBB+5ZmbU+SDA+NHi7L7aFqTu9IpO3nPHsayt+SibYUVplTjb+a7vc5uCWVzwKBrYHJ0ut5u4AfwtwCrdWeIYmf55Qa7uSW2tODknw9TnWhdFzDZXvDOVsgia1+1JIsHsiu5oPuJV9JXo9MQ7TrLUPLesKmOf5D6bmFjhu4bA5vHEPG0Rq+njaHhp1aaCgd1bTX0jdL3V5lpL5MWpcajR4uq0WyTSPjaA3flmiu3rDT1r4MQRbd1fTNszRvLSP1ZdK61MgCl51qdlnNtWA5VimRFrcsGcdoXcuxRxxfCbRjMeZpMdnU1WY5OsrSxxTHbQE2eTKjSuQBRWnjJhkhc8MXNRJqri0KrqLuOk9MvQXYo1FLHdpXLbvsVF/2datX4ooLN2Hd6pVL/s7v3bdr/ou19tT5/zm656t7cOfDT2HL+tW44eqL2TqpDOXDXafXNFnu27UPDzy2H1vWr8YVF24S26H9lfr4jQMv486Hn8LeQ0cW+d5w9cU4d8Oaoj4sfXz20FE88Nh+HDw8hwce27/YjsZPK1Mi2mduvLS+STIkvqnuutUrxXHLy1AZDh6ew6N7DmH72WfghqsvXizD6VIaP46vR2ctdWhfb73l5n2jssZrVDqgniMaFKEFmGjBLdK1Ul6yv/jWwapsuJL8eZbd0nHMJfmkOinQIwXLcDnSUnkpKKQ2H5qkXxpwc+3lW5dlHbYGpnD56nIDW6lvXKCNJ9sr5asFbXn6KOnQQqN72VtyZ1MoYU4eF1lvKGxpb8fBItNvCvG0woe1/uSupQSF5VyUWlZWS7sSlTL+5joo6YPjS/si7cu1vnny0ksycDBtzjaguVetOrTQ6F52S85z+pWsyeftIctsqoUnSl9hSyZTLSspnREt/ZdmUW1F4lkxlNrl+kFdY1zOfcntpZ1BQPvmCaHNr5f0wMlDx4Rz/0kZjy3Pvnc1BYzwZbfkAKNfydaglhLVBONoASsaD8pPAwXVBPVos2jrisHSLsdLO/0m1x3XN46vttqy5syrBWZJY8KBbLgch9Zn3yILpdFll7WclZUs91xAjdVa7Amr5IIZaGZXT1ilFgBSCuWsCbDR6njCVmkfNU9GDYTXoxfpecmvtUBUOZ1Z9OAJepIyHtcE4SQZ/vSP//PBUZ31Zv3K52W1/ZbnTDOLDFLb2rlwnrBKj9++ZD8ohZ5a6lj3wLSPnrRaLfBc7V7J5lK70uN0ZtWDtW0p4zEgr7ZKkOAuaalCCCsA7AKwP8b4SyGEjQA+B+ACAM8AeGeM8UULr15Qxx5QWA9U1QNZrYXwluSq6XOLnvJ6lj62QHhb9NKbPH2rlU2q36K7LmmpAHwQwB8DeHDh/48BuHHh7xsB3FbiocFle0AFxwBb5KgFklnDo8RrEtQCI55E2z0y0XKw2R598WTbpf9DCXE9yfLVCCFsBfBPANyVXX4HgHsX/r4XwK+YP0MLlJYe9+/e5606CJ+hKMn38YeeLMop9cXDo8RrEuRpu7ecFn4tbdK6+f89+qKNtdZ2iazL+N8F8CEAZ2TXzokxHgCAGOOBEMLZJSYpu2y+LAZkd4+UMZaSZRnVdCBeIz9uO5DrgQOySICe1oMlLEcU99CPBHCxlPWMvcX1RsE1uTvUIh9tSwJvXXv5Vhx46ag5m7Gk8yS/5ThxizswUfFlDyH8EoAXYoy7QwhvLXJcXv96ANcDwMZzt3UNTEjUavypIQ8/zrUk9dEC6Kk5hrkkd2/91OQnsOQPoORxvWlGt9q2cvlTGSsASeJXczSZhSwz+5UAfjmEcA2AUwGsDSH8EYDnQwibF2b1zQBe4CrHGHcC2AkAl77hsvjeq7arudVrZjBK3BdTmjG5Ly2dWTR+dDbSZiUKuNCgqol6GKZy+SXQUg08VMtrTvlyR2NT0oAn0li0GtI8qwqqI8tzYZGTjq025tLqIl3HipWrWOXC8LLHGG8CcBMALMzs/y7G+G+3wJgAACAASURBVK4Qwn8E8B4AH134/YUSr1UnnySCB1pmMEocL8uMKa0qNH50NrKsWhLgQoOqJqqBRWr6AOwHGpZme8pXm6VSuxIvgHeZSdDX9L/H9aaVtawqpEzFeVlp5m2FcCcqrcyGOhHmowDuCyG8D8BeANdZK3KzSGmmbN3D01mb2+ukv9PKQ4OS0r7ks7S0IpH2Yq2zd2kG1mYVrUwJuizVyWc0rs8euwznlirBqSlZ9svaPlk6ecfyXPSGcJdWA5rrzfWyxxi/AuArC38fBHCVp34iSw46OlO27uEttgAaxACUZyNuf1XKYU/3Yq2ztxfWap1NStBlKcAml4Xrs2dvTtuoOVLZsl/W9sm0rGW/3Hr0s0Sl1cAHf/jKnFR3onDZT/z+J3f8/YVvU09jobm/W7Nz1iRK8CSMaMlJbklWYDlJhcJLa05f4WCuKdGF59Qc2sckW65DCVJqgTBr8GHp7AEOfitBVbVTY1IfOH6SHixJTWqg21LZ0WSXffHIXDFwQ5spa76Qki2A+zp7yibyzMp0/1qTCslyr6UOzRjrnZ2kPTrVIbcKsngENNsLYE9tJkFVtVNjPBZwy/7e0+8aLxWlib7sG9aswnuv2q7ua3vDIj0WzxrraI0sGvzWKn+tvKU6GiTYQh5LuFU2S3vaOFr00Fs+D39P2y1eqtFkl22BFQ4NB+XaGRoOKkEma7K1erLjevramh14UjDZmvEbSrZ8bDw6s5aFApcdTXbZGgNVbyCIpx1P2zVyetyAPepoZS3bBcCfHXhSQKea8RtKNu1gjdZxKtFUE05yRhFPFljJ4NVq+KLEZYH1HHpYypSrZTJNMfupj5pxU4rj1oyIdAy0+HBa35IdWNNPzQGJmu5Khx96xs9jeLUQd+ikR2dWI+HoDnZM5DmcTytLD8+zHKJXcwjgUAdGWg6xHOqAwDEc8KhR6aDFHodwToJ6H/Qp8cFYDnZMrjcte4k0C3KrAMlNlx/KJ7mMar7cnjraSoS6YzRXDu2jpW8WuWhees0VKa2YLO6klpzzgLxa8Rz3THlde/nWZfn5PYdtSjrl8v57MhdZZNBcj6M62DG53oAyhFACZ+SuEclNZ8nlVQNkac255nHHJKJ9tPTNIhfd12uuSAnYVGM/8JIEfa3Zs3J55Tx9koizXbRkLtJk8GQ+pjQV15sG8KeBAwmayGVQ7eEq8wTEaPU8ARo17h4LRNMSyMO5/2gG1gR3TWVpjnlNfinYp5crVYKuapBbWpfrW87PC3OVYMO1mYssz4Un4+0iSev7IX7ont2yFynt2VrJYgvQ9tLp3tD7QQt/SW6LTUArW7MvntR4te7ZJz2OvUiSE2PZs2vZZaW9iHQWl5d6Zyel9bQ9mSSLZz/rsW5TuXPZjkWo59hxe+BSllmLDBbS9EL36EkWzTtBidODh5/HUyT1pcYmwJH0PIxmz57IE4JKw0F7tsm1S69Z9/yewIehEnRIctMTZzi7iQUSbAlgacnsqulFCkbx5MzX9GDhV+PH72ET4Oi4ORFmYvBAZ5uT4t+7bzWy1dg5hpbbAtmlv1vTdHn4eWxEEr+hn0ONTAkne1HKQQf8aFbpQc8eOoI7vvwknj10RLwntanV9bSZvrStfaJ8S/9rPP5yz0Hc8eUnceClo4v3qJwl/QxBpT5sXrd6EFk8Y3TgpaPLZJTqc9fTtXR8lUcGy/PsfV6BKbveNPIsdVsiimqXVS3LMc9ytQYmWXKZ1fajlzutR7Reqywlfq19tPSppk7LczdRA91n7r5rx2+8//0iIACoiz+nRhXt6B6LgY6SxaCoETXWaEZHyfBH48NzYIjUJy4XAK1jOXIoEQf5tIJpch1KcfKWmHIKLa3NcyDpmQK0vM+FdHyXlhPBU6f03I3GQJdy0CWyBiiU4s+pUYXy0AxQtZlpawA29CBKLXtoIi2Hm6VPNBeABFbi+FPSsrlIfU38avPgSXppPeJJ0rMl4xDta+ojd81iAPXUaclqNFXXm8XNY4HW0mARja+2qpBIm/0sbjQ6C3kgr9oM5nH3UX7SjJbrnUI/NWhtSd683dLMblll1Lgvtaw2ngw4tK+5PmjfJJehZbbmrkvQ2qSH4/JgR092TkvOuBaYoTb7eXLbAXY4r1SXk8sDnS3NaLneAR76ac3cI2V2vX83n1vds8qo2bNrIbmeDDi0T1Rn2irTM1tb8v5TPQyVXdZN+YkwgC8vtkYeNwyFc1rypFtcQvReCYZb6r8F/tnDjcO5nnK4KM2cK41fCWKs5a6vgdhassFK46i5ICn/Ky/atGhRl8qWXG95XzU4rmXsJfh0gv3u/cGRl9mKGKE1vjXAwbrP4sATFuCNtpesPQGEk08rw7VdSxo4CIA5v31pptX27DUAKk82WKmvWhu5HpLdRyKNryfXnScHHx2LtFo86ZQ1ayU5p2qN56h3QgMpZPSaS87B0bljS/aHnoQUJeL205a+Sftwi+XX0v9SOGVtQgrJfsKFeL552yaWr1SXk5+zy9SEvUpU66Wh17QsthJ02WJXkvh+/cG7n7vlI7/9CU7eqVrjOfJYxy2zCOXF7QtzHq0zJW2Hy43vrSf1qae/Ov2vzVKa/CX7iTZTemwvQ/mgpf5oZBkTbWYH7HYD6TmmfLWZfao56GqptEe3hBN6Mry2yOjlWwPJrOn/UPDNFphoC2TV2kZPsoyJRd8eW0CpHe1EmNGEuCbSDqK30lBZZ7Vsqr3k6yH7NPmPsayFJvHMTKIulBDXiWLjLZSWLtxB9F4eNXWtfFva0Or2kH2a/MdY1kKTeGamLdNUs8vmJMFkPXBUDYBDy2qZaCUjVg7LpUamkkw5tcAhLW1o8eeSEdJj1KuB2HoARBywxwJzrok3tzwzJR15DGkWflofSzQauKxGEky2pq7lyFst4EEyYuXuupIxrwY4lN+r6bfEn/aR4+8xfNVAbD0AIg7YI/XR0gcPrFWjGoNoi6uth2s1p6m+7FxOMAmUogEONEOGltvukacPLstBxvHTABwScTJZ8tR5KMmVjpaWQCR5HyW9ciAj2g4F03C6s+ghkUWWEqDF0pbHuGkBQ9G8d7VGw5rx94CulpG0mR/ix5I3XrpXmxusVM+T223onN+1/Er53vJ2Jb168u316k+NLEPT2PPUHTd54+mePd+TlAITPPtwzx7Yk9utFWxTk6dO2w/ToBO6R6cBLBxQxhK2agnUkPKle/TBBdFIxzC35nCT9KwFXFmexR7kyRtPabQnwni+WJav3TRnBAsNfXpJzekxvVctnjGoWVUNtToa0zPTIhPGOrPnVJr1LF87zYLfe0bwkATjrPmCW+ok2CkXtuqdKbg2uZVUWjF4vCjSCTnayswSbqzJWeq/xythGRONJHlb4L7HhTW+ZLW1WLcB2YLfE0rppRorrmSJtVr5rQEsNVZoCa7s9aJYTsjxJPOoSV3lgajWjKNGJehybzK97CGE9QDuAvA6ABHAvwTwBIDPAbgAwDMA3hljfLFWkJJlstVyPWkopaXtGpk8sNNekFiJTyvUuEZOre0ectZYz2v123uciiSt7/MfAPcC+PWFv1cBWA/gYwBuXLh2I4DbSnwuef0bRShs6bD5VqjqmOGttW0PWTevn8aL/h6qz71kmQZUdyi4t1U+KHv24sweQlgL4B8C+LWFj8McgLkQwjsAvHWh2L0AvgLgwxqvFM/emlUV4IEWGrUskSYVaeVte8i6eX0KQOqVeXVoWTz971W2JoOOhX+P58yyjH8NgO8A+HQI4fUAdgP4AIBzYowHACDGeCCEcHaJ0emnnIxLhUwk+ZKLO2iQy45Ss/zVDgKk2Vby31LdHkuv2qw2JeKyotRkwuH0Yc0IxNWR2tX0kIA211xyDgAsyZ6T+GoAqtLhj55x1ABgnq2DNBbcuJUAVBayvOwnA7gMwG/FGB8NIdyJ+WW7iUII1wO4HgA2nrutmFWVy8DBfdG8XzdPlk9p1uid7ZO2m/Pt8SXnjJ61mXBollrvEVe1+f8pVBnAsmdIypxrGXOuryUqGfWshkppLLTjyVP/tb5IZHnZ9wHYF2N8NMmI+Zf9+RDC5oVZfTOAF7jKMcadAHYCwKVvuCxKRzYnol9jDcaZyDNbaV9I+lWmcFALBFaTRSrLzd5SnjYLX6pLC0zUQ95cadbVADc2pTEp1bfCkks59Ep1Svw5Khn88hVJDWSbUvFljzF+O4TwbAjhJ2KMTwC4CsA3Fn7eA+CjC7+/UOLlyVTjObjPM1tpX8hSFlhPtk9OFk9Zi1uq1P/SoZW1ZNG3tiqQSBobbUws9XtkP9LqWMpyVMoqy2WtbTng1Opn/y0Anw0hrALwLQDvxfw5cfeFEN4HYC+A67yNW4JbuC85/QpbZn/K1zLTaDOC9FXX6kgBJZqbR1vhlAI0aDnPbJXa5/bdSZZ8Bk1lpRVOyS5BbSMtmWqk+pbnjavrWTn1JI6/xf4gkelljzF+DcCbmFtXuVoj5AGIcPtOesJKax52raxVdq0Ona2tOd2s2XClvWmtpRkAu++m++e8bGrDsorRZlPrUdnSvZ7ArJqx7kFcPkCL/UGiicJlP/H7n9zx9xe+jQ068CQpoPDY2vO+KHlgjxLMUkuGYYFBSoEZXOZVKaiFQo49gRtcFlh6kk2C4eYn26SyCZpLz9/TxtFzQk5JX9I1qodS4A53XqCkv6GCczyBXamOdiLMRF/2/3D77+34y1Mvw9pTV+KKCzdh3eof/b7nq3tw58NPLd5LlJeh187dsGbJ71ZFSzJwROVKdfceOoIHHtu/hEe6t2X9/FHEmpxUhtTOFx//Nh54bD8OHp5b5P/21/3YEhnu27UPDzy2H1vWry7qUKKcR+Kf9Jtk2H72Gbjh6osX/8/L3rdrfhY8eHgOj+45tCiLNo45n5LeS/qSrlE9aONFeVA9W2RoIU02aSxTncNf+9Kxm2/60Mc5vhPFxm9YswqSNX5wqKCBWmTwwDhrZLDYGnro0LIn1uwTHptID7lrPQ4lOXuMWS216HA02WUTXHbSEEIrvzHyrYFfWtqZVDbVXvqo6VMLdFWDZ1v4Dp39VvofK1Y+HseQXTbBZXtl2hwqc+eY+KY6nmy7lnZa+uhpu5c+avpUozuOVw3f3vot6TX9rx3sOJXjn7T4aotxhRorWvi1Znil7XiO5rUYduhxVbkRy6Mria92zJRkNNXIktWFy9YrGeY8xk0p9p3LwuPJ7EqPYdbyJ3j6xum81Lf0XFF9pPtf+swdz4zCQHf3XZ/a8enbblo04nAGjRpDXQs/rq7HoEXbSUY4C1+LYScZzI7OHVti8PLqSuJLjWI5z28ceNlteKJtc7KkNqgRjyOPcZO2lfqYDIqcDFLfcl5UV9RAXNs3qyy5POm5ovpI92+95eZ9UvKKqbjetEP50hfRkpNc+opa3CaWAwJ7z2ilOloZbhYpuWFqMrRo+eok/VvGKJclzZTc4Zo9MvZY9OtZvXlWZNJKrCWvHMdfWjGMJlONdmSzFPDAlZXqSNBSb8aXGhgkBUBYABfeMtpx0Zz8PTK0aPJZ9MTJAiw/VNMqp6dtKwCnRDUHLgL8sdRcHY8snjz8lEZzZDMFWFi+jFJ21Xymp8fi0vZa99bSLKqBPbSyNVlaKVly8lH5OVuDlNmVy1pbAkXlwJnzN65ZspeuOf1H28OXjqXWxsKjX+3ZoSCo1j18qd+jA9WkPTunVAqwoHtUbl9D904coEPad3r27NqeSgJAaGAPrWySVwJ7WEjbU1LSbA1Ud1T+fN9Ysh/kwJlnDx1dspf2gKIse3gqi2csPPrVnh0Kgmrdw5f6PTpQjUY0qMMSXkqJAyMceOnoskQAXJ3agBgJAJH/5gJWSscZS6GcNaGunlNOaMBNrjstzLKkx/x+aUxqT/9JRENduUAhady0UFfPiT6UvzfhB+Vr1e9oQDWevPFjPYVlUjLUnMLiqWM5jWVaY9LrpBmaw96Tp57Lfz/0s9SDP1py0E2KLFDMHnynQTUylGCzHrisBVLqgfsOPSat7VhWW966PeSyyj3Ysyp9BYb4STP70BlZOZKyk3qgjrXQTwleOSnYZWtmXgvfkrxa2aH674HuttAkMh9b9Yqxzew98qvVttmSlbTWXSe5wiwush664mLUW/hpsnl0N3T/PZmBWqhVvx4XJi3j6dNUjn+qOcLHA4XliEInNXePBHKoBWlIrjDPoYr00EaPGymXrQSUsZDmBvLoTopn94CiNNcZfc44V6HF9dYCrbXwsBz4aQFBrVu9cjygmkQ1R/jkZRLV5J7jspJSKuUG89Sh9/I2OZmkuilDiZZzX5ohuYwnLVQDENHy4dFsRB5QlLYq0J4zTZ+e/tK+3b97Hwt68eQb5MgD7JFoKi97co2kHOBXXrSpKqdbbkzx5gW35GVrOvjewEfrm+Tm4Vxy2r2e/bDILRGX205yYUluy9wVxvXZk9OOewYl4vK4e8u2PsdS3zxjMIqZ/ZGnDwLg9zqePG20Llc2Ucs+yUs1OcxoHS0nn3avZz8sckuk7Ws1qLE0A3N9TnwtOe0kWCtHHoiqVLb1ObasZEo0lZfd4t7w5PH25Fin7WizKgVjUHnzNrRsqtKqQgNuUHCLJUe51k663nuWl3RH+edjtHnd6sVrpRWZBOTR9KHllk8krYY4+bWzBiS+niywUl9qT5qRaCovu7SH5WZey5fMk2Odk4ESrc/NLoAvm6q1b9bgHkufuDot2UktJMmZj1HSh0UWya6hBe5oueUTWe0FuQxAOTjLs9os9UV6hrwBMImmml02J2qt9CRtkJIW1GQrBeRw1dyaLmV21Sy/VF7OIyBZ6C19Kukl7xsN4rBYt1vCSjVZqIXZkkiE60fSsyXRhcRH8zDQgKsaPdT0UfMeUH6aNX6ioJqt23/KDPVsgS1OAyLbC2ZZA5Pt0Y4mfw3/GmqFFbfI12OMerfD1StBgDEWUI0nu6y2RylZRz3WYu0EFIulXrL85ns2jzXXsu/WZJHqUlsDZ48AeOt20oOlHy02AW2/6z0nj9pYLOMote2R09rHXF6LziTvgedZn+jLrp31pllkJUuqtCfzWIs5K7HH9ypZfrlzurwn1lD5atFVmq2BOzVX8pBYrNItln9tv2vFX5QwCSV5PZ6dGvuH5bxAjiTvgedZH00gjIdq/LweXhaLZ8k6yvGtkdfSV6mMxetR+l26VyNviTwylPzW3vP8PPL3eg49em2xxk81b7wneKQm97fWRq+giJpgC08ZT/8nlZfeQ5MODMnvabrrIVdtcI+1bI2MGMueneag8wSP1IIJegQQePj3AuvU9L+mTz3AGhb+NTJpdVqfnR5y1Qb3SHLSsr3BUFPNQce5FqjbiB4maHGn5FQT1JLIEoQj5e/m+Hryp6WcZZ48bVLQjCXLLOdW6hEsouWNl9xHrZl5aZ+4vtVkztXcq5KbTnO9JRmkspZniY7NaFxvnkw1peuTIK3tVjdRqYzmArO2k3h4eHGuLEv9Hm6/oce41b3Wc2xa+yrJAmUZP5UQV4484IkeGUI5krKnWjLQSrzyGVLKb5/LX/raW8gTzktlyGWkWWA9M602e1uPmraQJ3e9JRsu1WEukwR48mQH5vjWPMcS6Gx0Ia4clUIjNXdMr/2mxI9ru+Ty4Fx6EgSWttcKi/SE81IZchkBFENwE5Vchhp8uCUXeu3+vuQ242SyQHcB3X3L8a15jj2BPIlG87JLpAVzaCGMJX4cgIEGXVDAiYe/FqBBs55qQBZPttNEmitHCp3lsqDSLLCe4KQaVxkFA+UgIKktC9DH4srTQkhpcFI+fo88fZDNyMs9Q1zZXAbPc+wJ0U00+pfd+lW2zgiWJALpy10DmrAEaABYtnKQgCyWoBnaVk3oLPAjHeYzjmUlwpHnhBxtVUD1Qakm8YPWJte39BzQ04osqyEK8KGrOyqL9TkebGYPIdwA4NcBRAB/DeC9ANYA+ByACwA8A+CdMcYXNT5zrx5bAi/0hKLm0MQSCMGSJ90yw3jgm4m0mUb6gue/6WzKzbwtemgBjUjhpta2S8SNUdKHFF7K9Yebpa0rEU3P+ewMYHFW1fLo0xk4L+vRkQRdpisIrFi5SmQiWe7SD4BzAewBsHrh//sA/BqAjwG4ceHajQBuK/GigTBDBR9oltReFtShvAc9g3x6W7mHbruXBbzFG2GRv8ZL0TtASvp/xRlnPhcbQTUnA1gdQngF8zP6cwBuAvDWhfv3AvgKgA9rTGggzFDQRMsezUMefq0QylL9acA5J9W2pmfPflazQ7TIUiOTVrZFn9Lv5hNhAHwAwGEA3wHw2YVr3yNlXizxoX52DxxwKHir5V4NvLVWzlKOeQ3y6albar9V/qHgydb2rGUkaK0nl33vvnHjOJG88SGEDQDeAWAbgO8B+O8hhHcVP0U/qn89gOsB4Pzzz19yzwMHHDofnHavBt5aK6fER4Nm0jKWuqX2W+Wn9XtDPz18a6C1nlz2vftmiVa09I2SZRn/CwD2xBi/AwAhhM8D+FkAz4cQNscYD4QQNgN4gascY9wJYCcAvOlNb4r5vdyYVTKqtC5Jawx0petcmVY3oMWIl4jmRpNi6j2uw17bEGrUzPm2GKY8hyuWZMrvJaNbMnhZs9bSwzBL7knNeEz7pLlBad8s42V52fcCeEsIYQ2AowCuArALwPcBvAfARxd+f8HAawnlbhMA7BcW4DOQeqkmh1vpulSmxQ0o8eFk0HKj0Zj6lnhrD0kgKC5jrEWu0izqGRsNHJXcn7nLrJS1lrrkShluuf5w17hngLpBad8sVHzZY4yPhhDuB/AYgFcB/BXmZ+rTAdwXQngf5j8I15V4UdebNFvl7g0KYPDOBoCcmUXi5zn62EOeWcnrlpSOHZb4t8htkdFi4KLycnyoK1ObpUvErS44EA3Xrga24mZgCnbhsgZbssrWZLiVyGSNjzHeAuAWcvkHmJ/lzURDXKWTStKsxAEY8nISaaAMSxaTofadnlnJuxIZMk98jT5qQTTS6oVmZmldtZSeA7rqLIGt6Ayc6lBAErcKoG1z+vBkuJVoogi60085GZcycENu/0K/mpZc4Im4rzL9enL8LLnc8uteskA7E1n2YqXTUmrsBiU7Qmll5lkFePKlW+wonpVXaVWUz8A0zz0HeaV86fNlccFpfdOeWytNNOrttt/9xI7nNv8cDh6ewwOP7cfeQ0fwwGP7sfbUlbjiwnnh7/nqHtz58FPYsn41brj6Ypy7YQ2uuHATvvj4t/HAY/uxZf3qxbIS3bdr32LZZw8dXVJv3eqVIj/adoo8kq57KZer1Ickp9ZWkovq0VKX8kh16f+5LPft2rfkHtcOV1+6XtL32lNX4u2v+zG2L562Nf2m54u2ffDwHB7dcwhb1q9elCH1P79H25GeL649aZy469pzm5MW9TaV5BU0qcC1l29dlvObhpVakiAk0k4t1RJI1OQo9yR20JJWSGQJFU0JPrhTcUtEQy41PUshurU50CV9aDJ4Tm2tSbqhnVdAE6po+pZCfnvn5ffkjZ9adlkadMFlOS3tb62BIK37TG6fCPgz0WplJLKEikoBNhaS9sScDJZTaiwejdJ+W5PBc2przdhoASbUcp/KcPw9gTaSLBzVZB9ONJqotxIM0FKnpp0WGYG6DKYtbWr7ud6w0BL/2n61jJvWR1qmZmxq2rH0w8N3sOdZgtYN8aOlpSpBB3vBFi2w1qEzrkrt9cqy66kzNPRTk88D3/XAWT3t9XiGPPq13GuBAGMs2WVryANv7QVr7QV9tVINnFOTyVNnaOinJp+ljRo4q6e9Hs9Q7dK8BzxbK0NpNDnoEnkyjVLjkiWXWSn7J1eGGgs5OTUDTMnwosktZUrVjHyeOtqhmMci1D7mffUYLC3ZY0t9SUYyS946SybadE8zxkoHc3L6pUbZxI/LN2gxRlOSnpnRGOgs5AGeUOOSN2uJtYznOKKaHOWa3FomGYk8dah+OQiz1EcuUKMmi01NXyzHMlvakzLqcMbY1CeLfumYSyAeqzGaUk3evtG97BKAwwL26AUPTaQd+scZdDhI5bWXlw92tIBtLIARC6AFWHp4pQWQI8FDOWNTjZHQQxzoxaoHD2nALK0d6ZnhoLWJPAAqaZzSfS1TzehedgkOqbl5aqCTFpK+zul/bkbkIJWAnies14GJlj0g4As4KsFDc3lLK6YeJLkBuTK5TF7i4LK0LUtQiyVfomUmL6140/0Va9adKfGY6J79E7//yR1/f+HbTPtaunfiABgS2MNDXM5vCrDQTgBJJO39OGAPpZY9moVPfn3DaSvVnPIaUIbuO3vl6/eQx+ZC71lO+OHy57952yYWyMI9dzV2JAuV+KT7X3/w7udu+chvf4LjMZWz3iz7WsusbfnKl0gLmvEEH5SyqWpf7l57NAugJQVqSDnlLSuo1vP3Wshrc8nJ4/XInym6itOeuxo7koVKfNL9k05Zs1biMVW4bG6FTDMO/XJZLKkWy6xkSda+4B7+GkmWegs8kvZVsxJb+JVmCM8ZapxeJBm4657+9yANhptWbRqUWzqHkIPUaitSSpqXg54zVxq3L33mjmd23PyRP+DamejLfvddn9rx6dtuWgwGyAMraMBKIi2owxocACwPGkm/uYCH1JaHv0ZSwEpN4Eaus28ceLkYxEKpFIxj0bcWnFQTCGPpfw/Sgmdo4FUesELl1J4LOk6W4Cnp2cz1Yh23W2+5ed8oXW8a9LOVj1SGWpKHhrtybbfATTWdDQUxbpGlpuwkqUVnWl1P3yww3y66kqB1Q/xc8vo3DpLltDXDqLVM78yrk8paWwPZ9JYZmlp0pWXmrcnE2+vZrHneShmFsWLl41F4/06q/0z4KRno7t+9r1g2GUF6le1RxiOTpb5Hpo8/9GR121o7vXQ3NLXoiqvrKVsji6VszfMmPUOpKF/cZwAABRtJREFUH6NxvSUDnQUm6XFZUBeZJb6YI8nFZIEx1hjHPC63FF+d97F07DDXd8nYZNGL5p4rUasxzuJOo8dda0YyyYXX4uLk+umBe2vtSDkAqLFwNK63PJ69RB6XhcVFVpOF1JPh1gLksMSOW/pWyiUu9QtYHkvuzZxr7SulVqCL1Z1Wgjlz/PLnq8XFycnjgXtr7Uh9oTH2o3G9WQJhEnmCJbTsIok8mUKoG8YTsOGZIS2kBV0kd6UnM4tnlrZkAir1g8vSUwqasawCLNl/WvRtmbUtAKQWNygnl8R3dKAaD3mCJSzH19ZmCrECIWpmSC9f6UhlwJ6ZxTNLWzMBacTxKGV27RXi2arvUputEG7P6rXEd3SgmpoQVwvVhK1yM5uUr86TN8wzK7XCTel+LUFge81kpdlEgxpr/EshxK37Zko1qxeOPDaXGhuRZyVD7UoJHKbN7KN72WuAFhTsYQGEpDJaJlMJDKPJa5HfwtfT7wTu2H72GUUAh8RD0xnNKsv1nQJ8NP50LCh/SxZeTW5K2ph4njcql0V3mnwl/XJyJpBOqkPBYX/3V1985eabPvRxjsdUlvGeHOJDtCW1p53pZc0bZgE/pPBJempILZXa5EJcPeGfNUAZjaieabhqzrcmD326poUoJz1YnjcaVttyShF9liT5aFi39HxRII52ZPNUXnZLaOCQbUnXtb2k1TJbY2FvDZIotakF+7Tw9wT7cPJwQSfUSyHtUS3PkBZu6tEDDXwB6k8p0k5mLYV1S89Zbsv54A9fmZNkmMrLPknoZM2sNHQChkm1Q9trhSX3Igs8lJatgd9aVl0tcNZe2WU98jaRBK0b4idll/XAAHtlPe0B+ewN4e0lQ43uamChQ0OYPVl2a2TI/+/5PLRmIe6h59HBZRN5YIAeGKO3zSF5DAUx7aW7Gljo0BBmC7TUAxvW+Pd8HlqgzCU5vTKMBi6brPGaO0aCONYAJIC6o5d6gElaobVS2VbdSbBLy1FOJTixtf9SGU+8ucZXGnMte5Alc670PGjZZWuATh6XnAcuO7pDIm5/6In44x9+MN7+0BPFshaq4ddbhhb+vctKZWraeecn/2xQPXnl6lEn9SnvWw9+NTqqGZMVZ5z5XDxeDomwGCVK2UMtR/72kkHK0qrJmxtrcleLVtYqS6ks5eepKx2nPaSxzzN+pb5ozw11h+ZuUXpks+Zqq9GRJBeXdVgqm/iPzvWmUU2mzdL9npBJ2kbNqSAa9NETSCHx9/TNU7cFTlxLrdBc630p4CjPQWeBwNboSJKLg9OWIMya6210CDoLlcIda45GrpXBAlWV9mCW/XdNMEhNbjuLrmqORK4hCwzXc1y0pa9awFHpiGlLG9oeXtKrlr1WCuc+rk6EsVBN9tChZEihhXm2Vq1sadb2BFJ4AEOlflio5kjkGrKAgLQgFImsJ8OUjhPXqPRscjqT9Kplr01lLTIlCjFGU8EeFEL4DoDvAxD3FU20YuWqFWvWnfnDIy99F8pyxklnQpPX02bvsnyZM7Fi5csD6GGxvWM/OPLySaesWZt+N7azXL95vwCwfRlmrHla2tbaZfI6eLA6E/RqKcv0/8djjGdxIkz0ZQeAEMKuGOObJtpoA83kHZZm8k6OpgKqmdGMZjR5mr3sM5rRCULTeNlZS+GIaSbvsDSTd0I08T37jGY0o+nQbBk/oxmdIDR72Wc0oxOEZi/7jGZ0gtDsZZ/RjE4Qmr3sM5rRCUL/D+NVQ7mKgl8IAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,\n",
      "        1, 1, 1], dtype=torch.int16)\n"
     ]
    }
   ],
   "source": [
    "#Plot Adj matrix\n",
    "\n",
    "W = data.W\n",
    "plt.spy(W,precision=0.01, markersize=1)\n",
    "plt.show()\n",
    "\n",
    "idx = np.argsort(data.rand_idx) \n",
    "W = data.W\n",
    "W2 = W[idx,:]\n",
    "W2 = W2[:,idx]\n",
    "plt.spy(W2,precision=0.01, markersize=1)\n",
    "plt.show()\n",
    "\n",
    "target = data.node_label\n",
    "target = target[idx]\n",
    "print(target)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'nb_clusters': 5, 'size_min': 5, 'size_max': 35, 'p': 0.5, 'q': 0.2, 'p_pattern': 0.5, 'q_pattern': 0.5, 'vocab_size': 3, 'size_subgraph': 20}\n",
      "pattern: 0\n",
      "pattern: 1\n",
      "pattern: 2\n",
      "pattern: 3\n",
      "pattern: 4\n",
      "pattern: 5\n",
      "pattern: 6\n",
      "pattern: 7\n",
      "pattern: 8\n",
      "pattern: 9\n",
      "pattern: 10\n",
      "pattern: 11\n",
      "pattern: 12\n",
      "pattern: 13\n",
      "pattern: 14\n",
      "pattern: 15\n",
      "pattern: 16\n",
      "pattern: 17\n",
      "pattern: 18\n",
      "pattern: 19\n",
      "pattern: 20\n",
      "pattern: 21\n",
      "pattern: 22\n",
      "pattern: 23\n",
      "pattern: 24\n",
      "pattern: 25\n",
      "pattern: 26\n",
      "pattern: 27\n",
      "pattern: 28\n",
      "pattern: 29\n",
      "pattern: 30\n",
      "pattern: 31\n",
      "pattern: 32\n",
      "pattern: 33\n",
      "pattern: 34\n",
      "pattern: 35\n",
      "pattern: 36\n",
      "pattern: 37\n",
      "pattern: 38\n",
      "pattern: 39\n",
      "pattern: 40\n",
      "pattern: 41\n",
      "pattern: 42\n",
      "pattern: 43\n",
      "pattern: 44\n",
      "pattern: 45\n",
      "pattern: 46\n",
      "pattern: 47\n",
      "pattern: 48\n",
      "pattern: 49\n",
      "pattern: 50\n",
      "pattern: 51\n",
      "pattern: 52\n",
      "pattern: 53\n",
      "pattern: 54\n",
      "pattern: 55\n",
      "pattern: 56\n",
      "pattern: 57\n",
      "pattern: 58\n",
      "pattern: 59\n",
      "pattern: 60\n",
      "pattern: 61\n",
      "pattern: 62\n",
      "pattern: 63\n",
      "pattern: 64\n",
      "pattern: 65\n",
      "pattern: 66\n",
      "pattern: 67\n",
      "pattern: 68\n",
      "pattern: 69\n",
      "pattern: 70\n",
      "pattern: 71\n",
      "pattern: 72\n",
      "pattern: 73\n",
      "pattern: 74\n",
      "pattern: 75\n",
      "pattern: 76\n",
      "pattern: 77\n",
      "pattern: 78\n",
      "pattern: 79\n",
      "pattern: 80\n",
      "pattern: 81\n",
      "pattern: 82\n",
      "pattern: 83\n",
      "pattern: 84\n",
      "pattern: 85\n",
      "pattern: 86\n",
      "pattern: 87\n",
      "pattern: 88\n",
      "pattern: 89\n",
      "pattern: 90\n",
      "pattern: 91\n",
      "pattern: 92\n",
      "pattern: 93\n",
      "pattern: 94\n",
      "pattern: 95\n",
      "pattern: 96\n",
      "pattern: 97\n",
      "pattern: 98\n",
      "pattern: 99\n",
      "10000 2000 2000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAS/UlEQVR4nO3df5BdZ33f8fcncjC/YizXkiMktTIdQWN3BkNVlZQ2IYjGbkiR/6hn1DYZNVVGMxk3kzBpsdT0R9KpWpF00nSmIRkNkCgJwaNQqNWkSRFKXaYzYLEGA5ZtjQUy0iJhLc6QBNJRkPPtH/e4XEv3aq927929++z7NaM55zznObvfI+1+9tFzfmyqCklSW75tuQuQJI2f4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXboOSX4lyb9a7jqk+cT73LWaJHkG+NGq+thy1yJNkiN3qZPkhuWuQRoXw12rRpLfAP4i8N+TfD3Ju5JUkj1JzgJ/0PX77SRfSfJHST6e5M6+j/FrSf5dt/6WJLNJfirJxSQXkvzIspycdAXDXatGVf0wcBb4e1X1SuBIt+t7ge8C7u62fw/YCqwHPg184Bof9juBVwEbgT3ALyVZO/7qpetjuEvwM1X1jar6vwBV9f6q+pOqugT8DPD6JK8acuw3gX9bVd+sqv8BfB143ZJULV2D4S7BuRdWkqxJcjDJF5L8MfBMt+vWIcc+V1WX+7b/FHjlZMqURme4a7UZdHtYf9s/BHYCb6M33bKla89ky5LGy3DXavMs8Jpr7P8O4BLwHPBy4N8vRVHSuBnuWm3+A/Avk3wN+PsD9v868CXgy8ATwCeXsDZpbHyISZIa5MhdkhpkuEtSgwx3SWqQ4S5JDZqKFyXdeuuttWXLluUuQ5JWlEcfffSrVbVu0L6pCPctW7YwMzOz3GVI0oqS5EvD9jktI0kNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDZqKJ1SlabZl3+8ObH/m4NuXuBJpdI7cJalBhrskNchpGWmBnK7RNHPkLkkNMtwlqUGGuyQ1aKRwT3Jzkg8leSrJk0m+O8ktSY4lebpbru3rvz/J6SSnktw9ufIlSYOMOnL/z8DvV9VfAV4PPAnsA45X1VbgeLdNkjuAXcCdwD3Ae5KsGXfhkqTh5g33JDcB3wO8D6Cq/qyqvgbsBA533Q4D93brO4EHq+pSVZ0BTgPbx124JGm4UUburwHmgF9N8pkk703yCuC2qroA0C3Xd/03Auf6jp/t2iRJS2SUcL8BeCPwy1X1BuAbdFMwQ2RAW13VKdmbZCbJzNzc3EjFSpJGM0q4zwKzVfVIt/0hemH/bJINAN3yYl//zX3HbwLOX/lBq+pQVW2rqm3r1q1baP2SpAHmDfeq+gpwLsnruqYdwBPAUWB317YbeKhbPwrsSnJjktuBrcCJsVYtSbqmUV8/8OPAB5K8BPgi8CP0fjAcSbIHOAvcB1BVJ5McofcD4DJwf1U9P/bKJUlDjRTuVfUYsG3Arh1D+h8ADiyiLmlifCeMVgOfUJWkBhnuktQgw12SGmS4S1KD/GUdWtG8OCoNZrhLS8QfRFpKTstIUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQSOGe5Jkkn0/yWJKZru2WJMeSPN0t1/b135/kdJJTSe6eVPGSpMGuZ+T+fVV1V1Vt67b3AceraitwvNsmyR3ALuBO4B7gPUnWjLFmSdI8FjMtsxM43K0fBu7ta3+wqi5V1RngNLB9EZ9HknSdRg33Aj6a5NEke7u226rqAkC3XN+1bwTO9R0727W9SJK9SWaSzMzNzS2seknSQDeM2O/NVXU+yXrgWJKnrtE3A9rqqoaqQ8AhgG3btl21X5K0cCON3KvqfLe8CHyE3jTLs0k2AHTLi133WWBz3+GbgPPjKliSNL95wz3JK5J8xwvrwPcDjwNHgd1dt93AQ936UWBXkhuT3A5sBU6Mu3BJ0nCjTMvcBnwkyQv9f6uqfj/Jp4AjSfYAZ4H7AKrqZJIjwBPAZeD+qnp+ItVLkgaaN9yr6ovA6we0PwfsGHLMAeDAoquTJC3IqBdUJS2xLft+d2D7MwffvsSVaCXy9QOS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQvyBbU8VfCi2NhyN3SWqQ4S5JDRo53JOsSfKZJL/Tbd+S5FiSp7vl2r6++5OcTnIqyd2TKFySNNz1zLn/BPAkcFO3vQ84XlUHk+zrth9IcgewC7gTeDXwsSSvrarnx1i3pCt4vUL9Rhq5J9kEvB14b1/zTuBwt34YuLev/cGqulRVZ4DTwPbxlCtJGsWo0zK/CLwL+PO+ttuq6gJAt1zftW8EzvX1m+3aXiTJ3iQzSWbm5uauu3BJ0nDzhnuSHwQuVtWjI37MDGirqxqqDlXVtqratm7duhE/tCRpFKPMub8ZeEeSHwBeCtyU5DeBZ5NsqKoLSTYAF7v+s8DmvuM3AefHWbQ0CcPmrKWVaN6Re1Xtr6pNVbWF3oXSP6iqHwKOAru7bruBh7r1o8CuJDcmuR3YCpwYe+WSpKEW84TqQeBIkj3AWeA+gKo6meQI8ARwGbjfO2UkaWldV7hX1cPAw936c8COIf0OAAcWWZskaYF8QlWSGmS4S1KDDHdJapDhLkkNMtwlqUH+sg41yQeStNo5cpekBjlyl1YxXxPcLkfuktQgw12SGmS4S1KDnHPXRI1rTte7X6TrY7hLY+YPIk0Dw11aZv4w0CQ45y5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ2aN9yTvDTJiSSfTXIyyc927bckOZbk6W65tu+Y/UlOJzmV5O5JnoAk6WqjjNwvAW+tqtcDdwH3JHkTsA84XlVbgePdNknuAHYBdwL3AO9JsmYSxUuSBpv33TJVVcDXu81v7/4UsBN4S9d+GHgYeKBrf7CqLgFnkpwGtgOfGGfh0mrlu2g0ipHm3JOsSfIYcBE4VlWPALdV1QWAbrm+674RONd3+GzXJklaIiOFe1U9X1V3AZuA7Un+6jW6Z9CHuKpTsjfJTJKZubm50aqVJI3kuu6Wqaqv0Zt+uQd4NskGgG55ses2C2zuO2wTcH7AxzpUVduqatu6desWULokaZhR7pZZl+Tmbv1lwNuAp4CjwO6u227goW79KLAryY1Jbge2AifGXbgkabhRflnHBuBwd8fLtwFHqup3knwCOJJkD3AWuA+gqk4mOQI8AVwG7q+q5ydTvpbauH5tnqTJGuVumc8BbxjQ/hywY8gxB4ADi65OkrQgPqEqSQ0y3CWpQYa7JDXIcJekBo1yt4w0dj5CL02WI3dJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg3yISWqcD4ytTo7cJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoO8FVLSyIbdVvnMwbcvcSWajyN3SWqQ4S5JDTLcJalB84Z7ks1J/leSJ5OcTPITXfstSY4lebpbru07Zn+S00lOJbl7kicgSbraKCP3y8BPVdV3AW8C7k9yB7APOF5VW4Hj3Tbdvl3AncA9wHuSrJlE8ZKkwea9W6aqLgAXuvU/SfIksBHYCbyl63YYeBh4oGt/sKouAWeSnAa2A58Yd/GaHr6cSpou1zXnnmQL8AbgEeC2Lvhf+AGwvuu2ETjXd9hs13blx9qbZCbJzNzc3PVXLkkaauRwT/JK4L8CP1lVf3ytrgPa6qqGqkNVta2qtq1bt27UMiRJIxgp3JN8O71g/0BVfbhrfjbJhm7/BuBi1z4LbO47fBNwfjzlSpJGMcrdMgHeBzxZVb/Qt+sosLtb3w081Ne+K8mNSW4HtgInxleyJGk+o7x+4M3ADwOfT/JY1/YvgIPAkSR7gLPAfQBVdTLJEeAJenfa3F9Vz4+9ckkT4wXylW+Uu2X+D4Pn0QF2DDnmAHBgEXVJkhbBJ1QlqUGGuyQ1yFf+rnK+wlVqkyN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkE+oaiDfCiitbI7cJalBjtwlTYzvLlo+jtwlqUGGuyQ1yHCXpAYZ7pLUIC+oNsYLWJLAcJc0Bj4XMX2clpGkBhnuktQgw12SGjRvuCd5f5KLSR7va7slybEkT3fLtX379ic5neRUkrsnVbgkabhRLqj+GvBfgF/va9sHHK+qg0n2ddsPJLkD2AXcCbwa+FiS11bV8+MtW1KLvNtrfOYduVfVx4E/vKJ5J3C4Wz8M3NvX/mBVXaqqM8BpYPuYapUkjWihc+63VdUFgG65vmvfCJzr6zfbtV0lyd4kM0lm5ubmFliGJGmQcV9QzYC2GtSxqg5V1baq2rZu3boxlyFJq9tCw/3ZJBsAuuXFrn0W2NzXbxNwfuHlSZIWYqFPqB4FdgMHu+VDfe2/leQX6F1Q3QqcWGyRktriE62TN2+4J/kg8Bbg1iSzwL+hF+pHkuwBzgL3AVTVySRHgCeAy8D93ikjSUtv3nCvqn8wZNeOIf0PAAcWU5TGz5GStLr4hKokNchwl6QGGe6S1CDf5z7FfBRb0kI5cpekBjlyX4G880XSfAx3SSuWU5fDOS0jSQ0y3CWpQYa7JDXIcJekBnlBdQl58UdaGO8Qu36O3CWpQYa7JDXIaZkp4H85JY2bI3dJapAjd0nN8eYFw13SKrKaQt9pGUlqkOEuSQ1yWmYCvPtF0nJz5C5JDXLkvgiO0CVNK8N9BIa41LaFfI9P+x02E5uWSXJPklNJTifZN6nPI0m62kRG7knWAL8E/B1gFvhUkqNV9cQkPt+4OEKXtFjTci/9pKZltgOnq+qLAEkeBHYCEwl3Q1nSUrve3Fnq0J9UuG8EzvVtzwJ/o79Dkr3A3m7z60lOTaiWUdwKfHUZP/84rPRzsP7lt9LPYUXWn3f//9WF1P+Xhu2YVLhnQFu9aKPqEHBoQp//uiSZqapty13HYqz0c7D+5bfSz8H6X2xSF1Rngc1925uA8xP6XJKkK0wq3D8FbE1ye5KXALuAoxP6XJKkK0xkWqaqLif5p8D/BNYA76+qk5P4XGMyFdNDi7TSz8H6l99KPwfr75Oqmr+XJGlF8d0yktQgw12SGrQqwz3JzUk+lOSpJE8m+e4ktyQ5luTpbrl2uescJsk7k5xM8niSDyZ56bTXn+T9SS4mebyvbWjNSfZ3r644leTu5an6W4bU//Pd19Dnknwkyc19+6a+/r59/yxJJbm1r22q6ofh55Dkx7s6Tyb5ub72qTqHIV9DdyX5ZJLHkswk2d63b3H1V9Wq+wMcBn60W38JcDPwc8C+rm0f8O7lrnNI7RuBM8DLuu0jwD+e9vqB7wHeCDze1zawZuAO4LPAjcDtwBeANVNY//cDN3Tr715p9Xftm+nd+PAl4NZprf8a/wbfB3wMuLHbXj+t5zCk/o8Cf7db/wHg4XHVv+pG7kluoveX/D6AqvqzqvoavdcjHO66HQbuXZ4KR3ID8LIkNwAvp/cMwVTXX1UfB/7wiuZhNe8EHqyqS1V1BjhN75UWy2ZQ/VX10aq63G1+kt7zHLBC6u/8J+BdvPghw6mrH4aew48BB6vqUtfnYtc+decwpP4CburWX8W3ngdadP2rLtyB1wBzwK8m+UyS9yZ5BXBbVV0A6Jbrl7PIYarqy8B/BM4CF4A/qqqPskLqv8Kwmge9vmLjEtd2vf4J8Hvd+oqoP8k7gC9X1Wev2LUi6u+8FvjbSR5J8r+T/PWufaWcw08CP5/kHL3v6/1d+6LrX43hfgO9/xr9clW9AfgGvSmBFaGbl95J779qrwZekeSHlreqsZv39RXTJMlPA5eBD7zQNKDbVNWf5OXATwP/etDuAW1TVX+fG4C1wJuAfw4cSRJWzjn8GPDOqtoMvJNuRoEx1L8aw30WmK2qR7rtD9EL+2eTbADolheHHL/c3gacqaq5qvom8GHgb7Jy6u83rOYV8/qKJLuBHwT+UXWTpayM+v8yvQHCZ5M8Q6/GTyf5TlZG/S+YBT5cPSeAP6f3Aq6Vcg676X0PA/w235p6WXT9qy7cq+orwLkkr+uadtB7FfFRen/RdMuHlqG8UZwF3pTk5d0IZQfwJCun/n7Daj4K7EpyY5Lbga3AiWWo75qS3AM8ALyjqv60b9fU119Vn6+q9VW1paq20AuTN3bfH1Nff5//BrwVIMlr6d0g8VVWzjmcB763W38r8HS3vvj6l/Pq8TJetb4LmAE+R++LYy3wF4Dj3V/uceCW5a7zGvX/LPAU8DjwG/SuqE91/cAH6V0j+Ca9INlzrZrpTRl8AThFdzfBFNZ/mt686GPdn19ZSfVfsf8ZurtlprH+a/wbvAT4ze574dPAW6f1HIbU/7eAR+ndGfMI8NfGVb+vH5CkBq26aRlJWg0Md0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSg/wft2E8DChbyNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEICAYAAABYoZ8gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAASN0lEQVR4nO3dfYylZXnH8e+vbEHBUKA7i8BiB81qi6ZVOqW+VEtdX1CJS5NisNVsFbOt8aWSWl1KoukfJPgSrU2rzUbRtQq4opZNrApug6ZJAQcUBVbrKrgMrOxQ0BpNkMWrf5yHeFxmdnbOy87M3t9PsjnnuZ/nzLmu3Znfufc+z3kmVYUkqR2/ttQFSJIOLYNfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr80gCRnJZlZ6jqkQRj8ktQYg1+SGmPwq2lJNie5ar+xDyT5pySvSbIzyU+SfD/JXy1VndIoGfxq3RXAS5McC5DkCOAVwOXAXuAc4FjgNcD7k5yxVIVKo2Lwq2lV9QPgZuDcbuj5wM+q6vqq+nxVfa96vgJcAzx3qWqVRsXgl3qz+1d29/+82ybJS5Jcn+T+JD8CXgqsXqIapZEx+CX4NHBWkrXAnwKXJzkK+AzwXuDEqjoO+A8gS1emNBoGv5pXVbPAdcBHgTuqaidwJHAUMAvsS/IS4EVLVqQ0Qga/1HM58ILulqr6CfBmYBvwAL0loO1LVp00QvEXsUhSW5zxS1JjDH5JaozBL0mNMfglqTGrFjogyWX0Pra+t6qett++twLvASaq6r5u7CLgAuBh4M1V9aWFnmP16tU1OTm5+OolqWE33XTTfVU1sdjHLRj8wMeAfwY+3j+Y5FTghcDuvrHTgfOBpwInA19O8uSqevhATzA5Ocn09PTiKpekxiX5wSCPW3Cpp6q+Ctw/x673A28D+s8H3QBcWVUPVtUdwC7gzEEKkySNx0Br/EleDtxdVbfst+sU4K6+7ZlubK6vsSnJdJLp2dnZQcqQJA1g0cGf5GjgYuAdc+2eY2zOT4hV1ZaqmqqqqYmJRS9RSZIGdDBr/Pt7EnAacEsSgLXAzUnOpDfDP7Xv2LXAPcMWKUkanUXP+KvqW1W1pqomq2qSXtifUVU/pHctk/OTHJXkNGAdcONIK5YkDWXB4E9yBfDfwFOSzCS5YL5jq+o2ehe1uh34IvCGhc7okSQdWgsu9VTVKxfYP7nf9iXAJcOVJUkaFz+5K0mNMfglqTGDnNUjCZjc/Pk5x++89GWHuBJpcZzxS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY1ZMPiTXJZkb5Jb+8bek+TbSb6Z5HNJjuvbd1GSXUm+k+TF4ypckjSYg5nxfww4e7+xa4GnVdXvAv8DXASQ5HTgfOCp3WM+mOSIkVUrSRragsFfVV8F7t9v7Jqq2tdtXg+s7e5vAK6sqger6g5gF3DmCOuVJA1pFGv8rwW+0N0/Bbirb99MN/YoSTYlmU4yPTs7O4IyJEkHY6jgT3IxsA/45CNDcxxWcz22qrZU1VRVTU1MTAxThiRpEVYN+sAkG4FzgPVV9Ui4zwCn9h22Frhn8PIkSaM2UPAnORt4O/DHVfWzvl3bgcuTvA84GVgH3Dh0ldJhYHLz5+ccv/PSlx3iStS6BYM/yRXAWcDqJDPAO+mdxXMUcG0SgOur6q+r6rYk24Db6S0BvaGqHh5X8ZKkxVsw+KvqlXMMf+QAx18CXDJMUZKk8fGTu5LUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDVm4F+9KGm8/I1dGhdn/JLUGINfkhpj8EtSYwx+SWqMwS9JjVkw+JNclmRvklv7xk5Icm2S73a3x/ftuyjJriTfSfLicRUuSRrMwcz4Pwacvd/YZmBHVa0DdnTbJDkdOB94aveYDyY5YmTVSpKGtmDwV9VXgfv3G94AbO3ubwXO7Ru/sqoerKo7gF3AmSOqVZI0AoOu8Z9YVXsAuts13fgpwF19x810Y4+SZFOS6STTs7OzA5YhSVqsUb+5mznGaq4Dq2pLVU1V1dTExMSIy5AkzWfQ4L83yUkA3e3ebnwGOLXvuLXAPYOXJ0katUGv1bMd2Ahc2t1e3Td+eZL3AScD64Abhy1SGiWvgaPWLRj8Sa4AzgJWJ5kB3kkv8LcluQDYDZwHUFW3JdkG3A7sA95QVQ+PqXZJ0gAWDP6qeuU8u9bPc/wlwCXDFCVJGh8/uStJjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktSYQa/HLy0LXltfWjyDXzrMzffiCL5AtsqlHklqjMEvSY1xqUfLimv20vg545ekxhj8ktQYl3qkJXags26kcXDGL0mNGWrGn+RC4HVAAd8CXgMcDXwKmATuBF5RVQ8MVaW0SM6ipfkNPONPcgrwZmCqqp4GHAGcD2wGdlTVOmBHty1JWiaGXeNfBTw2yUP0Zvr3ABcBZ3X7twLXAW8f8nmksfNUUrVi4Bl/Vd0NvBfYDewBflxV1wAnVtWe7pg9wJq5Hp9kU5LpJNOzs7ODliFJWqRhlnqOBzYApwEnA8ckedXBPr6qtlTVVFVNTUxMDFqGJGmRhjmr5wXAHVU1W1UPAZ8Fng3cm+QkgO527/BlSpJGZZjg3w08M8nRSQKsB3YC24GN3TEbgauHK1GSNEoDv7lbVTckuQq4GdgHfB3YAjwO2JbkAnovDueNolBJ0mgMdVZPVb0TeOd+ww/Sm/1LkpYhP7krSY0x+CWpMV6kTVqAl3/Q4cbgl1YYP2GsYbnUI0mNMfglqTEu9Uh6FJeTDm/O+CWpMQa/JDXG4Jekxhj8ktQYg1+SGuNZPRoJzwKRVg5n/JLUGINfkhpj8EtSYwx+SWqMb+5KhwkvH62D5YxfkhrjjF8aMWfeWu6c8UtSYwx+SWrMUMGf5LgkVyX5dpKdSZ6V5IQk1yb5bnd7/KiKlSQNb9gZ/weAL1bVbwO/B+wENgM7qmodsKPbliQtEwMHf5JjgecBHwGoqp9X1Y+ADcDW7rCtwLnDFilJGp1hZvxPBGaBjyb5epIPJzkGOLGq9gB0t2vmenCSTUmmk0zPzs4OUYYkaTGGCf5VwBnAh6rqGcBPWcSyTlVtqaqpqpqamJgYogxJ0mIME/wzwExV3dBtX0XvheDeJCcBdLd7hytRkjRKAwd/Vf0QuCvJU7qh9cDtwHZgYze2Ebh6qAolSSM17Cd33wR8MsmRwPeB19B7MdmW5AJgN3DekM8hSRqhoYK/qr4BTM2xa/0wX1eSND5eq0crgte/kUbHSzZIUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYz+PXong+vbTyGfwaq/leKO689GWHuBJJj3CpR5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxngev5aEHwSTlo7BLzXMF+A2udQjSY0x+CWpMUMv9SQ5ApgG7q6qc5KcAHwKmATuBF5RVQ8M+zw6tFwCkA5fo5jx/w2ws297M7CjqtYBO7ptSdIyMVTwJ1kLvAz4cN/wBmBrd38rcO4wzyFJGq1hZ/z/CLwN+EXf2IlVtQegu10z1wOTbEoynWR6dnZ2yDIkSQdr4OBPcg6wt6puGuTxVbWlqqaqampiYmLQMiRJizTMm7vPAV6e5KXAY4Bjk3wCuDfJSVW1J8lJwN5RFCpJGo2BZ/xVdVFVra2qSeB84D+r6lXAdmBjd9hG4Oqhq5Qkjcw4Prl7KbAtyQXAbuC8MTyHpGVksaf/+qs3l9ZIgr+qrgOu6+7/L7B+FF9XkjR6fnJXkhrjRdokHTQ/0X14cMYvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhrjZZkb52V2tRTm+77zN3MdGs74JakxBr8kNcbgl6TGGPyS1BiDX5IaM/BZPUlOBT4OPB74BbClqj6Q5ATgU8AkcCfwiqp6YPhSJR3uPNvn0Bhmxr8P+Nuq+h3gmcAbkpwObAZ2VNU6YEe3LUlaJgYO/qraU1U3d/d/AuwETgE2AFu7w7YC5w5bpCRpdEbyAa4kk8AzgBuAE6tqD/ReHJKsmecxm4BNAE94whNGUcZhx//2Sgfmz8hghn5zN8njgM8Ab6mq/zvYx1XVlqqaqqqpiYmJYcuQJB2koYI/ya/TC/1PVtVnu+F7k5zU7T8J2DtciZKkURo4+JME+Aiws6re17drO7Cxu78RuHrw8iRJozbMGv9zgFcD30ryjW7s74FLgW1JLgB2A+cNV6IkaZQGDv6q+i8g8+xeP+jX1XC82qakhXhZ5hXIcJc0DC/ZIEmNMfglqTEu9Uha9lzeHC1n/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5Ia4yUbDiF/P6h0aPizdmDO+CWpMQa/JDXGpZ5lwCsPSjqUnPFLUmOc8UtqxiD/u57vDeGV/AaywT+ElfwPL+ngHI5LsS71SFJjxjbjT3I28AHgCODDVXXpuJ5r3A7HV3xJ7RpL8Cc5AvgX4IXADPC1JNur6vZxPJ9LLpKWu+WUU+Na6jkT2FVV36+qnwNXAhvG9FySpEVIVY3+iyZ/BpxdVa/rtl8N/GFVvbHvmE3Apm7zKcB3Rl7IwVkN3LdEzz1K9rH8HC692Mfy0t/Hb1XVxGK/wLjW+DPH2K+8wlTVFmDLmJ7/oCWZrqqppa5jWPax/BwuvdjH8jKKPsa11DMDnNq3vRa4Z0zPJUlahHEF/9eAdUlOS3IkcD6wfUzPJUlahLEs9VTVviRvBL5E73TOy6rqtnE81wgs+XLTiNjH8nO49GIfy8vQfYzlzV1J0vLlJ3clqTEGvyQ1pqngT3JckquSfDvJziTPSnJCkmuTfLe7PX6p61xIkguT3Jbk1iRXJHnMSukjyWVJ9ia5tW9s3tqTXJRkV5LvJHnx0lT9aPP08Z7ue+ubST6X5Li+fSumj759b01SSVb3jS3LPmD+XpK8qav3tiTv7htflr3M87319CTXJ/lGkukkZ/btW3wfVdXMH2Ar8Lru/pHAccC7gc3d2GbgXUtd5wI9nALcATy2294G/OVK6QN4HnAGcGvf2Jy1A6cDtwBHAacB3wOOWOoeDtDHi4BV3f13rdQ+uvFT6Z2c8QNg9XLv4wD/Jn8CfBk4qttes9x7maePa4CXdPdfClw3TB/NzPiTHEvvL/QjAFX186r6Eb1LSWztDtsKnLs0FS7KKuCxSVYBR9P7jMSK6KOqvgrcv9/wfLVvAK6sqger6g5gF73LgSy5ufqoqmuqal+3eT29z6/ACuuj837gbfzqBy+XbR8wby+vBy6tqge7Y/Z248u2l3n6KODY7v5v8MvPRQ3URzPBDzwRmAU+muTrST6c5BjgxKraA9DdrlnKIhdSVXcD7wV2A3uAH1fVNaywPvYzX+2nAHf1HTfTja0ErwW+0N1fUX0keTlwd1Xdst+uFdVH58nAc5PckOQrSf6gG19pvbwFeE+Su+j9/F/UjQ/UR0vBv4ref58+VFXPAH5Kb1lhRenWvzfQ+2/dycAxSV61tFWNzYKX/liOklwM7AM++cjQHIctyz6SHA1cDLxjrt1zjC3LPvqsAo4Hngn8HbAtSVh5vbweuLCqTgUupFu5YMA+Wgr+GWCmqm7otq+i90Jwb5KTALrbvfM8frl4AXBHVc1W1UPAZ4Fns/L66Ddf7Svu0h9JNgLnAH9R3SIsK6uPJ9GbVNyS5E56td6c5PGsrD4eMQN8tnpuBH5B7yJnK62XjfR+1gE+zS+Xcwbqo5ngr6ofAncleUo3tB64nd6lJDZ2YxuBq5egvMXYDTwzydHdzGU9sJOV10e/+WrfDpyf5KgkpwHrgBuXoL6D0v3yobcDL6+qn/XtWjF9VNW3qmpNVU1W1SS9YDmj+/lZMX30+Xfg+QBJnkzvpI77WHm93AP8cXf/+cB3u/uD9bHU72Af4nfLnw5MA9+k9w1xPPCbwI7uL3IHcMJS13kQffwD8G3gVuDf6L2jvyL6AK6g997EQ/RC5YID1U5v2eF79C7b/ZKlrn+BPnbRW2/9RvfnX1diH/vtv5PurJ7l3McB/k2OBD7R/azcDDx/ufcyTx9/BNxE7wyeG4DfH6YPL9kgSY1pZqlHktRj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TG/D9lxU1dZ4KQzgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARdklEQVR4nO3da7BdZX3H8e+vpHKzFNIEGgkY7KDWOlO1R4qXKjWiKEh4USodcYLiZLReGasE7ZTpC6fxMm190dZJFZsqA6WUlnS8EWOt43QAA4ISAgNKCoFIDrVYLzNg5N8XezFs4znknL33ueznfD8zZ/bez1p7718g55fnPHutdVJVSJLa8ksLHUCSNHqWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5a4lKcnuJK8c8jUuSPL1UWWSRslyl6QGWe5acpJ8BjgR+PckP0ry/iSnJvmvJA8nuTXJaX37X5Dku0l+mOSeJG9I8pvAJ4AXda/x8AL9caQpxcsPaClKsht4S1V9OcnxwLeANwJfBNYCVwLPBn4C7AVeWFV3JlkFLK+qnUku6F7jpQvxZ5CejDN3Cc4HPl9Vn6+qx6pqG7ADeG23/THguUkOr6q9VbVzwZJKM2S5S/B04NxuSebhbonlpcCqqvox8HrgrcDeJJ9L8uyFDCvNhOWupap/PfI+4DNVdXTf15FVtQmgqr5UVacDq4A7gL+f4jWkRcVy11L1IPCM7v5ngdcleXWSQ5IcluS0JKuTHJfk7CRHAo8APwJ+1vcaq5M8Zf7jS0/OctdS9RfAn3ZLMK8H1gEfACbpzeTfR+/745eA9wIPAN8HXg78cfcaXwF2At9L8tC8ppcOwqNlJKlBztwlqUGWuyQ1yHKXpAZZ7pLUoGUH2yHJZcBZwL6qem439lHgdcCjwHeAN1XVw922S4AL6R0u9q6q+tLB3mPFihW1Zs2aQf8MkrQk3XTTTQ9V1cqpth30aJkkL6N3bO8/9pX7q4CvVNX+JB8GqKqLkzwHuAI4BXga8GXgmVX1s6lfvWdiYqJ27Ngxyz+WJC1tSW6qqompth10Waaqvkbv+N7+seuqan/38HpgdXd/HXBlVT1SVfcAd9MreknSPBrFmvubgS9094+ndwLI4/Z0Y5KkeTRUuSf5ILAfuPzxoSl2m3LdJ8mGJDuS7JicnBwmhiTpAAOXe5L19D5ofUM9sXC/Bzihb7fV9E7b/gVVtbmqJqpqYuXKKT8PkCQNaKByT3IGcDFwdlX9pG/TVuC8JIcmOQk4Gbhx+JiSpNmYyaGQVwCnASuS7AEuBS4BDgW2JQG4vqre2v12mquA2+kt17z9YEfKSJJGb1FcOMxDISVp9oY6FFKSNH4sd0lq0EHX3KWlYs3Gz005vnvTmfOcRBqeM3dJapDlLkkNcllGGjGXd7QYOHOXpAY5c5cWmDN9zQVn7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIa5OUH1KTpTukHT+vX0uDMXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktSgg5Z7ksuS7EtyW9/Y8iTbktzV3R7Tt+2SJHcnuTPJq+cquCRpejOZuf8DcMYBYxuB7VV1MrC9e0yS5wDnAb/VPedvkxwysrSSpBk5aLlX1deA7x8wvA7Y0t3fApzTN35lVT1SVfcAdwOnjCirJGmGBr22zHFVtRegqvYmObYbPx64vm+/Pd3YL0iyAdgAcOKJJw4YQ0vFdNeKmY/rxCzke0uDGvUHqplirKbasao2V9VEVU2sXLlyxDEkaWkbtNwfTLIKoLvd143vAU7o22818MDg8SRJgxi03LcC67v764Fr+8bPS3JokpOAk4Ebh4soSZqtg665J7kCOA1YkWQPcCmwCbgqyYXAvcC5AFW1M8lVwO3AfuDtVfWzOcouSZrGQcu9qv5omk1rp9n/Q8CHhgklSRqOZ6hKUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDBr3kr7TkTXcpYGkxcOYuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGuShkNI88dBJzSdn7pLUIMtdkhrksoyWHJdHtBQ4c5ekBlnuktQgy12SGuSauzRmpvvMYPemM+c5iRYzZ+6S1CDLXZIaZLlLUoOGKvckFyXZmeS2JFckOSzJ8iTbktzV3R4zqrCSpJkZuNyTHA+8C5ioqucChwDnARuB7VV1MrC9eyxJmkfDLsssAw5Psgw4AngAWAds6bZvAc4Z8j0kSbM0cLlX1f3Ax4B7gb3AD6rqOuC4qtrb7bMXOHaq5yfZkGRHkh2Tk5ODxpAkTWGYZZlj6M3STwKeBhyZ5PyZPr+qNlfVRFVNrFy5ctAYkqQpDLMs80rgnqqarKqfAtcALwYeTLIKoLvdN3xMSdJsDFPu9wKnJjkiSYC1wC5gK7C+22c9cO1wESVJszXw5Qeq6oYkVwM3A/uBbwKbgacCVyW5kN4/AOeOIqgkaeaGurZMVV0KXHrA8CP0ZvGSpAXiGaqS1CDLXZIaZLlLUoMsd0lqkOUuSQ3yNzFJi9R0v3FJmgln7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ3yqpBS457s6pK7N505j0k0n5y5S1KDnLlrQUw3m3QmKY2GM3dJapDlLkkNcllGI7FQyyz+Kjppas7cJalBlrskNWiock9ydJKrk9yRZFeSFyVZnmRbkru622NGFVaSNDPDztw/Dnyxqp4N/DawC9gIbK+qk4Ht3WNJ0jwauNyTHAW8DPgUQFU9WlUPA+uALd1uW4Bzhg0pSZqdYWbuzwAmgU8n+WaSTyY5EjiuqvYCdLfHTvXkJBuS7EiyY3JycogYkqQDDXMo5DLgBcA7q+qGJB9nFkswVbUZ2AwwMTFRQ+RQQzy0URqNYWbue4A9VXVD9/hqemX/YJJVAN3tvuEiSpJma+CZe1V9L8l9SZ5VVXcCa4Hbu6/1wKbu9tqRJJU0crP9Sclr/4yPYc9QfSdweZKnAN8F3kTvp4GrklwI3AucO+R7SJJmaahyr6pbgIkpNq0d5nUlScPxDFVJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGDXvJX+lJ+ZuVpIXhzF2SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yEMhNaXpDmHcvenMeU6imfKwU/Vz5i5JDbLcJalBlrskNchyl6QGWe6S1KChj5ZJcgiwA7i/qs5Kshz4J2ANsBv4w6r632HfR9Li5dFVi88oZu7vBnb1Pd4IbK+qk4Ht3WNJ0jwaqtyTrAbOBD7ZN7wO2NLd3wKcM8x7SJJmb9iZ+18D7wce6xs7rqr2AnS3x071xCQbkuxIsmNycnLIGJKkfgOXe5KzgH1VddMgz6+qzVU1UVUTK1euHDSGJGkKw3yg+hLg7CSvBQ4DjkryWeDBJKuqam+SVcC+UQTV4uAp7kub///Hx8Az96q6pKpWV9Ua4DzgK1V1PrAVWN/tth64duiUkqRZmYvj3DcBpye5Czi9eyxJmkcjuSpkVX0V+Gp3/3+AtaN4XUnSYDxDVZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1aCSX/NX48jfraC5N9/dr96Yz5znJ0uPMXZIa5Mx9iXCGLi0tztwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgzwUUtK88+SmuefMXZIaZLlLUoMsd0lq0MDlnuSEJP+RZFeSnUne3Y0vT7ItyV3d7TGjiytJmolhPlDdD7y3qm5O8ivATUm2ARcA26tqU5KNwEbg4uGjSmqdH7SOzsAz96raW1U3d/d/COwCjgfWAVu63bYA5wwbUpI0OyM5FDLJGuD5wA3AcVW1F3r/ACQ5dprnbAA2AJx44omjiCG8+qOknqE/UE3yVOBfgPdU1f/N9HlVtbmqJqpqYuXKlcPGkCT1Garck/wyvWK/vKqu6YYfTLKq274K2DdcREnSbA1ztEyATwG7quov+zZtBdZ399cD1w4eT5I0iGHW3F8CvBH4dpJburEPAJuAq5JcCNwLnDtcREnSbA1c7lX1dSDTbF476OtKkobnGaqS1CDLXZIaZLlLUoMsd0lqkOUuSQ3yNzEtYl5ESdKgnLlLUoMsd0lqkOUuSQ2y3CWpQX6gOoa8Zrukg3HmLkkNstwlqUGWuyQ1yDV3SYueJ/TNnjN3SWqQM/d55OxDGq3Zfk8tpe9BZ+6S1CDLXZIa5LLMEJbSj3iSxoszd0lqkDN3Sc3xEh3O3CWpSc7cFwFnGZJGzZm7JDXIcpekBrksMwOzXTZxmUVqx7ge8jxnM/ckZyS5M8ndSTbO1ftIkn7RnMzckxwC/A1wOrAH+EaSrVV1+1y836hmyov9X2JJi99s+2iuemeuZu6nAHdX1Xer6lHgSmDdHL2XJOkAqarRv2jyB8AZVfWW7vEbgd+tqnf07bMB2NA9fBZw58iDDGYF8NBChxiC+RfWOOcf5+ywNPM/vapWTrVhrj5QzRRjP/evSFVtBjbP0fsPLMmOqppY6ByDMv/CGuf845wdzH+guVqW2QOc0Pd4NfDAHL2XJOkAc1Xu3wBOTnJSkqcA5wFb5+i9JEkHmJNlmaran+QdwJeAQ4DLqmrnXLzXHFh0S0WzZP6FNc75xzk7mP/nzMkHqpKkheXlBySpQZa7JDVoSZd7kqOTXJ3kjiS7krwoyfIk25Lc1d0es9A5p5PkoiQ7k9yW5Iokhy3m/EkuS7IvyW19Y9PmTXJJd/mKO5O8emFSP2Ga/B/t/v58K8m/Jjm6b9uiz9+37U+SVJIVfWNjkT/JO7uMO5N8pG980edP8rwk1ye5JcmOJKf0bRsuf1Ut2S9gC/CW7v5TgKOBjwAbu7GNwIcXOuc02Y8H7gEO7x5fBVywmPMDLwNeANzWNzZlXuA5wK3AocBJwHeAQxZh/lcBy7r7Hx63/N34CfQOfvhvYMU45Qd+H/gycGj3+Ngxy38d8Jru/muBr44q/5KduSc5it5/7E8BVNWjVfUwvcskbOl22wKcszAJZ2QZcHiSZcAR9M4lWLT5q+prwPcPGJ4u7zrgyqp6pKruAe6md1mLBTNV/qq6rqr2dw+vp3dOB4xJ/s5fAe/n5080HJf8bwM2VdUj3T77uvFxyV/AUd39X+WJ84GGzr9kyx14BjAJfDrJN5N8MsmRwHFVtReguz12IUNOp6ruBz4G3AvsBX5QVdcxJvn7TJf3eOC+vv32dGOL2ZuBL3T3xyJ/krOB+6vq1gM2jUV+4JnA7yW5Icl/JnlhNz4u+d8DfDTJffS+ny/pxofOv5TLfRm9H5H+rqqeD/yY3rLAWOjWptfR+5HtacCRSc5f2FQjddBLWCwmST4I7Acuf3xoit0WVf4kRwAfBP5sqs1TjC2q/J1lwDHAqcD7gKuShPHJ/zbgoqo6AbiIbiWBEeRfyuW+B9hTVTd0j6+mV/YPJlkF0N3um+b5C+2VwD1VNVlVPwWuAV7M+OR/3HR5x+YSFknWA2cBb6huwZTxyP8b9CYHtybZTS/jzUl+nfHID72c11TPjcBj9C7ANS7519P73gX4Z55Yehk6/5It96r6HnBfkmd1Q2uB2+ldJmF9N7YeuHYB4s3EvcCpSY7oZiprgV2MT/7HTZd3K3BekkOTnAScDNy4APmeVJIzgIuBs6vqJ32bFn3+qvp2VR1bVWuqag29QnlB972x6PN3/g14BUCSZ9I7MOIhxif/A8DLu/uvAO7q7g+ffyE/PV7oL+B5wA7gW/T+khwD/BqwvfuPvB1YvtA5nyT/nwN3ALcBn6H3yfqizQ9cQe/zgZ/SK5ILnywvvSWD79C7HPRrFmn+u+mtjd7SfX1inPIfsH033dEy45KfXpl/tvseuBl4xZjlfylwE70jY24AfmdU+b38gCQ1aMkuy0hSyyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KD/By9h8tAJ79sWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time (sec): 158.54225397109985\n"
     ]
    }
   ],
   "source": [
    "# Generate and save SBM graphs\n",
    "\n",
    "class DotDict(dict):\n",
    "    def __init__(self, **kwds):\n",
    "        self.update(kwds)\n",
    "        self.__dict__ = self\n",
    "\n",
    "def plot_histo_graphs(dataset, title):\n",
    "    # histogram of graph sizes\n",
    "    graph_sizes = []\n",
    "    for graph in dataset:\n",
    "        graph_sizes.append(graph.nb_nodes)\n",
    "    plt.figure(1)\n",
    "    plt.hist(graph_sizes, bins=50)\n",
    "    plt.title(title)\n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "\n",
    "# configuration for 100 patterns 100/20 \n",
    "nb_pattern_instances = 100 # nb of patterns\n",
    "nb_train_graphs_per_pattern_instance = 100 # train per pattern\n",
    "nb_test_graphs_per_pattern_instance = 20 # test, val per pattern\n",
    "SBM_parameters = {}\n",
    "SBM_parameters['nb_clusters'] = 5 \n",
    "SBM_parameters['size_min'] = 5 \n",
    "SBM_parameters['size_max'] = 35 \n",
    "SBM_parameters['p'] = 0.5 \n",
    "SBM_parameters['q'] = 0.2 \n",
    "SBM_parameters['p_pattern'] = 0.5 \n",
    "SBM_parameters['q_pattern'] = 0.5  \n",
    "SBM_parameters['vocab_size'] = 3 \n",
    "SBM_parameters['size_subgraph'] = 20 \n",
    "print(SBM_parameters)\n",
    "    \n",
    "\n",
    "dataset_train = []\n",
    "dataset_val = []\n",
    "dataset_test = []\n",
    "for idx in range(nb_pattern_instances):\n",
    "    \n",
    "    print('pattern:',idx)\n",
    "    \n",
    "    SBM_parameters['W0'] = random_pattern(SBM_parameters['size_subgraph'],SBM_parameters['p'])\n",
    "    SBM_parameters['u0'] = np.random.randint(SBM_parameters['vocab_size'],size=SBM_parameters['size_subgraph'])\n",
    "\n",
    "    for _ in range(nb_train_graphs_per_pattern_instance):\n",
    "        data = generate_SBM_graph(SBM_parameters)\n",
    "        graph = DotDict()\n",
    "        graph.nb_nodes = data.nb_nodes\n",
    "        graph.W = data.W\n",
    "        graph.rand_idx = data.rand_idx\n",
    "        graph.node_feat = data.node_feat\n",
    "        graph.node_label = data.node_label\n",
    "        dataset_train.append(graph)\n",
    "\n",
    "    for _ in range(nb_test_graphs_per_pattern_instance):\n",
    "        data = generate_SBM_graph(SBM_parameters)\n",
    "        graph = DotDict()\n",
    "        graph.nb_nodes = data.nb_nodes\n",
    "        graph.W = data.W\n",
    "        graph.rand_idx = data.rand_idx\n",
    "        graph.node_feat = data.node_feat\n",
    "        graph.node_label = data.node_label\n",
    "        dataset_val.append(graph)\n",
    "\n",
    "    for _ in range(nb_test_graphs_per_pattern_instance):\n",
    "        data = generate_SBM_graph(SBM_parameters)\n",
    "        graph = DotDict()\n",
    "        graph.nb_nodes = data.nb_nodes\n",
    "        graph.W = data.W\n",
    "        graph.rand_idx = data.rand_idx\n",
    "        graph.node_feat = data.node_feat\n",
    "        graph.node_label = data.node_label\n",
    "        dataset_test.append(graph)\n",
    "\n",
    "\n",
    "print(len(dataset_train),len(dataset_val),len(dataset_test))\n",
    "\n",
    "\n",
    "plot_histo_graphs(dataset_train,'train')\n",
    "plot_histo_graphs(dataset_val,'val')\n",
    "plot_histo_graphs(dataset_test,'test')\n",
    "\n",
    "\n",
    "with open('SBM_PATTERN_train.pkl',\"wb\") as f:\n",
    "    pickle.dump(dataset_train,f)\n",
    "with open('SBM_PATTERN_val.pkl',\"wb\") as f:\n",
    "    pickle.dump(dataset_val,f)\n",
    "with open('SBM_PATTERN_test.pkl',\"wb\") as f:\n",
    "    pickle.dump(dataset_test,f)\n",
    "    \n",
    "    \n",
    "print('Time (sec):',time.time() - start) # 163s\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert to DGL format and save with pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/xbresson/Documents/Dropbox/06_NTU_2017_now/03_my_codes/34_benchmark20/GITHUB_benchmark_project/benchmarking-gnn\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "os.chdir('../../') # go to root folder of the project\n",
    "print(os.getcwd())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import pickle\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from data.SBMs import SBMsDatasetDGL \n",
    "\n",
    "from data.data import LoadData\n",
    "from torch.utils.data import DataLoader\n",
    "from data.SBMs import SBMsDataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[I] Loading data ...\n",
      "preparing 10000 graphs for the TRAIN set...\n",
      "preparing 2000 graphs for the TEST set...\n",
      "preparing 2000 graphs for the VAL set...\n",
      "[I] Finished loading.\n",
      "[I] Data load time: 5083.2611s\n",
      "Time (sec): 5083.27174282074\n"
     ]
    }
   ],
   "source": [
    "DATASET_NAME = 'SBM_PATTERN'\n",
    "dataset = SBMsDatasetDGL(DATASET_NAME) # 4424s = 73min\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000\n",
      "2000\n",
      "2000\n",
      "(DGLGraph(num_nodes=119, num_edges=4842,\n",
      "         ndata_schemes={'feat': Scheme(shape=(), dtype=torch.int64)}\n",
      "         edata_schemes={'feat': Scheme(shape=(1,), dtype=torch.float32)}), tensor([0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0,\n",
      "        0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0,\n",
      "        1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0],\n",
      "       dtype=torch.int16))\n",
      "(DGLGraph(num_nodes=134, num_edges=5956,\n",
      "         ndata_schemes={'feat': Scheme(shape=(), dtype=torch.int64)}\n",
      "         edata_schemes={'feat': Scheme(shape=(1,), dtype=torch.float32)}), tensor([0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1,\n",
      "        0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,\n",
      "        0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], dtype=torch.int16))\n",
      "(DGLGraph(num_nodes=118, num_edges=4686,\n",
      "         ndata_schemes={'feat': Scheme(shape=(), dtype=torch.int64)}\n",
      "         edata_schemes={'feat': Scheme(shape=(1,), dtype=torch.float32)}), tensor([1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,\n",
      "        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,\n",
      "        0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0,\n",
      "        0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0,\n",
      "        0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1],\n",
      "       dtype=torch.int16))\n"
     ]
    }
   ],
   "source": [
    "print(len(dataset.train))\n",
    "print(len(dataset.val))\n",
    "print(len(dataset.test))\n",
    "\n",
    "print(dataset.train[0])\n",
    "print(dataset.val[0])\n",
    "print(dataset.test[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time (sec): 21.121261835098267\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "\n",
    "with open('data/SBMs/SBM_PATTERN.pkl','wb') as f:\n",
    "        pickle.dump([dataset.train,dataset.val,dataset.test],f)\n",
    "        \n",
    "print('Time (sec):',time.time() - start) # 21s\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test load function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[I] Loading dataset SBM_PATTERN...\n",
      "train, test, val sizes : 10000 2000 2000\n",
      "[I] Finished loading.\n",
      "[I] Data load time: 30.4815s\n"
     ]
    }
   ],
   "source": [
    "DATASET_NAME = 'SBM_PATTERN'\n",
    "dataset = LoadData(DATASET_NAME) # 30s\n",
    "trainset, valset, testset = dataset.train, dataset.val, dataset.test\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time (sec): 0.0003371238708496094\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "\n",
    "batch_size = 10\n",
    "collate = SBMsDataset.collate\n",
    "train_loader = DataLoader(trainset, batch_size=batch_size, shuffle=True, collate_fn=collate)\n",
    "\n",
    "print('Time (sec):',time.time() - start) #0.0006"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}