code/hosgns_3Way_Stat_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/Stat/Aggtime%d.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 3WAY P(Stat)"
   ]
  },
  {
   "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": 10,
   "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.4505593478679657 (26.5536 seconds)\n",
      "step  200 - loss: 0.23769976198673248 (38.2648 seconds)\n",
      "step  400 - loss: 0.195443794131279 (39.0625 seconds)\n",
      "step  600 - loss: 0.1774454265832901 (38.5544 seconds)\n",
      "step  800 - loss: 0.16678138077259064 (38.7622 seconds)\n",
      "step 1000 - loss: 0.15928062796592712 (38.9199 seconds)\n",
      "step 1200 - loss: 0.15361125767230988 (36.0006 seconds)\n",
      "step 1400 - loss: 0.14907245337963104 (38.1222 seconds)\n",
      "step 1600 - loss: 0.14525440335273743 (39.3366 seconds)\n",
      "step 1800 - loss: 0.14200319349765778 (37.8488 seconds)\n",
      "step 2000 - loss: 0.13908013701438904 (38.8870 seconds)\n",
      "step 2200 - loss: 0.1364922672510147 (36.8351 seconds)\n",
      "step 2400 - loss: 0.13410311937332153 (36.0536 seconds)\n",
      "step 2600 - loss: 0.13187317550182343 (38.5140 seconds)\n",
      "step 2800 - loss: 0.1298067271709442 (39.5125 seconds)\n",
      "step 3000 - loss: 0.12786470353603363 (36.6575 seconds)\n",
      "step 3200 - loss: 0.12602567672729492 (37.5195 seconds)\n",
      "step 3400 - loss: 0.12427119165658951 (36.1811 seconds)\n",
      "step 3600 - loss: 0.1225828304886818 (41.9727 seconds)\n",
      "step 3800 - loss: 0.12098202109336853 (37.6604 seconds)\n",
      "step 4000 - loss: 0.11948291212320328 (36.8746 seconds)\n",
      "\n",
      "Total time: 788.0978 seconds.\n"
     ]
    }
   ],
   "source": [
    "# train 3way-hosgns \n",
    "print_every=200\n",
    "solver = hosgns_walks.HOSGNSSolver(network=supra_G.to_undirected(), \n",
    "                         window_size=1,\n",
    "                         emb_dim=128,\n",
    "                         order=3,\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 = solver.model.factors\n",
    "\n",
    "#save embeddings\n",
    "# pickle.dump(list(solver.model.factors), open( '../output/%s_Stat_Aggtime%dWindow%dWay%s.pkl' % (dataset, AGGR_TIME, 1, '3Redu'), \"wb\" ), protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load embeddings\n",
    "emb_W, emb_T, emb_C = pickle.load(open( '../output/%s_Stat_Aggtime%dWindow%dWay%s.pkl' % (dataset, AGGR_TIME, 1, '3Redu'), \"rb\") )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, '$ -\\\\sum_{ijk}[ p_{ijk} \\\\log \\\\sigma (m_{ijk}) + \\\\kappa p_ip_jp_k \\\\log \\\\sigma (-m_{ijk})]$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAE3CAYAAABrZHykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6N0lEQVR4nO3deZwU1bn/8c+XAVRWkWXYR1Q0briAuwngiho13mwuv2hMjDGJJuZmUWI0aOJNvNlMXILGNXEhi5qgUZGrjLsRcF9QETCACIIsjrLz/P441VA03TNdPd1dM9PP+/WqV3dXn6p6uqamn65Tp86RmeGcc84Vo13aATjnnGu9PIk455wrmicR55xzRfMk4pxzrmieRJxzzhXNk4hzzrmieRJxzjlXNE8izjnnitY+7QCqgaRuQJ8iFl1iZktLHY9zzpWKJ5HKOBX4QxHLXQz8rMSxOOdcyXh1VmWsTTsA55wrB3nfWeUnqSMwETgaMOBMM7s13aicc675PIlUiKRtgAeAkcB64FQz+2u6UTnnXPN4EqkgSV2AycCBhCquz5rZvelG5ZxzxfMkUmGStgUeAfYBVgPHm9nkVINyzrkieRJJgaRewKPAbsDHwBgzezzdqJxzLjlvnZUCM1sMHAHMBDoB/5K0X77ykm6RZNFUX4kYJY2T1KJ+YUi6SlKq1X+SvivpJUn+v5OQH1Ol09hxKOnHse+LeeWOxf8RUmJmC4DDgXeArsAkScMaWeQ94CDgmxUIr8WRtCPwdeDSlEMZT7hx9IyU43DN1IKOqWI0dhzeTPiuuL8SgXgSSZGZ/YeQSJYCPYDJkvrlKb7azJ4xs9cqFmDLcj7woplNSzMIM1sJ/An4fpLlJM2RNK4sQblinU8LOKaK0dhxaGbzzewZ4P1KxOJJJH2fBLaNnt8XnaEUTNKzkv4q6TJJb0taFZ3mHl7ySMP2xkh6WtJKScsl/UPSLlllTpE0I4rlZUknSKrProqT1FnSFZJmSloTOwXPTN+Lym0F/D/gjhzxvBhV931N0mtRXE9J2lFS96i6YqGkpZKulqQS7IYJwG6SDi7BunKS9KikybHXQ6N9eq+krrH5Rf39i1muXMdaIcdUVK5kx1Vjx1QzPkelj8WyH4cFMTOfUpqAE4F1hBsQ7wJq8pS7BZiTY357YCXwEaHF1wnAScBrwHKgVzNiGxcOj83mjSHc4zI52taphOs67wMDojJHAhuAfwDHEU63ZwHvAvWxdQmYBDQAY6PlLonW/zahu5ddo7Ijo300IiuejsAa4D/R9o4HTo/2x0PAVOAywvWn30frOLEEf7d20f69LMEyc4BxCcovA34ZPT8MWAL8EmjX3L9/McuV4lgr9pgqx3GV75hqxjFR8WOxqeOQ8L0xrxSfr9E4yr0Bn/IeAKOAVdHBNBno2EjZW8idRIZFyz9KLAHF/kFOakZ8uf7hpwFvAe1j84YQ7nn5TfT6KeAVopZ/0bx9o3jqY/O+GX0pHJm1jbujL5D48hdEZTtmlc2s9y9Z8/8Wzf9cbF5NFOfYEv39HgceyvOeCF+68WlO9CUSn5fvR8OQKP7TgG9EX4inl+rvX8xypTjWij2mynFc5TummnE8pHIsNnEc3kIFkohXZxVB0hE5TpFzTfV5lh9O6AZlK+DfhH/ANUWEMjx6/JGZrY/NnxE99oy210PSlOh5J0mPS6qRtK2kDwo5rZbUmfCP8hczW5eZb2azgSeBkZJqgBHAXRYdxVGZ54DZWas8E5hsW94jMwPoEV8e6A+syLGP9okeL8ma3xl4ycz+Hpu3DeGLe0n0eTbuk+h10v3yfhRXLiMJXxLxqY7QoWZ83sN5ls98rtOBawlnMH/KUS7x3z/hchv3SaHLRMsVdFwVckxF5cpxXOU8pprxv12SY7HEx2FFeC++xXmKcErclI+zZ0R1vQ8SWmS9ChxrZg1FxrEv8K6ZPZk1P3NQzQOw0J386Oj5x4TrMEgaAUzL+sLOpwfhF3auazbvEb4kewEdgEU5yizMPJFUS/hS+G6Ocv2i9cVtTbgxM9s+wH/M7I0c8/+cNW+v6PFF2HyfRK+T7peVhC+DXKYD2U22JwL3AdfH5n2YZ/m9o8eewJvAmZJ+Y2Ybssol/vsnXC6+TwpaJlLocVXIMQXlOa7yHVPF/m+X5Fgs8XFYEZ5EihD9oWc0WTCLpEGEqqtehOqNo8zsg2aEsi8wP8f8LxIO8sej7V4GrDOzyyRdBHQys4sIX3TPRmVGAr8GLjSz/8uxzqWE0/K+Od7rS/hVtZjwCzvX2Cm1hPpi2PTlsNmXR/SL8xjC9aG4JYQvnGz7AM9lraNvFM9zOcquB16Kym3cJ9HrpPtlu+jzbsHMPiRU08TjWkP4Ei6kJdA+wOuEL5MTgTsJ1wpuyyqX+O+fcLn4PilomUjO/Zdj2UKOKSjPcZXzmCr2f5sSHYulPA4rxauzKkThLvXJwCDCr6cjzezdZqyvHeEXTZ2k9rH5/Qn1wleb2UfR7OGEX8cQ/UqMnu8HTI8O4h8Dn86TQIjWNR34fPRPmdleHXAw8GhUzTEN+Gz8FDyqvhsSW92y6PETWZu5gPCPfV3W/BlAB0kDsz7/MHL/gwI8n2P+DAtNI2HzfQLJ98sQIPtXZ6nsDTwexfpXwrWAcZI6ZAoU+/dPuNwIYFrCZaDA46qQYyoqV47jaotjqlglPhZb0nFYmHJfdPHJIFRdTSP86loK7JVw+VvIurBO6DLFCL/CbiPcb3IGoQXKs8DWsbLvAf2i53OBQbHnMwgXfNtlrX8c+VvSPEBofXIKobrlfaB/VObIKK5/AMcS6vVnEX4dPhKVEeEfroHwJXQ0cA3hQue3cnz+7aN1/lds3i7RvOOyyv6IUE2U/XmeA/6Ua58k2S/R+9tGsZ6V4G84hwJaZxGqsAw4Jzbvc9G8r5fg759kubmEHz0FL9PY/iv2mCrHcZXrmGrG/3fJjsVSHod466y2MwG3RwdZ0umZ2MEwJ2ud/y8qswehrv1jwhnOVUDXWLmBhGoUCNUBC6PnfQnVE7cAv80R8xb/8NH8McDThLrY5cA/gV2yypxK+HW0mnDd5yTCr7F7YmUGE64TNLCpOuT4Rvbhv4GbY69Pjj5/v6xyfwOeyJrXIYrl/Ox9knS/RGVOI7Ss65ngGJhDYUnkiOhzHRibp2j/zSP6wi7m759wufg+KWiZpvZfc46pchxX2cdUM/6/S3Islvo4xJNI25mAekqQRIg1CwV+C8wtYNsnAvdGz48D/hU9Px64N1rng8APy/j5B0YH+8XNWMeXoy+YTiWIZ+M+KWa/EH41/znlYyrx3z/hcvF9UtAyre24KuUx1YzPUPT/Z77jkE1NzG+tRBLxC+uVcTKhNUhS8dYjdYSLi48S7jHZl83r9PMZzqa62uz61ufMbJ2kzwH1khaYWXZLkkQUBt/6DfB/hAt+OwA/JPwqvKEZq/5ztJ5vAr9qToxsvk8gwX6RtDehJc0ezYyhuYr5+ydZLr5PCl0GWtdxVcpjqlhF/X82cRxeBPw0ep6rMURppZWBfSp8ItTfjoimXQi/NJaT8BcYofni4WWOtSNwD6Gueg3hGtBEYI8SrPtA4JtliLng/UKoejkl5eOh2L9/wctl9kmx22otx1W5jqkiYynJcUhozpz5vhhW7rh9PJEqIGkooZXPs8A3bMv7DaqS75ct+T6pvNa+zz2JOOecK5rfJ1JBzezewDnnWhw/E2kBJB1BaH1xVNqxOOdcEn4mUkEK4zBcEj2/SNLl0VubdW8gaVqUWJxzrkXzJFJZze5+xDnnWhJPIpXVWJvwn0fPjzaz7F5snXOuRfIkUiFRR28bzGyBpD6EwXDmRj19tgOeIXQh0aqa9znnqpsnkcqJV2XtR9adqcBZwK6SfphCbM45VxRPIpXTZPcGhN5avyDpSynE55xziXkT3xRIeorQjUS+4VGdc65V8DORCpI0VNLzwMvAlLTjcc655vIzEeecc0XzMxHnnHNFq7rxRHr16mXbb799Uct+9NFHdO7cubQBlYDHlYzHlYzHlUxbjGv69OmLzax3zjfT7kO/0tPw4cOtWFOmTCl62XLyuJLxuJLxuJJpi3EB0yzPd6pXZznnnCuaJxHnnHNF8yTinHOuaJ5EnHPOFc2TiHPOuaKlnkQkjZH0hqSZki5spNx+ktZL+lxs3hxJL0t6QdK0fMs655wrj1TvE5FUA1wDHAnMA6ZKmmhmr+UodwUwKcdqRpvZ4rIH65xzbgtpn4nsD8w0s1lmtgaYAJyYo9x5wF3AokoGt9GGDXD55fR49tlUNu+ccy1V2klkADA39npeNG8jSQOAk4DxOZY34CFJ0yWdXbYo27WDX/2KXk8/XbZNOOdca5R2tyfKMS+7R8grgQvMbL20RfFDzOzdaKTAyZJmmNljW2wkJJizAWpra6mvr08c6IiePWk/f35Ry5ZbQ0ODx5WAx5WMx5VM1cWV71b2SkzAQcCk2OuxwNisMrOBOdHUQKjS+kyOdY0Dvt/UNovu9uT44+3DHXYobtkya4vdLJSTx5WMx5VMW4yLFtztyVRgqKQhkjoCJwMT4wXMbIiZbW9m2wN/B75pZv+Q1FlSVwBJnYGjgFfKFmldHVsvXFi21TvnXGuUanWWma2TdC6h1VUNcJOZvSrpnOj9XNdBMmqBe6IqrvbAHWb2YNmCrauj/UcfwfLl0L172TbjnHOtSdrXRDCz+4H7s+blTB5m9uXY81nAXmUNLm7w4PD4zjswbFjFNuuccy1Z2tVZrUddXXj8z3/SjcM551oQTyKFyiSRd95JNw7nnGtBPIkUqk8fNnTo4GcizjkX40mkUO3asapPHz8Tcc65GE8iCayurfUk4pxzMZ5EEljVp49XZznnXIwnkQRW9e0LCxbAmjVph+Kccy2CJ5EEVtfWghnMndt0YeecqwKeRBJY1adPeOJVWs45B3gSSWRVbW144hfXnXMO8CSSyOrevcMTTyLOOQd4EknEOnaEfv28Oss55yKeRJKqq/MzEeeci3gSSWrwYE8izjkX8SSSVF1daOK7YUPakTjnXOo8iSRVVwerV8OiRWlH4pxzqfMkklR8cCrnnKtyBScRSTWSOuWYf5ik30n6uaQhpQ2vBfLBqZxzbqMkZyK/Aj6QtHGAcUknA5OB84ALgGclDSptiC2MD07lnHMbJUkinwKmmNny2LyfAMuA04EfAtsC/12q4Fqk7t2hWzdPIs45R7IkMgiYmXkhaQdgF+AqM7vNzH4FPACMKW2ILVBdnVdnOeccyZJIN2BF7PUhgAEPxua9CgxMEoCkMZLekDRT0oWNlNtP0npJn0u6bMn5DYfOOQckSyILgPiF8yOAlcD02LwuwLpCVyipBrgGOAbYDThF0m55yl0BTEq6bFn4DYfOOQckSyLPACdI+rSkI4DPAY+Y2dpYmR2A+QnWuT8w08xmmdkaYAJwYo5y5wF3AYuKWLb06upg2TJYsaLJos4515a1T1D2fwhf0v+MXm8ALs+8KakbMIrwZV6oAUB8hKd5wAHxApIGACcBhwH7JVk2to6zgbMBamtrqa+vTxDiJg0NDdTX19O7oYHdgal33cVHQ9Jv1ZyJq6XxuJLxuJLxuJIpV1wFJxEze1nSAcAZ0ay/mNnUWJFhwEPAnQm2r1ybynp9JXCBma2XNiteyLJhptn1wPUAI0aMsFGjRiUIcZP6+npGjRoFW20FP/0p+/XpA0Wuq5Q2xtXCeFzJeFzJeFzJlCuuJGcimNnLwPfzvPcE8ETC7c8jtPrKGAi8m1VmBDAhSiC9gGMlrStw2fLwGw6dcw5ImERykdQB2AP42MzeSLj4VGBodKf7fOBk4NR4ATPbWF8k6RbgPjP7h6T2TS1bNn37QocOfnHdOVf1knR78gVJf5W0XWzejoRmvdOA1yTdHX25F8TM1gHnElpdvQ781cxelXSOpHOKWbbQbTdLu3YwaJAnEedc1UtyJvIVoL+ZfRCb92tgJ+ARoCfhwvuZwB8LXamZ3Q/cnzVvfJ6yX25q2YrxGw6dcy5RE9/dCNVPwMbWWMcSzgCOIDS5nUFIIm2f3yvinHOJkkhvwg2HGQcRzmQmAET3i0wGdixZdC1ZXR28+y6sWZN2JM45l5okSeRDoHvs9UhCk9p4i6xVQNcSxNXy1dWBGcxPcm+lc861LUmSyFvAMZK2ktQR+DzwkpktjpWpY/O7ytsuH5zKOecSJZHrCd2avEVoDbUDcFNWmQMIrbXaPr9XxDnnCk8iZnYr8AugE6Fa6+poAsIIh8D2wJTShthCDYruc/QzEedcFUt6x/qPgB/lefsJoAfwUXODahW23hpqaz2JOOeqWrPvWM+IetKtrqZKfq+Ic67KJU4ikgYThsPdhzAc7nLCmCK3mVl1/Syvq4MXX0w7CuecS02iJCLpa8DvgY5s3ovuZ4CLJX3HzK4rXXgt3ODBcO+9oamvcnUq7JxzbVuSvrMOB8YDqwnjiBwG7Bo9/oxwj8g1UbnqUFcHq1bB+++nHYlzzqUiyZnIDwg3HA43s7dj898A6iXdSqjW+gHwcOlCbMEyzXzfeQf69Ek3FuecS0GS+0T2J/ST9XauN6P5f4vKVQe/4dA5V+WSJJFtgMVNlHk/Klcd/IZD51yVS5JE3iFc/2jMaKB6vlG33Ra6dvUzEedc1UqSRO4B9pN0raRt429I6i7pd4SqrLtLGF/LJnmX8M65qpbkwvrPgROAc4DTJL1I6Bq+L7A3offeGVG56uE3HDrnqliSvrNWAAcTRi2sAQ4l9OT7yej1H4FDonLVw89EnHNVLGnfWcuBr0v6FvAJQkeMy4E3okGpqk9dHXzwATQ0QJcuaUfjnHMVleSayEZmts7MXjGzJ6PHtQCSrpG0pLQhtnDeQss5V8WKSiKN6EToT6tgksZIekPSTEkX5nj/REkvSXpB0jRJh8bemyPp5cx7zQ+/CH6viHOuipWsF99iSKoBrgGOBOYBUyVNNLPXYsUeBiaamUkaBvyVUJWWMTprdMXKit+17pxzVabUZyJJ7Q/MNLNZUVfyE4AT4wXMrMHMLHrZmTCue8vRrx+0b+/VWc65qpR2EhkAzI29nhfN24ykkyTNAP4FfCX2lgEPSZou6eyyRppPTQ0MHOhnIs65qqRNP/JLsDLpZuB0M6spsPzngaPN7Kzo9ZeA/c3svDzlPwVcYmZHRK/7m9m7kvoAk4HzzOyxHMudDZwNUFtbO3zChAlFfDpoaGigS44WWHuffz5av57nr7qqqPU2V7640uZxJeNxJeNxJdOcuEaPHj3dzEbkfNPMSjYBNwPrE5Q/CJgUez0WGNvEMrOBXjnmjwO+39Q2hw8fbsWaMmVK7jdOP91s0KCi19tceeNKmceVjMeVjMeVTHPiAqZZnu/URi+sS5qVMGH1Slh+KjBU0hBgPnAycGpWDDsBb5uZSdqXMCDWEkmdgXZm9mH0/CjgsoTbL43Bg2H+fFi7Fjp0SCUE55xLQ1Ots7YvYp0F14+Z2TpJ5wKTCHe932Rmr0o6J3p/PPBZ4HRJa4GVwBejhFIL3KMwomB74A4ze7CIeJuvrg42bAiJZPvtUwnBOefS0FQSGVLuAMzsfuD+rHnjY8+vAK7IsdwsYK9yx1eQ+A2HnkScc1Wk0SRiZt7kqBB+w6Fzrkql3cS3bfAk4pyrUp5ESmGbbcIY637DoXOuyngSKRXvEt45V4U8iZSKD07lnKtCnkRKpa4unImUsAcA55xr6TyJlMrgwbByJSxOr0Nh55yrNE8ipeKDUznnqpAnkVLxZr7OuSqUaFAqSZ2ArwJ7AwOBXB1FmZkd3vzQWhkfnMo5V4UKTiLRqIIPAb0BNVK0Oq8sb7cddO7s1VnOuaqSpDrrSkIC+QmhY8YOZtYux1TQWCJtjuT3ijjnqk6S6qwDgbvM7GflCqbVyzTzdc65KpHkTKQB8G/IxvgNh865KpMkiTwCHFCuQNqEwYPDfSIffZR2JM45VxFJksiPgF0lXahoJCiXxe8Vcc5VmYKviZjZLEmHAk8BX5P0ArA8d1H7aonia13iSWTXXdONxTnnKiBJE9+BwD+BHtGUb9RDI9xLUn38hkPnXJVJ0jrrSmBn4CbgVuBdYF0ZYmq9+veHmhpPIs65qpEkiRwGTDKzs8oVTKvXvj0MHOjXRJxzVSPJhfV2wMvlCqTN8BsOnXNVJEkSeQbYo9QBSBoj6Q1JMyVdmOP9EyW9JOkFSdOii/sFLZsKv+HQOVdFkiSRi4BRkk4u1cYl1QDXAMcAuwGnSNotq9jDwF5mtjfwFeCGBMtWXl0dzJ8P6/xykXOu7UtyTeQ4wg2Ht0s6B5hO/ia+Py1wnfsDM81sFoCkCcCJwGuxlTXEyndmUwePTS6bisGDYf16ePfdTa21nHOujUqSRMbFnn8qmnIxoNAkMgCYG3s9jxx3xUs6Cfg50IeQzApetuLiXcJ7EnHOtXFJksjoMmw/153vW3Qlb2b3APdI+hQhQR1R6LIAks4Gzgaora2lvr6+qGAbGhqaXLbTwoXsD7z24IMsWr++qO2UI640eFzJeFzJeFzJlC0uM0ttAg4iNBvOvB4LjG1imdlAr2KWNTOGDx9uxZoyZUrThRoazMDs8suL3k5SBcWVAo8rGY8rGY8rmebEBUyzPN+paQ+POxUYKmmIpI7AycDEeAFJO2X66pK0L9ARWFLIsqno3Bl69fIWWs65qpBoeNxSM7N1ks4FJgE1wE1m9mp04R4zGw98Fjhd0lpgJfDFKDPmXDaVD5Jt8GC/4dA5VxWSjrHeD/gxcDThwnbHHMXMzJJ07Hg/cH/WvPGx51cAVxS6bItQVwczZqQdhXPOlV2SDhgHAM8CtcCrwFaEQapWAztE63qB3M1+q0tdHUyaBGZh2FznnGujklwTuQToC4wxs72ieTeb2ScISWQSsA3wX6UNsRUaPBg+/hg++CDtSJxzrqySJJGjgQfN7P+y3zCzecDnCUnk0hLF1nrF7xVxzrk2LEkS6UuoxspYT0gawMY7yycT7hqvbp5EnHNVIkkSWcHmF9KXEi6uxy0Hejc3qFYvc6e6t9ByzrVxSZLIO8Cg2OsXgcMkdQKQ1A44itD9SHXr1Qu22cbPRJxzbV6SJPIwMFpSh+j1rUB/4ClJvwSeBHYH/lLaEFshKVRp+ZmIc66NS3KfyI2EKqxewAIzu03ScOA8YFhUZgJweWlDbKV8cCrnXBUo+EzEzN4ysyvMbEFs3neBfoR+rPqZ2almtqoMcbY+PjiVc64KJLnZ8HRgoZlNis83s/eB90sdWKtXVwfvvw8rV4brI8451wYluSZyEzCmXIG0Od5CyzlXBZIkkfcSlq9ufq+Ic64KJEkKDxJaZ3kiKUQmifiZiHOuDUuSEC4CugI3SupVpnjajv79oV07PxNxzrVpSZr43km4I/104GRJcwhVXNlD0pqZHV6a8FqxDh1gwABPIs65Ni1JEhkVe74VsEs0Zcs5znlV8hsOnXNtXJL7RNoVONWUM+BWxW84dM61cX6RvJzq6mDePFi/Pu1InHOuLDyJlNPgwbBuHSxY0HRZ55xrhZq8JiKpN2HckLlmZtG8kcDIHMWfN7N7SxtiKxa/V2TgwHRjcc65Mmg0iUjaDpgJTMtqcTUK+EmORVZIGmJmS0sXYisWTyKHHJJuLM45VwZNVWedRrg3ZGyO9wz4Wmz6IdAtWqZgksZIekPSTEkX5nj/NEkvRdNTkvaKvTdH0suSXpA0Lcl2K8K7PnHOtXFNVWeNAWaY2bO53jSzG+OvJZ0KHAtcXcjGJdUA1wBHEgazmippopm9Fis2GxhpZkslHQNcDxwQe3+0mS0uZHsV16ULbLedt9ByzrVZTZ2JDAOeSrC+54A9E5TfH5hpZrPMbA1hPJLNxmg3s6di1WPPAK3r4oJ3Ce+ca8OaSiK9CHelZ3sB+FOO+QuBngm2PwCYG3s9jy3HbY/7KvBA7LUBD0maLunsBNutnMGDvTrLOddmNVWdtR7YOnummf0T+GeO8h1Jdse6cszLubyk0YQkcmhs9iFm9q6kPsBkSTPM7LEcy54NnA1QW1tLfX19ghA3aWhoSLzsTjU19J01iyemTAnD5pZBMXFVgseVjMeVjMeVTNniMrO8EzALmNhYmazyE4HZCcofBEyKvR4LjM1RbhjwNrBzI+saB3y/qW0OHz7cijVlypTkC/3612ZgtmRJ0dttSlFxVYDHlYzHlYzHlUxz4iK00M35ndpUddZUQvfvPZpKRlGZw4CcF+EbWf9QSUMkdQROjhJRfL2DgbuBL5nZm7H5nSV1zTwHjgJeSbDtyvAWWs65NqypJHIn0Bm4VmqyLuYawk2JdxS6cTNbB5wLTAJeB/5qZq9KOkfSOVGxSwjXWa7NaspbCzwh6UVC4vqXmT1Y6LYrxgencs61YY1eEzGzf0h6EvgC0FfST4HHzWwtgKQOwKcIY42MBJ6ycL2kYGZ2P3B/1rzxsednAWflWG4WsFf2/BbHk4hzrg0rpCv4LxDOFEYSEsY6SUsIF8B7ResQoSrp82WKs/Xq3Ru23tqrs5xzbVKTHTCa2QLgQOBSQhPcDkBfoF/0fB7hovZBZparOXB1k7xLeOdcm1XQoFRm9jEhiVwqaRAhiQhYYGZzG13Ywe67w5Qp8MEH4Q5255xrIxJ3BW9mc81sqpk96wmkQOPGwbJlcNFFaUfinHMl5eOJVMKwYXDeeXDddTCt5fUT6ZxzxSo6iUj6jqRZpQymTRs3Dmpr4Zvf9JEOnXNtRnPORLYF6koUR9vXvTv8+tcwdSrceGPT5Z1zrhXw6qxKOuUUGDkSxo6FxS2z93rnnEvCk0glSXDNNbBiRUgkzjnXynkSqbTdd4fzz4cbboBnnkk7Gueca5bmJJF/AF8pURzV5ZJLoH9/+Na3/CK7c65VKzqJmNmLZnZrKYOpGl27wm9/C889F5r9OudcK1VUEpH0SUnflnRx9PjJUgfW5n3+83D44eEGxEWL0o7GOeeKUlC3JxmSDgFuAnbKzCIaiVDSW8BXzezJkkbYVklw9dXhRsQLLoCbb047IuecS6zgMxFJw4HJwFDgMeAy4BvR4+PAzoTxzvctQ5xt0yc+Ad/7HtxyCzzpudc51/okqc66nHDmcqKZjTazS83suuhxFHASYYz1y8sQZ9v14x/DoEHhTvZ169KOxjnnEkmSRA4G7jaze3O9GQ1GdU9UzhWqc2e48kp46SW49tq0o3HOuUSSJJENwMwmyrxFdI3EJXDSSXD00XDxxfCeD8ninGs9kiSRaTQ9HO1ehPHOXRISXHUVrFoFP/hB2tE451zBkiSRHwNHSvpGrjclfQs4HLi4FIFVnaFD4Yc/hNtug8ceSzsa55wrSJImvkcBjwBXSzqf0CJrIVALHEpotfUgcLSko2PLmZn9tDThtnFjx8Kf/xwusj//PHTokHZEzjnXqCRJZFzs+dBoynZMNMUZkDeJSBoD/A6oAW4ws19kvX8acEH0sgH4hpm9WMiyrU6nTvD738OJJ4bqrf/+77Qjcs65RiVJIqNLvXFJNcA1wJHAPGCqpIlm9lqs2GxgpJktlXQMcD1wQIHLtj7HHw/HHQc/+QmcfHLoY8s551qogpOImT1ahu3vD8w0s1kAkiYAJwIbE4GZPRUr/wwwsNBlWyUJfve70Nvv974Hd96ZdkTOOZdX2l3BDwDmxl7Pi+bl81XggSKXbT123DFcH5kwAR55JO1onHMuL5mld1uHpM8DR5vZWdHrLwH7m9l5OcqOBq4FDjWzJQmXPRs4G6C2tnb4hAkTioq3oaGBLl26FLVsUu1Wr2a/r3yFDR06MO2Pf8QaucheybiS8LiS8biS8biSaU5co0ePnm5mI3K+aWY5J8LNheuAnWOv1xcwrcu3zhzbOAiYFHs9Fhibo9ww4O1MLEmWzZ6GDx9uxZoyZUrRyxblvvvMwOyKKxotVvG4CuRxJeNxJeNxJdOcuIBpluc7tbFrIo8RWlZ9nPW6lKYCQyUNAeYDJwOnxgtIGgzcDXzJzN5Msmyrd9xxoaXWj38crpV873vQLu0aSOec2yRvErHQqWLe16VgZusknQtMIjTTvcnMXpV0TvT+eOASoCdwrSQIZzoj8i1b6hhTd/PN8LWvhRsR778f/vSn0GGjc861AInGEymEpK3MbHWh5c3sfuD+rHnjY8/PAs4qdNk2p0cP+NvfQnfx3/52GH9k/Hj44hfTjsw55xKNJ3JBAWU6EsZed6UkwZlnwgsvhDFITj4ZTj8dVqxIOzLnXJVLUsH+P5JOyfdmdPPf3wndo7hy2HFHePxxGDcO7rgD9toLnngi7aicc1UsSRJ5Crg5amq7mSiB/AX4NHB1iWJzubRvH+5mf/zxcJF95EiG3HgjrF2bdmTOuSqUJIkcD8wC7pa0R2amwtXuPwH/BfzRzL5T2hBdTgcdFKq3zjiDuttug0MOgTffbHIx55wrpYKTiJktA8YAK4EHJGXuDr8JOAX4s5l9veQRuvy6doWbbuKVcePg7bdhn33g+ushxRtInXPVJdFNB2b2H0Ivvd2ASZJuAM4A/gacWfrwXCEWjxwZhtc9+GD4+tfhM5+B999POyznXBVIfOeahW7YP0voCv5M4J/AqWa2ocSxuSQGDIBJk+A3v4EHH4Q994QHHmh6Oeeca4a894lIOr2JZR8CDiTcp3FadCMgAGb2p5JE55Jp1w6++104/HA47TQ49lj41rfC2O21tWlH55xrgxq72fAWcndzomh+JmuMz5pnhAvtLi3DhsHUqaEn4CuvDNdJPve5kFAOPjjcd+KccyXQWBLxaxyt2dZbw29/C9/4Blx7bbjj/c47Ye+9w/C7p54KnTunHaVzrpVrrO+sWysZiCuTnXcOZyOXXw633w7XXANnnx364jrzzJBkhuYa6dg555rmXcJWi86dQ/J44YVwo+KYMWEc9513Ds/vvRfWr087SudcK9OsJCLpBElXSvqdpM+WKihXRhIcemio2vrPf+DSS+Hll+GEE2CnneCKK2Dx4rSjdM61Eo0mEUnHS3pM0sgc790M3AN8GzgP+Kuku8oTpiuLfv3gkktgzpzQU/D228OFF8LAgXDGGfDss37jonOuUU2diZwA7Av8Oz5T0qcJNxl+DPwMuIDQJcpnGuuk0bVQHTqE1ltTpsArr8BXvwp33w0HHBDOTr7zHXjoIVhdcA//zrkq0VQS2R942sxWZc3/CqEp75lmdomZ/RL4JLAKOK30YbqK2X33cPF9/ny47jrYddfQRPjoo6FnTzjpJLjhBnj33bQjdc61AE0lkb6Esc2zfQpYBmysvjKz94B/AfuUKjiXom7dwoX4++6DJUvC4+mnw3PPhZEWBwyA4cNDddizz8IG77DAuWrUVBLpAXwQnxGNeb4d8EQ0gHvcbMJQtq4t6dQpjPd+7bXh+slLL8HPfw7bbBOaDh9wQLi+8uUvw9//DsuXpx2xc65Cmhoe90NgYNa84dHj83mWya76cm2JFPrl2nPPcBF+yZLQZ9d998HEiXDrrWHMk09+ksE77BDOUPbfH7p0STty51wZNJVEXgaOk9TFzBqieScRrofkGlJvCLCghPG5lq5nz3D3+6mnwrp18MwzIaH861/sMGUK3Hhj6NNrzz3DGCiZaaedvPsV59qAppLI7cB1wKOSbgV2Jlw4fw+YEi8YDU51KPB0GeJ0rUH79uEelEMPhV/8gicmTuTQDh3g6adDcrnjDhg/PpTt1QsOPHBTUtlvPz9bca4VaiqJ3EgYsfBoYG9CB4trge+YWfbtzYcTLsT/X5IAJI0BfgfUADeY2S+y3v8EcDOhqfFFZvar2HtzCFVu64F1ZjYiybZdea3r1g1GjYJjjgkz1q+H118PSSUz3XdfeK9du9Bx5EEHheSy//6hO5aamtTid841rdEkYmYbJB1HGLnwYGAJcLeZvZCjeC9CMphY6MajsdmvAY4E5gFTJU00s9dixT4g3ND4mTyrGW1mfot1a1BTA3vsEaavfS3M++AD+Pe/N52t3H47/OEP4b2ttw5NjvfcMySYYcPC8z590vsMzrnNNHUmQjTY1O3R1Fi5CcCEhNvfH5hpZrMAJE0ATgQ2JhEzWwQsipKZa2u22y6cqWSfrUyfHrpjeemlMMjWLbdsWqa2dlNiyTzutltIOs65imoyiZTZAGBu7PU84IAEyxvwkCQDrjOz60sZnEtB/GwlbtGikFQyieXll0OT41VRY8B27UJnkpmWY7vsEl4PHepd3jtXRtryVo8Kblz6PHC0mZ0Vvf4SsL+ZnZej7DigIeuaSH8ze1dSH2AycJ6ZPZZj2bOBswFqa2uHT5iQ9IQpaGhooEsLvPhbtXGtX8828+fTZfZsOr/9Np1nz6bLrFlsk3U3/arevVk5cCAfDxrEykGDWNKrFzZ0KKv79sVa0DWXqv07FsnjSqY5cY0ePXp6vmvOaZ+JzAMGxV4PBAruT8PM3o0eF0m6h1A9tkUSic5QrgcYMWKEjRo1qqhg6+vrKXbZcvK4snz8Mbz1Frz5Jrz5Jlu/8QZbv/kmPR57DJYtY6dMuQ4dQlPjnXcO0y67hDOXHXaA/v3D2U0F+d8xGY8rmXLFlXYSmQoMlTQEmA+cDJxayIKSOgPtzOzD6PlRwGVli9S1Hp06wV57hSnODBYv5rkJE9i3c+eQZN54Izw+8ACsWbOpbMeOUFcXejYeMiRM8ee9e/t9Ls6RchIxs3WSzgUmEZr43mRmr0o6J3p/vKS+wDSgG7BB0vnAboTWYPeE21NoD9xhZg+m8DFcayFB796s2HPP0PQ4bv16eOedkFDmzIHZs8M0Z07o0Th7jJVOnTYllXhyGTw4TJ5kXJVI+0wEM7sfuD9r3vjY8/fYsusVgBXAXjnmO5dcTU2oxtphh9zvf/hhSDLx5JJ5/vjjsGLF5uW32iqMyzJ4MAwaFKbs5926lf1jOVduqScR51qFrl1ztxrLWLo0JJS5c8OIkXPnbno+ZUroOj97+OFu3TZLLnVr18Lbb4frMZmpZ8+KX5txLglPIs6VQo8eYdp339zvr1sHCxbkTjJz58LUqQxZvBhuvnnz5Tp0CD0kxxNL//6hK/746+7dvfrMpcKTiHOV0L79pqqsgw/OWeTRhx5i5C67hLOWXNOMGfDII7Bs2ZYLb7019O3b9FRb6zdlupLyJOJcC2GZFmF1dY0X/PjjcFaTnWTeey9MM2fCE09s2RggY9ttt0wstbWhO5nsx222KfnndG2LJxHnWptOnWDHHcPUmLVr4f33NyWXXNP06SEhNTTkXkeXLhzQrVtIbPHkEn/eu3eYttsunHG5quJ/cefaqg4dNl0zacrKlaFrmUWLYOHCzR5XvPwy2wDMmhU6yly8OPdwyFK4LpRJKr17hy7/G3vtVWutnicR51yotspTlfZ6fT218ftq1q8PvS8vXBimxYvDGU9myrx+80148skw+mV2y7SMzp1DC7RevbZ8zDWvZ89wJuZaDE8izrlkamo2nUnka/Ict2FDaAyQnWQyz5cs2fQ4a1Z4zNV4IGObbaBnT0ZstVVIetttF5JLrsf48w4dSrUHXIwnEedcebVrt+kLfZddCltm3bpwtpOdZGKPq958ky6rV8Mrr4SyH3wQlsuna9ctE0yPHls+z57XqZM3n26EJxHnXMvTvv2mC/h5vJLdoaBZ6FlgyZJNSSXzPNfj3Lnh+dKljSefjh23TDKNTJ1mzw4devboURWt2zyJOOfaBin0AtCtW+jHrFBmoXVaJqHEH3PNmzcvjGezdOmW3d0QuhLfaKutNiWYbbfd/HnmdeZ59tS9e6to7dbyI3TOuXKSQlVX165N36OTbf36cP1m6dKNj68+8QS79+8f5mVP770XbhpdtixMuVq5xXXpsnmy6d492dSlS9m7zfEk4pxzxaqpCddYevbcOOv9Dh227CU6l8wZUCahxJNRrmnpUpg/H157DZYvD1Nj1XCw6eyse3f27t49jApaYp5EnHMuDfEzoEGDmi6fzSz0XpBJKE1MHy9ZwrYl/xCeRJxzrnWSwn02nTsXdEPpm/X1FHDbaWLex7RzzrmieRJxzjlXNE8izjnniuZJxDnnXNE8iTjnnCuaJxHnnHNF8yTinHOuaJ5EnHPOFU1mlnYMFSXpfeCdIhfvBeQZuDpVHlcyHlcyHlcybTGuOjPrneuNqksizSFpmpmNSDuObB5XMh5XMh5XMtUWl1dnOeecK5onEeecc0XzJJLM9WkHkIfHlYzHlYzHlUxVxeXXRJxzzhXNz0Scc84VzZNIASSNkfSGpJmSLkxh+3MkvSzpBUnTonnbSZos6a3osUes/Ngo1jckHV3COG6StEjSK7F5ieOQNDz6PDMl/V6SyhDXOEnzo332gqRjU4hrkKQpkl6X9Kqk70TzU91njcSV6j6TtLWkZyW9GMV1aTQ/7f2VL67Uj7FonTWSnpd0X/S6svvLzHxqZAJqgLeBHYCOwIvAbhWOYQ7QK2ve/wIXRs8vBK6Inu8WxbgVMCSKvaZEcXwK2Bd4pTlxAM8CBwECHgCOKUNc44Dv5yhbybj6AftGz7sCb0bbT3WfNRJXqvssWkeX6HkH4N/AgS1gf+WLK/VjLFrnfwN3APel8T/pZyJN2x+YaWazzGwNMAE4MeWYIMRwa/T8VuAzsfkTzGy1mc0GZhI+Q7OZ2WPAB82JQ1I/oJuZPW3h6P1TbJlSxpVPJeNaYGbPRc8/BF4HBpDyPmskrnwqFZeZWUP0skM0Genvr3xx5VOxY0zSQOA44Ias7Vdsf3kSadoAYG7s9Twa/4crBwMekjRd0tnRvFozWwDhSwHoE82vdLxJ4xgQPa9EfOdKekmhuitzSp9KXJK2B/Yh/IptMfssKy5IeZ9FVTMvAIuAyWbWIvZXnrgg/WPsSuCHwIbYvIruL08iTctVN1jpJm2HmNm+wDHAtyR9qpGyLSFeyB9HpeL7A7AjsDewAPh1WnFJ6gLcBZxvZisaK1rJ2HLElfo+M7P1ZrY3MJDwK3mPRoqnHVeq+0vSp4FFZja90EXKEZcnkabNAwbFXg8E3q1kAGb2bvS4CLiHUD21MDoNJXpcFBWvdLxJ45gXPS9rfGa2MPrH3wD8kU1VehWNS1IHwhf17WZ2dzQ79X2WK66Wss+iWJYB9cAYWsD+yhVXC9hfhwAnSJpDqGY/TNJtVHh/eRJp2lRgqKQhkjoCJwMTK7VxSZ0ldc08B44CXoliOCMqdgbwz+j5ROBkSVtJGgIMJVw0K5dEcUSn1x9KOjBqAXJ6bJmSyfwTRU4i7LOKxhWt50bgdTP7TeytVPdZvrjS3meSekvaNnq+DXAEMIP091fOuNLeX2Y21swGmtn2hO+lR8zs/1Hp/VXoFfhqnoBjCS1Y3gYuqvC2dyC0qHgReDWzfaAn8DDwVvS4XWyZi6JY36AErT9i672TcNq+lvDr5avFxAGMIPzDvQ1cTXTTa4nj+jPwMvBS9M/TL4W4DiVUC7wEvBBNx6a9zxqJK9V9BgwDno+2/wpwSbHHeoXiSv0Yi613FJtaZ1V0f/kd684554rm1VnOOeeK5knEOedc0TyJOOecK5onEeecc0XzJOKcc65onkRcVZFkkupLsJ56Sd600VU9TyKuoqIv8STTl9OOuTWRtK2kyxS6Jm+QtFqhu/JnJP1a0j5Z5W+J9vP2KYXsWrn2aQfgqs6lOeadD3QHfgcsy3rvhRJvf1fg4xKs53SgUwnWUzKS+gNPAtsDs4DbCb0bDwA+QdjPKwk3zjlXEp5EXEWZ2bjsedHZRnfgSjObU+btzyjRev5TivWU2GWEBHITcJZl3UkcddPRL8dyzhXNq7Nci5W57iCpo6RLFEZjWy3pluj97pJ+IOkRSfMkrZH0vqSJkg7Ms84trokojFBnkkZJ+pzCKHYfS/pA0gRJW3SLneuaSLS8RevbW9K/JC2L1vWopIPzxNRP0s0KozOujKqizoivr8Bdlln/VdkJBDYfRySzL9jUx9LsWBXinKz4tpP0c4WREFdKWi7pYUlH5fgsX85UQ0o6TtJTkj6StFTS3yUNzbFMraRfRX/fj6J99kZU1bZDgZ/dpcTPRFxrcBewH2HEtX+wqVfSXYHLgceAfwFLgcHACcAxko43swcTbOeb0bITgUeBA4AvAntJ2tvMVhe4nhGEMR6eJgwWNBj4LPBwtJ43MgUl9QGeIpxBPBY97wtcCzyUIHaAJdHjzhRWDXgpYfChvdi8KjHziKQ6Qq+12wOPAw8CnYFPAw9K+rqZ/THHuv+LMHTBPdHyexP2wWhJB2f2gaROhCq4HYHJwL2ErsnrCIMo/Z1QNedaqlJ0/uWTT82ZCMP/GrB91vx6NnUU2CvHct3zzM90Zf16jvcMqM+aNy6avwLYM+u9O6L3vpArtqx5o6KyBnw5672vR/OvzZp/YzT/iqz5ewGro/fGFbgfz419jv8l9Dbbs4llbsm177M+5wbg5Kz52xIS1UrCIEiZ+V+O7YNPZy3znWj+w7F5x0fzfptj2x2Brmkfnz41Pnl1lmsNLjazxdkzzWx5nvnzCL9gPyFpcILt/N7MXs6al/mVnWSI4SfN7JaseTcB6+LrURha4BRgOfCzeGEze5EwTGkS1wA/Jwzf+gPCL/vFkmZL+qOkvZKsLCo/ErjLzCZkxbcM+AmwNeEMI9sjZnZf1ryrCb3EHhad4cStzF6Bma2xMHyva8G8Osu1BnnHQ5F0COEX7kGEYUA7ZhUZABR6EXxajnmZ4UR75Hiv4PWY2VpJC7PWswuwDTAtz5flE8BZhW7Uws/3H0n6X+Bo4EBgX0K13FnAmZK+Ybmrn3I5KHrsnue6TO/ocdcc7z2aI771kp4gVF3tA7wTlZsPXChpX+B+QvXWC2a2vsA4XYo8ibjW4L1cMyWdRDjjWEX41f028BGh+mUU4Vf0Vgm2syzHvHXRY00z15NZV3w93aPHhXnK55vfqOgs4S/RlBnM7ELgx8BVkiaaWSHr7hk9HhlN+XTJMS/f+jN/y+5RrCuiRhCXEq5HHR29v1jStcDPzGxtAbG6lHgScS1e9As7l58Ca4ARZvZ6/A1J1xGSSEuWGW+9Ns/7+eYnYmYfARdLGkUYkOoQ4O7Gloksjx6/Y2a/T7jZfLH3zVp3pvrxq5IE7AYcBnwLuITQgvTihNt2FeTXRFxrthPwWo4E0o7wZdnSzSBcCximaAjkLKX+DJkqM8XmZaqMcp1pPRM9frKIbW2RwCXVsOkzbXHDowWvmtlVbDrz+UwR23YV5EnEtWZzgKHRndrAxvHDf0L4RduimdkaQpVTd0JV00bRRe3Tk6wvumdm9zzvHQqMJlSpPR17K9MseIsGCGY2jdCs978kfSXPeveMmilnO0zSp7PmnUu4HjLFzN6Jlt9DubtcyZzJlKJ3AVdGXp3lWrPfAuOB5yXdRRhj/RBCArmX0Hy0pbuQUH3zQ0kHEO4T6Qd8gXCR+TOEazyFOA34X0kzCGcRCwj3dOwebUPA98zs3dgyDxNacv1R0t+BBmCZmV0dvX8q8Ahwo6RvA/8mXPMZSBh7fA/CBfhFbO5e4B5J9wAzCU2WjyV0w/LNWLkjgN9IeopwZrYoWveJ0ef+ZYGf3aXEk4hrtczsOkmrCX1CnUGoGnocOJPQ7LTFJxEzWxjdyf4/hC/ZA4A3CF+0HxGSyIq8K9jcmcBxhIQxinD9QYTWT3cCfzCzJ7K2P0nS94CvAd8ltG57h9AcFzObJ2k4cB5hn55GqPp6D3gNuArIbhYN4ZrL9cBFUUxro3ljzezNWLlJwJXApwiJoxsh+U0GfmNmTxX42V1KlP+apXMuTZIuB34EjDGzSWnHUwiFftBuBs7Mca+Ma4P8mohzKYtf04nN2xP4NqH6Z4t7LpxrKbw6y7n0TZM0E3iFUIU1lFAF1A44x8xWpRmcc43xJOJc+q4jXPs4BehKuHA9CfiVmdWnFpVzBfBrIs4554rm10Scc84VzZOIc865onkScc45VzRPIs4554rmScQ551zRPIk455wr2v8H7gt8actZ0ckAAAAASUVORK5CYII=\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_{ijk}[ p_{ijk} \\log \\sigma (m_{ijk}) + \\kappa p_ip_jp_k \\log \\sigma (-m_{ijk})]$', fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Node Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F1_macro: 0.2734977768190634 +- 0.0786175050511122\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": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F1_macro: 0.9692040426786029 +- 0.001010137601193648\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]*emb_T.numpy()[emb3_train_idx]\n",
    "    X_test = emb_W.numpy()[emb1_test_idx]*emb_C.numpy()[emb2_test_idx]*emb_T.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/Stat/Aggtime%d.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.45056596398353577 (27.8326 seconds)\n",
      "step  200 - loss: 0.23527264595031738 (36.4922 seconds)\n",
      "step  400 - loss: 0.18689385056495667 (35.5936 seconds)\n",
      "step  600 - loss: 0.1660294234752655 (34.8449 seconds)\n",
      "step  800 - loss: 0.15345735847949982 (33.9806 seconds)\n",
      "step 1000 - loss: 0.14460784196853638 (34.3189 seconds)\n",
      "step 1200 - loss: 0.13793399930000305 (37.1656 seconds)\n",
      "step 1400 - loss: 0.1325712502002716 (38.9705 seconds)\n",
      "step 1600 - loss: 0.12813670933246613 (33.9316 seconds)\n",
      "step 1800 - loss: 0.1243031844496727 (36.0547 seconds)\n",
      "step 2000 - loss: 0.12091726064682007 (37.1183 seconds)\n",
      "step 2200 - loss: 0.11789632588624954 (35.7180 seconds)\n",
      "step 2400 - loss: 0.11513346433639526 (32.9361 seconds)\n",
      "step 2600 - loss: 0.11256448179483414 (36.7340 seconds)\n",
      "step 2800 - loss: 0.11018118262290955 (35.2570 seconds)\n",
      "step 3000 - loss: 0.10793272405862808 (33.4903 seconds)\n",
      "step 3200 - loss: 0.10580471158027649 (34.2925 seconds)\n",
      "step 3400 - loss: 0.10379299521446228 (37.6415 seconds)\n",
      "step 3600 - loss: 0.10186532139778137 (34.0595 seconds)\n",
      "step 3800 - loss: 0.10003628581762314 (32.9694 seconds)\n",
      "step 4000 - loss: 0.09830280393362045 (38.4831 seconds)\n",
      "\n",
      "Total time: 737.8897 seconds.\n"
     ]
    }
   ],
   "source": [
    "# train 3way-hosgns \n",
    "print_every=200\n",
    "solver = hosgns_walks.HOSGNSSolver(network=supra_G.to_undirected(), \n",
    "                         window_size=1,\n",
    "                         emb_dim=128,\n",
    "                         order=3,\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 = 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.8661268754945614 +- 0.005821062608457102\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]*emb_T.numpy()[emb3_train_idx]\n",
    "    X_test = emb_W.numpy()[emb1_test_idx]*emb_C.numpy()[emb2_test_idx]*emb_T.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": []
  },
  {
   "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
}