code/hosgns_4Way_Dyn_OpenABM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "%pylab inline\n",
    "import pickle\n",
    "\n",
    "import itertools\n",
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "import networkx as nx\n",
    "import umap\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.config.threading.set_intra_op_parallelism_threads(36)\n",
    "tf.config.threading.set_inter_op_parallelism_threads(36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %env CUDA_VISIBLE_DEVICES=0\n",
    "# gpus = tf.config.experimental.list_physical_devices('GPU')\n",
    "# if gpus:\n",
    "#     try:\n",
    "#     # Currently, memory growth needs to be the same across GPUs\n",
    "#         for gpu in gpus:\n",
    "#             tf.config.experimental.set_memory_growth(gpu, True)\n",
    "#         logical_gpus = tf.config.experimental.list_logical_devices('GPU')\n",
    "#         print(len(gpus), \"Physical GPUs,\", len(logical_gpus), \"Logical GPUs\")\n",
    "#     except RuntimeError as e:\n",
    "#     # Memory growth must be set before GPUs have been initialized\n",
    "#         print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hosgns_walks\n",
    "from utils import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocess dataset\n",
    "dataset = 'OpenABM-Covid19-Interactions-5k-20'\n",
    "AGGR_TIME = 1\n",
    "partial_times, s_temp_net, df_tnet = load_temp_data(dataset, AGGR_TIME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename nodes and time slices\n",
    "map_node_index = {node:index for index, node in \\\n",
    "                  enumerate(np.unique(s_temp_net.reset_index()[['i', 'j']].values.flatten()))}\n",
    "map_time_index = {time:index for index, time in \\\n",
    "                  enumerate(np.unique(s_temp_net.reset_index().tslice.values))}\n",
    "NR_NODES = len(map_node_index)\n",
    "NR_TIMES = len(map_time_index)\n",
    "\n",
    "df_tnet['tslice'] = df_tnet['tslice'].map(map_time_index) \n",
    "df_tnet['i'] = df_tnet['i'].map(map_node_index) \n",
    "df_tnet['j'] = df_tnet['j'].map(map_node_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load network and default list of active nodes\n",
    "supra_G = pickle.load(open('../preprocessed/SupraAdjacencyMatrix/%s/Dyn/Aggtime%dWeight1ReptimeNon.gpickle' % (dataset, AGGR_TIME), 'rb'))\n",
    "pat_active_time = pickle.load(open('../preprocessed/SupraAdjacencyMatrix/%s/PatActiveTimeAggtime%d.pkl'% (dataset,AGGR_TIME), 'rb'))\n",
    "node_active_time = np.array([str(map_node_index[int(n_node.split('-')[0])])+'-'+\\\n",
    "                                 str(map_time_index[int(n_node.split('-')[1])]) for n_node in pat_active_time])\n",
    "times = [int(x.split('-')[1]) for x in node_active_time]\n",
    "node_tuple_time = [(int(n_node.split('-')[0]), int(n_node.split('-')[1])) for n_node in node_active_time]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import f1_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.utils.extmath import cartesian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TENSOR FACTORIZATION of 4WAY P(Dyn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1280
    },
    "colab_type": "code",
    "id": "OAUpjCOqXCK1",
    "outputId": "5d16a443-9628-47f3-bab4-1a1d6c74236e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampling Random Walks...\n",
      "Warmup...\n",
      "WARNING:tensorflow:From /home/spiaggesi/anaconda3/envs/tf2/lib/python3.6/site-packages/tensorflow/python/ops/array_ops.py:5049: calling gather (from tensorflow.python.ops.array_ops) with validate_indices is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "The `validate_indices` argument has no effect. Indices are always validated on CPU and never validated on GPU.\n",
      "Training...\n",
      "step    0 - loss: 0.4505613148212433 (29.9376 seconds)\n",
      "step  200 - loss: 0.4140888750553131 (42.3284 seconds)\n",
      "step  400 - loss: 0.3824400007724762 (40.7112 seconds)\n",
      "step  600 - loss: 0.3675401210784912 (39.9300 seconds)\n",
      "step  800 - loss: 0.35868701338768005 (40.7223 seconds)\n",
      "step 1000 - loss: 0.3525395393371582 (41.5090 seconds)\n",
      "step 1200 - loss: 0.3481978178024292 (40.0792 seconds)\n",
      "step 1400 - loss: 0.3445282578468323 (39.4327 seconds)\n",
      "step 1600 - loss: 0.3416595458984375 (40.3869 seconds)\n",
      "step 1800 - loss: 0.3390478789806366 (39.2934 seconds)\n",
      "step 2000 - loss: 0.3369506299495697 (42.0722 seconds)\n",
      "step 2200 - loss: 0.33506134152412415 (39.0037 seconds)\n",
      "step 2400 - loss: 0.3333771228790283 (39.9336 seconds)\n",
      "step 2600 - loss: 0.33187270164489746 (41.5047 seconds)\n",
      "step 2800 - loss: 0.33052903413772583 (38.5487 seconds)\n",
      "step 3000 - loss: 0.32932376861572266 (39.5251 seconds)\n",
      "step 3200 - loss: 0.32812681794166565 (38.6274 seconds)\n",
      "step 3400 - loss: 0.32699280977249146 (38.8175 seconds)\n",
      "step 3600 - loss: 0.32599109411239624 (38.0392 seconds)\n",
      "step 3800 - loss: 0.32502201199531555 (39.1147 seconds)\n",
      "step 4000 - loss: 0.3241663873195648 (39.5100 seconds)\n",
      "\n",
      "Total time: 829.0316 seconds.\n"
     ]
    }
   ],
   "source": [
    "# train 4way-hosgns \n",
    "print_every=200\n",
    "solver = hosgns_walks.HOSGNSSolver(network=supra_G.to_undirected(), \n",
    "                         window_size=10,\n",
    "                         emb_dim=128,\n",
    "                         order=4,\n",
    "                         iters=4000,\n",
    "                         batch_size=20000, \n",
    "                         negative_samples=5,\n",
    "                         learning_rate=tf.keras.optimizers.schedules.PolynomialDecay(0.05, 4000),\n",
    "                         warmup_steps=1000,\n",
    "                         random_state=42)\n",
    "loss = solver.train(print_loss='sg', print_every=print_every)\n",
    "emb_W, emb_T, emb_C, emb_S = solver.model.factors\n",
    "\n",
    "#save embeddings\n",
    "# pickle.dump(list(solver.model.factors), open( '../output/%s_Dyn_Aggtime%dWeight1ReptimeNonWindow%dWay%s.pkl' % (dataset, AGGR_TIME, 10, '4'), \"wb\" ), protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load embeddings\n",
    "emb_W, emb_T, emb_C, emb_S = pickle.load(open( '../output/%s_Dyn_Aggtime%dWeight1ReptimeNonWindow%dWay%s.pkl' % (dataset, AGGR_TIME, 10, '4'), \"rb\" ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, '$ -\\\\sum_{ijkll}[ p_{ijkl} \\\\log \\\\sigma (m_{ijk}) + \\\\kappa p_ip_jp_kp_l \\\\log \\\\sigma (-m_{ijkl})]$')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAE3CAYAAABrZHykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9nUlEQVR4nO3deZgU1dXH8e+PEYiACgoCggIqLgQRHUSJG+gQcWOi0bglqHndYjSaxYVojJoYY9zjEmPUaBIVjSuuiMoYjYKAOwIRFQUX3JdBBcHz/nGrmaLpnunq6ema5Xyep57uvrWdqanpM3Vv1b0yM5xzzrlitEs7AOeccy2XJxHnnHNF8yTinHOuaJ5EnHPOFc2TiHPOuaJ5EnHOOVc0TyLOOeeK5knEOedc0VZLO4CWTNKawLpFrPqhmX1c6nicc67cPIk0zsHAX4pY7zfA70sci3POlZ1XZzXO12kH4JxzaZL3nVU8SR2AicBugAGHm9kN6UblnHPl40mkkSStDjwA7AwsBw42s1vTjco558rDk0gJSOoCTAa2I1Rxfd/M7kk3Kueca3qeREpEUlfgUWArYAmwt5lNTjUo55xrYp5ESkhSd+AxYBDwBTDGzB5PNyrnnGs6fndWCZnZB0AVMA/oBNwnaZuG1pN0vSSLppomDjNfDGdKajb/UUi6TFKqVYKSfi7pBUn+d1JGfi42jXzns6TTY98/C5Nu1/84SszM3gF2Bd4A1gAmSRpSwKrvAiOAY5swvBZB0kbA0cBZKYdyFeFh0kNTjsOlpBmdi6WQ73z+O+G75/5iNupJpAmY2ZuERPIx0A2YLKl3A6stMbOpZvZykwfY/J0IPG9mM9IMwsy+BP4B/CrJepLmSzqzSYJy5XYizeBcLIV857OZvWVmU4H3i9muJ5GmsyPQNXp/b3SFkpikpyXdKulsSa9K+iq6JN21ZJHm3/cYSU9J+lLSp5LukrRpjuUOkjQniu1FSWMl1WRXzUnqLOk8SfMkLY1dQmemX0rqCPwQuCnHfp6Pqv6OlPRyFNeTkjaStFZU7bBI0seSLpekEhyGCcAgSd8pwbbykvSYpMmxzwOjY3qPpDWisqLOhXKvV+ptRNtpVudisVI6h+NKfz6bmU8lnoBqYBnhAcTbgYoGlr8emJ+jfDXgS2Ax4c6vscA+wMvAp0D3EsZ8ZjgdVnweQ3juZXK034MJbT3vA31iy40GvgHuAvYkXCq/BrwN1MSWEzAJqAXGR+udEe3jVUI3MJsTnrcxYFhWfB2ApcCb0b72BsZFx+YhYDpwNqFN6s/RNqpLcFzaRcf67ATrzAfOTLifT4Dzo/e7AB8C5wPtGnMulHu9Eu27WZ+LjTiXUjmHCz2fCd9DCxNvs5QB+mQAI4GvohNgMtChgHWuJ3cSGRJt5zFiiSh2cu9Twriz/3BnAK8Aq8XKBhCeg7koVvYk8BLRnX5R2dZRfDWxsmOjP/DRWfu9I/oyyNwpeEq0XIes5TLbvCWr/N9R+X6xsooozvElOjaPAw/lmSfCl2V8mh99GcTL8v4jER1XAw4BfhJ9uY0rxblQ7vVKtO9mfS424jxK7Rwu5HymyCTi1VmApKocl7O5ppoGtlNJ6AalIzCN8EeytBGhVUavvzaz5bHyOdHrOtF+u0maEr3vJOlxSRWSukr6KNclcQPzOhNO+FvMbFmm3MxeB/5L+BJAUgUwDLjdorMwWu4Z4PWszR4OTLZVn52ZA3SLrb8e8FmO47ZV9HpGVnln4AUzuy1Wtjrhi/vD2M+U+BjFvB/FlcvOhD/2+NSP0MlmvOyReraf+dnGAVcSrmL+kbVM4nMh4XorjknC9eo7fo3eRnM8FxvxfVH0OZzv/I0+5zx+9RzX+s7nxLwX3+BJwuVrQ77INyOqn32QcEfWLGAPM6ttZFxbA2+b2X+zyjMnwEIAC93Kj4ref0Foj0HSMGBG/I8qZhjwTJ553Qj/Yedqx3mX8CUJ0B1oD7yXY7lFmTeSekb7+3mO5XpH28z4FuFhzWxbAW+a2dwc5f/MKtsyen0+U1DkMcr4kvBHnctMIPs27onAvcDVsbLP69n+0Oh1HeB/wOGSLjKzb2LLJD4XEq634pgkWY/we813/JJsoyWdi8V+XxR9Duc7fyP5fgf5jmt953NinkRY8UuZ0+CCeUhan1B11Z1QlfFdM/uoBKFtDbyVo/wAwgn6eLT/s4FlZna2pNOATmZ2GuHL7elomZ2BC4FTzezhaF6+O04+Jlxe98oxrxd1/+F/QPgvO9eYKj0Jdb9Q94e+0hdB9J/U7oR2o4wPCV8c2bYCnslav1cUzzM5ll0OvBBbtphjlLF29LOuwsw+J+s4SlpK+PIs9I6erYDZhC+GauBmQr3/v2LLJD4XEq4XPyYFr0f9xy/JNlrMudiI74uiz+F6zl/I8zsg/3HNez4Xw6uzGknhKfXJwPqE/3hGm9nbJdhuO8J/I/0krRYrX49Qp3u5mS2OiisJ/xFD9F9J9H4bYGZ0Ap4O7BX7ctyG0JC3imi7M4H9Y9UbSOoHfIdQx01URTED+H78kjmq1hsQ2+Qn0etmWbs6hfBH+tdY2RygvaS+WcdiCLn/0ACezVE+x8ItjRnFHKOMAUD2f4+lNBR4PIr3VkK9/pmS2kPx50LC9VYck4Tr5Tx+RWyjRZyLxSrBOZzv/IX853C+41ra87nUjTZtaSJUXc0g/Kf0MbBlkdu5nqyGdULXKUb4D+pfhOdODiXcPfI08K3Ysu8CvaP3C4D1Y+/nEBp522VtfwHQL/b5THLfEfMA4S6SgwhVLe8D68WWGx3FeRewB6Fe/zXCf3qPRsuI8MdTS/gC2Q24gtBo+dOsuPpH29s3VrZpVLZn1rK/JlQTZf9szwD/yCpLfIyi+V2jOI9I8PucT4F3ZxGqsAw4Jla2X1R2dCPPhSTrxY9J0vVWOX5FbKNFnIuN+K5o1DlMnvO3gd/BSse1ofMZvzur/BNwY3RiJJ2m5vjlzc8q+2G07GBC/foXhCudy4A1Ysv1JVSdQLiUXxS970WoSrgeuDhr272A97LKVvrDjcrGAE8R6lA/Be4GNs1xHA4m/GezhNAetA/hP6s7Y8tsQGgrqKWuKmPvPMd1GvD32OcDo2PRO2u5fwNPZJW1j+I4sTHHKLbuIYS77dZJcF7Mp/AkUhX9bNvFyhQdv4WEevnE50LCc2jFMUm4Xn3nWJJttJhzsRHfFUWfw/nO3/p+B7mOa0PnM55Eyj8BNZQwiRC7FRS4GFhQQAzVwD3R+z2B+6L3ewP3RNt8EDg5ts7ewP1NeFz6Rifqb4pc/7Doi6JTieJJfIxi6z4A/DPl8yzxuZBwvRXHJOF69Z1jSbbRZs7FImPIef7W9zvId1xznc/U3aZ+QzFJxBvWG+dAwn+KSeW646MfoWHwMcKzJltTVwdan0rq6lmz6/qfMbNlkvYDaiS9Y2b/pP6GzEQUBuW6CHiY0Fi3IXAy4T+8a4rc7D+jbRwLXFCCMIs5RkgaSrgjZnAJYmiMYs6FJOtl17EXul59xy/JNtrSuViMfOcv5PkdAANZ9YaPoeQ+n08Dfhe9z3UjRP3Syq4+rfSfQP/o5BhGqDsV4b+fRP89EW493LXMsXcA7iTUOy8ltA1NBAY3crvbAcc2QbwFHyNCFcpBKZ8bxZ4LBa8XPybF7q8UMZfgWLWoc7HIWIr+G893PhNubc58/wxJul0fT6QVkDSQcFfP08BPbOXnCxx+jHLxY9JyNOfflScR55xzRfPnRIpUTDca8XkNdBfhnHMtgjesF8ka2dVItFy+rh6cc65F8CuRIimMj3BG9P40SedEs1bqgkDSDElVsXkzcrx3zrkWyZNI8Rrb1Ujerh6cc66l8CRSvPqePTg3er+bmcV7BfUrEedcq+JJpAhRh2zfmNk7ktYlDFqzIOqRsx0wldC1wzexdXoBHc3sjfj7VH4A55wrEU8ixYlXZWVfXTwDHAFsLunk2Dp+FeKca3U8iRSnwW40CD2x/kDSj2LzPIk451oVf9iwkSQ9Sejeob7hT51zrlXyK5EiSRoo6VngRWBK2vE451wa/ErEOedc0fxKxDnnXNHaXLcn3bt3t/79+xe17uLFi+ncuXNpAyoBjysZjysZjyuZ1hjXzJkzPzCzHjlnpt0/frmnyspKK9aUKVOKXrcpeVzJeFzJeFzJtMa4CP0B5vxO9eos55xzRfMk4pxzrmieRJxzzhXNk4hzzrmieRJxzjlXNE8izjnniuZJxDnnXNE8iRTiq6/g/PPpNsM73nXOuThPIoXo0AHOP5+eDz2UdiTOOdeseBIpRLt2UFVFt5kzwTusdM65FTyJFKqqio4ffQSzZqUdiXPONRueRApVVRVeH3443Ticc64Z8SRSqA024Iu+fT2JOOdcjCeRBD6urISaGvj667RDcc65ZsGTSAIfV1bC4sUwdWraoTjnXLPgSSSBT7baKtyp5VVazjkHeBJJZFmXLrDNNp5EnHMu4kkkqaoqmDYNPv007Uiccy51nkSSqqqC5cvhscfSjsQ551LnSSSpESOgUyev0nLOOTyJJNexI+y0kycR55yjGSQRSWMkzZU0T9Kp9Sy3jaTlkvbLKq+Q9Kyke5s+2khVFcyeDQsXlm2XzjnXHKWaRCRVAFcAuwODgIMkDcqz3HnApBybOQGY3ZRxriLTBcojj5R1t84519ykfSUyHJhnZq+Z2VJgAlCdY7njgduB9+KFkvoCewLXNHWgK9liC1h3Xa/Scs61eWknkT7AgtjnhVHZCpL6APsAV+VY/xLgZOCbJoovt3btYNddQxLxruGdc23YainvXznKsr+VLwFOMbPlUt3ikvYC3jOzmZJG1rsT6SjgKICePXtSU1NTVLC1tbUr1u3Vty+bvfsu06+/nsUDBhS1vVKJx9WceFzJeFzJeFzJNFlcZpbaBIwAJsU+jwfGZy3zOjA/mmoJVVrfA84lXLnMB94FvgD+1dA+KysrrVhTpkyp+/DGG2ZgdvHFRW+vVFaKqxnxuJLxuJLxuJJpTFzADMvznZp2ddZ0YKCkAZI6AAcCE+MLmNkAM+tvZv2B24BjzewuMxtvZn2j8gOBR83sh2WLfIMNYJNNvF3EOdempZpEzGwZcBzhrqvZwK1mNkvSMZKOSTO2glRVha7hly5NOxLnnEtF2m0imNn9wP1ZZbka0TGzw/KU1wA1JQ6tYVVVcOWVoS+tHXcs++6dcy5taVdntWyjRnnX8M65Nq3gJBI9Gd4pR/kuki6VdK6kdG9TKreuXUPX8JMnpx2Jc86lIsmVyAXAR5LWyhRIOhCYTHgY8BTgaUnrlzbEZq6qCp5+2ruGd861SUmSyE7AFDOLf1v+FvgEGEd46K8r8ItSBdcijB7tXcM759qsJElkfWBe5oOkDYFNgcvM7F9mdgHwADCmtCE2c9tt513DO+farCRJZE3gs9jn7QlPlz8YK5sF9C1BXC1Hpmt4bxdxzrVBSZLIO0C84bwK+BKYGSvrAiwrQVwtS1UVzJnjXcM759qcJElkKjBW0l6SqoD9CE+Jfx1bZkPgrVIG2CKMHh1evWt451wbkySJ/CFa/m7CE+YdgHMyMyWtCYwEppUwvpZh8ODQNbxXaTnn2piCn1g3sxclbQscGhXdYmbTY4sMAR4Cbi5hfC1DdtfwytU5sXPOtT6Juj0xsxeBX+WZ9wTwRCmCapFGj4abb4ZZs8KViXPOtQGN7vZEUntJW0natBQBtVi77hpevUrLOdeGJOn25AeSbpW0dqxsI8JtvTOAlyXdISn1Th1T4V3DO+faoCRXIj8GNjOzj2JlFwIbA1OAFwjjox9euvBamKqq8OS6dw3vnGsjkiSRQYRBpIAVd2PtQRgDpAoYDsyhLSeR0aNh8eLQNbxzzrUBSZJID8IDhxkjCA3zEwCi50UmAxuVLLqWZuTIcKeWt4s459qIJEnkc2Ct2OedCd2exO/I+gpYowRxtUyZruG9XcQ510YkSSKvALtL6hiNh74/8IKZfRBbph/wXikDbHFGj/au4Z1zbUaSJHI1oVuTVwjjoW8IXJe1zLaEu7Xarqqq0DV8TU3akTjnXJMrOImY2Q3AH4FOhGqty6MJCCMcAv0Jd2q1Xd41vHOuDUn6xPqvgV/nmf0E0A1Y3NigWrSOHWHnnT2JOOfahEY/sZ5hZkvN7FMza3tdwWfzruGdc21E4iQiaQNJp0u6XdIj0VPqp0nqV0wAksZImitpnqRT61luG0nLJe0XfV5f0hRJsyXNknRCMftvElVV4dWvRpxzrVyiJCLpSGAucBawDzAK+B7wO2CupKMTbq8CuALYnfAw40GSBuVZ7jxCF/QZy4BfmtnmwHbAT3Otm4pM1/CeRJxzrVySvrN2Ba4ClhDGEdkF2Dx6/T3hGZErouUKNRyYZ2avmdlSwoOL1TmWOx64ndjtw2b2jpk9E73/nHDHWJ8E+2467dqFq5FM1/DOOddKJbkSOYnwwGGlmZ1hZjVmNjd6PQOoBGqj5QrVB1gQ+7yQrEQgqQ/hqueqfBuR1B/YiuY0IFZVFSxaBC+9lHYkzjnXZJLcnTWc0E/Wq7lmmtmrkv4NfD/BNnON3pT9r/slwClmtlw5BnuS1IVwlXKimX2WcyfSUcBRAD179qSmyGc4amtrC163Y+fOjADmXXUVC/ffv6j9NUVc5eRxJeNxJeNxJdNkcZlZQRPwJfCHBpb5A/Blgm2OACbFPo8Hxmct8zowP5pqCVVa34vmtSe0k/yi0H1WVlZasaZMmZJshU03Ndtjj6L3V6jEcZWJx5WMx5WMx5VMY+ICZlie79Qk1VlvENo/6jMKeDPBNqcDAyUNiLpSORCYGF/AzAaYWX8z6w/cBhxrZncpXJZcC8w2s4sS7LN8vGt451wrlySJ3AlsI+lKSV3jMyStJelSQpXXHYVu0MIzJccRriZmE6rLZkk6RtIxDay+PfAjYBdJz0XTHgl+nqZXVRW6hp86Ne1InHOuSSRpEzkXGAscAxwi6XlC1/C9gKGE3nvnRMsVzMzuB+7PKsvZiG5mh8XeP0HuNpXmI9M1/MMPw047pR2Nc86VXJK+sz4DvgP8DagAdiD05Ltj9PlvwPaWp3G7TeraFYYP9+dFnHOtVqKHDS10a3I00BUYQkggQ4BuZna0mX1c+hBbuKoq7xreOddqFdV3lpktM7OXzOy/0evXAJKukPRhaUNs4bxreOdcK1ayDhgjnQhXKS5jxAjvGt4512qVOom4bB06hAb2Bx7wLlCcc62OJ5Fy2HtvePVVePnltCNxzrmS8iRSDmPHhteJE+tfzjnnWhhPIuWw3nqwzTZw991pR+KccyXlSaRcqqth2jR45520I3HOuZKp94l1Sa8l3F73RsTSuo0dC6efDvfcA0cdlXY0zjlXEg11e9K/iG36LUi5DB4MAwaEdhFPIs65VqKhJDKgLFG0BVKo0vrLX6C2Frp0STsi55xrtHqTiJm9Ua5A2oSxY+GSS+Chh2DffdOOxjnnGs0b1stpxx2hWze/1dc512p4Eimn1VaDPfeEe++FZcvSjsY55xrNk0i5VVfDhx/Ck0+mHYlzzjWaJ5Fy22230J+WP3jonGsFPImU2xprwC67hCTiHTI651o4TyJpqK4OHTLOnp12JM451yieRNKw997h1au0nHMtnCeRNPTpA8OG+a2+zrkWr6En1lciqRPwf8BQoC/QPsdiZma7Nj60Vq66Gs44A959F3r1Sjsa55wrSsFXIpKGAK8BlwCHA6OBkXmmgkkaI2mupHmSTq1nuW0kLZe0X9J1m6Xq6tCwfs89aUfinHNFS1KddQnQA/gtoWPG9mbWLsdUUegGJVUAVwC7A4OAgyQNyrPcecCkpOs2W4MHQ//+3i7inGvRkiSR7YDbzez3ZvammS0vwf6HA/PM7DUzWwpMAKpzLHc8cDvwXhHrNk+ZDhkffhgWL047GuecK0qSJFILlLpDxj7AgtjnhVHZCpL6APsAVyVdt9mrroYlS0KHjM451wIlaVh/FNi2xPtXjrLsJ/AuAU4xs+XSSosXsm5YUDoKOAqgZ8+e1NTUJA4UoLa2tuh1c8a1fDnf6dKFD//6V+Z061b0dkodV6l4XMl4XMl4XMk0WVxmVtAEbAi8D5wKqND1GtjmCGBS7PN4YHzWMq8D86OpllCl9b1C1s01VVZWWrGmTJlS9Lp5HXKI2TrrmC1bVvQmmiSuEvC4kvG4kvG4kmlMXMAMy/OdWvCViJm9JmkH4EngSEnPAZ/mXtT+r8DNTgcGShoAvAUcCByctbEVA2NJuh6418zukrRaQ+u2CNXVcOONoUPGHXdMOxrnnEuk4CQiqS9wN9AtmvKNemiEZ0kaZGbLJB1HuOuqArjOzGZJOiaan90O0uC6hf48zcaYMXUdMnoScc61MEnaRC4BNgGuA24A3gYaPSiGmd0P3J9VljN5mNlhDa3b4qyxBowaFZLI+eeHu7acc66FSJJEdiG0QRzRVMG0WdXVcOyxMGcObL552tE451zBktzi2w54sakCadO8Q0bnXAuVJIlMBQY3VSBtWt++UFnpScQ51+IkSSKnASMlHdhUwbRp1dUwbVrokNE551qIJG0iexIeOLwxuntqJvlv8f1dKYJrUzK9+t57LxzhzU7OuZYhSRI5M/Z+p2jKxQBPIkltsUVdh4yeRJxzLUSSJDKqyaJw4dbesWPh6qtDh4ydO6cdkXPONSjJE+uPNWUgjlCl9ec/w+TJ8L3vpR2Nc841yIfHbU523BG6dvW7tJxzLYYnkeakfXvYY4/QuL68FMO1OOdc00qURCT1lnRFNBztl9FwtdlTo7tCadOqq+GDD+Cpp9KOxDnnGpRkjPU+wAzgaGAx0BF4E3gFWE4Y3+N54PHSh9mGjBkTrki8Sss51wIkuRI5A+gFjDGzLaOyv5vZZoSxRiYBqwP7ljbENmbNNWGXXUISsZxjbDnnXLORJInsBjxoZg9nzzCzhcD+hCRyVolia7vGjoVXXgkdMjrnXDOWJIn0AuLjdSwnJA0AzKwWmAxUlya0Nmzs2PA6cWK6cTjnXAOSJJHPgA6xzx8DfbKW+RTo0dig2jzvkNE510IkSSJvAOvHPj8P7CKpE4CkdsB3gYWlC68NGzsWpk6FRYvSjsQ55/JKkkQeAUZJah99vgFYD3hS0vnAf4FvA7eUNsQ2qro6NKzfe2/akTjnXF5Jksi1wHlAdwAz+xdwKWGMkV8C2xISyDkljrFtGjIE+vXzKi3nXLNWcBIxs1fM7DwzeydW9nOgNzAC6G1mB5vZV00QZ9sjhauRyZNDh4zOOdcMJXnYcJyk3bLLzex9M5tmZl55X2pjx8JXX4VE4pxzzVCS6qzrgDFNFYjLYaedQoeMfquvc66ZSpJE3k24vGusTIeMEyfCZ5+lHY1zzq0iSVJ4kHB3VkkTiaQxkuZGnTqemmN+taQXJD0naYakHWLzfi5plqSXJN0s6VuljK1ZOPFE+PhjOPnktCNxzrlVJEkIpwFrANdK6l6KnUuqAK4AdgcGAQdJGpS12CPAlmY2FPgxcE20bh/gZ8AwMxsMVAAHliKuZmWbbeDnP4e//hVqatKOxjnnVpJkeNybCU+kjwMOlDSfUMWV3UugmdmuBW5zODDPzF4DkDSB0G3Ky7GN1caW75y1v9WA1SV9DXQC3i74p2lJzj67buz1F16ATp3Sjsg554BkSWRk7H1HYNNoypak69k+wILY54WE501WImkf4FxgXWBPADN7S9IFhO7ovwQeMrOHEuy75ejUCa65BkaOhN/8Bi68MO2InHMOAFmK3Y1L2h/YzcyOiD7/CBhuZsfnWX4n4Awzq5LUDbgdOAD4BPg3cFv0EGT2ekcBRwH07NmzcsKECUXFW1tbS5cuXYpatxQGXnwx6917L89edhmfDaqr9Us7rnw8rmQ8rmQ8rmQaE9eoUaNmmtmwnDPNLLWJ8JDipNjn8cD4BtZ5nfDU/P7AtbHyccCVDe2zsrLSijVlypSi1y2JTz81W399s803N/vqqxXFqceVh8eVjMeVjMeVTGPiAmZYnu/UtG/ZnQ4MlDRAUgdCw/hKD0VI2liSovdbE3oS/pBQjbWdpE7R/F2B2WWNvtzWXBOuvhpmz4bf/z7taJxzruE2EUk9COOGLIgyEpJ2BnbOsfizZnZPoTs3s2WSjiOMilgBXGdmsyQdE82/Cvg+MC5qPP8SOCCKY5qk24BngGXAs8DVhe67xRozBsaNgz/+EfbbD7bcsuF1nHOuidSbRCStDcwjXMrE77gaCfw2xyqfSRpgZh8XGoCZ3Q/cn1V2Vez9eYSOH3Ot+9s8cbRuF18MkybB4YfDtGlpR+Oca8Maqs46hPBsyPgc8ww4MjadDKwZreOa0tprw5VXwrPPwgUXpB2Nc64Na6g6awwwx8yezjXTzK6Nf5Z0MLAHcHlpwnN57btvqM466yw6Xd36a/Gcc81TQ1ciQ4AnE2zvGWCL4sNxiVx+OXTuzKZ/+hMsX552NM65NqihJNKd8FR6tueAf+QoXwSs08iYXKF69oRLLmGtWbPgiivSjsY51wY1lESWA6t0amhmd5vZ4TmW70CyJ9ZdY/3wh3y47bYwfjy8/nra0Tjn2piGksh7wCYJtrdptI4rF4n//eIXUFEBRx4ZxmV3zrkyaSiJTCd0/96toQ1Fy+wC5GyEd01nybrrwvnnwyOPwLXXNryCc86VSENJ5GZCz7lXZp4ar8cVhIcSbypFYC6hI48MHTT+8pfw1ltpR+OcayPqTSJmdhfwX+AHwKOSdpHUPjNfUntJu0p6lNAR4pNmdndTBuzyaNcu9PT79dfwk594tZZzriwK6TvrB8AsQjcnk4FaSW9LeguoBR4iPME+i9ApokvLRhuFPrXuuQeK7KnYOeeSaDCJmNk7wHbAWYTxPtoDvYDe0fuFwJnACDPLdTuwK6cTToBtt4Xjj4f3/B4H51zTKqgXXzP7wszOMrN+QD/CwFHbAf3MrJ+ZnW1mi5syUFegigq47jr4/HP42c/SjsY518ol7grezBaY2XQze9rMFjS8hiu7QYPCCIi33BKG1XXOuSaS9ngirqmcckroJv4nP4G3W+fQ88659BWdRCSdIOm1UgbjSqh9e7j++lCt9Z3vwP/+l3ZEzrlWqDFXIl0J7SOuuRo6FGpq4IsvYPvtYfr0tCNyzrUyXp3V2lVWwn//C126wKhRMHly2hE551oRTyJtwcCBIZFstBHsuac/Q+KcKxlPIm3FeuvBY4/BiBFw8MFw2WVpR+ScawUak0TuAn5cojhcOXTtCg8+CNXV4RmS00/37lGcc41SdBIxs+fN7IZSBuPKYPXV4d//hiOOgHPOgaOOgmXL0o7KOddCNTTGek6SdgS2AtYCPgWeNbPHSxmYa0KrrQZXXx1GRjznHPjgA7jpppBgnHMugURXIpK2lzQXqAEuJvSndTFQI2mOpO2TBiBpjKS5kuZJOjXH/GpJL0h6TtIMSTvE5nWVdFu079mSRiTdf5slhc4a//xnuOsuGDMGPvkk7aiccy1MwVcikioJvfh+C3iMkEjeJXTGOArYCXhI0o5m9kyB26wgjEMymtCR43RJE83s5dhijwATzcwkDQFuBTaL5l0KPGhm+0nqAHQq9OdxkeOPhx49YNw42Hnn0GbSu3faUTnnWogk1VnnRMtXm9k9WfPOklQN3BYtt3uB2xwOzDOz1wAkTQCqgRVJxMxqY8t3JhrDXdKahMR1WLTcUmBpgp/HZRx4IKyzDuyzT3i6/aGHwm3BzjnXgCTVWd8B7siRQACIBqO6M1quUH2AeCeOC6OylUjaR9Ic4D7q7gjbEHgf+LukZyVdI6lzgn27uNGjYcoUqK0NT7fPnJl2RM65FkBW4C2ekj4BLjez0+tZ5hzgp2bWtcBt7g/sZmZHRJ9/BAw3s+PzLL8TcIaZVUkaBkwFtjezaZIuBT4zs9/kWO8o4CiAnj17Vk4o8mG72tpaunTpUtS6TamUca2+YAFbnnQSq332GbN+9zs+rqxsFnGVkseVjMeVTGuMa9SoUTPNbFjOmWZW0AQ8DNzTwDL3Ag8l2OYIYFLs83hgfAPrvA50J7TFzI+V7wjc19A+KysrrVhTpkwpet2mVPK43nrLbPBgs/btzS6+2GzZsqI202aOV4l4XMl4XMk0Ji5ghuX5Tk1SnXU6MFrST3LNlPRTYFdglSuBekwHBkoaEDWMHwhMzNruxpIUvd8a6AB8aGEUxQWSNo0W3ZVYW4prhPXWg//8J1Rx/fznsN128NxzaUflnGuGkjSsfxd4FLhc0onA48AioCewAzAQeBDYTdJusfXMzH6Xa4NmtkzSccAkoAK4zsxmSTommn8V8H1gnKSvgS+BA6LMCHA8cGOUgF4DDk/w87j6dOsG994Lt94anm4fNgx+8Qs480zo5DfBOeeCJEnkzNj7gdGUbXdWvTPLgJxJBMDM7gfuzyq7Kvb+POC8POs+B+Sup3ONJ8EBB4QrklNOgfPPh9tug7/8BXbbreH1nXOtXpIkMqrJonDN29prw9/+Bj/6UegmZcyY0InjxRfDuuumHZ1zLkUFJxEze6wpA3EtwE47wfPPw7nnwh/+AA88ABdeCIcdFq5anHNtjncF75Lp2DG0izz/PHz72/DjH8Muu/jwu861UZ5EXHE23zyMT3L11fDsszBkSOiLa6l3GuBcW5I3iUj6RtIySZvEPi8vYPJ+xduKdu3gyCNh9uwwRslvfgNbbQVPPpl2ZM65MqmvTeQ/hDurvsj67NzKeveGW24JnTgee2zoNuWYY1htjz3Sjsw518TyJhEzG1nfZ+dWseeeMGsWnHEGXHop291wAxx9dHjOZMCAtKNzzjWBkreJSOpY6m26FqRLF7joInj2WT7YYQe4/HLYeGPYf3+YOjXt6JxzJVZwEpF0SgHLdCCMve7auiFDmPPrX8P8+XDSSfDwwzBiROhq/rbbfEhe51qJJFcif5B0UL6Z0QBTtxG6R3Eu6NMH/vhHWLAALrsMFi0KVyUDB8Kll8Lnn6cdoXOuEZIkkScJY3es8uR6lEBuAfYCLi9RbK416dIFjjsuPE9yxx3Qty+ceGJ4PekkePPNtCN0zhUhSRLZm9DJ4R2SBmcKox52/wHsC/zNzE4obYiuVamoCCMoPv44TJsGe+wRuk/ZcMPQlcqMGWlH6JxLoOAkYmafAGMIPek+ICkzAuF1wEHAP83s6JJH6Fqv4cPh5pvh1VfDVcl998E224TuVW6/3R9cdK4FSHR3lpm9Seild01gkqRrgEOBf+PdsLti9esHF1wQ2k0uvjhUbe23XxjX5IQTfCwT55qxxLf4mtnzhDE+BhISx93AwWb2TYljc23NmmuGK5J58+D++2HXXeGqq8JT8EOHhob4999PO0rnXEzehw0ljWtg3YeA7QhjgRyiWC+uZvaPkkTn2qbVVoPddw/TRx+FKq/rrw8J5qSTYK+9Qs/Bu+8O7dunHKxzbVt93Z5cT+5uThSVZ7LGVVllRmhod67x1l4bfvrTML34ItxwA/zzn3DnnWEskx/+EA4/HAYPbnhbzrmSqy+JeBuHa1622CK0nZx7Ljz4IPz97/DnP4cn5CsrQzI56KCQeJxzZVFf31k3lDMQ5wrWvj3svXeY3n8fbropJJTjjgvjwO+9N4wdG6q7evRIO1rnWjUfT8S1bD161N3B9eyzcMwx8MQTcOih0LNn6GrlnHPCfPNOqJ0rtUYlEUljJV0i6VJJ3y9VUM4VJXMH19tvw/Tp8Nvfhj66Tj893OG1wQZw9NGs8+ST8MUXDW7OOdewepOIpL0l/UfSzjnm/R24E/gZcDxwq6TbmyZM5xJo1w6GDQtJZPp0eOcduPba8HDjTTexxWmnhXaTPfaAK64InUQ654rS0JXIWGBrYFq8UNJehIcMvwB+D5xC6BLle/V10piLpDGS5kqaJ+nUHPOrJb0g6TlJMyTtkDW/QtKzku5Nsl/XhvTqFcaCv/12+OADnr/gglDt9coroR1lwIBwd9cpp4TuWL7+Ou2InWsxGkoiw4GnzOyrrPIfE27lPdzMzjCz84Edga+AQwrdedRx4xWEp+AHAQdJGpS12CPAlmY2NNrvNVnzTwBmF7pP18Z17MjHlZVwySUhicydCxdeGNpPLroodLmyzjqhYf6yy8LQv96W4lxeDSWRXsCrOcp3Aj4BVlRfmdm7wH3AVgn2PxyYZ2avmdlSYAJQHV/AzGrNVvwVdyb27IqkvsCerJpYnCvMJpuEO7oeeQQ++CCMdXLwwWGExp/9DAYNCm0phx8e7gJbtCjtiJ1rVup7TgSgG/BRvEDSBsDawD2xL/eM1wlVYIXqAyyIfV4IbJu9kKR9gHOBdQlJI+MS4GRgjQT7dC63tdaC738/TACvvRYG05o8Ge6+Ozw1DzBkCIweHaYdd4ROnVIL2bm0adU8EJspfQjcZ2bjYmX7EK5AzjazM7OWPw842sy6FrRzaX9gNzM7Ivr8I2C4mR2fZ/mdgDPMrCpql9nDzI6VNBL4lZntlWe9o4CjAHr27Fk5YcKEQsJbRW1tLV26dClq3abkcSVTVFzLl7PGvHl0mzGDbjNnstZLL9Hu66/5pn17Ph08mI8rK/m4spLagQOxioryxVUGHlcyrTGuUaNGzTSzYTlnmlneCagBPgS6xMr+ASwHqnIsfyswu75tZi0/ApgU+zweGN/AOq8D3QlXJguB+cC7hEb+fzW0z8rKSivWlClTil63KXlcyZQkrsWLzR580OxXvzLbckuz0HJi1rmz2c47m518stntt5stXFjeuJqAx5VMa4wLmGF5vlMbqs66Efgr8JikG4BNCA3n7wJT4gtGg1PtADxVQGLLmA4MlDQAeAs4EDg4a7sbA6+amUnaGugAfGhm46OkQ+xK5IcJ9u1c8Tp1gt12CxOEtpJHH4WnngqDbV18cd1dXn36wLbbwnbbhdfKSujcOb3YnSuhhpLItYQRC3cDhhI6WPwaOMHMlmctuyuhIf7hQnduZsskHQdMAiqA68xslqRjovlXEbqdHyfpa8KAWAdEmdG55qNnz9Bv10HRHe5LloSn5KdNg6lTw+sdd4R5FRXhluJtt62bNt88tdCda4x6k4iZfSNpT8LIhd8hVG3dYWbP5Vi8O3ApMDFJAGZ2P6E7+XjZVbH35wHnNbCNGkLVm3PNQ8eOdQniZz8LZe+/D08/HRLKtGlwyy1w9dVh3pprsuVGG8Euu8DWW4cn7DfZJCQc55qxhq5EsDDY1I3RVN9yEwi36DrncunRA/bcM0wA33wD//vfiqRSMWUKXH55uIqBUGU2dGhdUtl663DLcYcOqf0IzmVrMIk455pIu3aw2WZhOvRQnqmpYeT228OcOaEzyWeeCdMNN4TkAiGBDB4cEkpm2mILv83YpcaTiHPNSfv2ISlssQWMi+6s/+YbePXVkFAyyeXOO+Ga6Bnbdu1g003r1hs8OEwDBnh1mGtynkSca+7atYOBA8N0wAGhzAwWLKhLKs89BzNmwK231q23+urw7W/XJZXBg0OS6d0bYsNZO9cYnkSca4mk0B3LBhtAdaynoNra0N/Xiy/CSy+F6cEH6562B+jWrS6hDB4cEs3mm0P37p5cXGKeRJxrTbp0gW22CVPcBx+E/sBeeqkuwdx4I3z6ad0y3bqF9plNN61rq9l0U7RsWXl/BteieBJxri3o3h123jlMGWawcGFIKHPnhmnOHJg0aaUrlx0rKmCjjVZKLCte11mn/D+La1Y8iTjXVkmw/vph2n33led9+umKxLJg0iT6ffVVSDAPPghLl9Yt1707bLxxmDbaKEyZ9z16ePVYG+BJxDm3qrXWCiNBDh/O6+uvT7+RI0P58uVhJMg5c+quXObNg//8J1SPxTuT6NKlLqFkJ5m+fcMNA67F8yTinCtcpmpro43qHprMWLIEXn893I48b154ffXV0AYzceLKI0Z26BBuQd5wQ+jfH/r1C6+Z9z17+lVMC+FJxDlXGh071rWbZFu+PLS/ZJJL5nX+/PDE/kcfrbz8t75Vl1iyE0z//mHIY9cseBJxzjW9ioqQAPr1g113XXX+55/DG2+EpDJ//srvn3km9DsW16EDw3v0CLcmx5NN5v1668Fq/vVWDn6UnXPpW2ONugcic1m8OCSWWHKpnT6dTp9/DvfdB+++u/LyFRXhhoHsq5nMa9++3gdZiXgScc41f507h84nBw1aUfRyTQ3rZhr8v/wS3nyzLsnEXx9+GN5+e+VGfyk8uZ95YHODDULSiX9eZx1vlymAJxHnXMu3+urhuZVNN809f+nS0E1MPLm8+WZd1zF3313Xe3J8m7mSy/rrh4HG+vQJV1BtPNF4EnHOtX4dOtTdVZaLWWh3ySSWN99cebr//lWrzCBcIfXpE9pg1lsP+vSh7xdfwHvv1ZX37h1uFGilPIk455wE664bpmHDci+zZAm89VZIKm+/Haa33qp7/9RT8PbbbLxkCVxxxcrrrrPOiiRD7951ySX+2qtXi2yn8STinHOF6NgxPNey4Yb5lzHjiXvuYYcBA1ZNMm+9FaYXXoBFi8Jtz9m6d8+fZHr3DommZ89mNX6MJxHnnCsViWVrrlk3tks+y5eH6rN33gkJJv6aeT9rVnifK9l07hySybrrNvzarVuT9g7gScQ558qtoiJcVfTqFYY+zuebb0IPzJmrmUWLwvTee3Wvr70GU6eGpPTNN6tuY7XVoEcPhnbvHq6CSsyTiHPONVft2tW11QwdWv+yy5eHJ//jSSb2/otFi+jaBCF6EnHOudagoiL0nNyjR87Z/6upYb0m2K13o+mcc65oqScRSWMkzZU0T9KpOeZXS3pB0nOSZkjaISpfX9IUSbMlzZJ0Qvmjd865ti3V6ixJFcAVwGhgITBd0kQzezm22CPARDMzSUOAW4HNgGXAL83sGUlrADMlTc5a1znnXBNK+0pkODDPzF4zs6XABKA6voCZ1Zqt6PSmM2BR+Ttm9kz0/nNgNtCnbJE755xDFu+UrNw7l/YDxpjZEdHnHwHbmtlxWcvtA5wLrAvsaWZPZc3vD/wHGGxmn+XYz1HAUQA9e/asnDBhQlHx1tbW0qVLl6LWbUoeVzIeVzIeVzKtMa5Ro0bNNLPcj/KbWWoTsD9wTezzj4DL6ll+J+DhrLIuwExg30L2WVlZacWaMmVK0es2JY8rGY8rGY8rmdYYFzDD8nynpl2dtRBYP/a5L/B2voXN7D/ARpK6A0hqD9wO3GhmdzRloM4551aVdhKZDgyUNEBSB+BAYGJ8AUkbS6GvZUlbAx2AD6Oya4HZZnZRmeN2zjlHym0iAJL2AC4BKoDrzOwcSccAmNlVkk4BxgFfA18CJ5nZE9Gtvo8DLwKZZ/1/bWb3N7C/94E3igy3O/BBkes2JY8rGY8rGY8rmdYYVz8zy/kUY+pJpCWRNMPyNS6lyONKxuNKxuNKpq3FlXZ1lnPOuRbMk4hzzrmieRJJ5uq0A8jD40rG40rG40qmTcXlbSLOOeeK5lcizjnniuZJpAAN9TRchv3Pl/RipifjqGxtSZMlvRK9dostPz6Kda6k3UoYx3WS3pP0UqwscRySKqOfZ56kP2eeAypxXGdKeis6Zs9Ft5KXO66cPU2nfczqiSvVYybpW5KelvR8FNdZUXnaxytfXKmfY9E2KyQ9K+ne6HN5j1e+R9l9WtGtSgXwKrAh4UHH54FBZY5hPtA9q+xPwKnR+1OB86L3g6IYOwIDotgrShTHTsDWwEuNiQN4GhgBCHgA2L0J4joT+FWOZcsZV29g6+j9GsD/ov2neszqiSvVYxZto0v0vj0wDdiuGRyvfHGlfo5F2/wFcBNwbxp/k34l0rAGexpOSTVwQ/T+BuB7sfIJZrbEzF4H5hF+hkaz0O3MR42JQ1JvYE0ze8rC2fuP2DqljCufcsaVr6fpVI9ZPXHlU664zMxqo4/to8lI/3jliyufsp1jkvoCewLXZO2/bMfLk0jD+gALYp8XUv4u5w14SNJMhR6JAXqa2TsQvhQIPRxD+eNNGkef6H054jtOYUCz62KX9KnEpdDT9FaE/2KbzTHLigtSPmZR1cxzwHvAZDNrFscrT1yQ/jl2CXAydb12QJmPlyeRhuWqGyz3LW3bm9nWwO7ATyXtVM+yzSFeyB9HueL7C7ARMBR4B7gwrbgkdSF0FHqi5RiqIL5oOWPLEVfqx8zMlpvZUEJnrMMlDa5n8bTjSvV4SdoLeM/MZha6SlPE5UmkYYl6Gm4KZvZ29PoecCehempRdBlK9PpetHi5400ax8LofZPGZ2aLoj/8b4C/UVelV9a4lLun6dSPWa64mssxi2L5BKgBxtAMjleuuJrB8doeGCtpPqGafRdJ/6LMx8uTSMMa7Gm4KUnqrDD8L5I6A98FXopiODRa7FDg7uj9ROBASR0lDQAGEhrNmkqiOKLL688lbRfdATIutk7JZP6IIvsQjllZ44q2k6un6VSPWb640j5mknpI6hq9Xx2oAuaQ/vHKGVfax8vMxptZXzPrT/heetTMfki5j1ehLfBteQL2INzB8ipwWpn3vSHhjorngVmZ/QPrEMaffyV6XTu2zmlRrHMpwd0fse3eTLhs/5rw38v/FRMHMIzwB/cqcDnRQ68ljuufhB6eX4j+eHqnENcOhGqBF4DnommPtI9ZPXGlesyAIcCz0f5fAs4o9lwvU1ypn2Ox7Y6k7u6ssh4vf2LdOedc0bw6yznnXNE8iTjnnCuaJxHnnHNF8yTinHOuaJ5EnHPOFc2TiGtTJJmkmhJsp0aS39ro2jxPIq6soi/xJNNhacfckkjqKulsha7JayUtUeiufKqkCyVtlbX89dFx7p9SyK6FWy3tAFybc1aOshOBtYBLgU+y5j1X4v1vDnxRgu2MAzqVYDslI2k94L9Af+A14EZC78Z9gM0Ix/lLwoNzzpWEJxFXVmZ2ZnZZdLWxFnCJmc1v4v3PKdF23izFdkrsbEICuQ44wrKeJI666eidYz3niubVWa7ZyrQ7SOog6QyF0diWSLo+mr+WpJMkPSppoaSlkt6XNFHSdnm2uUqbiMIIdSZppKT9FEax+0LSR5ImSFqlW+xcbSLR+hZtb6ik+yR9Em3rMUnfyRNTb0l/Vxid8cuoKurQ+PYKPGSZ7V+WnUBg5XFEMseCuj6WXo9VIc7Pim9tSecqjIT4paRPJT0i6bs5fpbDMtWQkvaU9KSkxZI+lnSbpIE51ukp6YLo97s4OmZzo6q2DQv82V1K/ErEtQS3A9sQRly7i7peSTcHzgH+A9wHfAxsAIwFdpe0t5k9mGA/x0brTgQeA7YFDgC2lDTUzJYUuJ1hhDEeniIMFrQB8H3gkWg7czMLSloXeJJwBfGf6H0v4ErgoQSxA3wYvW5CYdWAZxEGH9qSlasSM69I6kfotbY/8DjwINAZ2At4UNLRZva3HNvelzB0wZ3R+kMJx2CUpO9kjoGkToQquI2AycA9hK7J+xEGUbqNUDXnmqtSdP7lk0+NmQjD/xrQP6u8hrqOArvnWG+tPOWZrqxn55hnQE1W2ZlR+WfAFlnzborm/SBXbFllI6NlDTgsa97RUfmVWeXXRuXnZZVvCSyJ5p1Z4HE8LvZz/InQ2+w6Daxzfa5jn/VzfgMcmFXelZCoviQMgpQpPyx2DPbKWueEqPyRWNneUdnFOfbdAVgj7fPTp/onr85yLcFvzOyD7EIz+zRP+ULCf7CbSdogwX7+bGYvZpVl/stOMsTwf83s+qyy64Bl8e0oDC1wEPAp8Pv4wmb2PGGY0iSuAM4lDN96EuE/+w8kvS7pb5K2TLKxaPmdgdvNbEJWfJ8AvwW+RbjCyPaomd2bVXY5oZfYXaIrnLgvszdgZkstDN/rmjGvznItQd7xUCRtT/gPdwRhGNAOWYv0AQptBJ+RoywznGi3HPMK3o6ZfS1pUdZ2NgVWB2bk+bJ8Ajii0J1a+Pf915L+BOwGbAdsTaiWOwI4XNJPLHf1Uy4jote18rTL9IheN88x77Ec8S2X9ASh6mor4I1oubeAUyVtDdxPqN56zsyWFxinS5EnEdcSvJurUNI+hCuOrwj/db8KLCZUv4wk/BfdMcF+PslRtix6rWjkdjLbim9nreh1UZ7l85XXK7pKuCWaMoOZnQqcDlwmaaKZFbLtdaLX0dGUT5ccZfm2n/ldrhXF+ll0E8RZhPao3aL5H0i6Evi9mX1dQKwuJZ5EXLMX/Yedy++ApcAwM5sdnyHpr4Qk0pxlxlvvmWd+vvJEzGwx8BtJIwkDUm0P3FHfOpFPo9cTzOzPCXebL/ZeWdvOVD/+nyQBg4BdgJ8CZxDuIP1Nwn27MvI2EdeSbQy8nCOBtCN8WTZ3cwhtAUMUDYGcpdQ/Q6bKTLGyTJVRriutqdHrjkXsa5UELqmCup9plQceLZhlZpdRd+XzvSL27crIk4hryeYDA6MntYEV44f/lvAfbbNmZksJVU5rEaqaVogatccl2V70zMy388zbARhFqFJ7KjYrc1vwKjcgmNkMwm29+0r6cZ7tbhHdppxtF0l7ZZUdR2gPmWJmb0TrD1buLlcyVzKl6F3ANSGvznIt2cXAVcCzkm4njLG+PSGB3EO4fbS5O5VQfXOypG0Jz4n0Bn5AaGT+HqGNpxCHAH+SNIdwFfEO4ZmOb0f7EPBLM3s7ts4jhDu5/ibpNqAW+MTMLo/mHww8Clwr6WfANEKbT1/C2OODCQ3w77Gye4A7Jd0JzCPcsrwHoRuWY2PLVQEXSXqScGX2XrTt6ujnPr/An92lxJOIa7HM7K+SlhD6hDqUUDX0OHA44bbTZp9EzGxR9CT7HwhfstsCcwlftIsJSeSzvBtY2eHAnoSEMZLQ/iDC3U83A38xsyey9j9J0i+BI4GfE+5ue4NwOy5mtlBSJXA84ZgeQqj6ehd4GbgMyL4tGkKby9XAaVFMX0dl483sf7HlJgGXADsREseahOQ3GbjIzJ4s8Gd3KVH+NkvnXJoknQP8GhhjZpPSjqcQCv2g/R04PMezMq4V8jYR51IWb9OJlW0B/IxQ/bPKMxfONRdeneVc+mZImge8RKjCGkioAmoHHGNmX6UZnHP18STiXPr+Smj7OAhYg9BwPQm4wMxqUovKuQJ4m4hzzrmieZuIc865onkScc45VzRPIs4554rmScQ551zRPIk455wrmicR55xzRft/q9qvdBKlFkcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot hosgns loss function \n",
    "\n",
    "plt.plot(np.arange(len(loss['sg']))*print_every, np.array([t for t in loss['sg']]) ,'r')\n",
    "plt.grid()\n",
    "plt.xlabel(\"Training Steps\", fontsize=20)\n",
    "plt.ylabel(\"Skip-Gram Loss\", fontsize=20)\n",
    "plt.title(r'$ -\\sum_{ijkl}[ p_{ijkl} \\log \\sigma (m_{ijk}) + \\kappa p_ip_jp_kp_l \\log \\sigma (-m_{ijkl})]$', fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Node Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F1_macro: 0.609783365528168 +- 0.05558362140265357\n"
     ]
    }
   ],
   "source": [
    "y = get_infection_label('SIR', 0.125, 0.001, pat_active_time, dataset, AGGR_TIME)\n",
    "nc_res=[]\n",
    "for train_test_dict in make_train_test_splits_NC(y, n_splits=30, starting_test_size=0.3, \\\n",
    "                       node_active_list=node_tuple_time, random_state=100):\n",
    "    emb1_train_idx, emb2_train_idx, y_train_idx = train_test_dict['train']\n",
    "    emb1_test_idx, emb2_test_idx, y_test_idx = train_test_dict['test']\n",
    "    \n",
    "    X_train = emb_W.numpy()[emb1_train_idx]*emb_T.numpy()[emb2_train_idx]\n",
    "    X_test = emb_W.numpy()[emb1_test_idx]*emb_T.numpy()[emb2_test_idx]\n",
    "    y_train = y[y_train_idx]\n",
    "    y_test = y[y_test_idx]\n",
    "    \n",
    "    model_clf = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, n_jobs=6, random_state=100))\n",
    "    y_test_pred = model_clf.fit(X_train, y_train).predict(X_test)\n",
    "    nc_res.append(f1_score(y_test, y_test_pred, average='macro'))\n",
    "print('F1_macro:', np.mean(nc_res), '+-', np.std(nc_res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Link Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F1_macro: 0.8205713345549522 +- 0.013561788280833464\n"
     ]
    }
   ],
   "source": [
    "lr_res=[]\n",
    "for train_test_dict in make_train_test_splits_LR(n_splits=6, starting_test_size=0.3, \\\n",
    "                       node_active_list=node_tuple_time, df_events=df_tnet, random_state=100):\n",
    "    emb1_train_idx, emb2_train_idx, emb3_train_idx, y_train = train_test_dict['train']\n",
    "    emb1_test_idx, emb2_test_idx, emb3_test_idx, y_test = train_test_dict['test']\n",
    "    X_train = emb_W.numpy()[emb1_train_idx]*emb_C.numpy()[emb2_train_idx]\\\n",
    "                *emb_T.numpy()[emb3_train_idx]*emb_S.numpy()[emb3_train_idx]\n",
    "    X_test = emb_W.numpy()[emb1_test_idx]*emb_C.numpy()[emb2_test_idx]\\\n",
    "                *emb_T.numpy()[emb3_test_idx]*emb_S.numpy()[emb3_test_idx]\n",
    "    \n",
    "    model_clf = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, n_jobs=6, random_state=100))\n",
    "    y_test_pred = model_clf.fit(X_train, y_train).predict(X_test)\n",
    "    lr_res.append(f1_score(y_test, y_test_pred, average='macro'))\n",
    "print('F1_macro:', np.mean(lr_res), '+-', np.std(lr_res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Link Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load network and default list of active nodes\n",
    "supra_G = pickle.load(open('../preprocessed/SupraAdjacencyMatrix/%s/Dyn/Aggtime%dWeight1ReptimeNon.gpickle' % (dataset+'_7030_0', AGGR_TIME), 'rb'))\n",
    "pat_active_time = pickle.load(open('../preprocessed/SupraAdjacencyMatrix/%s/PatActiveTimeAggtime%d.pkl'% (dataset ,AGGR_TIME), 'rb'))\n",
    "node_active_time = np.array([str(map_node_index[int(n_node.split('-')[0])])+'-'+\\\n",
    "                                 str(map_time_index[int(n_node.split('-')[1])]) for n_node in pat_active_time])\n",
    "times = [int(x.split('-')[1]) for x in node_active_time]\n",
    "node_tuple_time = [(int(n_node.split('-')[0]), int(n_node.split('-')[1])) for n_node in node_active_time]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1280
    },
    "colab_type": "code",
    "id": "OAUpjCOqXCK1",
    "outputId": "5d16a443-9628-47f3-bab4-1a1d6c74236e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampling Random Walks...\n",
      "Warmup...\n",
      "WARNING:tensorflow:From /home/spiaggesi/anaconda3/envs/tf2/lib/python3.6/site-packages/tensorflow/python/ops/array_ops.py:5049: calling gather (from tensorflow.python.ops.array_ops) with validate_indices is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "The `validate_indices` argument has no effect. Indices are always validated on CPU and never validated on GPU.\n",
      "Training...\n",
      "step    0 - loss: 0.4505610466003418 (27.9205 seconds)\n",
      "step  200 - loss: 0.4138459861278534 (38.8677 seconds)\n",
      "step  400 - loss: 0.37898558378219604 (36.5471 seconds)\n",
      "step  600 - loss: 0.36319392919540405 (35.1069 seconds)\n",
      "step  800 - loss: 0.35400718450546265 (38.5533 seconds)\n",
      "step 1000 - loss: 0.34787166118621826 (38.3080 seconds)\n",
      "step 1200 - loss: 0.3432738780975342 (35.9054 seconds)\n",
      "step 1400 - loss: 0.33967161178588867 (35.6171 seconds)\n",
      "step 1600 - loss: 0.33674106001853943 (35.7738 seconds)\n",
      "step 1800 - loss: 0.3343159258365631 (37.0103 seconds)\n",
      "step 2000 - loss: 0.3322913348674774 (34.2305 seconds)\n",
      "step 2200 - loss: 0.33043172955513 (36.8444 seconds)\n",
      "step 2400 - loss: 0.32875069975852966 (36.0171 seconds)\n",
      "step 2600 - loss: 0.32729747891426086 (35.1550 seconds)\n",
      "step 2800 - loss: 0.32592958211898804 (35.5750 seconds)\n",
      "step 3000 - loss: 0.3247249722480774 (37.0568 seconds)\n",
      "step 3200 - loss: 0.32358676195144653 (36.8186 seconds)\n",
      "step 3400 - loss: 0.3224763870239258 (32.8196 seconds)\n",
      "step 3600 - loss: 0.32144683599472046 (38.2842 seconds)\n",
      "step 3800 - loss: 0.3204757273197174 (36.8876 seconds)\n",
      "step 4000 - loss: 0.31956747174263 (35.3905 seconds)\n",
      "\n",
      "Total time: 754.6942 seconds.\n"
     ]
    }
   ],
   "source": [
    "# train 4way-hosgns \n",
    "print_every=200\n",
    "solver = hosgns_walks.HOSGNSSolver(network=supra_G.to_undirected(), \n",
    "                         window_size=10,\n",
    "                         emb_dim=128,\n",
    "                         order=4,\n",
    "                         iters=4000,\n",
    "                         batch_size=20000, \n",
    "                         negative_samples=5,\n",
    "                         learning_rate=tf.keras.optimizers.schedules.PolynomialDecay(0.05, 4000),\n",
    "                         warmup_steps=1000,\n",
    "                         random_state=42)\n",
    "loss = solver.train(print_loss='sg', print_every=print_every)\n",
    "emb_W, emb_T, emb_C, emb_S = solver.model.factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, _, df_tnet = load_temp_data(dataset, aggr_time=AGGR_TIME)\n",
    "_, _, df_samp = load_modified_temp_data(dataset, aggr_time=AGGR_TIME)\n",
    "df_del = df_tnet.loc[df_tnet.merge(df_samp, on=['tslice','i','j'], how='left', suffixes=('','_')).weight_.isnull()].reset_index(drop=True)\n",
    "\n",
    "df_samp['tslice'] = df_samp['tslice'].map(map_time_index) \n",
    "df_samp['i'] = df_samp['i'].map(map_node_index) \n",
    "df_samp['j'] = df_samp['j'].map(map_node_index)\n",
    "\n",
    "df_del['tslice'] = df_del['tslice'].map(map_time_index) \n",
    "df_del['i'] = df_del['i'].map(map_node_index) \n",
    "df_del['j'] = df_del['j'].map(map_node_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F1_macro: 0.8181839871297502 +- 0.009654671998388967\n"
     ]
    }
   ],
   "source": [
    "df_events = (df_samp, df_del)\n",
    "\n",
    "lp_res=[]\n",
    "for train_test_dict in make_train_test_splits_LP(n_splits=6, starting_test_size=0.3, \\\n",
    "                       node_active_list=node_tuple_time, df_events=df_events, random_state=100):\n",
    "    emb1_train_idx, emb2_train_idx, emb3_train_idx, y_train = train_test_dict['train']\n",
    "    emb1_test_idx, emb2_test_idx, emb3_test_idx, y_test = train_test_dict['test']\n",
    "    X_train = emb_W.numpy()[emb1_train_idx]*emb_C.numpy()[emb2_train_idx]\\\n",
    "                *emb_T.numpy()[emb3_train_idx]*emb_S.numpy()[emb3_train_idx]\n",
    "    X_test = emb_W.numpy()[emb1_test_idx]*emb_C.numpy()[emb2_test_idx]\\\n",
    "                *emb_T.numpy()[emb3_test_idx]*emb_S.numpy()[emb3_test_idx]\n",
    "    \n",
    "    model_clf = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, n_jobs=6, random_state=100))\n",
    "    y_test_pred = model_clf.fit(X_train, y_train).predict(X_test)\n",
    "    lp_res.append(f1_score(y_test, y_test_pred, average='macro'))\n",
    "print('F1_macro:', np.mean(lp_res), '+-', np.std(lp_res))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Copy of medium_matrix_factorization_tensorflow.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}