main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import random\n",
    "import pandas as pd\n",
    "import math\n",
    "from copy import copy\n",
    "from collections import defaultdict\n",
    "from matplotlib import pyplot as plt\n",
    "from itertools import combinations\n",
    "from tqdm import tqdm\n",
    "from joblib import Parallel, delayed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "seed = 123456\n",
    "random.seed(seed)\n",
    "np.random.seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def sample_graph_from_infection(g):\n",
    "    g_bone = g.copy()\n",
    "    g = nx.Graph()\n",
    "    g.add_nodes_from(g_bone)\n",
    "    \n",
    "    # sample active edges\n",
    "    rands = np.random.random(g_bone.number_of_edges())\n",
    "    \n",
    "    probas = np.array([g_bone[i][j]['p'] for i, j in g_bone.edges_iter()])\n",
    "    all_edges = g_bone.edges()\n",
    "    active_edges = [all_edges[i] for i in (rands < probas).nonzero()[0]]\n",
    "\n",
    "    g.add_edges_from(active_edges)\n",
    "    \n",
    "    for i, j in g.edges_iter():\n",
    "        g[i][j]['tried'] = False\n",
    "        g[i][j]['d'] = g_bone[i][j]['d']\n",
    "        g[i][j]['p'] = g_bone[i][j]['p']        \n",
    "        \n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def simulate_IC(g, s=None, debug=False):\n",
    "    \"\"\"return dict of node to infected times\n",
    "    edges in g should contain the follow information:\n",
    "\n",
    "    - d: the transmission delay\n",
    "    - p: infection probability\n",
    "    \"\"\"\n",
    "    infection_time = {n: float('inf') for n in g.nodes_iter()}\n",
    "    bp = {n: None for n in g.nodes_iter()}  # back tracker\n",
    "        \n",
    "    g = sample_graph_from_infection(g)\n",
    "    \n",
    "    if debug:\n",
    "        print(g.edges())\n",
    "    if s is None:\n",
    "        s = random.choice(g.nodes())\n",
    "        \n",
    "    queue = [s]\n",
    "    t = 0\n",
    "    infection_time[s] = t\n",
    "    infection_paths = []\n",
    "    while len(queue) > 0:\n",
    "        outbreak_nodes = [u for u in queue if infection_time[u] <= t]\n",
    "        queue = list(set(queue) - set(outbreak_nodes))\n",
    "        infected_nodes = set()\n",
    "        if debug:\n",
    "            print('outbreak_nodes {}'.format(outbreak_nodes))\n",
    "        for u in outbreak_nodes:\n",
    "            for v in g.neighbors(u):\n",
    "                if not g[u][v]['tried']:\n",
    "                    if debug:\n",
    "                        print('{} infects {}'.format(u, v))                    \n",
    "                    g[u][v]['tried'] = True\n",
    "                    if np.isinf(infection_time[v]):  # not infected yet\n",
    "                        queue.append(v)\n",
    "\n",
    "                    # update infection time\n",
    "                    if (infection_time[u] + g[u][v]['d'] < infection_time[v]):\n",
    "                        infection_time[v] = infection_time[u] + g[u][v]['d']\n",
    "                        infected_nodes.add(v)\n",
    "                        bp[v] = u\n",
    "        if len(infected_nodes) == 0:\n",
    "            break\n",
    "        n = min(infected_nodes, key=lambda n: infection_time[n])\n",
    "        t = infection_time[n]\n",
    "    infection_tree = nx.Graph()\n",
    "    infection_tree.add_nodes_from(g.nodes())\n",
    "    del bp[s]\n",
    "    infection_tree.add_edges_from(bp.items())\n",
    "    return infection_time, infection_tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-4173535b8d6d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtree\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msimulate_IC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'b'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'c'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mtree\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0medges\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'c'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'b'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'c'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "g = nx.Graph()\n",
    "a, b, c = 'a', 'b', 'c'\n",
    "g.add_nodes_from([a, b, c])\n",
    "g.add_edge(a, b, {'d': 3, 'p': 1})\n",
    "g.add_edge(a, c, {'d': 1, 'p': 1})\n",
    "g.add_edge(b, c, {'d': 1, 'p': 1})\n",
    "t, tree = simulate_IC(g, s=a)\n",
    "assert t == {'a': 0, 'b': 2, 'c': 1}\n",
    "assert tree.edges() == [('a', 'c'), ('b', 'c')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFaRJREFUeJztnV2MXOV5x/8PzFrBNrNQlPARCyJUEdKgxF6qJA3YGQWi\nIKImkdJKlFYQctGLQkFUakPIBXtVRZWqln7coBLyUfBWdmkBqUmAWnNRVfkSuBAMhNZqwKXZqHU9\nW1SpM9CnFzt2t5tzdnfeeee8s8/8ftKRZ9bH8ztvlvznnTnvOX9zdwEAQBzOKn0AAACQF4IdACAY\nBDsAQDAIdgCAYBDsAADBINgBAIKx5WA3swfNbNnMnlvzs/PN7Ekze9nMvmVm85M5TAAA2CqjzNgf\nkvTxdT+7R9LT7v5uSUckfSHXgQEAQBo2ygVKZnaZpCfc/X3D5y9J+oi7L5vZRZK67n7lZA4VAAC2\nwrjfsb/D3Zclyd1/LOnt4x8SAACMAydPAQCC0Rrz3y+b2YVrvor5Sd2OZsZNaQAAEnB3G2X/UWfs\nNtxO87ikzw4f3yrpsY3+sbuH3e67777ix8D4GBvji7elMMpyx0ck/b2kK8zsVTO7TdKXJH3MzF6W\ndP3wOQAAFGTLX8W4+801f3V9pmMBAIAMcPI0E51Op/QhTJTI44s8NonxzSIjrWMfS2TmTbkAAKJg\nZvIJnzwFAIAph2AHAAgGwQ4AEAyCHQAgGAQ7AEAwCPYMrKys6MSJE1pZWcG3jVz48E27L5kGL4v1\nSPT7fT948KDv37vXd7ZafsnOnb6z1fL9e/f6wYMHvd/v45tCFz580+5bzzA7R8vbUf9B6hYp2E+e\nPOn7Fxb82t27/bDkfcl9+Odhya/dvdv3Lyz4yZMn8U2RCx++afdVQbA3QL/f9/0LC37Hjh3+1vCX\nvH57S/Lbd+zw/QsLY7+bR/ZFHhs+fLkg2Bvg4MGDfu2uXbW/6LW/8Gt27fKlpSV8U+DCh2/afXUQ\n7A2wf+9eP7zJL/r0dkjyA/v24ZsCFz580+6rg2CfML1ez3e2Wme+Z9ts60u+c27Oe70evhkaGz58\nOUkJ9nEblGaKlZUVnbdjh+befHNL+89JmhsMND8/n+w8b/g6EX2Rx4Zv9nzzc3NaWVlRu91OduaC\ndewj0G63darf12CL+w8kDebm1Ov1kj7l9Ho99VutkL7IY8M3m77eYDAVoS4R7CPRbrd19VVX6fEt\n7v+YpJ+/6qrkX3ZkX+Sx4cNXnJR3s5RNAb5jd49/Zp5VMfjwlfHVoYTv2An2ETm9tvX2LaxtPXD1\n1dnW0kb0RR4bPny5INgb4vTVaNfs2uWH9P+vRjs0fPc+cPXV2a9+i+iLPDZ8+HJAsDdIv9/3paUl\nP7Bvn++cm/N5rS53OrBvny8tLU3kfhVRfZHHhg/fuKQEO52nGVhZWdH8/Lx6vV4jJ08i+yKPDR++\nFFI6Twn2TAz/x8e3zVz48G0TH2XWAACzDMEOABAMgh0AIBgEewZO12Q1Wc8V1Rd5bPjwNQXBnshg\nMNDS0pIO7Nuniy+4QOdJuviCC3Rg3z4tLS1pMNjqXSbwRR4bPnxFGHV9ZNUm6W5JP5D0nKSHJe2o\n2Cfv4s6CRK/nohoPH74yvipU4gIlSZdIOn46zCX9haRbKvab2MCbJHo9F9V4+PCV8dVRMth/JOl8\nSS1JT0i6vmK/iQy6aaLfiIibgOHDV8ZXR5FgX/XqTkn/KWlZ0tdr9pnIoJsmej0X1Xj48JXx1VFq\nxn6epL+V9DOSzpb0V5JurthvIoNukuj1XFTj4cNXxrcRKcGeoxrveknH3f2kJJnZo5I+LOmR9Tsu\nLi6eedzpdNTpdDLom4NqvLy+yGPDN3u+XNV43W5X3W53rNfIMWP/gKTnJb1Nkkn6iqTbK/bL/k7W\nNNFnDczY8eEr49sIJczYx17H7u7flXRY0rOS/mEY7g+M+7rTSPR6Lqrx8OEr48vOqO8EqZsCzNjd\n45+ZZ1UMPnxlfHUoYcZOsI9I9HouqvHw4Svjq4Ngb4jo9VxU4+HDV8ZXBcHeINHruajGw4evjG89\nKcFOg1IGotdzUY2HD18Zn0Q1XlFmpJ4rnAsfvm3ioxoPAGCWIdgBAIJBsGcgeosLDUr48JXxpUKw\nJxK9xYUGJXz4yviyMOoymtRNgZY7Rm9xoUEJH74yvirEOvbJE73FhQYlfPjK+Oog2Bsg+v0quFcM\nPnxlfHUQ7A0QvcWFBiV8+Mr46iDYJ0z0e0JzP3Z8+Mr4NiIl2HM0KM0MNCjl9UUeG77Z8+VqUMoB\nyx1HoN1u61S/r60ubhpIGszNqdfrJX3K6fV66rdaIX2Rx4ZvNn29wWAqQl0i2EcieosLDUr48JXx\nZSfl3SxlU4Dv2N3jn5lnVQw+fGV8dSjhO3aCfUSit7jQoIQPXxlfHQR7Q0RvcaFBCR++Mr4qCPYG\nid7iQoMSPnxlfOtJCXaKNjIQvcWFBiV8+Mr4JBqUijIjLS7hXPjwbRMfDUoAALMMwQ4AEAyCHQAg\nGAR7BqLXc1GNhw9fGV8qBHsi0eu5qMbDh6+MLwujro9M3RRoHXv0ei6q8fDhK+OrQqUuUJI0L+mQ\npBclvSDpgxX7TGzgTRK9notqPHz4yvjqKBnsX5F02/BxS1K7Yp+JDLppot+IiJuA4cNXxldHkWCX\ndK6kf9rCfhMZdNNEr+eiGg8fvjK+OkoF+/slfUfSQ5KekfSApHMq9pvIoJskej0X1Xj48JXxbURK\nsOeoxmtJWpB0u7t/38z+UNI9ku5bv+Pi4uKZx51OR51OJ4O+OajGy+uLPDZ8s+fLVY3X7XbV7XbH\neo0cM/YLJR1f8/xaSU9U7Jf9naxpos8amLHjw1fGtxFKmLGPvY7d3ZclvWZmVwx/dJ2kY+O+7jQS\nvZ6Lajx8+Mr4sjPqO0HVptXv2b8n6aikRyXNV+yT/Z2sBNHPzLMqBh++Mr46lDBjzxLsWxIFCfbo\n9VxU4+HDV8ZXB8HeENHruajGw4evjK8Kgr1BotdzUY2HD18Z33pSgp0GpQxEr+eiGg8fvjI+iWq8\nosxIPVc4Fz5828RHNR4AwCxDsAMABINgBwAIBsGegej1XFTj4cNXxpcKwZ5I9HouqvHw4Svjy8Ko\n6yNTNwVaxx69notqPHz4yviqEBcoTZ7o9VxU4+HDV8ZXB8HeANFvRMRNwPDhK+Org2BvgOj1XFTj\n4cNXxlcHwT5hot/sn6INfPjK+DYiJdhzVOPNDFTj5fVFHhu+2fPlqsbLAcsdR6DdbutUv6+tLm4a\nSBrMzanX6yV9yun1euq3WiF9kceGbzZ9vcFgKkJdIthHIno9F9V4+PCV8WUn5d0sZVOA79jd45+Z\nZ1UMPnxlfHUo4Tt2gn1EotdzUY2HD18ZXx0Ee0NEr+eiGg8fvjK+Kgj2Bolez0U1Hj58ZXzrSQl2\nGpQyEL2ei2o8fPjK+CSq8YoyI/Vc4Vz48G0TH9V4AACzDMEOABAMgj0D0VtcaFDCh6+MLxWCPZHo\nLS40KOHDV8aXhVGX0aRuCrTcMXqLCw1K+PCV8VWhkuvYtTr7f0bS4zV/P7GBN0n0FhcalPDhK+Or\no3Sw3y3pz6MHe/T7VXCvGHz4yvjqKBbskvZIekpSJ3qwR29xoUEJH74yvjpKBvshSXslfSRysEdv\ncaFBCR++Mr6NSAn2sRuUzOwTkpbd/aiZdSTVXiG1uLh45nGn01Gn0xlX3yg0KOX1RR4bvtnz5WpQ\n6na76na7Y71Gjtn670p6VdJxSf8q6Q1JX6vYL/s7WdNEnzUwY8eHr4xvI5QwYx97Hbu73+vul7r7\n5ZJuknTE3W8Z93WnkegtLjQo4cNXxpedUd8JNtoU/Dt29/hn5lkVgw9fGV8dSpixZw32DUVBgj16\niwsNSvjwlfHVQbA3RPQWFxqU8OEr46uCYG+Q6C0uNCjhw1fGt56UYKdoIwPRW1xoUMKHr4xPokGp\nKDPS4hLOhQ/fNvHRoAQAMMsQ7AAAwSDYAQCCQbBnIHo9F9V4+PCV8aVCsCcSvZ6Lajx8+Mr4sjDq\n+sjUTYHWsUev56IaDx++Mr4qxAVKkyd6PRfVePjwlfHVQbA3QPQbEXETMHz4yvjqINgbIHo9F9V4\n+PCV8dVBsE+Y6Df7p2gDH74yvo1ICfaxq/FmCarx8voijw3f7PlyVePlgOWOI9But3Wq39dWFzcN\nJA3m5tTr9ZI+5fR6PfVbrZC+yGPDN5u+3mAwFaEuEewjEb2ei2o8fPjK+LKT8m6WsinAd+zu8c/M\nsyoGH74yvjqU8B07wT4i0eu5qMbDh6+Mrw6CvSGi13NRjYcPXxlfFQR7g0Sv56IaDx++Mr71pAQ7\nDUoZiF7PRTUePnxlfBLVeEWZkXqucC58+LaJj2o8AIBZhmAHAAgGwZ6B6C0uNCjhw1fGlwrBnkj0\nFhcalPDhK+PLwqjLaFI3BVruGL3FhQYlfPjK+KpQiXXskvZIOiLpmKTnJd1Zs9/EBt4k0VtcaFDC\nh6+Mr45SwX6RpL3Dx7slvSzpyor9JjLopol+vwruFYMPXxlfHUWC/adeUPprSddV/Hwig26a6C0u\nNCjhw1fGV0fxYJf0Lkn/LGl3xd9NZNBNEr3FhQYlfPjK+DYiJdizNSiZ2W5JhyXd5e5vVO2zuLh4\n5nGn01Gn08mlbwQalPL6Io8N3+z5cjUodbtddbvdsV4j10y9JembWg31un2yv5M1TfRZAzN2fPjK\n+DZCCTP2XOvYvyzpmLvfn+n1ppLoLS40KOHDV8aXnVHfCdZvkq6R9Jako5KelfSMpBsq9sv+TlaC\n6GfmWRWDD18ZXx1KmLGPHexbFgUJ9ugtLjQo4cNXxlcHwd4Q0VtcaFDCh6+MrwqCvUGit7jQoIQP\nXxnfelKCnaKNDERvcaFBCR++Mj6JBqWizEiLSzgXPnzbxEeDEgDALEOwAwAEg2AHAAgGwZ6B6PVc\nVOPhw1fGlwrBnkj0ei6q8fDhK+PLwqjrI1M3BVrHHr2ei2o8fPjK+KoQFyhNnuj1XFTj4cNXxlcH\nwd4A0W9ExE3A8OEr46uDYG+A6PVcVOPhw1fGVwfBPmGi3+yfog18+Mr4NiIl2LNV480CVOPl9UUe\nG77Z8+WqxssByx1HoN1u61S/r60ubhpIGszNqdfrJX3K6fV66rdaIX2Rx4ZvNn29wWAqQl0i2Eci\nej0X1Xj48JXxZSfl3SxlU4Dv2N3jn5lnVQw+fGV8dSjhO3aCfUSi13NRjYcPXxlfHQR7Q0Sv56Ia\nDx++Mr4qCPYGiV7PRTUePnxlfOtJCXYalDIQvZ6Lajx8+Mr4JKrxijIj9VzhXPjwbRMf1XgAALMM\nwQ4AEAyCHQAgGAR7BqLXc1GNhw9fGV8qBHsi0eu5qMbDh6+MLwujro+s2iTdIOklST+U9PmafbKu\n7SxJ9HouqvHw4Svjq0IlLlDS6qz/HyVdptW7Vx6VdGXFfhMbeJNEr+eiGg8fvjK+OkoF+4ckfWPN\n83uqZu1Rgj36jYi4CRg+fGV8dZQK9s9IemDN81+T9EcV+01k0E0TvZ6Lajx8+Mr46igV7L9UEez3\nV+w3kUE3SfR6Lqrx8OEr49uIlGDPUY13QtKla57vkfR61Y6Li4tnHnc6HXU6nQz65qAaL68v8tjw\nzZ4vVzVet9tVt9sd6zVyzNjP1v+dPN2h1ZOn76nYL/s7WdNEnzUwY8eHr4xvI5QwYx97Hbu7vyXp\nDklPSnpB0pK7vzju604j0eu5qMbDh6+MLzujvhOkbgowY3ePf2aeVTH48JXx1aGEGTvBPiLR67mo\nxsOHr4yvDoK9IaLXc1GNhw9fGV8VBHuDrK/LunjnzkbruSL5Io8NH75xSQl2GpQysLKycmaZU1P1\nXFF9kceGD18KVOMBAASDajwAACDYAQCiQbADAASDYAcACAbBDgAQDIIdACAYBDsAQDAIdgCAYBDs\nAADBINgBAIJBsAMABINgBwAIBsEOABAMgh0AIBgEOwBAMAh2AIBgEOwAAMEg2AEAgkGwAwAEg2AH\nAAgGwQ4AEAyCHQAgGAQ7AEAwxgp2M/s9M3vRzI6a2V+aWTvXgQEAQBrjztiflPRed98r6RVJXxj/\nkLYn3W639CFMlMjjizw2ifHNImMFu7s/7e7/M3z6bUl7xj+k7Un0/7gijy/y2CTGN4vk/I79c5K+\nkfH1AAAggdZmO5jZU5IuXPsjSS7pi+7+xHCfL0oauPsjEzlKAADYMubu472A2a2Sfl3SR939vzfY\nbzwRAMCM4u42yv6bztg3wsxukPQ7kg5sFOopBwYAAGmMNWM3s1ck7ZD078MffdvdfyPHgQEAQBpj\nfxUDAADTRaNXnka8oMnMbjCzl8zsh2b2+dLHkxMz22NmR8zsmJk9b2Z3lj6mSWBmZ5nZM2b2eOlj\nyY2ZzZvZoeH/714wsw+WPqZcmNndZvYDM3vOzB42sx2lj2lczOxBM1s2s+fW/Ox8M3vSzF42s2+Z\n2fxmr9P0LQVCXdBkZmdJ+hNJH5f0Xkm/YmZXlj2qrLwp6bfc/eck/YKk24ON7zR3STpW+iAmxP2S\n/sbd3yPp/ZJeLHw8WTCzSyT9pqQFd3+fVs8X3lT2qLLwkFbzZC33SHra3d8t6Yi2kJuNBnvAC5o+\nIOkVd/+Ruw8kLUn6VOFjyoa7/9jdjw4fv6HVUHhn2aPKi5ntkXSjpD8rfSy5MbNzJe1394ckyd3f\ndPeVwoeVk7Ml7TKzlqSdkl4vfDxj4+5/J+k/1v34U5K+Onz8VUmf3ux1St4ELMIFTe+U9Nqa5ycU\nLPhOY2bvkrRX0nfKHkl2/kDSb2v12oxoXC7p38zsoeFXTQ+Y2TmlDyoH7v66pN+X9Kqkf5F0yt2f\nLntUE+Md7r4srU62JL19s3+QPdjN7Knhd16nt+eHf/7imn2iXNBUtYQzXECY2W5JhyXdNZy5h8DM\nPiFpefipxFT9+9zOtCQtSPpTd1+Q9F9a/Vi/7TGz87Q6k71M0iWSdpvZzWWPanoYax17Fe7+sY3+\nfnhB042SPprbXYATki5d83yPAnwcXMvwY+5hSV9398dKH09mrpH0STO7UdI5ks41s6+5+y2FjysX\nJyS95u7fHz4/LCnKCf7rJR1395OSZGaPSvqwpO0+Waxi2cwudPdlM7tI0k82+wdNr4o5fUHTJze7\noGmb8D1JP2tmlw3PyN8kKdrKii9LOubu95c+kNy4+73ufqm7X67V392RQKGu4cf318zsiuGPrlOc\nk8SvSvqQmb3NzEyrYwtxYlg//enxcUmfHT6+VdKmE6zsM/ZN+GOtXtD01OrvYntf0OTub5nZHVpd\n7XOWpAfdPcp/XDKzayT9qqTnzexZrX7NdK+7f7PskcEI3CnpYTObk3Rc0m2FjycL7v5dMzss6VlJ\ng+GfD5Q9qvExs0ckdSRdYGavSrpP0pckHTKzz2n1De2XN30dLlACAIgF1XgAAMEg2AEAgkGwAwAE\ng2AHAAgGwQ4AEAyCHQAgGAQ7AEAwCHYAgGD8L5cNzbyALiJJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa4a5ef9be0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 10\n",
    "g = nx.grid_2d_graph(n, n)\n",
    "for i, j in g.edges_iter():\n",
    "    g[i][j]['delay'] = 1\n",
    "pos = {n: np.array(n) for n in g.nodes()}    \n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, node_size=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def prepare_standard_infection_graph(g, p):\n",
    "    for i, j in g.edges_iter():\n",
    "        g[i][j]['p'] = p\n",
    "        g[i][j]['d'] = 1\n",
    "    return g\n",
    "\n",
    "\n",
    "def get_display_kwargs(g, infected_times, obs_nodes=set(), max_node_size=1000):\n",
    "    diff = 1\n",
    "    times = np.array(list(infected_times.values()))\n",
    "    times = times[np.invert(np.isinf(times))]\n",
    "    upper_bound = times.max() + diff\n",
    "    node_colors = [((upper_bound - infected_times[n] + diff) if not np.isinf(infected_times[n]) else 0)\n",
    "                   for n in g.nodes()]\n",
    "    def node_size(n):\n",
    "        if n in obs_nodes:\n",
    "            return max_node_size\n",
    "        elif not np.isinf(infected_times[n]):\n",
    "            return max_node_size / 4\n",
    "        else:\n",
    "            return max_node_size / 10\n",
    "    node_sizes = list(map(node_size, g.nodes()))\n",
    "    return {'node_size': node_sizes,\n",
    "            'node_color': node_colors,\n",
    "            'cmap': 'OrRd'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def add_colorbar(cvalues, cmap='OrRd'):\n",
    "    eps = np.maximum(0.0000000001, np.min(cvalues)/1000.)\n",
    "    vmin = np.min(cvalues) - eps\n",
    "    vmax = np.max(cvalues)\n",
    "    norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)\n",
    "    scm = mpl.cm.ScalarMappable(norm, cmap)\n",
    "    scm.set_array(cvalues)\n",
    "    plt.colorbar(scm)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAEACAYAAABmohcVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXV4VMf3h9+JK+7uRYq7B3cLBHf34lK8SIsXtyLBIUiR\noi1uwV0CFAgJxEM8m2R3fn8E+qNk5e4mQPh23+fZB7Iz93MnN7vnzj1z5hwhpcSMGTNmzHxdLL72\nAMyYMWPGjNkYmzFjxkyqwGyMzZgxYyYVYDbGZsyYMZMKMBtjM2bMmEkFmI2xGTNmzKQCFBtjIcR6\nIYS/EOLuR++lF0KcEEI8EUIcF0Kk/TzDNGPGjJn/bYyZGW8EGn3y3gTgTynld8ApYGJKDcyMGTNm\n/ksIYzZ9CCHyAoeklKXe//wYqC2l9BdCZAPOSCmLfp6hmjFjxsz/Lsn1GWeRUvoDSCn9gMzJH5IZ\nM2bM/PcwL+CZMWPGTCrAKpnH+wshsn7kpgjQ1VEIYU6CYcaMGcVIKUVyjk8nhAxT3v2VlDJfcs6X\nXIydGYv3rw8cBHq+/38P4IC+g6WU39Rr2rRpX30M/8vjNY/ZPF5dr5QgDJih8AXkTZGTJgNjQtu2\nA5eAIkIIbyFEL+AXoIEQ4glQ//3PZsyYMZMqEApfqQHFbgopZWcdTfVTaCxmzJgxk6J8S4tiyfUZ\n/0/j4uLytYdgFN/aeME85i/BtzbelMTyaw/ACIyKM07WiYSQX+pcZsyY+bYRQiCTuYAnhJBzFfYd\nT/IXDJPLtzSLN2PGjBmjMNVnLITIJYQ4JYR4KIS4J4QY/v59RSkghBA9hBBe7/t1VzJWszE2Y8bM\n/ywWCl9aSABGSSmLA1WBIUKIoihIASGESA9MBSoClYFpSvL2mI2xGTNm/mcx1RhLKf2klLff/z8S\neATkAloB7u+7uQOttRzeCDghpQyTUr4DTgCNDY3VvIBnxowCVCoV3t7eREZG4uTkRO7cubGzs0u2\nrpSS169fExoairW1NdmyZSNDhgwpMGIICQnh7du3JCQkkCFDBnLlyoUQqSWQ68uQEr+tECIfUAa4\nAmSVH6WAEEJoSwGRE3j90c++79/Ti9kYmzGjh2fPnrFy6VLcN27E0ULgYGlJtFpNpEZDj569GDx8\nOIULFzZaNzw8nC1btrBy4UICA/zJaGdLvJQERMdQp1Ytho4ZS7169bCwMO7hVaPRcOLECZYvXMC5\nCxfJ5uSApRAEx8SSNVs2howeQ9euXXFycjJ6zN8iuq7eU+CZguOFEE7AHuAHKWWkwp3E2u4BBo8z\nR1OYMaMFtVrNqGHD2Oq+iaY2lrS0EuS0+v9AqTcJag4mSI7EaejYtStLVq7E0lJZINXhw4fp0akT\nFR1taW8nqOhg88+MNVqj4UhYDLtVYJctO4dP/kmOHDkU6b5+/ZoWDRugDgmiSzpbWmRJi4NlojmS\nUnLpXRTbQmK5GhbN1l27aNKkiZFX5cuRUtEUyxX2HUrSaAohhBVwGDgqpVzy/r1HgIv8/xQQp6WU\nxT45ruP7PgPf/7z6fb9desdrNsb/RkrJnTt3ePPmDenSpaNSpUpYWaXMA4SPjw8PHjzA0tKS8uXL\nkz59+hTRjY6OxtPTk9jYWAoVKmTSTE0bUkpu3ryJv78/GTJkoGLFiooNjiG8vb15+PAh1tbWVKhQ\ngbRpU6YuQWRkJFevXiUuLo4iRYpQoEABozU0Gg1d3dx4+tefzLIFZz2z0yiNZIpKkqe2Czv37zc4\nk921axfD+/bh18wOlLS30dlPSsmGd7Hsk9Zcun6DnDn1P+V6e3tTvWIFuqe1pk/2dHrdETfCohn8\nPJDVm9xxdXXVq6sNjUbDtWvXCA4OJkuWLJQvXz7F3R8pZYxXKuw7GK3GeDMQJKUc9dF7c4EQKeVc\nIcR4IL2UcsInx6UHrgPlSJycXwfKv/cf6x6v2Rj/Pzt37GDOz3OIjoqiUMEC+PkHEBIayvBhwxk1\nerTRj4wfePDgAZN+/JHzFy5Qrkwp4uPjuXPvAa5t2jDn55/JmjWrSbpRUVFMmzqVTe7ufFekEE6O\njty9/4CiRYvy008zqVmzpkm6AJvd3fll7i+oExLIny8vb976EREZycgRIxk2fLjJX77bt28zedIk\nrnheoVyZ0qhUKu7ef0iH9u2ZPWcOGTNmNEk3PDycKZMns3XbNooVLYKjgwO3796jVMlSzJo9m8qV\nKyvWmj1jBrsXL+RXOwtsFfyeKikZHauh5dDhTJ81S2e/O3fuUK9GddZkceQ7O2tFY1kfGsPptJm5\n+eChzs9fQkICZYsXoyXR9M2h7AZ/LyKG3l4BnL1yhRIlSig6RkrJyhUrWLR4EfZ2duTKmYOXr7wR\nFhaMGzuOnr16pZhRTiljvEph30H82xgLIaoD54B7JLoYJPAjcBXYDeQGvAE3KeU7IUR5YICUsv/7\n43sCk94fN0tKudngeM3GOJEF8+ezZs1qVi2dT706tf75UN26fZcfxkwif8FCbNrkbvSH7ebNmzRp\n0oRJ40bQp2cXHB0dAQgICGT+4uXsP3SUCxcuki1bNqN0o6KiqF+/HgXz5WbW9B/JlzcPAHFxcezZ\nd5CR46ewfv0GmjdvbpQuwIzp09m1ayerls6jVo1q//zO167fZNjoHylTphyrVq82+lpcunSJ1q1b\n89OUcfTo2hF7e3sA/Pz8mTPvV06ePsf58xfIlCmTUbrh4eG4uNSmbKkSTJ80jty5E2eRKpWKnR77\nGfvjdLZv30H9+oZ37sfGxpIrSxZW2EJuK+VPAW8S1PSP1eDjH4CDg4PWPt06tCfLuZP0yuCoWFdK\nSeeAGOZt2abTrXDgwAF+GtAXj++yGPU3WeETzLtq9flt8xZF4xg+bBhXrlxixeK5VKxQ9oPB5MLF\nKwz6YRytW7dh1uzZis+vj5QyxmsU9h3A19/0YTbGwL1792jYsAHXL/xJzpzZk7THxMRQs34LRo8Z\nR6dOnRTrajQaihUryqypE3Br20prnykzfubx05d47Nlj1JgnjB+P98vnbNuk3ShevXaTpm068eLF\nC5ydnRXrenp60q5dW25c/JMsWZIuFEdGRlLVpSmzZs+hVSvtv5M2EhISKFiwACt/nUuzJg219hk1\nbgrBYZG4uxucRPyLYUOHEh35jt9W/ar1Wpy7cIl2nfvw6tWrf24AutiyZQurRwxnvm4Pgk7Gx0Hv\nBYvo1atXkragoCAK5cnD4dxpSWeEkQfYHxrFlaJl+ePUKa3tDWrWoGnQS1pnTWeUbmBcAg3vvOal\nry/p0uk/9o8//mDM6FF4njtOmjRJP09BQcFUqNGArVu3UaNGDaPGoY2UMsbrFPbtx9c3xuY4Y2Dl\nihUM7t9LqyEGsLe3Z/L4UaxYoXQ5IJG//voLB3s72rm21Nln3Kih/HXqFL6+vop1Y2JiWL9hAzOn\nTdA5E6pUsRwutaqxbetWo8a8YsVyRgzpr9UQAzg5OTFxzHBWLDfuWhw+fJhcOXLoNMQAkyeM4uDB\nQwQFBSnWjYiIYNv27cyYPF7ntahVoxrly5bCw8PDoN76ZctopolTfP6PaaGJZ/2ypVrbdu3aRa00\n9kYbYoDGae25eOUyAQFJ04X7+vpy4+ZNmmROY7RuZhsramZ0VnRdVixfzoQxw7UaYoBMmTIycugA\nVq5cYfQ4PiffUtY2szEGjp84jpur/lle86YNuXHjJlFRUcp1jx3DzbWl3kdHZ2dnGtSrzSkdsx5t\n3Lp1i7x5clGwQH69/dzatOT48eOKdQFOnDhp8Fq0bdOCU6dPo1arleseP46bawu9fTJkSE+NapU5\ne/asYl1PT0++L16UXLn0Rxy4ubbkhIJr8drHh/wmGEyAfFYWvPZ9o7XN++ULCsgEk3TtLSzI7mDP\nmzdJtX18fMiTxglbE9czClpKvF+90ttHSsnJP//ETc+kAt5f4xMnTRrH5yIZO/C+OOY4YxL9rI46\n/Hwf+BBRYUx8po2NDfNmTzPYT0pJ9+7d6d5d0RZ2AEqVNLzo4ujowKnTp4zyI9ra2uLoqP9a2Ngk\nPsMbE2Xi4GDPr/MN+xPVajXt2rVTrAtQpXIFg30cHRw4eOigwWthD1hnUu7W+RhrIfAPCNB6Ditg\nRFbjZ68fiA4Pp2zZslrbCjvYAllM0rUWMGf2bIO+XiGEQRePo6MjKpXKpHF8Lr6lrG2p5abwVSlS\npAie127o7fPE6ylp06ZFrVYrrlawcMECPK/d1KsrpeThIy/OnTunWNfb2xvfN2+Jjo7Wq331+k16\n9+ptVIWFMmVKG7wWd+7eJ2fOHMZVm5g6zeC10Gg0PPF6xvXr1xXrenl58eKlN3Fx+l0LV6/fYtjQ\nYQb1ChUsQKjGtLWNUI2GAnnzatWdMmMGIRrTHoillKjtHfDy8kqie+/ePRLs9BtJfYRoYM4vvxi8\nLvnz5+P6jVt6ta5ev0mRIikTVplSfEsz49Qyjq9K//4DWLbqN73lXpauXEef3r2NCm/r0rUrR0/8\nha/vW519Tp05j1qjMWrRI3fu3FSpXIWtO3T7+qKjo1m3cSv9+vdXrAuJ12LpSv3XYsmKtfTvZ5xu\n9x492LP/EIGBuv3Bfxw9Qbr06SlXrpxi3cKFC1OsWDF27/1dZ5/w8Ajct+2kT9++BvUat2rNKY3i\n0/+LU2po3FL7o3zjxo05qdKgNmER+25MPNaOjhQsWDBJW7FixUiwsuZeRIzRuglSciw0hkaNGhns\n279ff5asWKuzXUpp0ufic2M2xt8Yrq6uqOITGDVuChpN0m/iug2bOXTkBMOGDzdKN3369IwZPZoW\n7boQEBCYpP3Bw8f06DeUX37+xegwsRk//cTkGT9z+uz5JG1RUVG4delD48aNKV68uFG6nTp1wj8w\niEnTZicxyFJKlq5Yy9kLlxk4aJBRutmyZWPwoEG0dOtGSEhokvZbt+/Sf+ho5sz52ehrMWvWbEaN\nn8qly1eTtEVERODasScd2ndQtAFk0NChHFepiTXSaKqk5IhKzWAdn5GKFSuSKUcOLkYa/xi/O0bN\n4BEjtU4ELC0tGTB0GNuD9T8laeNUcAS58+enTJkyBvv2HzAAz+u3WLRkpdbPxfRZc/H2eUPXbt2M\nHsfn5FtawDOHtr0nJCQEN7d2vHr5kr49u1KkcEHe+vmzaetOIqOiOXDgIEWKFDFaN/ERfSrLV6yg\no1tr6tSqQXx8PAf+OMaJP8+wbOlSkz/Ap06dolOnTpQtXZKObq1xdnLi2o1bbNyyg+bNmrNq9ep/\n/LvGEBAQgKtrG4ICA+nbswsF8ufDx/cNG7fsQK2RHDhwkPz59S8eakOj0TBh/HjWb9hA147tqFm9\nCrGqWPYfPMLpsxdZu2YN7dzcjNYFOHr0KN27d6dKpfK4ubbEwd4ez2s32bR1B27t3Fi6bJliH3fT\nunXJdfMqPe2VbcwA2BITz/NSZTlxLunN8QMbN25kyZiRbMjiiI2FMhPwJDaePn4RPPd+rXNDjL+/\nP98VKMDu4tkp7KgseVGsWkPHJ/6MW7yULl26KDrm1atXtGrVEqnR0KdHZ3LlzMGLl6/4bdM20qVP\nz/79vxsdL6+LlApt26mwb0e+fmib2Rh/hJQST09P3Ddtws/Pj+MnjuPhsYfGjRsnexuwr68vv61b\nx7179zhx4gSzZ8+me48eyd4GHBsby+7duzl29Ch//PEHAwYMoG+/fibdOD5GSsnFixfZsnkz/v7+\nnDh5kv3799OgQQOTdyJ+wNvbm9/WrePBgwecOHmCeXPn0bVbN6PiobURHR3Nzp07OXH8OEeOHmHI\n4CH07ddP6+O9Pnx8fKhcpgw91bE0U7BT7lhsPOssbLhy8xZ58+ouMqzRaGjbogUx1y4xJ6ODQYP8\nUhXPgIBoFqxZazC+3X3TJiaPGM6W77KRV882awCVRsMPz4NIX7k62/fsNervqdFo+Ouvv9i5Ywd7\n9uyhQ4cOdO3WjZo1a6boluiUMsZ6k0F8RAfMxjhV82GH0X9d93Nqp1bdx48f08jFhQqxUbhZC/Jo\nCXfzSVDjES/xtHXg6KlTirYVx8TE0LldW7yvXGaAkyVVHG2x+MSIRao1HA6PYW1EPD8v/lWRrxtg\n9cqVTJswnmHZ0tAqSxqcPhmzRkrOhUSyLCCSQlVrsGXXbmxtbRVpa+MLfN6SbYwNR1An4obZGKdq\nUquh+NK6n1M7NesGBASwcN48NqxbRyErC8rGxeJgIYjWSO7Y2PEkXk2vvn0YM36CUflF1Go1GzZs\nYNn8eUQEBNDQBjJaChIkvMSCk+Gx1K1Th7FTplClShWjxnzp0iXmz5rJmbNnaZo5DfmtEkPuAuM1\nHA5XkT5rNoaPG0+PHj2S/YTzLRjjfQr7umI2xqma1GwovqTu59T+FnRVKhV79+7l5vXrhIUEs8F9\nM+5bttCuXbtkJZiXUnLlyhUOHzpEsL8fGzZsZPa8eXTp0kVx2kxd+Pj4sH37dnxfvWTNqtWMmTiR\nFi1bUqlSpZRO5pMiWjq0k22MdcfY/JvWmI1xquZbMBRfQvdzan9rup9T26ybRDvZxvigwr4t+frG\n2LwDz4wZA8TGxibOjK9dIzw0FEtg69atKTIzvnz5Mn8cPkywnx/WwLz3M2ND+YsN4ePjw7Zt2/B9\n+QIbAZMnTaJ5ixZUrlw5WTPjmJgYPDw8uHPzJrYWgjGjRlK+YiVcXV2T5X/+XKSWsDUlmGfGevjW\nZhTm2WDK6gYEBLBg7lw2rFtHQQHfR0XhAMQA9x0deSYlPXv3YezEiSb5jJfOnUtEYAAuMoGMQAIS\nb0trTqviqePiwvip00zyGc/7aQZnz52nYRo78qHBUkCQRnJCBemyZuWHCRON9hm/ffuWBb/8gvum\njZR2dqCiDThYWhClVnNJBV7RKvr078+oseOMToGqjZSaGR9WGELYXCO/+szYbIz1kJoNxZfU/Zza\nqVX38ePHNKxdm7LhYbRSJ5BTy2zyjZQcsLTihnMajp0+rTiaopOrKy8vX6KnhZoKNlZJoimiNJLj\nsXG4qwU/L/6Vvv36KRrz6pUrmT5+HP2drWme1h6HT4ytRkouRalYHZlA4Rq12LrbQ9Fs9u7duzSt\nX4+GDpb0yJqGPFpC555Hq9jgH46n2oqjp04nO7QypYzxH5bKbjjN1JqvbozNO/A+4sNj46CBA2nd\nqhX29nYcPnzYqOxkuvD19WXG9Om0dXXF2cmJpUuWEBYWlmzdmJgY3N3d6dSxI87OTowdM4YnT54k\nW1dKyfnz5xnQv//7a2HPsWPHtO5QNBZvb2+mTJ6Ma5s2ODk5sXLFCiIiIpKtGx0dzfr16+nYoQNp\nnJ2ZOGECz58/N1rn9evX1K1RA7eQYAZr1FoNMUAOIRikUdPxXQj1a9bklYHsZ2q1ms5t2xJ55SK/\n2gkq2VonMcQAjhYCVwdbltlbMnXkSLZv325wzO6bNjF7wng2ZnOifXrHJIYYwEIIajjZsT6LI+8u\nX6Bbh/YG/55///03DV1cGJfJjin5Mmk1xAAFHWyZnT8zfZwE9WvV1Jph7mtgaWmh6JUaSJFRCCFG\nCiHuCyHuCiG2CSFMSM39dQkJCaFevbp069qFfLmy0r2TKwt+/omfZkyjRInieHl5maQrpWTqlCmU\nLFmSAL/XdGzbnNXLFnDp4lny58/P1i2Gqyzo4tSpU+TLl4+d27fQuH5NNq5ZhpVIoFatmvTu1ctg\n8hxdBAQEULNmDfr17UPh/Dnp3smVebOnMXHCOMqUKc2LFy9M0tVoNIwbO5ayZcsSHhpIZ7eWrF46\nn1N/HSdfvnzsUZBXVxdHjx4lb968HNi/h6YNa7NhzVISVFFUqVKZwYMGkZCgPH3lgO7daRgRTkOF\nj7j1haBpVCT9unbV22/z5s38fekik20tsFHgt81jZclcOwsG9+1LcHCwzn7+/v78MHgwy7M4kNvG\n8DKQrYXgl4z2PD5/lh07dujt26drF/pktKdFFmWbkzplS0drBwuG9O2jqP/nRlgIRS+txwqxXgjh\nL4S4+9F7O4UQN9+/XgghtGa/EkK8FELcEULcEkIk3aev7ZjkPiIKIXIAF4CiUso4IcQu4A/5Sc2n\n1OymiIuLo1atmlStWI6Fc39K4ktbt2EzM39ZxNWr14ze7jl71iz27dvD0d93JknY/uDhYxq1bM+K\n5Sto1bq1Ubo3btygSZPG7Nqyjjq1/13rLioqivZd+5IlWw42btxklG5MTAxVq1ahWaN6zJr+478W\ne6SULFu5jl9XrOXatetG16v7ceJEzpz5i8N7t5Mhw79rtd26fZembTqxaZO7osQ1H3Px4kXatGnN\n77s2U61qpX+1RURE0KZDT74rWpwVKw2Xp3zx4gXlihdnkzoeOyMWuuKkpIelNVfv3aNQoUJJ2qWU\nlC1alK6BvlSzVb7NGmC2SuIy4UfGjhunvX3mTG6t+JXpGY3L3vZXeAw7MuTkyp27WtsfPXpE7UoV\nOV82DzZG+JcjE9TUuvWKe0+8yJUrl1Fj+kBKuSlOGNiN+IGGMXFJzieEqAFEApullKW06C8A3kkp\nkxQ+FEL8TWIR0qSJWHSQUvNzS8BRJJa2dgBSxzOKQvbt24ettRWL5s3UuqjRr3d3WjRtyLKl2qs4\n6CI0NJQFCxdy0GOr1soZJYoXxX3dciZMnGC0f3Pa1KnMmjYxiSGGxLyyHtvWc/z4cR4+fGiU7o4d\nO8iWJXMSQwyJX5DhQ/pTu0ZVVq9SWuoxET8/P1auWsWB3VuSGGKAsmVKsXb5QiaacC0mT57Eork/\nJTHEkJi8f9/OTezavYu///7boNaq5cuphzTKEAPYCEEDJCt1fEauXbtG8Ns3VFYwc/2UVkLNysWL\ntboU1Go1q5Ytpb2D8dv1azvb4f3iBbduaU+NuXLJEtpndjbKEAM4WVnSMkta1ii4+X1ukjMzllJe\nAPQZ0/aArkcLgZH2NdnGWEr5BlhIYqVUXxLvFH8mV/dLsnbtGoYN6qs35Gf44H6s37DBKJ/ptq1b\nadKwns5yTgB1XWpiaWHBhQsXFOu+fv2aK55X6NpJd1IdBwcH+vbswrq1utMeamPt2jUMH6z/Wvww\npD9r1xmnu9ndnXZtWpA5s+6V9mZNGhL27h03b+rPe/wxT58+5dGjR7Rvq/vJIk0aZ3p06cj6334z\nqHf099+paYRL42NqJiRw9MABrW3Hjh2jNhosTQgrK2FtSVxkpFb/96NHj7COj6e4whngx1gJQX17\nS06cOKG1/ejhQzQzonjqxzRNa8exg0q3XHw+hBCKXibo1gT8pJS6FiUkcFwIcU0IoWgFNtnGWAiR\nDmgF5AVyAE5CiM7J1f2SeHl5Ublieb19vitSmNjYWN69e2ekrv7cvEIIKlUoa9Si2/PnzylRvKjO\nKsQfqFShnNGLeV5eTw1ei9KlvsfPz5/Y2FgjdA1fCwsLCypWKGuUf/7p06eULV3KYHa6Sgp134WH\nk87Eh+N0AkLDw7W2hQQEkB7TFj+FEKS3sSYkJCRJW2hoKBltjHN7fEwGJMGBSdO7AryLiCCjCTN5\ngEw2VoS+S/4CdXJJzszYAJ3QPSsGqCalrAA0BYa8d3noJSU2fdQH/pZShgAIIfYB1YAkS8DTp0//\n5/8uLi64uLikwOmTj42NDVEGqmYkJCQQFRVllJ/UysqK6ZO1+/k+xs8/gH79+tFPYQgTQKGChnPz\nRkVFc+HiBaPu/Pb2dkRFRZMxYwadfeLi4oiPjzdYhudjbG1t+b7oFIP93r7xo3PnznTurPx+Xup7\nwzmbo6KjOXrsmKKyS/HWpmXoiwfC3r3Teg4LoJ+T6RtEgt+F6Yw5zmNtCZlMS10Zp5EsWbiQ+QsX\nJmmztxAkmLjOE6eR2Fgrv0mcOXOGM2fOmHQufej6e99OUHM7wbQoKSGEJYnpLHTOLqSUfu//DRRC\n7Acqkbi2ppOU8Bl7A1WEEHYi8TevBzzS1nH69On/vFKLIQZo1LARHvu0P15+4PCRE1SuXElxOSAp\nJUeOHGHP/kN6faARERFcuXoDHx8fxboxMTG8Cwvn+d/6oxo89h9k3tx5Ro25Xbt2Bq/F3v2HqF+v\nnlG6O3fuxGPfIb26wcEh3H3wiMDAQMW64eHhvPZ9i4+P/mUKj30HWblypUG9ihUq8MpEA+QtJeXK\nlNGqu3T5cl5Zm7ZDLVZKom1s8Pf3T6Lr4+NDuJUNKhNLRb20tGbt2rVax1yyRAmeRplW0+5ZtIrc\nuXMr7u/i4vIv+5BS6AplK29rTR9Hu39eetCWf74B8Oi9izbpAUI4CCGc3v/fEWgI3Dc01pTwGV8F\n9gC3gDskDtw4h+JXZvCQIaxcu1FneaSYmBhm/rKQIUOGGqVbr149omNi2bNP9w75eYuWU69uXaO2\nv9rZ2dGnd2+mzPhFp6G/eu0mZ85doouBcKtPGTJkKL+uWKu1Mgkk3jzmzF/CkKHGXYvmzZvj8+YN\nfxzV7p8EmD13MS1btjBqB5ezszNdOndm2qy5Oq/FuQuXuHHrLm4KEtf3HT6cEw6m+UmP2zvSV0el\njw4dOnApVkWYCXHaf8bEUb1qFbJkSVp0NGfOnFQoV5YT4caXXQpKUHMpIkbndek1eAg7Q5W7oj5m\n5zsVvY38vnwOkhnath24BBQRQngLIXq9b+rAJy4KIUR2IcTh9z9mBS4IIW4BV4BDUkrdH/z3pEg0\nhZRyhpSymJSylJSyh5QyPiV0vxQlS5Zk9KjR1G7Ykj9Pnf3Xl/rW7bs0atGekqVK07FjR6N0LSws\n2LFjJ0NGTmDJ8jVERkb+0xYQEMiYCVPZ4bGfZcuXGz3mqdOm8cLbh269B/Hylfc/78fFxbFthwct\n2nVh06ZNRidsr1y5Mn379KV2w1acOXfhX9fi2vWbNGjuRs2atWipo9abLqysrNixYyc9+w9n1doN\n/yqm+vatH8NGTuDYn6dZtGixUboAs+fM4fa9h/QdNILXr33/eV+lUrFpyw7ade7N9u3bFblV3Nzc\neC7B18jZsZ+UPNZo6NChg9b2TJky0bJ5C/5QGbc4KKXkgIUNw8bqdncNHTsODxVGR6HsC1fR1rUt\n6dKl09qyvzjSAAAgAElEQVTetWtXroRG8ibWuK/zsygVz2PjadOmjVHHfQ6Ss4AnpewspcwhpbSV\nUuaRUm58/34vKeXaT/q+lVI2f///F1LKMlLKslLKklLKXxSN9UvF/qbmOOMP7Nq5k9lzZhMdFUWh\nggV44vUUtUbDD8N/YOSoUSbnf33w4AGTJ03i3PnzlCtTiqCgIF56+9DW1ZXZc+YYldfgY6Kiopg2\ndSqb3N0pUrgganUC3q99KVasGDNm/ETNmknD3pSy2d2dufPmkhAfT/58eXn85AkWllaMHDGSocOG\nmZxs5vbt20yZPJkrnlcoW7oU/v7+ePu8oUP79syeM8fo2OUPhIeHM2XyZLZu20axokWIi1PxytuH\n0qVKM3PWLCpXrqxYa/aMGeycP5/ZcbHYKvg946Rkso0dbUaMYPqsJCGn/3Dnzh3qVa/GQjtLCiv0\nS2+NTeBS1pzcfPhQ5+cvISGBMsWK0TQ6hO7plPnxH8TEMTgwmnNXPPVu4548fjynNq9nQ+Es2Cr4\n/EerNXR54k+XUWMZO2GCorFoI6XijC9kVbZZpYZ/2FffDm02xp8gpeTOnTu8ffuWpk2bkpCQkOyS\nSx/w8fHh4cOHNGrUiJCQENKnTxpvawrR0dFcvXqVOnXq4OXlReHCKVMuXUrJrVu38PPzo1mzZil6\nLby9vXn06BGNGzfm3bt3yS4/9YHIyEiuXbtG3bp1ef78uaIipJ+i0Wjo6ubG4+PHmRQXi5Megxwl\nJT/b2pG/bj127t9v8Ia9a9cuhvXpzWxbC0roiVSQUrJVpeaIrSOXbtww6Mby9vameoUKdLBW0yOt\nnd6b5e1oFaOCYljtvhlXV1e9umq1mvatWhF805NlBTIlqR7yMe/iExj4PIji9RqyYeu2ZGWHSylj\nfDGb9ln/p1T3e2c2xqmZ1JrE5kvrfk7t1KqrVqsZPXw4WzZupL7U0ESdQPaPjIuflBy1tOKkEHTs\n2o2lq1YpvlEdPnyY7p06Uc7GilYynnI2Vv8YrmiN5KQqngPCGoccOTh88k/FieZfv35Ni4YNSAj0\nx80GmnyULEgjJZ5RKnarJDej49i224PGjRsr0k1ISGDYgAHs3b0Lt8zOdMriTE67/w8lfBmjYltg\nJPsDwundrx9zFy5KqSoiyTbGl3Iom/BUexNqNsapmdRqKL607ufUTu26z58/Z+XSpWzasAF7wMHS\nkoCwMCydnOjRsyeDhg836UkkIiKCLVu2sHzBfAL9A0hvY0NoeBgqaxvq1K7NsLFjqVu3rtFGTaPR\ncPLkSZYvmM/ZCxfIYm9PdEQE8ba2ZMuWnaFjx9KlSxecnJyMHvOTJ09YtWwpWzZvxtnKioToKCwd\nHIjWSHr36cvAoUNNqhqujZQyxpdz6g7R/JiqviFmY5yaSe2G4kvpfk7tb0VXpVLx+vVrIiIiKFeu\nHDExMclKLP+BD+FpISEhlClThuDgYDJkUGZADBESEoKfnx8lSpTA29ubXLlypUjJpdjYWHx8fChc\nuDDPnj0jV65cKZ5YPqWM8dU8yiJzKnkHmY1xakVKiYWFxTdhKD637ufU/tZ0P9fnwvx506qdbGN8\nLV/SnDDaqPgy8Ksb49SRyDMVkJCQwO+//07D+nVIlzYNlpaWWFtZUqlCWdavX/+vUCwz/x3i4uLY\nsWMH1cqVw8nODktLSyyBKmXKsH37dlQq0zZFSCnx9PSkW/v2ZHR2xtrKCgugQPbszJw+HT8/P5PH\n/ObNG6ZPmUL+7Nmxs7ZGAJnSpKFnp05cv37dZN3Y2Fi2bNlCtTJlcH4fJpjGwYHalSqye/duk1O2\nfk4+V26Kz4F5ZgwcPHiQIYMHkjtrOga2r02TWqVI6+RAVIyKize9WL37HJdvezFl6jR++GFEsv94\n3/BM5T+lu2XzZkYNG0ZOKakTEUFxErdLxwIPgdNOTrwWgnm//krv3r0V6z569Igurq4EvH5NkzgV\ntZGkA9TAawnHbWw5p9HQrm1bVvz2m2J3SExMDIP79GHf/v3Ut7ehuYWGvFaJN49gjeTPODUH1JAj\nX3627d1rVDWOdWvW8OO4sRS1s8bNFio42OJgIYjSSK5EqditglfxGhYtX07nLl0U6+oipWbGNwom\n3SijjfLPA776zPg/b4w3bFjP5B/Hs3PBIGqU/05nv2ev/Gk7YjkNm7ZmwcJFJhlkjUZDaGgomTJl\nQqPRpOgdOSIigjRp0hAfH4+VVcrVmVWpVISFhZE1a9YUNZpqtZrQ0FAyZ86c4tciPDyctGnTJisU\nb97PP7N41iyGR0eTT08/b2CpgwNDxo/nx6lTDereuHGDxnXq0CUmmkYCrZU+ACKlZLmNLfHFinHi\n7DmDSaGioqJoULMmzi+eM8ZG4KhjV5laSg7EJrAZS06cOUuZMmUMjvmnqVPZuGQxizM5UMhOd76J\nBzFxjA6KYcy06YwYPdqgrj5SyhjfLKQshr/cM/+vboz/026KP//8k0kTx3Nq43i9hhigUN6snHWf\nwIkjv7N06RKjzhMaGsrMn34ib948FCpUCFsbG0qW/J6VK1Yk69FOSsmOHTuoXr0aWbNmxcHBgWzZ\nsjFm9Gh8fHxM1gW4desW3bp1JWPGjBQvXgwrKytc27Th/PnzydINCgpi2tSp5M6diyJFimBjY0PZ\nsmVYt3atUdU4PkWj0bB582YqV65E9uzZcXBwIHv27EwYP97oR/5dO3eyeNYsJhswxAB5gEnR0Syf\nO9dg1RYfHx+aNWjAoJgomlgInYYYwEkIxsWpsH/0iC7t2unVlVLSsXVrMr54xhRb3YYYwFIIXO2t\nGS7UNK1fj7dvtacA+MDGDRvYsGQxm7I56TXEACXsbdiUzZF506exd+9evX2/FJ8xa1uK8582xpMm\njmPF5K4Uyac73/DHpEvjyJ7FQ5g18ydiYpTlAnj79i1Vq1bh6ZMHHN67jdC3z4h558vSBbPYu2c3\nzZs1MyoV5Qc0Gg19evdm/rxfGPvDIMIDXhAV7I3nuWNo4qOpXLky9+8bzE2ilT0eHjRq1JAyJYrg\n7XWbIB8vQt48pWHd6nTs2MHoxPIf8Pb2pnLlSvi+fsHJw3sIefOUmFAf5s2awpbNm3Bt08akm5Na\nraZLl84sX7aEKeNH/HMtLvx1iKjwECpWrKA4LadGo2HciBH0j45GaUxDemBQdDQTRo3SWy9x0bx5\n1IiOprrCcDULIfghTsXVc+e4du2azn6XL1/mrucVxtpa6DXwH1PXzprq6jiWLFqks098fDyTxo1l\nXkZ7MurZ7PEx2aytmJXRnokjR3w2l5kxWAih6JUa+M8a4+vXrxPg/5YWdfTn2P2UwvmyUalUQXbt\n2qWof+fOnejk1prN61dSutT3QOIjWF2XWpw47EG6NI6M11FORx8rli/n8eMHXPjrMK1bNv3HNVGw\nQH4WzZvFvFlTaNmyBfHxxuUVeP78OYMGD+bEIQ9GjxjyT1UOZ2dnBvbrxYW/DjNz1kw8PT2N0pVS\n0q5dWwb17cFvq36lRPGiQGL+job16/DX0X1ITTzTFDzqf8r8efPwf/uGcycP0rxpo39cE0UKF2LZ\n4l+YOnE0rVq1VFQY4MSJE9hGR6P/OSkphQBnlYqjR49qbY+JicF940aaa4yb/VsLQeM4Fcv1GM1l\nCxbQykKDlZFGxdVKsH7tWp2LkIcOHSKnwOjE9RUdbLCIjOD06dNGHfc5sLC0UPRKDaSOUXwF1qxe\nRf92tU2qDDuwfW3WrDKc3OfWrVv8/fffTJ6g3X9maWnJkgWz2bJ1K+E6kpJrQ6PR8OuSX/l13myd\nvsQundzInTMHBw/qzhinjZUrVtCnR2fKlC6ptT1/vryMGznUaFfNpUuXCHv3jlE/DNbabm1tzdKF\nc1j3229GRa7Ex8ezbPlyliyYrXOhq2+vbjg62HP8+HGDeisWLcIlIiJJzkQl1I6IYPmCBVrb9u7d\nSxEL8a9dfEppiGT/779r/YyEhIRw5OhRmtoav06Q18qSglYW/P679oocqxYtpJ2t8bNbIQTt7ASr\nFuu+gXwpzG6Kb4DHD+9TtWzSwpFKqFamMI8eG37s3btnD107ttO7iJQ9ezaqV63EsWPHFJ//1q1b\n2FhbU7FCWb39enbriMfu3Yp1Afbs3UOPrvqz0/Xo2pG9e/cZVYJqj4cH3bu017ujLH++vJT6vjin\nTp1SrHvlyhWyZc1MST0J5oUQ9OzaUVH16UcPHmBqZo8iwJPHj7W2PXzwgCJRUSbpphOCLDY2vHr1\nKknbixcvyOlgTxoTtx8Xi4vl0SOt6cd5/MSLMiaUcwIoY2/NYyPrL34OvqXQtv+sMY6KjsbBzrQP\nmqO9LRERkQb/wAsWLiCrlkKkn+Lk6ECHDh0Uf3AqVKiAra2NwQ9R1iyZOXbsqGJdIQR+fv4Gx5wu\nXVrU6sRIBaW6a9etVXQtHOztaNGihWLdWrVqYWdr+O+YNUtm9uzZY1DP580bTN1LZgu88ffXqjv/\nl19M1gVQhYdTqlQprZ+F2DDTyxvZAzNnzNA65oDgIOxNnDXaCwuiFK6rfE7MM+NvgLRp0xAWYdqH\n5V1EFOnSpTFYNWLqlKl4PdVVr/D/8X3jx+HDhxVXt3j06BHBIaF6F4sAvJ4+p23bdop1pZQULFjA\n4Ji9vX1wcnJGo9Eo1h05YqSia+EXEMipU6cU6964cYOAwGCDi0Vez57TvXt3g3r58uTB1O090UCu\nHDm06k6cMsVkXQBLZ2cePHiQRPfWrVtYOacxWTdSCGbMnKl1zNmzZCHSxAoiERoNzo7G579Iacwz\n42+ASlWqc/T8PZOOPXLuDlUqJy0L/yldunZlh8d+ovQ8nj589ISnz1/QoEEDxecvWrQo2bNn58ix\nkzr7SClZs34zPXr2VKwL0L1bd9asd9fbZ816d3p0727Uh7hb9+64b9uld8fa9Ru3CAoOoVatWop1\ny5Yti62dHafO6A65U6vVrNu4VdG1qFKjBndMfOS/bWFBlWrVtOtWqcJNR9MqiPhKSSRQsGDBJG1F\nihQhKD6etybUc5NS4mlhrTPXc+UqVTlvYtml8zEJVElGPu2UQlhYKHqlBlLHKL4CAwcOYvOBC0TH\nGP9hW737LIOH/mCwX968eWnRojm9+g/XGtUQGvqO7n2HMG7sWIPVjT9l0qTJDBv947+qfHxASsm4\nH6eTKXNmoxPM9+3Xj+N/nsFjr/Y6eKfOnGPdxq0M01FeSBffffcdtWrWov+QUVrjiQMDg+g14Acm\nTpho1EYNIQSTfpzEgKGjtZbN0mg0/DD6R74r8h0VKlQwqDds1ChO29kZXcdZA5y2s2P4mDFa2xs1\nakSkrR1PTJhpHrW0plefPlqT8Tg4ONCtR3cOJBivey9eTYKjE/Xq1dPaPnT0aHbHaIwOUYuXkr2R\ncQwZMcLoMaU0ZjfFN0D+/PmpUqUK6/ecNeq4054PCQyNVpwLdvXqNcSoEqhQvT7rN23l7xcveeL1\nlHkLl1K6sgsuLnUZMXKk0eNv1aoVY0aPoVLNhkyaNpt79x/y8pU3HnsP4NKwFWcuXGH//t+NfgTL\nlCkTR44cYcS4KXTq3p/TZ8/z6tVrLl2+St9BI+jYfQC7d+/WOkszxMZNm3jjF0SV2o3ZvG0XL16+\n4tFjL+bMW0yZKnVo3boN/QcMMFq3U+fODBgwkPLV6zF91lwePHzMi5ev2Ll7HzXqNef2/UfsVBiK\nWL58eXLmzctlI8dwFcicMyeVKml/YrK0tGTwiBHstbE1yri9k5I/gUHDhunsM3j4D/wRpzaqvp6U\nkh0awRA9FWxq1KiBXYaMnIowLg7+SFg0BQsXplSpUkYd9zkQlhaKXlqPFWK9EMJfCHH3o/emCSF8\nhBA337+0GgIhRGMhxGMhhJcQYryisf6Xt0M/evQIl9o12TS7N41qGP7gPHruS/0+89mwaQtNmjRR\nfB6NRsPx48dZtXIl9+7fIygoCFdXV4YMGarzy6uUx48fs2rlSo4dP4avry+VK1emX7/+uLq6Gj3b\n/ph3797hvmkTW7Zuwd/fn7CwMCaMn0Cfvn1NLhMFiS6DI0eOsHrVKh48fEBwcDDt3dozZOhQypUz\nLub7U+7fv8+qlSs5cfIEb9+8pVq1avQfMIBWrVphbUTZeE9PT5rUrcuI6GhFkRXPgUUODhw+eZJq\nOtwUkFiFpFq5cpT1fkUnjdrgjTJKSqba2NFqyBBmz52rt+/YkSP5c+MG5hnYgQeJhvi3WDV3c+bh\n3NWrOOpxn5w7dw7Xpk1ZmcVBUbzxzfdVRI6fOUv58uUN9tdFSm2HflpVWf6Nwpe9kpxPCFEDiAQ2\nSylLvX9vGhAhpdQZtyeEsAC8gHrAG+Aa0FFKqT3U5sNx/2VjDHDx4kVc27RiQp8m9G1XG0eHpLGq\nCQlq9v95neFztrFg0RK6deuWrHOm5uQ4X1o7teoePXqULu3a0S46mhqANlMeT2Lp4N0ODrjv3EmL\nFi0M6r59+5a61atTyO8tnRLiyajFIEspeSJhuZ0d9Tt0ZOW6dQYNt0ajYUCvXlz4fT9jrCXFrLXH\nHQeqNayPk7zMko2T588rurHu27eP/j26MzqNDY3T2mOtZSwqjeRAWDQrIxLYvncvDRs2NKirj5Qy\nxs+qKdu+U+jSE63nE0LkJbG688fGOFJKuVDPeasA06SUTd7/PAGQUkq9d9T/vDGGxBny+LGjuHjp\nEl1bVKdxje9J5+xAZLSKy7efsW7POfLkzcusOXOpU6dOss+XWg3Q19BOzbo3btxgzNCh3Llzh9oJ\nCRSLj8ceiAEeW1lx1tqa70uWZP6yZUY94YSFhTF2xAh27dpFeUtLasZGkw6BGvCWkpMOjsQ4ODBh\n2jQGDByo2NUkpWTl8uXMmzUL5zgVzTVx5LGyxAoI1mj4U1hzIzaOTp06MXfxYtKkUR6FcfnyZSb8\nMJzHDx/RxtmG8jYWOFpYEKnR4Bmn4UCEinJlyzF36VLKltUf/66EFDPGNYoq6lvowmNjjHEPIBy4\nDoyWUoZ9ckxboJGUsv/7n7sClaSUehdazMb4I16+fMmaNau5duUS4RHheD3xolOXbgwYOEhRdiul\npGYD9KW1vwVdLy8vVi5dys3LlwmPiODp06f0GTSIwcOHU7Sosi+7NsLCwti8eTMHd+4gJDiEJ0+e\n0KR5c/oMGULDhg1NriOnVqs5fvw4G1etwtf7FXfv3aN6teq07NSJbt26GWWEP+Xhw4esXrqUOzeu\nc/PmTcqXL0+5KlUZPHw4hQqZtolKGylljJ/X0r4Z6Mq7KDzf/X+U01IdlT60GOPMQJCUUgohZgHZ\npZR9PjmmHdDwE2NcUUqpd9XfbIz18C0Yii+h+zm1vzXdz6lt1k2inWxj/HftEor6Fjj7QJExVtL2\n3k0xXUrZ+P3PitwUKZf49n+Ev//+m2vXrv2TB+D27dspMiuOiEhMnBIUFAQkJqSpU6eOUQtL2tBo\nNFy8eJGXL18CiTkQatWqRebMysrN6OPp06fcvHnzn2tx//59vv/++2TrhoWFcfr0aYKDg4HEVKYu\nLi7JzsOsVqs5f/483t6J4X779++ndu3aKVJT7vHjx9y+fZuIiAgg0bVVrFixZOuGhoZy5swZQkJC\nADh9+jS1atUyOQ/zB9RqNWfPnv0nlervv/+Oi4sL6dIpK12vjwcPHnD3bmKAwc6dOylXrpxRieq/\nJMIq2QFj4v0r8QchskkpP+RkdQW0pUa8BhR6b6zfAh2BTgbPpHSnU3JfiadKnajVanno0CHZpGF9\nmSlDWunaoJLs1aaWbFK9hMybM6usVL6s3LRpk4yJiTFa+8GDB3LwoAEyfdo0sn6V72WP5tVk65rf\ny6qli8jsWTPLKVMmSx8fH6N1Q0JC5IIFC2Sh/HlkyUK5ZecGFaRbrRKyWfWSMl0aJ9m1Uwd5+fJl\no3UTEhLk/v37ZQOXWjJLhrSybY2Sslf9srJx6XwyV5YMslrF8nLr1q1SpVIZrX3nzh3Zv09vmT6N\nk2xYupDsUbOEbFkqj6xYOLfMlS2LnDF9mnz79q3RuoGBgfKXOXNk7ixZZCFnZ1nbyUlWBlkxTRrp\nbGcnu3XsKK9du2a0bnx8vPTw8JDVK1SQGeztZSVnZ1nTwUEWA5nB3l5WLVdO7tq1S8bFxRmtfePG\nDdmjc2fpbGcnK6RJI+s5OsryIIs4O8ucmTLJObNmyYCAAKN1/f395ayffpI5M2WSxdKkkQ3SOMsa\nQsgqadPKNHZ2snfXrvLWrVtG66pUKrl9+3ZZrWwZmdXRUTbIlEHWs7OWDTJlkJkcHKRL5cpy7969\nMj4+3mhtbby3F8m2OS/rl1b00nY+YDuJ0RAqEusI9AI2A3eB28DvQNb3fbMDhz86tjHwBHgKTFAy\n3v+8myI6OpounTrwwus+P3SsRfsG5bD/KGeFWq3h6MUHLN19npBo+OPoCbJnV5b/eMWK5fw0bSoD\nWlaib/NK5Mry71nJgxd+rD7giceZe2zdvlPxCvTt27dp0awxtYrnYnCzslQplvtfizwh4dG4n7zD\n0kPX6NClB7/MnafI/xgREUF719YEeT9jeMMStKtcGNuPVuUT1BoO3/ybZScfEGOdhkNHjyuagUsp\nWbhgPgt+nsPgmoXpU7UQ2dL+O9vcHZ9g1lx6zu93fdm9L3FGq4SrV6/SolEjvlepqBMTw6fRz+HA\nOQsLTtrZMXDECGbMmqVoQezdu3e0aNSIgIcPqR4ZSSn+/RipJvEbedHJiQzffccfJ0+SPn16g7pS\nSmZOn87yBQuoHxtLLY2GtJ/0+ZvEDSR3bW05cPQoVatWNagLcOHCBVybNaOiOoGmKhWFPglxC5WS\n48KCw1bWjJo4kQmTJim6FsHBwbRo0ADVi79pSzw1bK3/la4zTkrOxsazR1iRtURJ9h05kiy/NKSc\nm+JVQ2VPtXlP3E72+ZLLf9oYx8fH06xJQzLZxrJhamdsdIQDQeKXaPb642w7eZ9LV66SMWNGvdpL\nly5h2cJfODq/NwVy6O974e7ftJuyje27PKhfv77evvfv36euSy2WDWyAWy39LoPg8GhazdhNxTpN\nWLJUf8rP2NhYGtZ1obCjilW9XLDSk1pUSsnk3Vc4/CiY85c9DT76/jxnNltXL+OPAbXJnUF/voK/\nHr+hy+ZL7Dt42ODuwRs3btCgdm16R0VhKKI1jMRY4DaDBjFXR5rLD0RFRVGjUiXSP3tG67g4vTuj\nNMBBGxsC8ufn0vXrODnp//0mT5jAzmXLGBUdjSHTfQv4zcGB42fOULFiRb19L1++TLP69RkTF0t5\nAzfeICmZbmtHt1GjmTpjht6+4eHh1KhQgVKBfgw0kLw+QUp+VWl4k7cApy9fxv590VJTSClj7N1Y\nWWRHnmO3zMb4azLpx4ncvniM/fP7YKWwksH4pQfwCrbgwKEjOvtcv36dFk0bcXHFIPJlV+avPHf7\nOW7TtvPY65lOQ5+QkEDRwgWZ2r4iXeuVVqT7LjKGGmPcmTH3V9zc3HT2G/nDMF5fP8XOoY2wULA9\nVErJ8M3nCHHMyw4P3SV2zp8/T+e2rbg4siE50ynLzXDioS89d1zl2YtXODs7a+2jUqnInzMnHYKD\nURpUFg7MdHBgze7dNGvWTGe/Pt2788jDg06xsYryGktgl50dBVu3xn3HDp39jh8/Ti9XV6ZGRyeZ\nDeviBrA1fXpe+PrqNG5RUVEUyJWToZERVFIYgREsJWNsbNl68BB169bV2a9z27bEnDrJGFsLRbNo\njZTMVEkKtWvPinXrFI1FGylmjJso23iS5+iNr26MU2Q7tBAirRDCQwjxSAjxQAihPfNIKiImJoa1\na9fw66g2ig0xwIwBTbl06RLPn+vOQLZ08UJGt6+h2BAD1CpTkMaVi7Jx40adfQ4fPkzWNDaKDTFA\nOid7ZnWvza8LdC/kRkRE4O7uzqIu1RUZYkj8ssxpX5XjJ07orbe3ZME8JtQrptgQAzQsnpPqBTKz\nRU9NuT179pBVpVJsiAHSAK2io5k/c6bOPsHBwez28KCFQkMMias7LWJj2bd/P4GBgTr7zZ85k1ZG\nGGKA8kCu+Hh268lLvW3bNgqr1YoNMUBGIWivUrFo9mydfd6+fcsfR44wyEaZIYbEMkdDrBPHFJaM\n1J4pxX8xN8US4IiUshhQGtCerToV4eHhQflieSiY27ioAztba3q2qMKqlSu0tgcFBXHo8B/0aGI4\nKc2nDGpVidUrlulM2r5y6WIGNlFuiD/QvHIRvF+95M6dO1rbt27dikuJPOTKqH0Wqgtnexs6VS/K\n2jWrtbb7+vpy6vRpulQyPo/FwKoFWLV0sc7QqSVz51InMtJo3crA3bt3ddbE27hhAyWFwLgrAY5A\naQsLftMxG3z+/Dk3btygipG6AC6RkSzRsR1aSsnyefNoEmN8gs46Ai5cuvRP9MmnrFu9mrr2NjgZ\naawyWVpQyc4Gd3f92f++BP8pYyyEcAZqSik3AkgpE6SUymsIfSU8dm6jRzPT9s73bFEJj93aE88c\nOnSIBpW+I2Na49MlVi6eB2sLDbdu3UrSFhoaymXPa7SroSxu8mOsLC3pVvd7du3aqbXdY/sWelQ3\nLWC/Z83v8Ni5XWvb/v37aVk6L84Gqgpro+532Yl4F8pjLZUz3rx5g5eXF6ZksrAGquqZaW7fsIHy\nJiZFrxATw45Nm7S27d27l8pqNaZkCykLvHzxQqvRfP78Of5v31LOhJy8dkJQQwj27duntX3Xls00\nEcbmr0ukMQns3LDBpGNTEmFpqeiVGkiJmXEBIEgIsfF9FqO1QgjTPfdfCH9/f/Jm17+wpou82TIS\nEBSsWzeLaSvJQgjyZMtAQEBAkrbAwECyZkyLrY1psbi5M6ch0F97yfqAgADyZjJtzHkyORMYFKK1\nzd/fnzxpTKtvIYQgVwZnrY/9AQEBZLS1xdSvUIaEBPx9fbW2BQYHK64KnUT3/fHa8PP1JYORxWE/\nYG+M0MEAACAASURBVAFktrXV+rkICAggi421yRWOM8Wp8H+bNPUoQGBICNlMLNaZzdKCwCDdLpsv\nxbc0M06JTR9WQDlgiJTyuhDiV2ACMO3TjtOnT//n/y4uLri4uKTA6b8OCfHxOv1oP7RXnhz9U8Le\nvaNp06Za27IY4XfVxrZt2/htw6Yk7zubWOfsA1FRkTqvxeQmpm+YCQsL0xnilslk1URWrV7N8tVJ\n3SvJKY0EiQVCtV0LS6B9MnTDwsJ0RlTkBjDxJg0wd948fpk3L8n7tgCZkxeippQzZ85w5syZFNdN\nLVU8lJASM2Mf4LWU8vr7n/eA9ifI6dOn//P62oY4a9asvHqrfRZjiFd+weTInk1r4PaGDRt4E2Ja\ngR0pJWGxGq5fv55ENyQkhOg4Dao440q9f+B1YDgjRo3WOuYK5cvzKsg0z5J3UAT58ubRqrts2TJe\nR8SZpCulJDxB8PDhwyS6vr6+xNvaYnxti0RCrKyYrqPUUPGiRdE+z1egC3xXuLBW3Tlz5xJi4m5L\nDRBnZ8erV6+S6D59+hSVgwMaEyOVgmxsWbR4sdYxF8yfDz+1aW4KP7WGzJmUr8e4uLj8yz6kFP+p\nSh9SSn/gtRDiw37IesDXLwtrALeOXXD/44ZJx246dBW39h20trVo0YKTV58QHGZ8JWDPh97Eayy0\nZr1Knz49VStXZM+FB0brJqjVbDl1nw4dtFd9duvcDfeLz4zWBdh0/gluHTtrbWvTpg0H77wiItb4\nx/NTT96SJn0GrYl4cuTIQZEiRbhptGpi2svL1ta0b699ntq5d29umBgfe93enk46Sju1bdsWT0tL\nTLk13QLy5c9Pnjx5krQVLFiQrNmzc9MEYxwrJRekxNXVVWt7h27dOSpNMxHHsKJj794mHZuSfEtu\nipS6JQwHtgkhbpMYTTEnhXQ/G+3bt+fGI2+evzbOrxUTG8emQ1cYNHiI1vZMmTLRonkz3I9e19qu\nj1UHrjJwyDCdu+UGDx/J6qPaIyL0cdjTizx581G6tPZIjK5du3LmgTc+wRFG6UbExLHj4mP6Dxio\ntT1nzpzUrVOHbVcNFyL9lNWX/2bQ8JG6XUHjx3PawAYLbXgCpUqV0plLoVfv3tyTEuOuBEQBdzQa\n+vbrp7W9YMGCVKhQgStG6gKccXLih/Hai0UIIRg6bhzH7I13YZ2WUKNaNa1GHqDfwIGciokzuihp\nkFrD1dg4evToYfSYUhwhlL1SASlijKWUd6SUFaWUZaSUrvKT/J6pETs7u/9j76yjo7q6KP57E1cC\nwQnuDsEpEixocC1WLLgUtyKltLgTILi7BQ/S4gkWrFCChCRIPMSTyczc748QykfGEwpddK/1VuXe\nOXO5vDnvvHPP2Rt398GMWXIYhQFijrM8T1G3bl2tskOjfhzH4n1XCHijfxrk0t3nnPb9i379+mmc\n06ZNG0Jj5ew4r79Djo5LYvq2i4wZr1n5xc7Ojr59+zJ251VUev7whBBM3Xed5q6uODk5aZw3evxE\n5p1/zOt3+r8peD96zdUX4VpJ/Dt37kyohQU39Laa1vRx1NqaCT/9pHGOo6MjXbt0wcvSEn1dkACO\nWVrSsUMHre3h46dP56i1NYb8OG4Br7RE8gA9e/bE38SEGwbILkUKwT4LC8ZOm6ZxTr58+WjdqhVr\n5Prr4KmEYFVq2pqyZTOkovrz4FuMjP+VmDX7Z1LNHPhh9k7kqdpzsUIIftlwGq+rT9m0RXMzAkD1\n6tWZMn0GLSZs4oUeDvnK/Rd0nbmL3Xv3a22zNjU15ejxk4zfeIH9l9SRRf0/ImMTcZu1l2ZtOmjt\nvgP4bf5CwoQdgzf9jkJHnlAIwbR9PlwKiGPtBu3lS/Xr12fYmHG0XPMHwVG664LP//WGPjuuc/CI\nl8buOwALCwuOe3uzzcYGfZJNMcAia2t6Dh2qtfsOYMWaNaQUL85hc3OdwqQq4Ki5ObGFC+vsOGve\nvDn9R45kkbU10Xqs2Q/YbG3NsTNntLYW29jYcOTkSZaaW3JbD4ccIdLknIaMG6+1+w5g7ebNPM+V\nlzUpKp15aYUQLElREVekGItWrNC5jn8CMhMTva6vAd90OzT8P1HQqG716eZaLQNR0MmrD1m57yrR\nSXD85BmDiYLc3WoyyC0jUdDDF29Ze/QGBy4aTxQ0tFVV6pT7f6KgyNhEtp69y8pjt+je6wd+mzff\nIKKg8KCnjGpWns61SmFp/v9EQcduP2fl2T9JsXDA6+RpvYmClixexMLf5jK0XhpRUL5PiILuBqcR\nBR19YBxRUPmUFBprIAq6KJNxztKSoT/+yKw5c/QmCmrbogUhf/7Jd/HxVEYzUZBjmTIc9/bWmyjo\nl9mzWblwIU2Sk2moiSjIyor75uZ4nT5N7dr6tYroIgqKek8UdMLUjHFTpzJp6lSDiYI6kkp9NURB\nfySnclAyJW+FShw8ceKrIQqK6Kte+fpT5Nx6/ou3Q3/zzhjSOIFPnTqFx8rl3Lh5g/rOpXGws+RN\naAR/BUaSN19+ho0cQ9euXbG0zKiRpw2PHj3CY/VKdu3chXOZQhTIace7mFjCYuUEhcUw0H0wgwcP\nIX/+/AbZjY6OZvPmzaxZtRwrE0GFonlQpCSRqJRx9UEAbm3aMGzkaL1/yOlQKpUcP34cj+VL8bvr\nR4PyhbG3NOVNeDR/vo2lcJFiDBv9I507dzZY8PT+/ft4rFzBvn17qVE0L3ntLIiOieNtsiA0Xs6g\nIcMY5O5O3rx5DbIbGRnJxg0bWLVkCRZJSRQQgqT4eIS9PX/J5XTo0IFR48YZLJCpUCg4evQoy+bP\n59HDh5QwNcVCqSQ8MZEwKyvKlCvHmEmTaN++vcG81H5+fqxYvJhDBw9Sytwce6WSdwkJxNnZkWRp\nyfAxYxgwaJDBvNRhYWFs8PTEY/ly7OVynBAkxcWjsLfnkTyFLp27MHLcOI3nB5ogl8s5dOgQqxYu\n4OkTfypbWaCKj0OyteNuYjIVK1Vi5MSJuLm5ZZqXGrLOGUf2a6bXXMfNZ/9zxl8bXrx4wa1bt4iN\njWXQoEH4+fllCbl8fHw8Fy5cIDIykv79++Pt7Y2Li0uWkcsHBgbSu3fvD+TyOXNmthI3jVzez8+P\nmJgY3N3defDgQZaQy8fGxn4glx8wYADnzp2jYcOGmf4Rq1SqD+Tyffr04dChQ1lKLn/v3j1iY2Nx\nd3fn0aNHWU4uP3DgQC5cuJCl5PKvX7+mT58+HD58mEaNGmVJHvfPP//kwYMH9OjRgz179lC1atUs\nJ5fPKmcc1V+/t80cm7z/c8ZfM/5tcjX/SQ19fruf0/Z/djPYzrQzjh7YXK+52Tec+eLO+Js+wPsU\nL1++ZMrkSTRpWI8aVSuRzcaSIYMHcffu3UzZTUhIwNPTkzYtmlG7ehWy2VjQub0bhw4dQqEwrokD\n0vKP58+fp1f3LtSvXZ1s1uY0b9yAhQsXfJA0MhbPnz9nwrixNKlXhxqVK5DNyowRQ4fw8KHug0Nt\niIuLw8PDg5ZNXKhVpSLZLM3o2r4tXl5eKJXGtnGkRcVnzpyha9u21K5YEWugaZ06LFu6lHfv3mVq\nzU+ePGHU8OHUdXamYqlSWAAjhw1Ty5thCGJiYlixYgWN69alSpkyWAKd2rTh1KlTGsmi9IFSqeTE\niRN0at2amhUqYAU0++47Vq1a9UFCy1g8evSIkYMHU9/ZGWugQbVqjB05kqdPn2bK7mfDv6i07b/I\nmDQ9s4njfuTa9Wv0bl6dFjVLks3GkoRkOdceBrL++E0KFirM3HkLadSokd52ExMTmTZ1Mtu2buW7\nCoXp5VKOAo52KFWCZ2+i2HT2AS/DYhk9dhzjxo03SA1427at/PrzLMxUctwblqFi4ZyYm8oIj03i\nwM2XHLv9jA7t2rNgyVKD8o73799n0tjR3Lp1mz61itGsdF7sLc2IT1Fw5UU4G64/o0TJUvy6cAn1\n6tXT225cXBxTJk5g186duJTIQ/eyechrZ4lCJXgaEcem+28JSVQwduIkRo0eY5A8/Yb16/l15kzM\nExJonphAIUnCFHiH4IqlFTcVSjp37syiFSv0OmRLx+3btxk7fDj37t+ncmoqhRUKLHivwWNqyj0z\nMypUqMDiVauoWVN/Ms93794xfvRo9u3bR1mZjEqJidiTdigYCtywtSXVxoYpM2YwZOhQg/bCY9Uq\n5s2Zg01yMvXj4shPWsQVA/ja2PBQpaJHjx4sWLrUoEO269evM2HECJ48fkxzlZIKSgVWSCQiuGdq\nxllJokpVZxatWqW2aclQZFVk/G5wS73mOqw79cUj42/eGV+9epWO7dsy5fsGDHSrhbVlxkMphULJ\n4csPGbXci4VLltGnj+5i9ujoaFq6NqWIveC3HxpROI96NYx7z0MYseYMhctWZdvO3TrzpkIIpk6e\nxKHd21jXvyH1yxRQ+2MNj01kntdtjv8ZytkLFylSpIjONZ8/f54eXToxs3l5+tYugZUavoNUpYoD\ndwIYe9iPVevWa61//bCW8HCaN3ahglUqsxuVpsAnlRTpuPM6ihEn/6Ri/cZs3LJN58NJCMHoYcM4\nsX07I1KSKSup5yKIFoLdJqb8lScv569c0VoXnY6TJ0/yfZcuNEhMzCC5lA4l8AD4w9qa7Xv24Obm\nptPumzdvaPzdd+R584bmcjnq3KEAAoED1ta06NGDNevX63TIKpUK9379uHjgAH0TEzNUlaQjGjho\nYUFIwYJcuHKFPHny6FzzoUOHGNS7N/1TkmkokzBTsxa5EJwTsMPCkt2HDuldGaQJWeaMh6rnefkU\nDmtOZvg+SZI2Am2AUPFeAVqSpAWAG2nP5OdAP3UslZIkvSTtGagCUoUQOp/W37Qzfvz4MS4N6rFl\nciea18rYdpth/stQmo7dwKatO2jZUvMTVy6X49qkERVzwbLBzXX+kJJSUukwZz/FqtRjraf2WtWF\nCxawfe1yzk9th6Od7rbd5af8WHflJddu3NJ6kHX37l1cG7uwp29dGpbSXbp371UkLdf8zu4Dh7XW\nqiYlJdGoXl3qZ0vllybldO5FglxBu12+VG/ZkSXLtdeqzpw+nf3LlzEnJRlbPaLHvUj4FizI9Tt+\nWqNCX19fmjduTOfERHS7bXhNmuM8cfYsdevW1TgvPj6e2s7OFA0IwFWh0ElenwR4WlvTZcQIftXA\nZ5yO8WPGcGr9esYnJqKr3kcAB01NeVGyJFdu3sTGRnP33qVLl+jQsiWz5SmU1KM54qFKMNfcgrOX\nLhlcvfIxssoZxwzTzxln81DrjOsB8cC2j5xxU+CCEEIlSdI80oRMp6j57hdANSGEPiXlwDeeM54w\ndjRTezbQyxEDlC2Shx3TuzFsyCCt+c0dO3YgS4xgqburXq+YVhZm7JvSkVPHj3DrluY26tDQUH79\nZQ7HxrXWyxEDjG5ZlbpF7Fm0MCMr18cYO3IYv7SqqJcjBqjs5Mj67rUY7j5Q6yHOhg0byJ4SrZcj\nBrAxN2V/txrs2raVR480U5wEBgayfPFiftLTEQN0FSryvH3LimXLtM4b0q8fTfR0xAAFgKaJibj3\n7at1L1auWIFlcLBejhjACuifmIjHihW8ePFC47wnT56wad06RuvhiCFNmaSTQoH5y5esU8Nclw4h\nBEP69mWEPFkvRwxQQSbRPyWZEQMG6DX/s0Mm6XepgRDiCvx/f44Q4pwQIj2h7wMabxMJA/3rN+uM\nAwIC8PH1ZUAbwxSiGjmXIJedJadPn1Y7LoRg9YqljOtY06AcsL2NJUNaVsVjleZocMP69XSqVZKC\nOQ3ToZjQuiqbNqwnJSVF7fjjx495/OgRvWsZpsjRqoITJqlJGqkPhRCsWb6UcXWKGkRlmN3KnP7O\nhfBYqXkv1q5eTWMJchhgV5IkOsnlrFmxQuPB6a1bt3gVGIihFP7lgPA3b7hxQ32DtlKpxGPZMhoZ\nIOcEYAfUUKnwWLlS45zVy5fTQKEwSJ1EAlomJbFq8WKNh4VXrlwhMSKCOgYecDWSSbx46s/9+/cN\n+tznwGc+v+sPnNIwJoAzkiTdlCRJPWHJJ/hmnfHaNR70aV5NbY5YF4a0rY7HSvXR1c2bN4mOCKN5\nNcOVM/o1r8rhw4eJjs74ZqNUKlm3ZhVDmpQz2G7p/Dmo4OSoUdHBY9VK+tcpjrkBWoCQ5tyG1C2K\nx/KlascvXrwIKQnUL2JY4wLAwOpF2LVrJ/FqpJXkcjkb162jlcJwNrgSMgnH1FROnDihdnzl0qVU\nTk42+IchAZWTk1muQXn6zJkzWCQnU9hAuwB15HI2b9ig9mGamJjItq1baWxEVU4pQIqL4/z582rH\nVy5cSMvkJIM5gU0lieYKBauWqr8v/lFo8L6XX0fyq6//h8tws9I00nLB6mVuoK4QojrQChj+PuWh\nFd+sM75x/SotaxlXqN6qTll8fG+qHfP19aVF9eIGRcXpyO1gS/mi+dRGFCEhIaSmpOBcVPeBizq0\nrJgfn2tX1Y7duHaZlmUN6wD8YLe8E743fNWO+fr60qKYo1EE307ZrCmcw15t+VhgYCAWKhUFjQxp\nqibE43Ptmtqx61euUMLIsrISKhXXNdj18fGhZJyhXHBpyA1YSZJaEVx/f39ymJpi+OMu7QFSMTER\nX18Nf38+PtQ0MmqsrlLhc+mScR/OSsjUX/UL5WRq3dIfLkMgSVJf0pyseu5YQAgR8v6f4cBh0K2d\n+80645jYWLLZGsdb62BrRawGMcyYmBiyWRnfVZfNxlJtLWhMTAzZbIxXs7K3Mic2Rv1ZQkxsnNFr\nzmZlTmyceka2d++iyWZufCdZNiszjXthk4kONRvgnYY67Lj4eL3yrupgCcQnqN+L6IgIMqNFZi2T\nqd2L2NhYbDLxnm2lUvEuSj2dflxiIurrXnTDVoI4DXvxT0KSJL0ubSbeX+n2WgATgbZCCLV5P0mS\nrCVJsn3/7zaAK6CzQD8rZJf+lbCxtiYx2TgVioRkOTINZVQAg1sbrgydjjeh4bRt21btWA5bY90E\nJKaksnvPXrZsz/hWZW9pRqLccNVpgIQUBanyZI17Ma6+YVHHx3gTFkmTJuqJXnKA0VJDycDWjRtZ\nu3FjhjFz0gjojUEqEBsTo3EvtHPFaUd4TAx16tRRO5Y7E3aTgcVLl7JYTUrBEkg2M+6hlyzA2kiS\n/qxEZugxJUnaBbgAjpIkBZEmJTeVtNvk7Pu/Zx8hxDBJkvIB64UQbYA8wGFJkgRpPnanEMJb1/d9\ns5FxmXIVuP4w0KjPXnvwEucqldVK1Zw6dYo7AcZ1vyUmywmOTCQ4ODiD3eTkZDA1JzDcuA6q6y8i\nWbRkmdo1u7Vrx/UXGcUu9YFPQBi1a9ZSa/fgwYPcDDNOaTkmWU54ioqwsLAMduPi4lBYWRFlZKnk\nU2sbNmzapHbNzZo1I9goqxAMNGzYUK3d7du388oIMnyAOEBhaUlMTEwGu5GRkSRZWKCbnFQ9Xtra\nsmfPHrVrrlunDo8NJJZPxyMEZcsZfr6R5chcNcX3Qoj8QggLIUQhIcRmIURJIURhIYTz+2vY+7lv\n3ztihBABIo3bvaoQoqIQYp5eS82yP/S/DEOGDcfz+E2URmh8rT12i8HDRqoda9asGeFxcm75q1cf\n1oY9Fx9St24dtU0JFhYW9OrVm/UXDJddCn2XwJl7AfTq1UvtuPuwkay7HmAUz8Da6y8ZPHK02rE2\nbdrwNCqBx2GGaw3s8AukebNmarsHbW1t6dy5M2ckw2/fUCF4qFRqbFYZMXYs921t9SaW/xj37ewY\nNX682rFOnToRBEQYYddXJqND+/Zqa6Nz5MhBq1atuGxEquINECxJtG/fXu348PHjOa2lBlkThBCc\ntrJh2LhxBn82yyHpeX0F+GadcbVq1cidNz/Hrhrm3J4Gh3PjUSDduqnXwDMxMWHwsOGs9DJMX0+p\nVOFxwo9hI8donDNsxEg2XXxEgoGacp4XHtK5c2ccHNR3AdavXx9TGzvOPn5jkN0Hr6N4Ehar8cds\nbm7OQPchrLrx0iC7CqWKdXdeMWzMWI1zRo4dy2lTM+QGPkBOyEzo3aePxkYHV1dXVDY2BkfHr4AU\nCwuNzUBWVlb07dePKway9KUCPlZWjNLi2EaNH88Fa2sMrac4Z27OwMGDsbBQr4ndtm1bQk1MeWpg\ndHxfCCR7e4OoAz4XsiBn/I/hm3XGAHPnLWT4Mi/8g/TTwXsXl0Snn3YwfcYMrcoLgwcP4Zp/KBtP\n6yeZKYRg3PqzOOQpqLWNtHTp0rRu246ea87qVONIx9n7gXicf8SkKZrldSRJ4teFSxi0x5eXeurg\nRcQn023LVWb/8qtWXuMRo0ZxJuAde+4F6WVXCMHIE/cpXr4S9evX1zivSpUqNGjalKXmFij1dMg+\nKhV/WFoxfkqGhqkPkMlkLFi2jGPW1uibEIoDvKytmb90qVbqy/GTJvHA1hZ9RbNUwD5LS2o3aED1\n6prPIerUqUPlOnXYbGmpU5kkHT7APVtbRo/V/MAzNTXlt8WLmW9uwTs99zhcCJaYWzJ/xYqvw8n9\nFxn/O9C0aVPmzltA4zGeXL6nucMJ0iLiBiPX0rxtJ0aN0hy9QpqS88kzZ5m95xoL918lVYvGXnxS\nCoNXnOTK8xgOHT2msyRuzbr1KLM50W7JSSJiNedjVSrBzit/0WvNOQ4cPkqJEtrrnt3c3Jgy82dc\nlp/DN0B7/vjx23c0WH6Wbj8MZJC7u9a5uXPn5oT3WSZd8Ge1z3OUWsrGYpNT6X/Ej0cKG/YcPKzz\nx7xt715UFSoyz9yCOC3OQiUEZ1SClZbWHPP21ijAmY7u3bszfvp0dlhb81brzDRin+3W1oyeNElj\nGigdBQoU4OTZsxyys+OqJGl1nInAdisrZBUqsOvAAa12JUli35EjJJYpwzpLSxK1zFUBZyWJXfb2\nnDx/XqdqzQ/9+jFw7FgmmFnwUkeE7K8STDCzYPLPP2tUnP6n8W/SwPumuSnS4eXlxYhhQyjgaMsQ\nt+q0rF0GB1sr4pPkXHv4kjVeN/F5+JKfZs5g9GjNisWfIigoiD49u/PU/wmDWlShV+NKFHC0R6lS\n8exNFBu977Hzwj1aNG/Bug2btGq+fYzU1FTGjx3D1q1baVu9BIMblaVioZxYmJoQHpvEfp+nrP39\nERa2Dmzbtccgcvx9+/bx48jhFM1hw9A6RXAt5/SetS2Vy89CWXM9gDtBEcyZ+ytDhg7T2+7z58/p\n3a0Lr4NeMtC5EN0rOv0fa9tGv2D23Q+iXdu2eKzfiLW1fkVVKSkpjBw8mD179lDfxITm8hQKSbxn\nbYOLSJw2t8CxQAF2HjpE+fL699Zt27aNsSNHklOlolJ8PEXhA2vbS9Iiy3BJYuGyZfQ3QJb+r7/+\nonuHDoQGB1M7KQlnlQo7/mZt87W0xE8IOnfpgsf69XqryyQlJTGkf38OHz5MHUmiYXIy+fh7L67L\nZPxuaUmBokXZffgwJUuW1HvNnuvWMXXCBIoBLZMSqChJWJP20PBTCU5a2/DWxISlHh706NFDb7ua\nkFXcFInT9XsoWP9y6Iuztqk9Rf0cV9pXfb1ITU0VR44cEa5NXIRDNjshk8mEuamJqFmtiti4caNI\nSEgw2va9e/fE4EEDRZ5cjsLU1ESYyCRRpGABMXXKZBEYGGi03fDwcDF//jxRpkRRYW1lKSQJ4Zg9\nm+jRpZO4fPmyUKlURtmVy+XiwIEDonH9uiKbnY2QJEmYmchEnWpVxdatW0VSUpLRa759+7YY0LeP\nyJ3DQZiamAhTmSSKFcwvZvw0Xbx69cpouyEhIWLunDmieP78wtLMTEggctjaip6dOwsfHx+j9yIl\nJUXs3r1b1K5aVdhYWgpJkoQJiJqVK4tdu3aJ5ORko+yqVCrh6+srenbtKnLY2QkTmUzIQBTOm1f8\nPHu2ePv2rVF2hRDi9evXYuZPP4nCefIIC1PTtL2wsxN9evQQN2/eNNpuUlKS2L59u6hVqZKws7QU\nEgg7KytRr1o1sW/fPiGXy422/Sne+4tM+5yknzrqdWXF92X2+i8y1gAhBDKZzKgKgy9hFz6f8sJ/\ne/E3/tuLz2/3I9uZjoyTZnbSa67V7INfPDL+Zps+NCE1NZWgoKAP3U6JiYl6vzLrwuvXrz8ocERE\nRGSJTh2kkZW/fZuW3QwMDKRQoUJZcngil8sJCgoi7n0bb1JSktaDS30hhODVq1dEve/8ioyMxNHR\nMdN2PydSUlIICgr6wJWRnJxssDitOgghCA4O/sBHEhUVlSWafem20u+L4OBgnJycsuS+SE5OJjg4\nrd7k2bNnFCxYUGNFxpfG13CGqC++6QO8jxEcHMz06dMo5JSfpi716NezM3kc7XAqkI+hQwbz4MED\no+wmJiayefNmajhXpkrF8vTs3JYCubJRolgRWjRrbLTckBCCixcv0q1zBwoXdKJDq6YUzuNAnRrO\nlClZjKWZkBt6+fIlkydNoGD+vLg2rMsPXdqSO5s1TvnyMnK48XJD6fJTzuXLUq1COXq2aU5+W0uK\nFy5IG9emmZYb+hx49uwZY0eNIn/OnDR2dqZ7o0Y4APlz5uTHTMgNxcbGsnr1asoVLUrVMmXo2KAB\n2YHC+fPj1qwZZ8+eNWovVCoVp0+fpnXjxmm26tYlF1CldGkqFi/O2rVr1ZIv6YMnT54wathQ8ufM\nSdPq1cknk9G4mjMFcuViwtgftdJ8fjH8J7uk5ou+0jSFEIKff57FimXL+L5lDYZ0qkvZon+fML8O\ne8fGI9dYf+Q6TZq6sn7jZr2jgEuXLtG1c0eql3FiiFsNmtcsjYlJ2vMvOSWV/b/fY+2xW7xLVnHs\nxGmdFQ/pCA8Pp2O7NkSFvmZIiyr0alqZbDaWH/481x4FsebkXU7f9Gf1mnV6H6ioVCqmTp7EBs91\n9GlQlsGNy1My398yRcERcaz//U82/vGINu3a47HWU291a29vb3p260qdAjkYVD4vLoVzIXv/ZSHO\nzAAAIABJREFUI0hMVXDwr9d4/hmCwsqeY2e8KVzYGH6zv5HZV2ilUsnYkSPZsWULTYWKlkoF+T76\n0YYIwWkTU7wlGd179WLFmjV6KzofP36c3t27Uw5olJBAOf6urkoGrgG/29pi6+TEyfPnyZ9fPxKn\n4OBgWjdpgjwkhJaJCbjIJCzfr1kIwV0hOGllzUOVYOf+/VoFEj6GQqFghLs7B/fupbWZjLbmMvKa\n/B3HvVIoOaqA0ykK+rkPYsGSpUYRZX2MrEpTpPzcWa+5FjMOfPE0xTftjIUQjBk9iqsXTuK1ZBB5\nc2qWMU9KltNn5k4SJHuOnTil0wmdPXuWnj26smN6N5pW184Ot+6oD3O2/8Efl67olDyPiIigXp1a\ndKhRiDl9XLTe9PdfhNB29j5mzPmVgQO1U6oKIXAf2J9HPn9weExLctprTkckJKfSfZU3ZnmKceCI\nl06pKC8vLwb16cX2VpX5rqD2dMTqOwGsfBDCZR9fvaSiPkV8fDw3btygSZMmPH/+nGLFihlsQ6VS\n0bNLF/y9zzBVB3l9ghD8Zm5JkcaN2XvkiE4ntHfvXob368fopCS0PXoFcNzUlCs5cnD9zh0KFCig\n1W5QUBB1qlWjTVwsHVRKremIRyrBXHNzPLfv0FmCplQq6dLWjdBrV5llLmGjpQwsVqViagpUbNmK\nLbt2ZSolkmXO+Jcues21mL7/izvjbzpNsWHDes6fPor36mFaHTGAlaU5u+f2wVQeydgf1bf/piMg\nIICePbqxf3ZPnY4YYHC72szq24hWLVzTOCg0QAhB5w5taV+9EHN/aKzzh1+pWF68537P9MkTuaaB\n2jEdK5Yv486lc5ya2EarIwawsTTj4JgWxL9+yrQpk7XOffLkCQP69GJ/26o6HTHAcOeijKqYj9bN\nmpKaqn+nYWxsLKNHjaJw4cLM+GkaNWs4U6dObZo1baqRIlITfpszh4dnzjBTDxURG0lihjwZ/wsX\n+HnGDK1z7927x9D+/RmvwxFDWqTsplBQLyqKlo0ba01ZKBQKWjZuRJvYGDoKlU4nWE4mMUsuZ2Dv\n3vz5p/YO1BlTphB89Sq/Wmh3xAD2MhkLLSRunzzBonl60TF8dvyb6oyzzBlLkiSTJOmOJEleWWXz\nc0KlUjHv119YP607Dnb6HdCZmpqwecb37Nixg4gIzSwDq1au4IcWztSvrH9UNtCtFsXz2rN//36N\nc65du8aboAB+6euit91STjmZ1asBC3+bq3FOamoqC+fPY5N7I2z1JNs3NzVh6+DGrPf01Cr/vnTh\nAtwrFaRaPv1VmYc5F8FBlYyXl363UmxsLC4uDYmPjeKuz+9cOX8c30veBPnfpVf3Dri5teHcuXN6\n2UpOTmbp4sX8KE/GQs/IzlySGCtPZsWyZSQmam65WPDLL7Q0kGC+jUJBytu3nDlzRuOcEydOIAsL\np4PQP8dcUibRVpHK4t9+0zgnPj6e1atXM9UctSKk6mAlk5hqDosWzEcuN44VMUvxL8oZZ2VkPBrQ\nLFr2lcHb25vsNubUrFDEoM85OtjSzqUymzZtUjuelJTEtq1bGNzWMDkngGHtamhUEAHwWLmcoa2q\nGpyP69m4EpcuX/5wAv4pvLy8KJbTloqFDKMoz5fdlqaVCrNt2za147Gxsezdu5d+FfRVk/sbg8rn\nZfXSxXrNnTZ1KlUrlWfDmmUULPj367yFhQV9e3XnwK5NfP/99yQl6WaQ279/P8UlKGDgDzSvJFFG\nJmPv3r1qxyMiIjh2/DgNDTyUkwCXuDiWL9CsYbhi/nxaJiYYnBZoIVQcPHhQ40Hvjh07cLY0I4+J\nYfdbEVMTCksShw8fNuhznwOSTL/ra0CWLEOSJCfSmO83ZIW9fwJbNq1nYPtaRuW13DvUYcsm9SrO\nx44dw7l0QYrmN7xUq1XtsrwKDlJbrRAfH8+x4yfo01T/brp02Fia092lokanuWX9WgY2NE71xN2l\nDFvWqxe1PHjwIA2K5CafnuKpH6Ndqfw8fPCAwEDtNKdxcXHs3LWL2dMnafy7bFCvLtWqVtL61pGO\nDStW4JpoHCl686QENqxQr9u3d+9eqkqSQTp16agDXLt+nbCwjG3qr1+/5rafH/WNeNXOLklUMzXV\nuC+bVq3CTRgu5wTQRiVnoxbdvn8M32BkvBSYAEYxD34RBAUGUr6YfkrIn6Jc0bwEv1LPcBYUFET5\nwsbVD5uYyChdOK/aCDYsLAzHbLZkN8KxAZQrmIPgwAC1Y8HBwZR3Mm7N5ZwcCX6teS/K2BuuMQhg\nbiKjWC4HXr16pXWer68vFcqVwclJe8VBl45t8dbyqp+O4FevKGzkj7OQJPHqtXrq1MCAAPLpEZmr\ngwWQ28KCN28y7vOrV68oYGGBuZFrLpCUSJCGB17wm9cUNVAXMR1FTU00von9o/gXEQVluulDkqTW\nQKgQ4q4kSS5o+aPNmjXrw7+7uLjg4uKS2a83GnK5HHMjVQzMzUxJTNIs1Di+h4vR60pMiNPI3JYv\nh3Hk5ACmJjK2b9vGuvUZ1S3sLM0xNzU8rQJgZiIjJiZW415Mr1fGKLsACbGx1KunU8eR2rV0K6vY\nWFvjdcxL55uQJWBm5H1hBoSGh6v9DhOgu1FW05AQG0vVqlXVjhUCo1VPTIFf587ll7kZzxTMAdNc\nGTmU9bIrQUqq/jnjP/74Q6PKeGbwVTDH6YmsiIy/A9pKkvQC2A00kiRJ7fvwrFmzPlxf0hFDGrNa\nWJRxApFhUXHkcsyhtr/cw8ODyHjjDy7kKlOuXr2awW5ISAjJCoxuioiISWDIsOFq11yxYgXCYrVx\nfWlGeGwS+fPlVWt30aJFhKcY3tCSDoWZJX5+flr7+f39/Ql4GaTzsOjGLT9Gjhipkx+gRLFivDPy\n/e6dgKKFCqm1O2P2bGKNrL0VgLCxwd/fP4PdBw8eIOyMf0jHmpry67x5atfslD8f0UYqfUSrBNmz\nqefPVgcXF5f/8w9ZhkwofUiStFGSpFBJku5/9P+yS5LkLUnSE0mSzkiSpLYMS5KkvpIk+b+f10ev\npRr1B/wIQoipIk2SpBhpD/8LQgi9vvxLonkrN/ad05dZ9v+x7+wdWrRooXbM1dWVY1f+JCnFcCW1\np8HhBIVG4ezsnGEsd+7cFCpUkHN+hnc5CSHYd/UpLVupV2FzbeXGPl/1KQxd2Of7lOYamgeaN2+O\n17NQUo1QU3kYFkOcQkWFChW0zitZsiRly5Zl38EjGufExsaxdeceBgwcqPN7W7Zvz2UdddOacNnU\nlJbt2qkda9GiBbesrPTmG/4YzwALW1uKFy+eYaxs2bKoLCwNJoAHUArBVRNTmjdvrna8ZRs3LiiM\nc8a/C4mWGkQH/lGYyPS71GMz8OnmTAbOCSFKAxeADOTYkiRlB2YANYBawExNTvtjfCXniP88+vXr\nx/HLD4iINqw1VKVSsfbQNYaNGKV2vHjx4lSvXo39vxvu6Nd5+dKv/wC1nAeSJDFs5BjWnLxrsN1r\nj4JIVslo3Lix2nH3wYPZd/0vYhLVit1qRKpCyfrfH2nciwoVKlCiVGmOPwsxeM2eD18zaOgwnQ0l\nAL/8Mpexk2Zw7fqNDGNxcXF07P4D3bp206sBZOiIEZxHItnABqUUITiLxLBR6veiRo0a5CpQQG9i\n+Y/xu5UVI8aOVVtFY2JiwtBRozhlBDeErxAUKlZMI8Xq8DFjOJ6qItXAvUhQCc4mpTJk+HCD15Tl\nyEQ5hRDiCvCppHo7YOv7f98KqHviNAe8hRAxQoh3gDegPnr7CFnqjIUQF4UQ6qWNvzI4OjrSvn17\n5m7SfajzMTZ7+eCYKw81atTQOGf4qB9ZsPsycYmaGzg+RcCbSLaduaOVI7hnz55c/TMQn8f6H4wo\nlEpm77rC0BGjNZbE5cuXj+auriw4pp8ySTrWnXtAiZKlqFixosY5I8aOZ+HtQBJT9T+V94+M4/CT\nNwxyH6zX/O+++46tW7fRrmsf3Dr1ZNvOvRw45MWEKbMoVq46pUqXZbmGKodPUbRoUerUrs1hA+ud\njkoyqlevrrGlXZIkxkyezDEbG4PUp4OAO0D/AQM0zhno7s4VpYpAA5xmihDss7Ri9GTNTTtly5al\nfMUKHEwxrKJid4oSF5eGarUc/3FkfTVFbiFEKIAQIgRQVw9aAP5Puev1+/+nFd9sZAywaPFSTt8I\nYPF2/RoCjl26z/Q1J9m6XXurZ+vWrfmuUTM6/bSTeD2izVdh72g1aQuzfp5D0aJFNc6zsbFhy7Yd\ndPrlAPdf6I42FUol7stPYOpQgOE6opTlq9ew93YQa7z1i932+/jz24l7bNy6Q+u8Ll26UKFuQ/qe\nekBSqu788ct3CXT08mPR0uV6czIAtGzZksDAQDp06spJ74v0HzIKUwsbfHx88VizRq8IOx2e27fj\nbWePt56v/ueE4KSNLet3aN+Lvn37UqpBA9ZaWurlkN8CS6ytWbtxo1ZWuzx58rBizRpmmVnwRg+H\nLBeCBeYWVGzSVCdvyaZdu9knmXFOT93Fo8mpnDW3wmOj+jr8fxyZyBlnAuoM6vyL+aa5KSCt/Kpl\n82ZULe7IhN6NqVgy4wMs6G0Uaw5cYdvJWxw9doKaNWvqtKtQKHAfNIA71y8ys29j2tQt94EkKB2J\nyXL2nPdj9pYLjBo7gQkTJuq15j179jBq+FAmd6lD32ZVMpS7CSG4eP8lc/dewyx7AfYfOqKXisjz\n589p6dqE+kWzM7ZVZcoWyOgAAsJiWOV9n/23Ajl+6oxeKiJyuZwfevXkie9lptUoStOiuTH55AcQ\nL1ew7/Frfrvxgqmz5zBipPrXfX2RWaKgv/76C9eGDakSG0M7pQInNQ/fN0Jw1MSU23b2nP79d71U\nRJKSkujesSNPL1+mXUICFcgYESUCVyWJo5ZpWnID9ch1A6zx8GDGhAn0kCfTSJKw/mTNKiG4LQR7\nLK0o06gxO/fv14v06v79+7Rq0oR6qcl0MpMooKbcLVChZG+q4L6NHWf+uGiQiog6ZBU3hcJDvQLL\nH/5vuej/t6jWnJN31X6fJEmFgWNCiErv//sx4CKECJUkKS/wuxCi7Cef6f5+zpD3/732/Tz1HUHp\nn/vWnTGk8QEvX74Mz7VrKFYgBy1ql3ovu5TCtYfBXPF7Rs+evRg/YaJO/bSPIYRgz549rFi6iNev\ngunRuBIFctqhVAmevYlm74W71Kldm7ETJhuspOvn58ei+b9x8tQp2tUtS8VCObAwMyU8JpH9V/2R\nzK0YOmI0gwcP1ptZDdK4hZcuWczG9Z6UyeeAa/l82FuZE5+cyqWn4fg+fUufH/oybvxEneQ1H0Ol\nUrF9+3ZWLlpAZGgInUvmIo+VOUqV4GlsCof931C/fn3GT5mmVzmbLmQF8XlYWBiLFyxgk6cnRSWo\nmJCAFZAEPLSx4bmAfgMGMH7yZPLkyaO3XaVSyaZNm1g+bx7vQkOpmZiIvRAogRALC25IEo0bNWLS\njBnUrl3boDVfu3aN+bNnc/HSJepLkF8uxxSIliQuW1iSI39+xkyZQt++fQ3q5Hz79i2L5s1jy6ZN\nlDU3pZI8GWuZRIJKcMfMkgCVYODgwYydMCFLeLqzzBmv1Zze+RimQzZqcsZFSHPGFd//93wgSggx\nX5KkSUB2IcTkTz6THbgFOJP2rL0FVHufP9a83v+c8d9ITU3Fy8sLX18fYmPesXHjJtau86Rbt27Y\n2hpfPgRw584djhw5QmREGBs2bGDmrJ/p0aOH1rSEPggLC2PXrl28DHjBurVr+HHceFxdm9OwYcNM\n1VjK5XIOHz7M7Vs3iXkXzaZNm1m/cRNdu3bNNNn+zZs38fI6SmRoKBs3bmT23F/5/vvvDXrQ6UJW\nqlCkpKRw8OBB7ty6RUxkJJu2bWPr9u107tw5UwTzQgh8fHw4fuwYESEhbNy8md8WLKBnz54GpWjU\n4dWrV+zatYtXAQGsXbuWiVOn4ta2LTVr1szUfZGUlMT+/fu55+fHqmXLGPXjjzjXqEHHjh2zlGA+\ny5zxOv3eKkwHb8jwfZIk7QJcAEfS5AlnAkeA/UBB0lL6XYQQ7yRJqgYMFkK4v//sD8A00tITvwgh\n1Le/fvx9/zljzfi3ydX8AzI437zdz2n7P7sZbGfeGW/Q7xDYdOC6L06h+Z/s0kdIj4x9fK4TGxOD\nqakJGzZsoHv37lkWGUdEhGFubsqvv/6aJZFxaGhoWmT8MgBLC3OmTJmCq6srLi4uWRIZ37p5g9iY\nd5iayNiyZUumI2MhRFpkfPQokWGhmMsk5s2blyWR8du3b9m5cydBL15gIZOYOmUKrVq35rvvvsvU\nXiQnJ3Pw4EH8bt0i7l00JqSR6GRFZHz9+nVOHD9OxNu3mAIL3kfGhqSA1OHVq1fs3LmTVwEBmAPT\np02jjZsbtWoZx8fyr0UmSe7/SfwXGZOWM162bCnrPddR3CknLb4rRzY7axISU7h67wVXbj/h+549\nGT9+okEKFEIIdu/ezYplS3j79jU92tSnQN4cKJUqngWGsOf4ZWrXqsXY8RM11gBrgp+fHwsXzOPU\nqdO0b9mAiqWLYG5uRnjkOw6cuISQTBg2fKRROeMlixexccN6yhXKSbPKBbG3tiQ+KYXLj9/i8ziY\n3n36Mn6CcTnjFQsXEBX6lk6FHMljZY5CpeJ5QiqHX4RQr149xk+dRv369Q3aixs3brDwlzmcO3+e\n1k45KWlpgrlMIjxFwbHweCyz52D42PEMGjRIbzUOSM8Zz2fThg2Ut7WkljnYmpqQoFDhKxc8jE/m\nh379M5UzjgsLo15yEg7vc8avLCy4ohI0cnFh0syZmcoZN7Y0o6AyFRMkooDfMcUhb16jcsbq8K+I\njDdrLhX9GKb9PL54ZPzNO+P0agrn0nkZ368FFUsVzDAn+G0kHrvPs+2YD0eOHqNWLd08DgqFAnf3\ngdzxvcqsMd1p3ah6BkeQmJTM3uNXmLV8D6PGjDWsmmLkCCYP78EPXVuR3eH/KyWEEFzyucucZdsx\ns87G/gOHDKimaEr9MrkZ26EWZdVQar4MiWbVsZvsu+rPsROnNfIlfAy5XE7fnt/jf/0yUyoVoknB\nnB8kl9KRkKrgwLO3/HY/iKmzf2bkKO0E/unYsnkzE8eMYkzx3HQpmBM7s/9/2RNCcCU8hsUvIshd\noQp7Dh3WK7L/66+/aNG4EQ3MBD/ksaOodcZ86MukFLaExvFHCpw8d17/aooOHQi4eoUeyUlUkaQM\ne5EoBOeFYI+5JfOWLWPgIO0qLelY6+HBzIkT6WMiaG5phrUsYzXFDbmCLSoZpRu4sEPPagpN+Fc4\n4y36NZ6Y/rD6P2f8JREZGUnd2jUZ2KE24/rp1gM79vsd3Gdt4/c/LlGuXDmN84QQuLsPJODJPQ55\nTMLWRjvT2qu3Ebj2mcWIMeMZMWKE1rnHjx9n0MD+nN6xgErltOtFKBQK3Ccu4k1UCsdPnNQaIYeG\nhlKnZnXGta3KUDfNDS3p2H/pIWPWX+Dy1etatfuEEPTq3pWoO75sblgWKx0sYC9jE2nvfY/p8xbq\nbF/et28fP7oPZE/N4pS01+5g5SoVo+8FoSxZgcMnTmqNkIODg6lTzZlRjhZ0yaObX+FwWAyLwhK5\nfvuO1jcnpVJJpzZteHfpEuNSU3Qyrb0SgulmFixev57vv/9e69ytW7YwdeQIllrKcNKxxylC8HOK\nilwNGrHn8GGDImSVSsW5c+fYs3s3Bw8epFu3bvTs1YsGDRpkafojy5zxtpF6zTXts/I/Z/wl0e+H\nPtiLCJZO1n6jf4yNBy+y7tBNbt7203jzHT9+nIljR3L94HzsbPXLrwYEh1CrwwRu3rqjMY+ckJBA\n4UIFObblN2pX0x2FQZpDbtV7Eq079GDMmDEa53Xv0pGiptHM7ad/umTVUV8O3ovk4pXrGufs27eP\nuT+O5HSLKljryYb29F08rifu8PCJv8aqgujoaIoXKsi+WiWo4GCjl125SkUX3+cMmPkL7u7uGue1\nbtqEMgGPGVEwh152Ada8juJegZKc+eOixjmbN29m8aiR/JaSrDflZYBKMMnMnBfBwRobP0JDQylV\ntCgeNqZ6U16mCMHIZMGU1R707NlTr88EBgbSrl1bEIL+fXrgVCA/LwOD2LBlJ9kcHDh06DD58hlH\nS/spsswZb9evXt2094ov7oz/PdntLEZkZCRHjhxh2mA3gz7Xr0N9oiLDuHnzpsY5q1cuZ+LgDno7\nYoCiBfPSp0Nj1q5do3HOzp07qVezkt6OGMDU1JSZY/visXqlRsa3t2/fcsb7LBO71tXbLsDg1tV5\n9tSfBw8eaJyzavFCxpXPr7cjBijpYEv7YnlZv26dxjlbNm/GJW92vR0xgLlMxo/FcrFq8UKNr9cB\nAQH4+PgwIL/+jGMA/fM5cPvWLZ49e6Z2XAjBst9+o1tykkHcw0VlErVMZGzamJH6NB0bPD1xsTQz\niHvYQpLoJSlYMV8/rbro6GiaNGlMnx6d8fO5wKjh7nRs34axo4fx550ruDauj6trMxISjCPm/1z4\nJjXw/m3YvHkzbVycyZndMO0FmUzGkC4N8Vilnuvg+fPn3Lp9m66tDW9cGPx9czZv2qRWlFQIgcfq\nlQztYzj1R93qFbGyMOXChQtqxz0919G1YQWy2RhWGWBmasKgFlXxWKVe0eHBgwc88/enVWH9D7fS\nMaBUXtav8VArSqpSqfBYtpS+BQxzmAANcmcjITKS69fVR/NrV6+iQy57rAyUGrKQyeiU2541GtQt\nbt68SeSbN1Qz4lW+VXIyq5csUfswVSqVrFmxnPYyw/ng6lqYEfQiAD8/P51zPdeto3YNZ8aOHpbh\njVCSJGZNn0ThggXYsX27wev4rDAx0e/6CvDNOuMzp47RtYVuUnJ16NqyFqc1qEZ4e3vj1qQmlhaG\nK1yULJqfQvlzcedORsKesLAwgoKCaVrf8DVLkkTXNg05deqk+jWfPE5XI0ngu9Yvx5nTGux6e9Om\ncC7MDHRsABUc7bE1kdSqFwcGBpIQG0NNR8NFjGSShFtuW06dVL/m015etM5unJpK6+zWnD6mXkT1\n9OnT1JWnYGKEMy4jgTw+nufPn2cYe/z4MaZyOaXNDK9SNZUkGpqk/T3pgud6T0YP15zakSSJ0cPd\n8VzvafA6Piu+Qdmlfx2io6PJncM4FYPcjvZERavvbMyMXYBcjtmIiopSazdXzuxGlyPlcnQgOjqj\n3XTbuR2Mq6PO5WBD9Dv16tBRUZHkMjP+FstpbUF09KcMhmnrzWltZfSBkaOZCdER4WrH3sXE4GiE\nY0uza0p0jIa9CAvDwUhhAEmSyG5mqvG+yG5A6eKncFApiQxXvxfpUCgUvHwZSPVq2itnalZ3xt//\nqdFr+SyQyfS7vgJ8s00f5mbmyA2gdfwY8lQlCKHRGQzvo57EXR+ER0Ti5qY+j50zh05+ao1ITVWy\ne9duNm/emmHM3toCucI4RY5UhZKUZM0SVOOq6uYQ1oSIqGiN9df5rIzT1gNQCME6T09Wr8sYxdmY\nyEgtbNw+pwpBTEy02r2QAX2MeENIR+S7GI01x04mMjBSHilVCBYuXszCxdqVuGUyGUlJSVrLAhMS\nErK0JTpL8JVEvfrg63gkfAEUKlyYP59pF7vUhEfPX1OieFG1UjV79+7laaD2SEMTlEoloZFxPH78\nOIPduLg4UlKVRL8zTirq0dNApv80Q+2aGzRowJ+BGZWH9bIbFE65smXU2t20aRP+ScY5eblSRZRK\n4uXLlxnsRkZGkoCMBCMfIP7JShYuWap2zTWqV8c/wTCS/XQ8S0jBuUpVtXZXrFrFK0vj0h/JQpBk\nYUFoaGgGu69evSLe3JwUIyuVgs0t8fT01ClH1dzVlf2H1Kdg0rH/kBeurs2MWsdnw39piq8fP/Qf\nxPoDl40qWvfcf4kf+qsvxHdzc+POw2cEBBuubnHyj9s4ORWkTJmM+VtbW1vc2rRm6/5TBttNSExi\nz9Hz9OmjXg3rh4GD2eB9X+2YLnievscPA9X3/3fq1InLryJ4m6A/yX46jgWEUKFCRbV1uzly5KBR\nw4YcCo4w2O47uYLTr6M01u32Hz6CvTHGOeM975LpP1x9nXi3bt24oVQSa8T9dlElqFunDrlz584w\nVqBAAapVdeYPPfmGP0aUUsWNZDldunTROXf4iBHMW7SC2Fj1wUBERCRLV61j2LCvQN3jY/yL0hRf\nxyq+AFxdXXkXL+fG/YyHItoQ+S6Oo+dv0b+/ep5UKysr+vTty7pdhimIAHjsOK1Rwghg2PCRrNnu\nZbAo6c5DZ2nQoD4FC2bsLgRo27YtL0Ji9CKs/xhvI+M4d+eZRidvb29Pt27d2OqfUWJeFzY8D2f4\n2HEax4ePHcfW1zEGP0z3BkXQskVztY4N0sjw/4xPJsBACargJDl+sYl069ZN7XjOnDlxa9MGbwN1\n4YUQnLKxYdREzd2ZIydO5IjMzOC9OCZX0qlTJxwcdFeltGrViqZNm9GkVUd8b9z+8F1CCC5fuY5L\n8/b06tkrS+hPsxSZ08D7R/F1rOILQCaTMXnqNAbN3MK7WP1qIxUKJf2mbaJXr15aOVtHjBjJloMX\nuHwjYyWAJqzfc4bnweFao5S6deuSv0BBps3foLdd/+dBzFy8mQkTNcvrmJmZMWHyFAYsO0F8kn5O\nSJ6qoO+SYwwa5I69veZ85Y8TJrLhrzfcCdNK5fp/WPMwiBhTS9q21VzG16RJE6zy5GPZU/0fIH/F\nJLDqRRjjp03XOMfS0pIfx41jYmAUyXoKqaaoVIwPjGLU6DFac6oTpk/nkJk5LwwQD90vM8E8bz6N\noqGQpiyjzJGTvQYocf+VquCAQjBWi+zSx5AkiRUrV9Kv3wC+7zeEitUbULthc8pWqYv7yPGMGzee\nX+bO1fv7/zFkQgPvn8Y33YEnhGDMmFFcvXAar9WjyZtLc4SQlCyn92RPEoUtx06c0kmd1wxOAAAg\nAElEQVS+c/bsWXp+350dS36kaT3tahjrdp1mzqr9/HHxMqVKldI6NyIignrf1aW9a23mThqotbri\n/qNnuP0wlRmzftapFiGEwH1gfx7dvMzhnzqTM5vmZor4pBR6zD+CmWNhDhw+qlPSyMvLi0F9erGl\nYTnq5tPc1SaEYM2jYDz8w7jse4MiRYpotfv27Vu+q1Gdzjks+LFkXq3VFXej4+l/O4CFqzzo2auX\nVrsqlYre3bsRcPl31hTPhb2WZopYhZKRzyPIX7seuw4e1FntsnfvXkb078d0uZwyWpoNhBDsk5lw\nLpsD1+/c0UnKFBQURN1q1eigSKa7hYnWvXggVzA9RcW6bdvp2LGjVrvqoFKpuHnzJrVr1+bWrVs4\nOztnORNcVnXgKU/M0GuuSeufv3gH3jftjCHtpv/559msWL6MHq1qM6RbI8qV+PvGfx0axYYDF9lw\n8BJNmjRj/cbNep8YX7p0ia5dOlOtYnGG9HClRUPnD5wIySly9p24wtqdZ4hJTOX4iVNqpdjVISIi\ngo4d2hERFsKQ3m706dyCbPa2H/48V28+YM22o5z+wxcPj7V0795dL7sqlYppUyezft06ejepxOBW\nzpRy+vsNICjsHetP+bHxjB9u7drjsdZTb0a4s2fP8n2XztTKm50BxXPh4vQ3WVBiqpKDz9+w6XkE\nClt7jp321psdLyQkhLbNXUkIfUuf/PZ0KpgL2/fdfkIIrkfEsvV1DFfC3rFp+w7atWunl12lUsm4\n0aPZuXULHXPZ8X0uOwp9VMERnCRnd3gcB8Ji6d6rF8tWe+jNCHf8+HH6dO9OJZlEy6REKkvSB2eW\nJAS/CzhlZY2tkxPHz53Tm2g+ODgYt2bNkIeG0FYlp5mlOVbvHb5KCG7LFRyRzLgnV7Br/35atNAp\nWKwV/waiIOXJmXrNNWk1+4s7Y52nqFl1pX3V14ugoCAxffo0kTdPLlHYKa+oVLaYyO1oL7I72Iuh\nQwaLBw8eGGU3MTFRbN68WdSoVlU45nAQ5UsXE/nz5BDZ7O1EC9emwsvLSygUCoPtqlQqcfHiRdGt\na2dhb28nShUvIgoVyCPy5c0tSpUsLpYuXSqio6ONWnNAQICYPGmiyJ0zhyiSP7eoVLKQyO1gI3I4\n2IuRI4aJx48fG2U3Pj5eeHp6isplSotc9raiXP7cIp+1hchmbSXaNGsqTp06JZRKpcF2VSqVOHfu\nnOjYupWwt7IUJXM5igJW5iK3va0oV7yoWLVqlYiJiTFqzc+ePRPjRo8WjvZ2olD2bKJcbkfhaGYi\nHO3txNhRI4W/v79RdmNjY8Xq1atF2SJFhKOVlShhby9ygLC1sBBuzZqJs2fPGrUXSqVSnD59WrRu\n0kTYWViIotmyiTwySThaW4nyxYqJtWvXiri4OKPW/Ck+52/6ve1M+xzl6dl6XVnxfZm9vvnI+FOk\npqYSHBxMbGwsVatWJSEhIdMyQ+l4/fo1kZGRVK5cmfDw8CzRCoM0PuaQkBDKli3Ly5cvKVSoUJa8\nNsrl8g974ezsTGJiIlZWxpVnfQwhxIe9qFKlChEREVrVjw1BdHQ0ISEhlCtXjsDAQAoWLJgle5GS\nkkJwcDBxcXE4OzuTlJSUKWL5dIj35WlRUVFUqVKFyMhIcuTQn6BIG6KioggJCaF8+fIEBQXh5OT0\nOZjVssyeGtuZj4zP/KzXXJPmM754ZPyfM9YAIQQymSzLb7bPZRf+fVJD/8a9+O+++Px2P7KdeWfs\nPUevuSauP2X4PkmSSgF7SdOxk4BiwE9CiBUfzWkIHAVevP9fh4QQvxiz3q/jGPErgEKh4MiRI7g2\na4qDQzZMTEwwMzOlZo3qbNy4kcTERKNt37t3j8GD3cmbJzfm5uaYmMgoUrgQU6dOISgoyGi7ERER\nzJ8/n7JlSmNtbY0kSeTM6UiP7t24fNm4GurPjdu3bzOgT29yOThgZmqKiQTFCuRnxvRpvH792mi7\noaGhzJ0zh5IFnbAyN0cCcmazp1fXLvj4+Bi9F3K5nN27d1O/ejXsrKzS7gvgO+eq7Nq1i5QU42qS\nhRD4+vrSp1s3HO3t0vYCKJY/P3NmzSIkxPA69XS8efOGWTNnUKKQE5bv9yJ3Dgf69e7JrVu3jLb7\nr0Qm6oyFEP5CiKpCCGegGpAAHFYz9ZIQwvn9ZZQjhv+cMZB22l+0SGEW/PYzvTs25tnNY6S8vUXY\nkz+YOa4fR/bvoGBBJ5YtW2bQjzooKIiGDerRumVz8turuHZ0FfFPvYl94o3XpjnEhz2napXK9Oje\nlbg4/TvrUlNTGT16FCVKFOfx3etsXDadsL8ukPzmJvf+2EutSkUY2L8vlStV5O7du8ZsSZbj+fPn\n1HGuSvtmTchz9won6xTnRbua+LethWfZ3AQf2EGF0qXo072bQQ++lJQUhvTvR6miRXi40YOFuSy4\nVasEj+uXxatcPgr7XaVHqxZUK19OLemQNmzfto2CeXKzcvQIOoa+5HTh7Nwpk48/Sueja3gQq38c\nScE8udm8aZNBdh8/fky18uXp3Kwp9meOs8FSxvlcdpzJnY3p8jjurFhK6aJFGdC79//aO8+oqo6u\nAT9DExDFTqyoWGKNsZckYu+9RKNir9ijJmp8ozExxig27B17wZLYsGLsvStiVwRBeq93vh+In8Lt\noEBynrXO4nJn7j57HQ777jOzi9oKfpqIiYlhQJ/eVCpXFp+D29nQqAwvRzYn6Pt2/PNtbUr73KJL\nq+bU+bIa3t7eBumcbcm4DLymwGMp5Ut1Z8kQVf/ryxRr167lp6mT2b76D76qW13jvEdPXtC53/c0\nb9mGuXPn6Vx7e/DgAU0aN2L0gE6MG9xNY/hXVHQMY3925dq9Z5w4eUpnAH58fDydOnZAJkSxYclM\nCuTPq3aelJItuw4ybto8du/eY3BfuRTi4uIICwvDzs7OaO/y9u3bNHNsyCj7vPQrZaexcllkQhKT\n7/rwMncBjp46rbNVVExMDK2bNiHH80fMts+PrYaaySop2RUQxlzfcA4cPUbt2rV16vzn7N9ZNGsW\nLgWsqKClDsaD2ATGB0YzbMIkJk+bplPu1atXadm4EQNFEm0tzdO0XEohQqVibjzElCnHkVOndO5b\nREVF0bxRQ4rGBrOwcQVy51Af5ZKkkqy59Zw/rrzg8PETVKumPexSE7GxsVhZWREbG/tR6lFk2DLF\nid/1mmvaeLLW8wkh1gBXpZRLU73fENgF+AC+wEQp5T2j9P0vG+Njx47Rp3cvPPetolyZkjrnh4aF\n8027gQwaMoIxYzT3aAsJCaFmjepMdv6WgT3b6pQrpWTsz4u589ifY8dPao1VHThwAIF+z9i1bq7O\n+F6Ao57n6TNiGmfPntPaHik1169fx8VlHnv27MXSMgdhYeG0a9uWcePHG2TYAwICqFG1ClPs89Cp\nuO4NSyklE+/4EFy8LAeOHtP6pdejUyeir5zDpXQBvUpTHg+M4CffcC7duKm1E/X2bduYMGQw6z+z\nwU6PovhvEpLo6x/JLNel9O7TR+M8Hx8falWtymgZR0NL3YWOVFLyW5zEsl4D9mgo+QnJ16xD61bY\n+j9iadNKGg38++x+4Mvk88+4fOOmQd05zp49y3wXFw4eOoSJiUBK6NixA+PGjadmTeNK0qojw4zx\nSfXF8z1vPObUjSfvfv9lw3GN5xNCmJNsaCtKKd+kGrMBVFLKaCFEK2ChlFJ7soAG0r1MIYQoJoQ4\nIYS4J4S4LYTQr89JFmDqlMks/XOyXoYYII9tbtzXzeXXX2cSExOjcd6KFStoUKOCXoYYkm88l5+d\nCQ0K0Fpb9sGDBxzYv5/Ny2fpZYgBmjnWY0T/bvyhZ0cHgF07d9KiRXOqVSrHC+8bBPp4E+z7kOaN\nG9Cjx7csX6a5G0lqXBcupLFtDr0MMSRfi9mVivHo1g1Onz6tcd7169c5feI4c0rl17tGcJMCueho\na8m82ZqvhUql4sfx4/gtv6VehhigoLkps/NZMWXC9yQlac6Cc5kzB0dVgl6GGJJrL/+YQ3D5zGmt\nnWXOnz/PvetXWNy4ol6GGKBz+SK0tc/DogXz9ZoPsHrVKrp160rjb+ry+tk9IgNf4PPoFjW/qEib\nNq3Zvm2b3rI+HULt4VitDD/3a/7u0EErkr3iNBXApJSRUsrot68PAeZCCKPCYdLtGQshPgM+k1Le\nePstcRXoIKX0SjUvS3nGV65coVvXTjy89JdB7dsB2n03mq49+9GvX780Y0lJSZRxKM3O5f+jRtXy\nBsldu+0A+07c4q/9B9SOjxkzhpym0fw2Vb8miyn4BwRRoX4nnj59pnMZ5PHjx9StW5ej+3dS7Ysq\nacafPnvOV03asnv3Hp1dsuPj47EvUpht1UtQXkfD0NSsefyaW/aV2b53n9rxwf36YnvmCM7FDAuJ\nexUbT/u7vrx87U/OnGmzDA8fPszE3t+xpZDh9ZJ7+0fz63o32rZN+yUcExNDsUKFWG4pKGpAeySA\nTTEJRLdqx/qtW9WO9+relapB3oysaVi5Uu+gCFruvsYLv9c6lxquXr1K27ZtOH3sb8o4pD3P7Tv3\naNK6M6dPn6F8ecPue3VkmGfs+Ydec00df9DmGW8FDksp09SfFULYSSn9376uDeyQUpY0Rt90e8ZS\nytdSyhtvX0cC9wHtuZtZgBUrljPEqYvBhhhgWL+urFi+VO3Y0aNHKZgvt8GGGKBHhyacO38eH5+0\npT3j4uLYtGkjQ5y6GCzXrlB+WjZuwKZNm3TOXbpkCQP7fqfWEAOUKmnPpHEjWbRooU5Z+/fvp7S1\nhcGGGKBbiQIcOXqUN2oKn0dGRrJz5y66FzK8hm9RSwtq5rFhx44daseXL5hPVwvNtaq10TWHZJmL\n+rrA7u7uVLAwM9gQA7TJYcqevXsJD09buD44OJgDBw/xXaViBsstlz8XlQrmYu/evTrnLlq0kO9H\nD1driAGqVK7IkAFOLF2yxGA9Pirp3MATQliRvHm3+733hgohUtqedBVC3BFCXAcWAOorRelBhkZT\nCCFKAtWAixkp92Pgdf8e9Wt9YdRn69euxn2vB2rH7t+/T70aFY2Sa21lSZUKDmp3ul+/fk1Oayvs\ni+uXGpuaujUr4+V1X+e8Xe676Ntbe/p03949cHffrbN63P3796lhY1wXitzmZpTKm5snT56kGXv5\n8iUFrXJQ0MI42dXMJfc1RFbcv3ePatbGFa6vZmWBl4b74t7du1SMN7yUKEBeExMKWVry/PnzNGNP\nnz6lVH5b8hlZbL92AWvu39d+X0gp2bXLnX59emqd169PD3a57zJKj49GOktoSiljpJQFpZQR7723\nQkq58u3rJVLKym9D4OpLKY22fRnW6ePtEsUuYMxbDzkN06dPf/fa0dERR0fHjDq9wURFRWFtZVwG\nVU5rKyIiIjR6TxNHqK+Vqw8J8XE0adJE7VjhzwoaLTentRXr1y1myRL1Hn0KFhYW2BXSfp48eWxJ\nSkrU+VQhgIkV1Zft1IeEyAiN3S2KWxrfasja1ITfXVzUdrewEgJLB/XlNXVhZWKC7+vXau8LM2B4\nLuMz9uIiwqlatarasQoFDO8FmEJOc1Nm/DKDGTNmaJ0nhCB/fu1LoXaFChKmoe2ULjw9PfH09DTq\ns9rJGoXj9SFDPGMhhBnJhnijlFL9Ih/JxjjlyExDDGBra0tYhNrvDJ2EhoWTJ08etfnlixYtIjLa\nOA8IAGGGp6dnGrk+Pj6kZwktLDySgYMG68yPd3AojfdD7TWeX7zwwcYmFyqVSqus32fPJsK4tm8A\nmOS04fLly2nkent7Y2KluaqcLiISVYyfNFGtzqXsSxBpQInLD+QmqShetIhauZOnTSMyHVsmpjY2\n3L17N43c69evY2JpfLp+WHwS03+ZqfO+yJ8/P0+ePtMqy/vhYwoX/swoPRwdHT+wDxmG+v27tEcW\nIKOWKdYC96SUuhcSswi169Tj4LEzRn324LEz1NWweVWnTh0On7xkcAF4gIDAEO56P1brAX322WdY\nWFhw7abupQZ1HDp+TqOX+T5OfZxYsSZtn7z3WbFmA32dnHSuq9apU4cTwdFGxSf7RsfxMjxKbdcT\ne3t7YqTkkZHtkU7Fqqhbr77asToNvuJ0dLxRck9Hx1Onvnq5devW5ZKZcfG4LxOTiFBJtVX9ypUr\nh194NM/DDM8QlVJy5GWozo1YgL5OTixftV7rnBVr3Ojr1NdgPT4q/6W2S0KIBkAvoLEQ4roQ4poQ\nIn21+T4Bw4YNw23b30RHaw5R08SydbsYoaG9Tq1atciXvwAenpcMlrt22wE6depM3rxpEzlMTU0Z\nMnQYy9bvNFiu18On3PF6rFft2kGDB+NxzJOd7uofcE54/sOqdZsYNVp3BGPDhg0ROXNxPtDwR9eN\nLwL5rlcvbGzSdq22sLBg4JAhbAk0/MnmbkQM/kmCNm3UN411HjeOXVGJJBn4BaKSkp3RSYzU0J2k\nRYsWRFjk4L4RTXD3JUr6DxqoNuLB2tqaPk5OrL1teD/HC6+CiTW10Lgs9j4jR41iw+YdHD3uqXZ8\n9979/H3wCIOHDFE7nmn8l4yxlPKslNJUSlktJY9bSnk4I5T7mJQqVYq6deuyetNu3ZPf4+TpS7wJ\nCtVYC1YIwQjnUcxdvl1rzGlqwiOiWL7xb5w1GHmAQYMG4f73MV6+MqxuwZ+uGxgwQP0/c2oKFCjA\nwYMHGTtpGj2dhnDy1GmeP3/JufOXGDR8LD2chrJjxw69ai8LIRgxbjzLngcZ5B2Hxiey5WUQI0Zr\nTqwZOsKZvQFhvInXv/eblJJV/hEMHTlSY5x2jRo1KFzCnsPhhn1JH42IJX/hwhqz+0xNTRk+Zgxb\nkgwrrhOSpOJwXBLDR2oOZxw+chRud14SFKO/Ry+lZMH1FwwbPUZnQXyAkiVL4u7uTu8BIxgwdDRn\nz13k+fOXeP5zht79hzNy/GT279+PnZ2d3jp8GrLPOsV/ujbFn3PnMWvBOjxOnNVr/n3vJ/QaNpWl\ny5Zr3bzq3bs3KlNLxk131esfLyYmju7DptOqTVtq1KihcZ6dnR1Tpkylbc/RBAaF6KXzguWbOHf5\nDhMmTNBrPkC1atW4e/cudet/zcSpM/mqaVtaduhO6TKfc/v2bYPW+wcNGkR4XjtmefnqdS2iE5MY\ncO0Zvfr2o2JFzVEp9vb2jBk/niGPAgnXs0v0slchPM2Zh9Fjx2qd57pmDXNC47mp53LFnZh4fg+N\nY8nadVqXbkaOHk1gocKsj0vS61pEqSST4yXDR4+idGnNMcTly5en3+Ah9Nh/g4g43V9OUkp+OfcQ\nf3Nbhg4dpnN+Cl9//TW3b9+m3OeVcR4/mYrVG/D95BnUqFWXO3fuaL13M41s5Bn/p9OhITm9s3Pn\nTvw4uh+D+3QmZ860myGJiYnsOXCCUT/+wdx5LvTRkvKaQkhICK1btaBk4Tz8PnkI9sXUb2zcvPsI\n56kLsHcoj9vGzToz66SUTJkyGfed21k1/398XU99y5s3gcH8vnAt+4+c5eix4zpbGOkiPeUS37x5\nQ3PHhpSPD2dS2c8oYq3eQ78ZEsmP93yp1qQ5a9w26vTYpJSMdXbGY8dWfi2el+q51SdqvIlPxNU3\nlMsmVhz95zTFiumOyT106BB9undjdC5z2tpaY6GmRVK8SnIwPJr5YQms37qVdu3a6ZTr5+dHkwYN\n+DzkDf3MTSigphmmlJJ7CUnMTRQ06f4tS1au1Lk+r1KpGDZoIJeOHmRRo8+pUVh9zRK/iBhmXnjM\nzVhzPE56psuTzRYlNC8u0j0RMK0zWqlnnBW4f/8+P0yayNlzZ+ndrS2tmtQnj20uIqNiOHfpJqs2\n7qZECXt+/W0WjRo10ltudHQ0P/00lQ3r19OgVhV6d25G0c8KkKRS8eiZD2u2HuL5qwDGjhvP+PHj\n9XpcTMHNzY1Zv/2KmalgaN/OVKlQlhw5LHgTGMzOv47xt8cpOnXsyJw/51KwoPEhcSmk9x8vMjKS\nH7//ni1bNlO/UB46FbKhkKU5iSrJk8hYtryO4E0STPhxMqPGjNE76UJKyZrVq5k9YzqW8bH0yGNJ\nGWsLzIUgKCGRA+HxnAqKoFvXLsxZsFDterwmrl69yg+jR3Pz5g065spBbQsTcpqYEKVScTlexd6I\nOCpXqcLshYv0Kj6UQlhYGJPGjmX79u3UtrSgsUwgn4kJicCzxCQOmFgQbWnFj//7H0OGDTPoWixd\nsoS5v88inzn0/7wQZfPZYG5ign9ULDseBeL5LICePXsye+48rY1k9SFbGONLi/Waa1p7lGKMsxLP\nnj1jxYoVXL50gfDwcLy9ven5XS+GDh1mdHUrSI5p3rJlC/v2uhMUGMR9r/s0bdqMXr2daNeund51\nJlIjpeTkyZOsW7uG58+fcevWLerWrUvTZi3o379/hnXPgIz7x4uIiGDTpk38vWM7QUGBeN33onmr\nVjgNHkLr1q2NyoiEZM/w2LFjrF22lJfPnnH79m3q1a9H6y7d6Nu3r17t6DXh7e3NssWLuX7xAuER\nETzy9qb/sGEMHzVabbSHvoSFheHm5sbf27YRFBzMgwdetG7bjgEjRtC8eXODvpzfJykpCQ8PDzas\nWsErn5fcvnWb+g0a0LZLN/r06ZNuI5xC9jDGrnrNNa09UjHGWZns1iEhu3X6yI5yP6ZsRW4a2ek3\nxpf1S882reWc6cY4wzLwFBQym6SkJE6fPv2ue8qePXto2LBhhvSU8/Ly4saNG++aANy/f58KFSqk\nW25ISAienp4EBwcDcPLkSb755hujnxBSSEpK4tSpU+/qnOzduxdHR8d0PSFkS7LI5pw+KJ6xFrKb\nR/Ff9QYDAwNZs3o1yxcvJLdMomwuK+KjooixsOJKQDAdO3TEedw4g+vtprTicp37J/fv3aN2vlzk\nFPA6LAKvRCj3+eeMnDCRTp06YW5uWHr2tWvXWLJgPrt376ZGobwUtDAlNDwSP2FGiEowbNRoBg0e\nbPB6f0BAAKtWrmDlElcKmJtQ1jYnsdFRRJlacMUviK5duuI8dmy6lt1SyBae8VXt6f8pmNYYkeme\nsWKMtZCVDdCnlPsxZadX7qVLl+jQuhWO+XLSt1gequX7sE5DYFwC218EsvZFEIOcRzF95ky9NsRC\nQ0Pp1LoVEU8e0i+fFc0K5MLivTXcBJXkaFA4G4JisbQvxb7DHnptDkop+XXGDJYtmE9/+3z0LF6A\nAqnqG98MiWSDTwgnAiPZs/8A9erV0+tanDlzhq4d2tOiaD4GlrPjiwK2H4wHRMex0duXlV6+jJ30\nAz9MnmxUdboUsoUxvqZf7W3T6sMz3RgryxSpCA0Nxd3dHV9fXyC5O4M+oVC6SEpK4vDhw9y8eROA\nc+fOUa9evXT9M6Rw584dDh9OzrPZvn07HTp0yJA28kFBQbi7u+Pv7w8kh2UZ0hVCE4mJiRw4cIA7\nd+4AyQbVkGiEFK5evUrrZk1xqVSU5kXUL0UUyGGOc9nCdC+enz5rVhAfF8fvf/6pVW5UVBTNvvma\nSuGB/FRefZsocxNB64K2tCyQm1kvXtLkqwb8c/GS2ozB9/nf1CnsW7uKQw3KYaeh0toXeW1wyWvD\ncb8Q2rVswaFjx6lVq5ZWuefPn6dTm9as/Ko8jYur96YLWefg+2ql6Fm2MN1c55OQEM+0n6drlauO\n169fs3t3crLUihUr6NKlCwUK6Nc84JMjsk8qRfbR9COTmJjIpIkTKVWqFB6H/iY2MoSe3TtTtWpV\nevbsobaWrL7s3buX0qVLMfOX6YQH+zNq+CD69XWievUvuXr1qtFynz59SqNGjrRo0ZwXT7yZNH4U\na1atoESJEumqKxsfH8+Y0aMpU6YMJ495EBsZQvcuHalYsSJ9+zoRFRVltOxtW7dSsmRJ/vxjNpGh\nb3AeNpAe335LrVo1uX37tt5y4uLi6NC6FX9ULKLREL9PQUsLNtUsxfa1qzlwQH3x/hTGOY+gRGgA\n/7PPp7OLiIkQTC2RjzKRwYwaNlTrXA8PD9yWL2NLzVIaDfH7NCmcl3mVitKhdSutnWWioqLo3K4t\nSxuU02iI36dITkvcm1Rh1aKFnDhxQuf8FKKjoxnQvz8VKlTg4rl/+GHCaP45eYyyZcsy0tmZ+Hjj\nanp8VLJR0ofiGZP86Ni/fz/8/V7hdeMcdnb/X0IxMjKSiVOm07x5M06cOKmzMWRqdu3cyegxo9m5\naQ0N6v9/QRaVSsWOXXtp1aoVHh4efPnllwbJ9fHxoWHDbxjrPIRRf237YM3yvpc3nXv0IyoqiomT\nJhkkV6VS0bNnD+Jjo3l4+yIFCvx/eFx4eARjJkyhdetWHDly1OBGlG4bNvDTtJ/Ys20dtWr+f/NX\n1TwVG7fsoGnTppw4cYJKlSrplLVr1y7KWJrRpqj+4XsFcpjzvUNB5v8+S2NtiqCgILbv3MnxL0ro\n/dQihODHYnlotHsPc9680bjOO//3WXzvUDDNsoQ2WhTJh5tfODt27KBvX/VFeDZv3kz1/DY0L6F/\n6c/PcloyoXIxFsyZTePGjXXOj4+Pp127thSxK8Azr2vY2v5/eFxQUDCDho+jW9eu7N6zJ92bjxlK\nFjG0+qB4xiTvYF+5fJl9Ozd+YIgBbGxsWLrwTwoVyMfKFSsMkhsXF4fzyJH8tXPjB4YYwMTEhB7d\nO/PHr9MYM8bwtoHTpv2E03fdGT9mRJrNowqfl+Po/l38Pns2fn5+Bsk9cOAAjx89ZNeWtR8YYoDc\nuXOxZvlCLMxM2bB+vUFyo6KiGDd+PIf2bvvAEEPytejbuwfTfhzP9+PH6yVvyby5OBUxPF62bdH8\n3Lp1S2Or+nVr19KkQG7yWxjmp+Q1N6NFQVvWrF6ldvzx48dcuXqVdga2iQLoWzg3S+bNVTsmpWTp\nfBcGljG8BnPXMoU5e/bcu+gTbWzevJmkhHjWr3L9wBAD5M+fj+2bVvHq1Uv27fNWxFMAACAASURB\nVNNYQVdBB4oxJrnV0BjnwVhZWakdF0Lww/hRLFu+zKANC3d3d76oUomaNTR7vb17duPx48fv1k/1\nITg4mL179zF2pObH4mLFivBt1w6sXqXeOGhi6ZIljB89XKPXa2JiwqTxI1m2XP+mpABbNm/m6wZ1\nqVRRc5LEoP69uXb9Go8fa6+n7Ovri5e3N80LGx6ylsPUhE5F82psu7TDbQOd8xi33t4lryU73NzU\njrm7u9OuSF6sjPAamxTOy9OnT9UazcePHxPw2o9GxQxfs81pbka7Unbv1n+1sXTpEiaNH6nR67Ww\nsGD8qGEsW6pf9MInI/1tl54JIW6+rUipthSjEGKREOKhEOKGEMLoMBXFGAOXLl+iRTPtj2r169XG\n19fPoLXjSxcv0qKpo9Y55ubmNG3UUGv339TcuXOHihXKp/FcU9OyWWMuXTKslOflK1do0VR7yneT\nRt9w585dg9YIL1++rPNaWFpa0vDrBly5ckXrvICAAArnyomZmnoR+lDEwpQ3Gp4Y3gQGUtTILiJF\nLc15ExSkdszf15diFsY9vpsKQdHcNgQEBKQZCwgIoJitjd5doVNT3NIU/9faqwCqVCquXbtOcx33\nRcvmTbis42/3yUn/mrEKcHxbkTLNLrMQohXgIKUsCwwFlhurqmKMDSAhPp48efIghNDrWKpnS/vI\nyEgGDBigt9yGDRsSGhqql+x/Tv+jt1whBFGR+tUIVqlU5MiRQ2+5mzbrboYKEB4eTo8ePbTK+vLL\nLwkPD9NLniaWLVumVrY6g2cIQUFBauXq08BVG6FhodSqVSuN3AYNGhAalr5r8ccff2i93qamph8t\nhO2jI0z0O7RIQLud7AC4Acjk/ne2Qgijqi8pxhioXas2Hke17yqfO38J+5L2OlsNvX+sX78ej2Oe\nWuUmJCRw4fJVbt++rbfcoKAgfP38CQxU74WlcPjoCSZ8P0FvuVJKHB0d8Th2Uqvc4yf/oWrVKgbJ\nXbhgoc5rERsby/Wbt3n06JFWWa9evSJKmJFoZHsk3/gk/vfLL2plV6pYkVex+tdIfp9XsQl8Xq6c\nWrm/zf4Dn3j961u/T5KURAkznj9/nkbuw4cPicQUlZHG8mVsEi7z5+v8+9WoUZ0jOu6Lw0eOU8vA\nxJqPTvo9Ywl4CCEuCyEGqxkvCrx87/dXb98zGMUYAyOcnVm4ZBXR0epb10gp+cNlMcOHDTcoLrhL\nly7cvH2XK1eva5yzaetOHBwcqFy5st5y8+XLR8eOHVjgqnlD0cfHl+279jFosLr7RzMjnJ1xWbSM\nuDj1LY1UKhVzXFwZPmy4QXK/69WL02cvcOeu5rZRq9dtovqX1XUWri9SpAiflyvHEb9gg3QAiEtS\nsedVCN27d1c73t2pL7tDjeth6B4SS3cnJ7VjXbp04W/fEGIMaDiQwnG/EEqVKkWJEiXSjDk4OFDo\ns8Kc9Ak0WG5UQiJ/P/XXqwPMiBHOzHFx1dgwIT4+HpfFyxk+YoTBenxU0u8Z15dS1gRaA85CiK9S\nn0HNZ4z6ZlSMMdCoUSNq1qpFx+5O+Pt/+JgaERHB8NETCQgMZshQ7XGkqcmRIwdLXF1p360PZ85e\n+GBMpVKxdbs7P/w0k4UL9au5+j4zZ/6K25YduCxcSkLCh57cfS9vmrXtyuQffzQ4SaNNmzaUKVuO\nrt8NSON5h4WFM3DYGBKSVPTt188guTlz5mS+iwutO/Xk0uVrH4wlJSWxfuNWZs52wWX+fL3kOX8/\nATdfw2O/978KomrVqpQrV07teP8BAzgeGE5QvGHtkUISEvF4E8bAQeq//BwcHKhZswZ/+2h/mlHH\nBr9wnL9X3xxAiORuKmsevTFY7q5HfjRoUF+tkU9Nr169MDW3oN/gkYSGfrgsEhQUzLe9B1O0aHE6\ndOhgsB4fF/WdPTwveTHDdc+7QxNSytdvf74B9gCp1419gPdboBcDfI3SVEmHTiYxMZEpkyezavVq\nmjVpSFmH0jx5+gyPY560bNmC5ctXGF16cO/evYwdOwa7QgVp9E0DoqOjOXTkBDa5crF69RqjOyQ8\ne/aM/v374e3tTef2bbC2tuLajdvcunOP6T//bLSXEh8fz8QJE3DbuJGWzRpTupQ9Dx8/4ejxU3To\n0J4lS5aSM6dx3Zm3bd3KhIkTKVmiOF83qEN4RASHPE6Qv0B+1q5dR5UqVfSSExcXh0OJ4sxwyE8b\nPZI+IDk1uv25hyx226QxzhhgSL++BB4/zNzSBfR6EpJS8uOzIHJ+3YR1mzZrnOfh4cGQnt+yv35Z\nCuoZa3zYN5gpD9/w+MVLjdE+UVFRlLEvwaLapfSONX4dFUuzg9dZv9NdrzhjSE76GDVyJLv37KF9\nm5YULlyIFy9ecejIcXp99x0u8+djYaF/DLU2Miwd+r5+exWmFXqnOZ8QwhowkVJGCiFyAkeAGVLK\nI+/NaQ04SynbCCHqAguklLo7/6rTVzHGH5KSDu3n58e0adN4+fJlhqZD37p1iylTpnD27NkMS4e+\ne/cuhw8fZsKECWzbto2OHTsanJChjuDgYNzd3Xn9+jX/+9//8PX1zdB06Lt37zJ16lQuXrxoVDr0\ntWvXaNG4EXMrFqGFDoP8Jjae3lee0rpPf73SoRvVr0fFsDf8pCMLTyUls14Ec906r17p0NOmTGbf\n2lVsrKE7C++4Xwhj7vjonQ7dvmULrenQKfhGxdL12G16DHNOVzq0s7Pzu3TojKydDRlojL226DXX\n9PPv1BnjUiR7w5LkBLnNUsrZQoihgJRSrnw7zxVoCUQB/aWUHz766auvYow1k1WL43xquR9TdkYV\nCmqYz5p+xfKqLRS07UUg614EMXjkaH7+5Re9vgBDQ0Pp1KY14Y+96ZfPiuZaCgVZlSzN3kOH9S4U\n9Nsvv7B0vgv9SuTjuxLaCwXtPXCQunX1c7R0FQryj45jo/crVnn5Me6HH5j043+gUJDXVr3mmn7e\nM9MLBSnGWAtZ1QB9arkfU3ZGyA0KCmLN6tUsW7SQ3DKRMjZWxEdHEfu2hGanjp1wHjfO4OWgxMRE\n9u3bh+vcP7l39y618+XCWoBfWDjeiYLyFSowcsJEOnbsaHAJzevXr7NkwXzc3d2TS2iamxIcEYE/\n5oRiwrBRoxk4aJBRJTRXr1rJCldXCpgLytjmJC46imhTCy77BdGtWzecx4zliy++MEiuOrKFMX6g\nPrknNabluyvGOCuTlQ3Qp5T7MWVnpFyVSvWuuLyTkxO7d+/O0OLyN2/eJDw8nCFDhnDv3r0MLy4/\naNAgTpw4kaHF5V+9eoWTkxN79uyhUaNG2Nra6v6wnmQLY+y9U6+5puW6KcY4K5MdDNCnkPsxZWc3\nuR9TtiI3jez0G+OHu/Saa1q2a6YbYyW0TeFfgUqlwsPDg55dOtOgxpfktjCjRcOvWbBgvt7Zipp4\n8OAB40aNpFGd2tSsVIHc5qaMHemMl5dXuuSGhYWxaNEiWjp+Q+2qlcltYca3nTpw6NAhVCqV0XKT\nkpI4cOAA3Tu2p0HN5GvRqrEjrq6u6SoFmz1RH9qW9sh8FM9YC9nNo/gveoNSSlavWsUfM38hp0yg\nX+mClM9jg4WpICAmnj0vgzn24g1du3Rlzvz5em2ypXD16lV+GDOaW7du0bN4fhrky0kuc1MiE5M4\nFxzFlpdBVKpUmdkLFhoUDRIaGsoP349nx44dNC5ekM7F82FnbUGCSuIdEsn6p4GES1MmTv2JocOG\n6b3JJqVk6RJX5v7+O/nNod/nhSiXLxdmJgL/qFh2Pg7k5NMAevbsyey589LdJTpbeMaPNMcQv49p\nmU6Z7hkrxlgLWdUAfWq5H1N2euRKKRk7ciRH3XewsG4ZatnlUWu4AqLjmHPzOWejJEc8T+kVqnjw\n4EGcenzLlLJ2dClRkBymaR8i41Uq9rwIZKb3a9Zt3kK7du10yvX19aW54zfUtpT8WK0kdtZpQxCl\nlFwNCGPsxUd8064jrstX6DTIKpWKYYMGcunoQRY1+pwahdV/6fhFxDDzwmNuxJhxxPMUdnZGlVEA\nsokxfrxXr7mmDh0z3RgryxTvER4ejuvixdSuVYvixYuRK5cNM6ZPf9eCyVhUKhUHDx6kfbt2lCxp\nj03OnPTt68TFixfTrfP9+/cZPWoU5cqVJae1NY0bN2Lr1q3p7roQEhLCgvnzqVmzxrtr8duvv75r\nwWQsSUlJ/PXXX7Ru1Qp7+xLY2ORk4IABRnU8mfHz/zizbxeHWlWj9md5NRqsQtY5mFuvHF0LWNCq\naWOdj+oXL17Eqce3bKhRiu9K2ak1xAAWJiZ8W7IQm2qWYkDv7zh37pxWuZGRkbRu1oT2ec1xqVdO\nrSGGZENU0y4PB1t+waVDfzFtyhStcgF+nDiBO54eeHStpdEQAxTOZcWSppVomd+Ets2bGty1JSAg\ngN9nzeKLL6pibW1NjRrVcZk3711366xH9lmmyBBjLIRoKYTwEkJ4CyF+yAiZn5p79+5RuXIlTnke\nY9aMyZw5tp+j+9157fucqlWr4uHhYZTc2NhYOrRvz5TJP9CxbTNOHNrNlXPHqFy+NN26deX78eON\n9i5cFy+mYcNvyJ3Tgt1b13Hn2hlGDHJi+VJXGjSoT2Cg4fUKAG7cuEHlypW5dOEsf/42jTPH9uPx\n106ePfGmSpUqeHp6GiU3MjKSli1a8MuMn+nRtR2eHnu5fOYYZUoVpX37dkz76Se9r8Xz589Z5OLC\n1saVyZNDv7Cy8V+UpLyIZ9HCBVrnjRoyiF8+L0yN/Lm0zkuhWr5czKpQBOdBA7Tq77p4MfaJ0Uyq\nVlKvpYfcFuZsbVSJ5UsW8+TJE43zHjx4wPrVq9jWtho2ehTFF0IwrX5Z7BLCWbFC/4qPp0+fpkqV\nKjzyvsfS+bO5d+0sLrOnc+3KRSpVqpSuFmIfC2FiqteRFUj3MoUQwgTwBpqQnJN9GeghpfRKNS/L\nLlOEhYVRpUplfv15Mk69vk0zfvbcRTp+25dTp05RsWJFg2T369eX6IgwNq9fniYWNSQklGZtu9Kr\nVx/G6dnhIoV9+/YxZsxoPD32UtL+w9oCUkomTZnOhSs3+Oef0wYF9gcGBlK1alUW/vkr3bqkrTNw\nwvMfejgN5fz58zoL+qSmW9euWFuasWb5QszMPjQab94E0rhVZ0aOHMXQYcN0ypry4w8Ee+zmjzpl\nDdLhZmAYvc884onPqzQ6AFy5coVOzZtyoVFFnf3v3kclJfU977P94GHq1KmTZjwpKQmHEsVYV6cU\n1QvlMUjnqZceYd2oHXPmzVM7PnbUSMyuHWfGV+UNknveJ4gRZ5/j9eQZJiba/bJnz55Ru3ZtNq9b\nRrMmjmnG9+w7wIixP3Djxo10LX2kkFHLFKqn2vsdpmBSqs2/YpmiNvBQSvlcSpkAbCO5xme2wW3D\nBurVrqnWEAM0qF+HsSOHsEDPIjYpPH/+nL//3s+6lYvUJgXkzZsHt9VLmPPnnwYvK/z2268snjcr\njSGG5Bt5zqzpBAUGcvr0aYPkrl61ihZNHdUaYoDGjt8wuH9vFi8yrLiRl5cX/5z+h5VLXNQawYIF\nC7BuxUJ+n/27xspgKcTHx7N21SoGljM8NfuLArYUtjTV2JR02aKF9CmW1yBDDMmNSZ2K5WXJfBe1\n4x4eHuQ3wWBDDDCgfGHWr12jtpJedHQ0G93cGFCluJpPaqdu0XxYJsVz/PhxnXNdFy+mb6/uag0x\nQKcObWjXujmrVq40WI+PSvqrtn0yMkKL1PU8fTCynmdmscFtA0MHqS99mMKgfr3Ztn17mgpp2ti8\naRM9u3XSWlSnYoXylHUoxdGjR/WW6+XlhZ+fH61bNtM4RwjBkAF9DO5V57bRjaED1Te+TGHowL5s\ncHMzaHll08aN9O31rdaaGTVrfEmB/Pn4559/tMp6/vw51qaCcnm114HQRLNCNlw4f17t2MWzZ2lS\nyLjEiCZ2tlzUIPfChQs0tTNOXwfbnOTJYa62HZW3tzeFc1tjb2tYo1xIvkeaF8+j197FBjc3hg3u\np3XO0IFObHDbYLAeH5f/1ppxhtXzzCx8fX0pV6aM1jl2doUwMzMzKE7T19eXcmV1P8qXK+tg0Cah\nr68vDqVK6szSMlRusmw/nTqXKFGMmJgYYmP1r/ur77Uo61Bap85hYWHkNqDDcmpyW5gRHqJ+wyks\nMpLc5satIeY2NyVcQ6eUsOAgchvY5PR9bC0t1N574eHh5M5h/LWwtTAlLCRE65zExESCg4MpXaqk\n1nnlyjrg56e9hdMnJ/3F5T8Zxt8d/48P8P6zssZ6ntOnT3/32tHREUdHxww4ffrJkycP/gEBFCtW\nROOcqKgoIiIiKFBA/8aPJiYmTBjnrHPew0ePWbN+E0OGDNFbdtEihZFSal0P9g94w+kzhq0ZW1tb\n4x/whnz5NO/Ih4aGkZiYiLW1/t6Yubk5pUp8pnOe96PH9O7dm969e2ud95mGSAR9iE5IYsWqVSxd\nmbZZq42ZKdFJxj3YRSeqiAgPU3u9BTCllmHr2+/jFxhMvXr11I6VNMIrTiEyLpG5Li7MdVG/vJKC\nmZkZQUHBWvsu+ge8SdM5Wl88PT2N3hjWShZZgtCHjND0MlBGCGEvhLAAegB/qZs4ffr0d0dWMcQA\nXbt0ZcOm7VrnbNnuTts2bQxqNXT16lW27tijdQ3Uz+81t+96ERYWprfcpKQkclhacvmK5g4iAOs3\nbmPDBjeDdB42dCgbNm3TKnfDpm18+213g+SePHmSjVt2as0se/rsOc9f+BAVFaVVVkREBLEmZvhH\nq+9GoovLoXGsWrNWreymzZpyJSjCKLlXgiJo2NBRrVy3jRu5ZGQHkcCYOCIxUXuPBAUFEZIIwTHG\nhTJeDo5h27ZtOv9+333Xk/UbtVdAW79xG127dDVKD0dHxw/sQ8bxH1qmkFImASNJLrx8F9gmpdTc\nWycLMmToULbu3MPFS+pDc3x8fJk524UxY8caJLdatWo4ODjw62z1u+BJSUmMmTCVPr17G5QNZWJi\nwtgxYxk7aarGVlGbt+7k5Stf2rdvb5DOI5ydWbNhCzdu3lY7/vTZc+bMd2X06DEGya1fvz558ubF\nZaH6Vu4JCQmM/n4KgwcN0ulx29jY0LVrVzZ6vzJIB4AXETFc8AvS2HZp6OixuPmGGxVu6OYXzvCx\n49SOdenShesBYTwLV//30sYmbz86deyo9h7Jly8fbVq3YstdH4PlegdFcPdNBB07dtQ5d/ToMcxb\ntIxHj9WH2N2+c4+Va90Y4az7SfCTYmKq35EFyBAfXkp5WEpZXkpZVko5OyNkfkqKFSvG+vXradul\nF7P/XMibN8nxuVFRUaxa60b9xq0ZN3acUd78li1b2bpzL04DR7wzcFJKjp/8h+ZtuxEaHsUfc+YY\nLNd55Eg+/7wSDRq3Yc++AyQmJrcJevzkKeMmTmXSTzP566+/DS7t6ODgwPJly2jerhtz57sSHJy8\nnhgREcGylWv5qklbpv00TW34ljaEEOza5c6y1RsYOGzMu154KpUKj6MnaNKqM8LEnBm//KKXPOcx\nY1nn/ZrYRMN6yq3xekXvPk4aN1WbN29OhKkZlwz0jq8FRxCoErRq1UrtuJWVFU79+rHyvmFfIHFJ\nSax7+BpnDUYewHnseFbd9SUhybB6FituvWTA4MF6NSKoUaMGM3+ZyTfN2rNk+RrCw5OvT0hIKPMX\nLaNpmy4sXrSY8uUNC6/7+GQfz1hJh36P27dvM9/FhZ27dpEjhwVRUdE0b9aMsePG0ahRI6PlhoSE\nsMTVlRUrVxAZGUViYiIlS9ozYvgIBg4aZHSrGikl27ZtY8kSV65du46VlSVCmNC/Xz/GjB2brg4l\n169fZ/58F3bv3oOlZQ6io2No1bIlY8eN4+uvvzZabmBgIK6LF7Ny1UpiY+OIj4+nTBkHnEc403/A\nALVhb5ro0aUzCfeuserrCpia6P6HOvTMn3FXnnHx2nWtfd+2b9vGxBHD2FevDIWtdBsq/5h42p9/\nyKxFrvTSstb96tUralf7gj+ql6BdKd3r5yopGXbGi0SHyrj/9bfGeVJKOrRuha3/I5Y2rYSJHnsE\nux/4MuXCMy5dv2lQ95Zz584x38WFAwcPYm1tRUxMLJ06dWTcuPFGtw9TR4bFGb86pddck6INMz3O\nWDHGaoiPjycsLAwbGxuNfceMQaVSERoaiqmpKblz586QlkspREREEBcXR548eQwyaLqIi4tLjl7I\nnRtLS8sMk6tSqQgJCcHMzMzoaxEbG0vrZk3IGeCDa/1yGjPxVFKy2duXX248Z7/HEb2K+vw5ezau\nc2azuro9VfJoDkm7FxbFgKvPGDr+eyb/NE2n3KtXr9K6WVOmVClOn/JFNX6JhMYlMOa8NyF5C3Po\n+AmdSzdRUVE0b9SQorHBLGxcgdwarkWSSrLm1nP+uPKCw8dPUK1aNZ06qyM2Npbw8HBsbW0zpMVX\najLOGOsXZ29S9GvFGCsopIe4uDjGjBjOtm3baO9QmP5l7CiX1wYLExPexMSx+4k/ax++Jo9dYTZu\n30GlSpX0lr3RzY0JY8dQLpclfQvn5qtCebAxS67advZNGG6+4dwPj+aPeS70HzBAb7leXl707t6V\nID9f+pexo4vDZxSyykG8SsXD0CjWP3rNvsev6datG4uWLdf7SzAmJgbnIYPZs2cPnSsUZUDFIpTN\nZ4O5qQn+UbHs8PJjzZ1XfFa8BJu276RsWeOjOz42GWaMfc/oNdekyFcfnE8IUQxwAz4DkoBVUsoP\nMp2EEA2BfUDKQvpuKeWvRitsyI54eo7kUykofBxev34tf/t1pixrX1xa5bCQpiYmsoBtbtmnR3d5\n4cIFqVKpjJIbFxcnt27dKr+pXVPmtraSQgiZ29pKflWzutyyZYuMjY01Sq5KpZIXL16Ufb/rKQvm\nsZWmJibSKoeFLFO8mJz5ywzp5+dnlFwppXz16pWc/vP/pEPxotLSwkKamZrKgnltZf8+veTly5eN\nlvspeWsv0m1zVH5n9TpSn++tEa729rUN8AD4PNWchsBf6dUz5VA8Y4V/JVJqj8H+r8j92LI/Fhnm\nGfupz4hMjUnhelrPJ4TYCyyWUh5/772GwAQppe7aqfrokBFCFBSyGh/L+GQ3uR9bdpYnA2pTCCFK\nAtUAdXnjdYUQ14UQB4QQhlURS0XG7fQoKCgoZDU0fBF5nr2K5zndJT+FEDbALmCMlDJ1rvtVwF5K\nGS2EaAXsBcoZraqyTKGgoJDVyLBlCv/Les01sauV5nxCCDNgP3BISrlQj/M9BWpIKY2qtK8sUygo\nKPx7SV+hoLXAPU2GWAhh997r2iQ7t0a3PFGWKRQUFP7FGOdvCiEaAL2A20KI6yRXopwC2JMcebES\n6CqEGA4kADGA+oLo+p5TWaZQUFDIamTYMsUb7cW0UjAp+GWmJ30onrGCgsK/mOwTSaIYYwUFhX8v\n2aiesWKMFRQU/r1koxhrxRgrKCj8i1GMsYKCgkLmo3jGCgoKClkBZc1YQUFBIfNRPGMFBQWFrIDi\nGSsoKChkOtmpYp1ijBUUFP7FKMZYQUFBIfNRPGMFBQWFrIBijBUUFBQyHyUdWkFBQSEroHjGCgoK\nCpmP4hkrKCgoZAGUDTwFBQWFrED2Mcbp8uGFEHOEEPeFEDeEEO5CiNwZpZiCgoJCuklHDzwhREsh\nhJcQwlsI8YOacQshxDYhxEMhxHkhRIn0qJreBZUjQCUpZTXgITA5nfKyFJ6enpmtgkFkN31B0flT\nkN30zViEnkeqTwlhArgCLYBKQE8hxOeppg0EgqWUZYEFwJz0aJouYyylPCalVL399QJQLD3yshrZ\n7SbObvqCovOnILvpm6EIE/2OtNQGHkopn0spE4BtQIdUczoAG96+3gU0SY+qGbnVOAA4lIHyFBQU\nFNKJcZ4xUBR4+d7vPm/fUztHSpkEhAoh8hmrqc4NPCHEUcDu/bdIbls9VUr599s5U4EEKeUWYxVR\nUFBQyHCMD21TZ6FTt7dPPUeomaM/Usp0HUBf4CyQQ8c8qRzKoRzKoe+RAbbpmQHne53qs3WBw+/9\n/iPwQ6o5h4A6b1+bAgHp0TddoW1CiJbAJOAbKWWctrlSyuwTY6KgoJDtkVKWTMfHLwNlhBD2gB/Q\nA+iZas7fJDujF4FuwIl0nA/x1qob92EhHgIWQNDbty5IKUekRyEFBQWFrMBbZ3MhyXtra6SUs4UQ\nM4DLUsr9QogcwEbgS5JtYA8p5TOjz5ceY6ygoKCgkDF80sTt7JIkoivYO6shhCgmhDghhLgnhLgt\nhBid2TrpgxDCRAhxTQjxV2brog9CCFshxM639/BdIUSdzNZJF0KIcUKIO0KIW0KIzUIIi8zWKTVC\niDVCCH8hxK333ssrhDgihHgghPAQQthmpo6fgk9dRSPLJ4noGeyd1UgExkspKwL1AOdsoDPAGOBe\nZithAAuBg1LKCsAXwP1M1kcrQogiwCigupSyKsnRUz0yVyu1rCP5/+19fgSOSSnLk7wWm+VsRUbz\nSY1xNkkS0SfYO0shpXwtpbzx9nUkyUYidUxklkIIUQxoDazObF30QQiRC/haSrkOQEqZKKUMz2S1\n9MEUyCmEMAOsAd9M1icNUsozQEiqt99PqNgAdPykSmUCmVlfLqsmiegT7J1lEUKUBKqRvMOblZkP\nTCQ5rCg7UBoIFEKse7u0slIIYZXZSmlDSukLzANeAK+AUCnlsczVSm8KSSn9IdnZAApmsj4fnQw3\nxkKIo2/Xp1KO229/tntvTlZOEtEn2DtLIoSwITktc8xbDzlLIoRoA/i/9eY1pkBlMcyA6sASKWV1\nIJrkR+ksixAiD8kepj1QBLARQnyXuVopaCLDS2hKKZtpGxdC9CX58bRxRp87g/AB3q++VIws+GiX\nmrePobuAjVLKfZmtjw4aAO2FEK0BKyCXEMJNSumUyXppwwd4KaW88vb3yCZTRgAAAUlJREFUXUBW\n39xtCjyRUgYDCCF2A/WBrOgEpcZfCGEnpfQXQnwGBGS2Qh+bTx1NkZIk0l5Xkkgm8i7Y++3Ocw8g\nO+z2rwXuSSkXZrYiupBSTpFSlpBSlib5+p7I4oaYt4/ML4UQ5d6+1YSsv/n4AqgrhLAUQgiSdc6q\nm46pn5D+Avq9fd0XyOoORrr51MXlF5OcJHI0+d7IekkiUsokIcRIkiM/UoK9s+oNDIAQogHQC7gt\nhLhO8rLKFCnl4czV7F/HaGCzEMIceAL0z2R9tCKlvCSE2AVcBxLe/lyZuVqlRQixBXAE8gshXgA/\nA7OBnUKIASR/qXTLPA0/DUrSh4KCgkIWIPt061NQUFD4F6MYYwUFBYUsgGKMFRQUFLIAijFWUFBQ\nyAIoxlhBQUEhC6AYYwUFBYUsgGKMFRQUFLIAijFWUFBQyAL8H9wvZhHC0iCTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa4a5ea4828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "infp = 0.53\n",
    "g = prepare_standard_infection_graph(g, infp)\n",
    "\n",
    "infection_times, tree = simulate_IC(g)\n",
    "\n",
    "kwargs =get_display_kwargs(g, infection_times)\n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, **kwargs)\n",
    "add_colorbar(kwargs['node_color'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def make_input(g, infp, fraction, sampling_method='uniform'):\n",
    "    while True:\n",
    "        infection_times, tree = simulate_IC(g)\n",
    "\n",
    "        cascade_size = np.count_nonzero(np.invert(np.isinf(list(infection_times.values()))))\n",
    "\n",
    "        sample_size = int(cascade_size * fraction)\n",
    "        infected_nodes = [n for n in g.nodes_iter() if not np.isinf(infection_times[n])]\n",
    "        \n",
    "        if sample_size > 1 and len(infected_nodes) > sample_size:\n",
    "            break    \n",
    "    if sampling_method == 'uniform':\n",
    "        obs_nodes = set(random.sample(infected_nodes, sample_size))\n",
    "    elif sampling_method == 'late_nodes':\n",
    "        obs_nodes = set(sorted(infected_nodes, key=lambda n: -infection_times[n])[:sample_size])\n",
    "    else:\n",
    "        raise ValueError('unknown sampling methods')\n",
    "\n",
    "    source = min(infection_times, key=lambda n: infection_times[n])\n",
    "    if source in obs_nodes:\n",
    "        obs_nodes.remove(source)\n",
    "    orig_obs_nodes = copy(obs_nodes)    \n",
    "    obs_infection_times = {n: infection_times[n] for n in obs_nodes}\n",
    "    print(\"{}/{} samples\".format(len(obs_nodes), cascade_size))\n",
    "    return orig_obs_nodes, infection_times, tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'obs_nodes' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-6e5f307c1148>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_display_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minfection_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobs_nodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mnx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw_networkx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpos\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwith_labels\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'obs_nodes' is not defined"
     ]
    }
   ],
   "source": [
    "kwargs = get_display_kwargs(g, infection_times, obs_nodes)\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def simulate_on_node_as_source(g, s, K, order_constraint_n, obs_nodes, infection_times):\n",
    "    success_cnt = 0\n",
    "    if order_constraint_n:\n",
    "        ordered_pairs = random.sample(list(combinations(obs_nodes, 2)), order_constraint_n)\n",
    "        ground_truth = np.array([infection_times[u] > infection_times[v] for u, v in ordered_pairs])\n",
    "\n",
    "    for i in range(K):\n",
    "        times, _ = simulate_IC(g, s=s)\n",
    "        inf_nodes = filter(lambda n: not np.isinf(times[n]), times)\n",
    "        if obs_nodes.issubset(set(inf_nodes)):\n",
    "            if order_constraint_n:\n",
    "                reality = np.array([times[u] > times[v] for u, v in ordered_pairs])\n",
    "                if (ground_truth == reality).all():\n",
    "                    success_cnt += 1\n",
    "            else:\n",
    "                success_cnt += 1\n",
    "    return s, success_cnt\n",
    "\n",
    "def estimate_source_likelihood(g, K, M, obs_nodes, infection_times):\n",
    "    result = Parallel(n_jobs=-1)(delayed(simulate_on_node_as_source)(g, s, K, M, obs_nodes, infection_times)\n",
    "                                 for s in g.nodes())\n",
    "    df = pd.DataFrame(result, columns=['node', 'cnt'])\n",
    "    df['mu'] = df['cnt'] / K\n",
    "\n",
    "    # stored here\n",
    "    mu = {r['node']: r['mu'] for i, r in df.iterrows()}\n",
    "    return mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'obs_nodes' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-73f24fe0f3ea>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_display_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minfection_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobs_nodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m \u001b[0mnx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw_networkx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpos\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwith_labels\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m# add_colorbar(kwargs['node_color'])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'obs_nodes' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAJPCAYAAACgtar/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHZBJREFUeJzt3X+M7Xdd5/HXu+2yxB/oarOa9Acai0L6T4Fslz82OG7J\n3qLdvf+g22tYWdMgu5vqRkNS/Ku38S8wG4zBXTEhLIjmKpogCsSy4hAJFBpLhb25bS+u1l7akIAQ\nsiaSpn3vH3Mok9npnXPvvOfOnDOPR3KTc8587jmfT+f0w3M+8+W0ujsAAMy56rAnAACwbgQWAMAw\ngQUAMExgAQAME1gAAMMEFgDAsD0Dq6reVVVfqqrPXWTMr1fV+ap6uKpumZ0iAMBqWeYE691JTjzf\nF6vqtUl+qLtfkuRNSX5zaG4AACtpz8Dq7k8k+epFhpxM8t7F2E8n+a6q+r6Z6QEArJ6Ja7CuS/LE\ntvtfXDwGAHAsTQRW7fKY//4OAHBsXTPwHBeS3LDt/vVJntxtYFUJLziGunu3H8RWiv0LjqfL3b+W\nPcGq7H5SlSQfTPIzSVJVr0ryte7+0vM9UXevxZ9777330OdgLeu5jnVbyzo57H+W3l/Wsgp/1mUd\n3fvbv/Y8waqq302ykeR7q+rvktyb5AVbe03/Vnd/uKp+vKq+kOQfkvzsvmYEALDi9gys7v7pJcbc\nPTMdAIDV55PcL9PGxsZhT2HMuqxlXdaRrNdaOHrW6f1lLUfPuqxjv2q/v2O8pBer6iv5esDhq6r0\nmlzkbv+C42U/+5cTLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsA\nYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsA\nYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsA\nYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGLRVYVXV7VT1SVY9V1T27fP2GqvpYVT1UVQ9X\n1WvnpwoAsBqquy8+oOqqJI8luS3Jk0keTHJndz+ybcw7kzzU3e+sqpcl+XB3/+Auz9V7vR6wXqoq\n3V2HPY/9sn/B8bOf/WuZE6xbk5zv7se7++kkZ5Kc3DHm2SQvWtz+7iRfvJzJAACsg2uWGHNdkie2\n3b+Qreja7r4k91fVLyT5tiSvmZkeAMDqWeYEa7ejsZ3n5KeSvLu7b0jyE0net9+JAQCsqmVOsC4k\nuXHb/euzdS3WdnclOZEk3f1AVb2wqq7t7i/vfLLTp08/d3tjYyMbGxuXOGXgKNvc3Mzm5uZhT+NA\n2L9gvU3uX8tc5H51kkezdZH7U0k+k+RUd5/bNuZDSX6/u9+zuMj9o919/S7P5SJROGZc5A6sqgO9\nyL27n0lyd5L7k5xNcqa7z1XVfVV1x2LYm5O8saoeTvI7Sd5wOZMBAFgHe55gjb6YnwDh2HGCBayq\ng/6YBgAALoHAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGLZUYFXV7VX1SFU9VlX3PM+Yn6qqs1X1+ap6\n3+w0AQBWR3X3xQdUXZXksSS3JXkyyYNJ7uzuR7aNuSnJ7yX5se7+elVd291f3uW5eq/XA9ZLVaW7\n67DnsV/2Lzh+9rN/LXOCdWuS8939eHc/neRMkpM7xrwxyW9099eTZLe4AgA4LpYJrOuSPLHt/oXF\nY9v9cJIfqapPVNUnq+rE1AQBAFbNNUuM2e1obOc5+TVJbkry6iQ3JvmLqrr5myda250+ffq52xsb\nG9nY2Fh2rsAK2NzczObm5mFP40DYv2C9Te5fy1yD9aokp7v79sX9tyTp7n7rtjH/I8mnuvu9i/v/\nK8k93f2XO57LNQxwzLgGC1hVB30N1oNJbqqqF1fVC5LcmeSDO8Z8IMm/Xkzm2iQvSfJ/LmdCAACr\nbs/A6u5nktyd5P4kZ5Oc6e5zVXVfVd2xGPOnSb5SVWeT/FmSN3f3Vw9w3gAAR9aevyIcfTFH7HDs\n+BUhsKoO+leEAABcAoEFADBMYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExg\nAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExg\nAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExg\nAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwLClAquqbq+qR6rqsaq65yLjXldVz1bVK+am\nCACwWvYMrKq6Ksk7kpxIcnOSU1X10l3GfUeSn0/ywPQkAQBWyTInWLcmOd/dj3f300nOJDm5y7hf\nSfLWJN8YnB8AwMpZJrCuS/LEtvsXFo89p6puSXJ9d394cG4AACvpmiXG1C6P9XNfrKokb0/yhj3+\nDgDAsbBMYF1IcuO2+9cneXLb/e/M1rVZm4vY+v4kf1RV/667H9r5ZKdPn37u9sbGRjY2Ni591sCR\ntbm5mc3NzcOexoGwf8F6m9y/qrsvPqDq6iSPJrktyVNJPpPkVHefe57xf57kl7r7s7t8rfd6PWC9\nVFW6e+VPte1fcPzsZ//a8xqs7n4myd1J7k9yNsmZ7j5XVfdV1R27/ZX4FSEAcIzteYI1+mJ+AoRj\nxwkWsKoO9AQLAIBLI7AAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgks\nAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGLRVYVXV7VT1SVY9V1T27fP0Xq+ps\nVT1cVR+tqhvmpwoAsBr2DKyquirJO5KcSHJzklNV9dIdwx5K8sruviXJHyb51emJAgCsimVOsG5N\ncr67H+/up5OcSXJy+4Du/nh3/+Pi7gNJrpudJgDA6lgmsK5L8sS2+xdy8YC6K8lH9jMpAIBVds0S\nY2qXx3rXgVWvT/LKJD+6n0kBAKyyZQLrQpIbt92/PsmTOwdV1WuS/HKSVy9+lbir06dPP3d7Y2Mj\nGxsbS04VWAWbm5vZ3Nw87GkcCPsXrLfJ/au6dz2M+taAqquTPJrktiRPJflMklPdfW7bmJcneX+S\nE9391xd5rt7r9YD1UlXp7t1OwleK/QuOn/3sX3teg9XdzyS5O8n9Sc4mOdPd56rqvqq6YzHsbUm+\nPcn7q+qzVfWBy5kMAMA62PMEa/TF/AQIx44TLGBVHegJFgAAl0ZgAQAME1gAAMMEFgDAMIEFADBM\nYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBM\nYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBM\nYAEADBNYAADDBBYAwDCBBQAwTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBM\nYAEADFsqsKrq9qp6pKoeq6p7dvn6C6rqTFWdr6pPVdWN81M9WjY3Nw97CmPWZS3rso5kvdbC0bNO\n7y9rOXrWZR37tWdgVdVVSd6R5ESSm5OcqqqX7hh2V5K/7+6XJPm1JG+bnuhRs05voHVZy7qsI1mv\ntXD0rNP7y1qOnnVZx34tc4J1a5Lz3f14dz+d5EySkzvGnEzynsXtP0hy29wUAQBWyzKBdV2SJ7bd\nv7B4bNcx3f1Mkq9V1feMzBAAYMVUd198QNXrkvyb7v65xf3XJ/kX3f1ft43534sxTy7uf2Ex5qs7\nnuviLwaspe6uw57Dftm/4Hi63P3rmiXGXEiy/aL165M8uWPME0luSPJkVV2d5EU742o/kwQ4bPYv\n4FIs8yvCB5PcVFUvrqoXJLkzyQd3jPnjJG9Y3P7JJB+bmyIAwGrZ8wSru5+pqruT3J+tIHtXd5+r\nqvuSPNjdf5LkXUl+u6rOJ/lKtiIMAOBY2vMaLAAALo1PcgcAGCawAACGCSwAgGECCwBgmMACABgm\nsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgm\nsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgm\nsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgm\nsAAAhu0ZWFX1rqr6UlV97iJjfr2qzlfVw1V1y+wUAQBWyzInWO9OcuL5vlhVr03yQ939kiRvSvKb\nQ3MDAFhJewZWd38iyVcvMuRkkvcuxn46yXdV1ffNTA8AYPVMXIN1XZIntt3/4uIxAIBj6ZqB56hd\nHutdB1bt+jiw3rp7t31ipdi/4Hi63P1r4gTrQpIbtt2/PsmTzze4u9fiz7333nvoc7CW9VzHuq1l\nnRz2P0vvL2tZhT/rso7u/e1fywZWZfeTqiT5YJKfSZKqelWSr3X3l/Y1KwCAFbbnrwir6neTbCT5\n3qr6uyT3JnlBku7u3+ruD1fVj1fVF5L8Q5KfPcgJAwAcdXsGVnf/9BJj7p6ZzurY2Ng47CmMWZe1\nrMs6kvVaC0fPOr2/rOXoWZd17Fft93eMl/RiVX0lXw84fFWVXpOL3O1fcLzsZ//yn8oBABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYNhSgVVVt1fVI1X1WFXds8vXb6iqj1XVQ1X1cFW9dn6qAACrobr74gOqrkryWJLb\nkjyZ5MEkd3b3I9vGvDPJQ939zqp6WZIPd/cP7vJcvdfrAeulqtLdddjz2C/7Fxw/+9m/ljnBujXJ\n+e5+vLufTnImyckdY55N8qLF7e9O8sXLmQwAwDq4Zokx1yV5Ytv9C9mKru3uS3J/Vf1Ckm9L8pqZ\n6QEArJ5lTrB2OxrbeU5+Ksm7u/uGJD+R5H37nRgAwKpa5gTrQpIbt92/PlvXYm13V5ITSdLdD1TV\nC6vq2u7+8s4nO3369HO3NzY2srGxcYlTBo6yzc3NbG5uHvY0DoT9C9bb5P61zEXuVyd5NFsXuT+V\n5DNJTnX3uW1jPpTk97v7PYuL3D/a3dfv8lwuEoVjxkXuwKo60Ivcu/uZJHcnuT/J2SRnuvtcVd1X\nVXcshr05yRur6uEkv5PkDZczGQCAdbDnCdboi/kJEI4dJ1jAqjroj2kAAOASCCwAgGECCwBgmMAC\nABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMAC\nABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMAC\nABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMAC\nABgmsAAAhgksAIBhSwVWVd1eVY9U1WNVdc/zjPmpqjpbVZ+vqvfNThMAYHVUd198QNVVSR5LcluS\nJ5M8mOTO7n5k25ibkvxekh/r7q9X1bXd/eVdnqv3ej1gvVRVursOex77Zf+C42c/+9cyJ1i3Jjnf\n3Y9399NJziQ5uWPMG5P8Rnd/PUl2iysAgONimcC6LskT2+5fWDy23Q8n+ZGq+kRVfbKqTkxNEABg\n1VyzxJjdjsZ2npNfk+SmJK9OcmOSv6iqm795ogUAcJwsE1gXshVN33R9tq7F2jnmU939bJK/rapH\nk7wkyV/ufLLTp08/d3tjYyMbGxuXNmPgSNvc3Mzm5uZhT+NA2L9gvU3uX8tc5H51kkezdZH7U0k+\nk+RUd5/bNubE4rH/WFXXZiusbunur+54LheJwjHjIndgVR3oRe7d/UySu5Pcn+RskjPdfa6q7quq\nOxZj/jTJV6rqbJI/S/LmnXEFAHBc7HmCNfpifgKEY8cJFrCqDvpjGgAAuAQCCwBgmMACABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYJjAAgAYJrAAAIYJLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAA\nhgksAIBhAgsAYNhSgVVVt1fVI1X1WFXdc5Fxr6uqZ6vqFXNTBABYLXsGVlVdleQdSU4kuTnJqap6\n6S7jviPJzyd5YHqSAACrZJkTrFuTnO/ux7v76SRnkpzcZdyvJHlrkm8Mzg8AYOUsE1jXJXli2/0L\ni8eeU1W3JLm+uz88ODcAgJV0zRJjapfH+rkvVlWStyd5wx5/J0ly+vTp525vbGxkY2NjiSkAq2Jz\nczObm5uHPY0DYf+C9Ta5f1V3X3xA1auSnO7u2xf335Kku/uti/svSvKFJP83W2H1/Um+kuTfdfdD\nO56r93o9YL1UVbr7eX/oWhX2Lzh+9rN/LRNYVyd5NMltSZ5K8pkkp7r73POM//Mkv9Tdn93lazYo\nOGYEFrCq9rN/7XkNVnc/k+TuJPcnOZvkTHefq6r7quqO3f5KLvIrQgCAdbfnCdboi/kJEI4dJ1jA\nqjrQEywAAC6NwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJ\nLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJ\nLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBgmMACABgmsAAAhgksAIBhAgsAYJjAAgAYJrAAAIYJ\nLACAYQILAGCYwAIAGCawAACGCSwAgGECCwBg2FKBVVW3V9UjVfVYVd2zy9d/sarOVtXDVfXRqrph\nfqoAAKthz8CqqquSvCPJiSQ3JzlVVS/dMeyhJK/s7luS/GGSX52eKADAqljmBOvWJOe7+/HufjrJ\nmSQntw/o7o939z8u7j6Q5LrZaQIArI5lAuu6JE9su38hFw+ou5J8ZD+TAgBYZdcsMaZ2eax3HVj1\n+iSvTPKj+5kUAMAqWyawLiS5cdv965M8uXNQVb0myS8nefXiV4m7On369HO3NzY2srGxseRUgVWw\nubmZzc3Nw57GgbB/wXqb3L+qe9fDqG8NqLo6yaNJbkvyVJLPJDnV3ee2jXl5kvcnOdHdf32R5+q9\nXg9YL1WV7t7tJHyl2L/g+NnP/rXnNVjd/UySu5Pcn+RskjPdfa6q7quqOxbD3pbk25O8v6o+W1Uf\nuJzJAACsgz1PsEZfzE+AcOw4wQJW1YGeYAEAcGkEFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAw\nTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAw\nTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAw\nTGABAAwTWAAAwwQWAMAwgQUAMExgAQAME1gAAMMEFgDAMIEFADBMYAEADBNYAADDBBYAwLClAquq\nbq+qR6rqsaq6Z5evv6CqzlTV+ar6VFXdOD9VAIDVsGdgVdVVSd6R5ESSm5OcqqqX7hh2V5K/7+6X\nJPm1JG+bnuhRs7m5edhTGLMua1mXdSTrtRaOnnV6f1nL0bMu69ivZU6wbk1yvrsf7+6nk5xJcnLH\nmJNJ3rO4/QdJbpub4tG0Tm+gdVnLuqwjWa+1cPSs0/vLWo6edVnHfi0TWNcleWLb/QuLx3Yd093P\nJPlaVX3PyAwBAFbMMoFVuzzWe4ypXcYAABwL1X3xDqqqVyU53d23L+6/JUl391u3jfnIYsynq+rq\nJE919z/f5blEFxxD3b3bD2orxf4Fx9Pl7l/XLDHmwSQ3VdWLkzyV5M4kp3aM+eMkb0jy6SQ/meRj\nk5MEOGz2L+BS7BlY3f1MVd2d5P5s/UrxXd19rqruS/Jgd/9Jkncl+e2qOp/kK9mKMACAY2nPXxEC\nAHBpDuST3Nflg0mXWMcvVtXZqnq4qj5aVTccxjyXsddato17XVU9W1WvuJLzuxTLrKWqfmrxvfl8\nVb3vSs9xWUu8x26oqo9V1UOL99lrD2Oee6mqd1XVl6rqcxcZ8+uLf+cfrqpbruT8LsW67F+JPexK\nzm9Z9q+j58D2r+4e/ZOtaPtCkhcn+SdJHk7y0h1j/nOS/764/e+TnJmexxVax48meeHi9n86iutY\ndi2Lcd+R5ONJPpnkFYc97318X25K8pdJXrS4f+1hz3sfa3lnkjctbr8syd8c9ryfZy3/KsktST73\nPF9/bZIPLW7/yyQPHPac9/E9OfL71yWsxR52xNZh/zqUtRzI/nUQJ1jr8sGke66juz/e3f+4uPtA\n/v/PBzsqlvmeJMmvJHlrkm9cycldomXW8sYkv9HdX0+S7v7yFZ7jspZZy7NJXrS4/d1JvngF57e0\n7v5Ekq9eZMjJJO9djP10ku+qqu+7EnO7ROuyfyX2sKPI/nUEHdT+dRCBtS4fTLrMOra7K8lHDnRG\nl2/PtSyOPK/v7g9fyYldhmW+Lz+c5Eeq6hNV9cmqOnHFZndpllnLfUn+Q1U9keRPkvz8FZrbtJ1r\n/WKO5v+Yr8v+ldjDjiL712q6rP1rmY9puFTr8sGky6xja2DV65O8MlvH7UfRRddSVZXk7dn6qI2L\n/Z2jYJnvyzXZOmZ/dZIbk/xFVd38zZ8Ij5Bl1nIqybu7++2Lz6R7X7b+m6CrZul/nw7ZuuxfiT3s\nKO5h9q9jtH8dxAnWhWy9Kb7p+iRP7hjzRJIbkqS2Ppj0Rd19seO5w7DMOlJVr0nyy0n+7eKY9Cja\nay3fma03/WZV/U2SVyX5oyN6kegy35cLSf6ou5/t7r9N8miSl1yZ6V2SZdZyV5LfT5LufiDJC6vq\n2iszvVEXsvh3fmHXf5+OgHXZvxJ72FHcw+xfx2n/OoCLxa7Oty58e0G2Lnx72Y4x/yXfukj0zhzB\nCyuXXMfLF2N+6LDnu9+17Bj/50leftjz3sf35USS/7m4fW2Sx5P8s8Oe+2Wu5UNJ3rC4/bIkFw57\n3hdZzw8k+fzzfO3H862LRF+Vo3uR+1rsX5ewFnvYEVuH/evQ1jO+fx3URG/PVnWfT/KWxWP3Jblj\ncfufZqtqz2frwsofOOx/uJe5jo9m69PtH0ry2SQfOOw5X+5adoz9WI7g/wPnUtaS5L8lOZvkr5L8\n5GHPeR/vsZcl+cRi83ooyW2HPefnWcfvZusnum8k+bskP5vkTUl+btuYdyw25L9a5ffXquxfS67F\nHnYE12H/uuLrOJD9yweNAgAMO5APGgUAOM4EFgDAMIEFADBMYAEADBNYAADDBBYAwDCBBQAwTGAB\nAAz7fwLd5806nddmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa4a5e60160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# source likelihood given observation\n",
    "# using different number of constraints\n",
    "K = 100\n",
    "M = 5  # number of pairwise order constraints\n",
    "fig, ax = plt.subplots(2, 2, num=1,\n",
    "                       sharex=True, sharey=True,\n",
    "                       figsize=(10, 10))\n",
    "\n",
    "\n",
    "kwargs = get_display_kwargs(g, infection_times, obs_nodes)\n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[0, 0], **kwargs)\n",
    "# add_colorbar(kwargs['node_color'])\n",
    "ax[0, 0].set_title('infection time')\n",
    "\n",
    "def run_with_n_constraints(M, i, j):\n",
    "    kwargs = get_display_kwargs(g, infection_times, obs_nodes=set())\n",
    "    mu = estimate_source_likelihood(g, K, M, obs_nodes, infection_times)\n",
    "    kwargs['node_color'] = [mu.get(n, 0) for n in g.nodes()]\n",
    "    nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[i, j], **kwargs)\n",
    "    ax[i, j].set_title('estimated source likeihood (with {} constraint)'.format(M))\n",
    "# add_colorbar(kwargs['node_color'])\n",
    "run_with_n_constraints(0, 0, 1)\n",
    "run_with_n_constraints(5, 1, 0)\n",
    "run_with_n_constraints(10, 1, 1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def simulate_given_passing_edge(g, q, u, s, K, M, max_attempts, infection_times, obs_nodes=set()):\n",
    "    \"\"\"\n",
    "    q: query node\n",
    "    u: direction node\n",
    "    s: source node\n",
    "    K: valid cascade number\n",
    "    M: number of order pairs\"\"\"\n",
    "    attempts = 0\n",
    "    valid_cnt = 0  # valid if q and u is infected\n",
    "    success_cnt = 0\n",
    "    if M:\n",
    "        ordered_pairs = random.sample(list(combinations(obs_nodes, 2)), M)\n",
    "        ground_truth = np.array([infection_times[u] > infection_times[v] for u, v in ordered_pairs])\n",
    "\n",
    "    \n",
    "    while valid_cnt < K and attempts < max_attempts:\n",
    "        attempts += 1\n",
    "        time, tree = simulate_IC(g, s=s)\n",
    "        if not np.isinf(time[q]) and not np.isinf(time[u]):\n",
    "            valid_cnt +=1 \n",
    "            if tree.has_edge(q, u) and time[q] - time[u] == g[q][u]['d']:\n",
    "                inf_nodes = set(filter(lambda n: not np.isinf(time[n]), time))\n",
    "                if obs_nodes.issubset(inf_nodes):\n",
    "                    if M:\n",
    "                        reality = np.array([time[u] > time[v] for u, v in ordered_pairs])\n",
    "                        if (ground_truth == reality).all():\n",
    "                            success_cnt += 1\n",
    "                    else:\n",
    "                        success_cnt += 1\n",
    "    return s, success_cnt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def estimate_source_likelihood_given_edge(g, q, u, K, M,\n",
    "                                          max_attempts, infection_times, obs_nodes):\n",
    "    cands = set(g.nodes()) - {q} - set(g.neighbors(q)) - obs_nodes\n",
    "    result = Parallel(n_jobs=-1)(delayed(simulate_given_passing_edge)(g, q, u, s=s, \n",
    "                                                                        K=K, M=M,\n",
    "                                                                        max_attempts=max_attempts,\n",
    "                                                                        infection_times=infection_times,\n",
    "                                                                        obs_nodes=obs_nodes)\n",
    "                                    for s in cands)\n",
    "    df = pd.DataFrame(result, columns=['node', 'cnt'])\n",
    "    df['mu'] = df['cnt'] / K\n",
    "\n",
    "    # stored here\n",
    "    mu = {r['node']: r['mu'] for i, r in df.iterrows()}\n",
    "    return mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# fixing query and edge\n",
    "# check the likelihood change with different number of node order constraints.\n",
    "K = 50\n",
    "M=5\n",
    "max_attempts=200\n",
    "q = random.choice(list(set(tree.nodes()) - obs_nodes))\n",
    "print(q)\n",
    "u = [u for u in g.neighbors(q)\n",
    "     if infection_times[q] - infection_times[u] == g[q][u]['d']]\n",
    "u = u[0]\n",
    "print(u)\n",
    "\n",
    "fig, ax = plt.subplots(3, 2, num=1, sharex=True, sharey=True,\n",
    "                       figsize=(10, 15))\n",
    "\n",
    "\n",
    "kwargs = get_display_kwargs(g, infection_times, obs_nodes)\n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[0, 0], **kwargs)\n",
    "# add_colorbar(kwargs['node_color'])\n",
    "ax[0, 0].set_title('infection time (+ observation)')\n",
    "\n",
    "kwargs = get_display_kwargs(g, infection_times, {q, u})\n",
    "mu = estimate_source_likelihood(g, K, M, obs_nodes, infection_times)\n",
    "kwargs['node_color'] = [mu.get(n) for n in g.nodes()]\n",
    "nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[0, 1], **kwargs)\n",
    "# add_colorbar(kwargs['node_color'])\n",
    "ax[0, 1].set_title('infection time (+ query and response)')\n",
    "\n",
    "def run_with_n_constraints(M, i, j):\n",
    "    kwargs = get_display_kwargs(g, infection_times, obs_nodes=set())\n",
    "    mu = estimate_source_likelihood_given_edge(g, q, u, K, M, max_attempts, infection_times, obs_nodes)\n",
    "    kwargs['node_color'] = [mu.get(n, 0) for n in g.nodes()]\n",
    "    nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[i, j], **kwargs)\n",
    "    ax[i, j].set_title('estimated source likeihood (with {} constraint)'.format(M))\n",
    "\n",
    "run_with_n_constraints(0, 1, 0)\n",
    "run_with_n_constraints(5, 1, 1)\n",
    "run_with_n_constraints(10, 2, 0)\n",
    "run_with_n_constraints(15, 2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "print(infection_times[q])\n",
    "print(infection_times[u])\n",
    "q, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def get_nodes_on_and_off_direction(g, start, end):\n",
    "    pos_nodes = set()\n",
    "    neg_nodes = set()\n",
    "    for t, path in nx.shortest_path(g, source=start).items():\n",
    "        if len(path) == 1:\n",
    "            pos_nodes.add(t)\n",
    "        elif path[1] == end:\n",
    "            # print(t, path)\n",
    "            pos_nodes.add(t)\n",
    "        else:\n",
    "            neg_nodes.add(t)\n",
    "    return pos_nodes, neg_nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "q = (9, 7)\n",
    "u = (8, 7)\n",
    "pos_nods, neg_nodes = get_nodes_on_and_off_direction(g, q, u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "node_color = [1 if n == q else (3 if n == u else 0)\n",
    "              for n in g.nodes() ]\n",
    "nx.draw_networkx(g, pos=pos, \n",
    "                 node_color=node_color,\n",
    "                 node_size=[400 if n in pos_nodes else 200 for n in g.nodes()],\n",
    "                 with_labels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# multiplicative weight update method\n",
    "# initial mu\n",
    "def mwu(g, orig_obs_nodes,\n",
    "        infection_times, \n",
    "        p_on_mu, \n",
    "        mu, \n",
    "        max_iter, \n",
    "        max_plot_iter=11,\n",
    "        debug=False, plot=False):\n",
    "    use_sampling = False\n",
    "    epsilon = 1e-5  # smoothing parameter\n",
    "    obs_nodes = copy(orig_obs_nodes)\n",
    "    K = 25\n",
    "    M = len(obs_nodes)\n",
    "    \n",
    "    source = min(infection_times, key=lambda n: infection_times[n])\n",
    "    \n",
    "\n",
    "    K = 50\n",
    "    max_attempts = 500\n",
    "\n",
    "    sp_len = nx.shortest_path_length(g, weight='d')\n",
    "\n",
    "    iter_i = 0\n",
    "\n",
    "    if plot:\n",
    "        ncols = 2\n",
    "        nrows = math.ceil((max_plot_iter + 1) / ncols)\n",
    "        width = 12\n",
    "        height = int(width * nrows / ncols)\n",
    "        fig, ax = plt.subplots(nrows, ncols,\n",
    "                               sharex=True, sharey=True,\n",
    "                               figsize=(width, height))\n",
    "        kwargs = get_display_kwargs(g, infection_times, obs_nodes=obs_nodes)\n",
    "        nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[0, 0], **kwargs)\n",
    "        ax[0, 0].set_title('ground truth')\n",
    "\n",
    "    cands = set(g.nodes())\n",
    "    if debug:\n",
    "        print('correct source is {}'.format(source))\n",
    "    all_nodes = set(g.nodes())\n",
    "    queried_nodes = set()\n",
    "    while iter_i < max_iter:\n",
    "        def phi(q, g, mu, sp_len):\n",
    "            return np.sum([mu[v] * sp_len[q][v] for v in g.nodes_iter()])\n",
    "        q = min(all_nodes, key=lambda q: phi(q, g, mu, sp_len))\n",
    "        if debug:\n",
    "            print('query node: {}'.format(q))\n",
    "\n",
    "        if plot and iter_i+1 <= max_plot_iter:\n",
    "            kwargs = get_display_kwargs(g, mu, obs_nodes={q})\n",
    "            likelihood_colors = [mu.get(n, 0) for n in g.nodes()]\n",
    "\n",
    "            i, j = int((iter_i + 1) / ncols), (iter_i + 1) % ncols\n",
    "            kwargs['node_color'] = likelihood_colors\n",
    "            nx.draw_networkx(g, pos=pos, with_labels=False, ax=ax[i, j], **kwargs)\n",
    "            ax[i, j].set_title('iter_i={}'.format(iter_i))\n",
    "\n",
    "        if q not in obs_nodes:\n",
    "            obs_nodes.add(q)\n",
    "            queried_nodes.add(q)\n",
    "\n",
    "        if np.isinf(infection_times[q]):  # q is uninfected        \n",
    "            if debug:\n",
    "                print('query is uninfected')\n",
    "            mu[q] = 0\n",
    "        else:        \n",
    "            if debug:\n",
    "                print('query is infected')\n",
    "            found_source = True\n",
    "            for u in g.neighbors(q):                \n",
    "                if u not in obs_nodes:  # need to query it\n",
    "                    obs_nodes.add(u)\n",
    "                    queried_nodes.add(u)\n",
    "\n",
    "                if infection_times[q] > infection_times[u]:\n",
    "                    # q later than u                    \n",
    "                    found_source = False\n",
    "                    # found the direction\n",
    "                    if infection_times[q] - infection_times[u] == g[q][u]['d']:\n",
    "                        mu[q] = 0\n",
    "                        if q in cands:\n",
    "                            cands.remove(q)\n",
    "                        break\n",
    "                else:\n",
    "                    # q earlier than u\n",
    "                    mu[u] = 0\n",
    "                    if u in cands:\n",
    "                        cands.remove(u)\n",
    "            \n",
    "            if found_source:\n",
    "                source = q\n",
    "                expected = min(infection_times, key=lambda n: infection_times[n])\n",
    "                print('**Found** source: {}'.format(source))\n",
    "                assert source == expected\n",
    "                print('query count {}'.format(len(queried_nodes)))\n",
    "                break\n",
    "            else:\n",
    "                assert infection_times[q] - infection_times[u] == g[q][u]['d']\n",
    "                if use_sampling:\n",
    "                    # update likelihood                      \n",
    "                    if debug:\n",
    "                        print('M={}'.format(len(obs_nodes)))                    \n",
    "                    result = Parallel(n_jobs=-1)(delayed(simulate_given_passing_edge)(g, q=q, u=u, s=s,\n",
    "                                                                                      obs_nodes=obs_nodes,\n",
    "                                                                                      K=K, M=len(obs_nodes),\n",
    "                                                                                      max_attempts=max_attempts,\n",
    "                                                                                      infection_times=infection_times)\n",
    "                                                 for s in cands)\n",
    "                    df = pd.DataFrame(result, columns=['node', 'cnt'])\n",
    "                    df['p'] = df['cnt'] / K\n",
    "                    # update mu\n",
    "                    for i, r in df.iterrows():\n",
    "                        mu[r['node']] *= r['p']\n",
    "                        mu[r['node']] += epsilon  # smoothing\n",
    "                else:\n",
    "                    pos_nodes, neg_nodes = get_nodes_on_and_off_direction(g, q, u)\n",
    "                    for n in pos_nodes:\n",
    "                        mu[n] *= p_on_mu ** sp_len[u][n]\n",
    "                    for n in neg_nodes:\n",
    "                        mu[n] *= (1 - p_on_mu) ** sp_len[u][n]\n",
    "        total = sum(mu.values())\n",
    "        mu = {n: (v / total) for n, v in mu.items()}  # normalize it\n",
    "        # print(mu)\n",
    "        if debug:\n",
    "            print('mu[{}] = {}'.format(source, mu[source]))\n",
    "        likelihood_colors = [mu.get(n, 0) for n in g.nodes()]\n",
    "        iter_i += 1\n",
    "    return len(queried_nodes)\n",
    "# def multiplicated_weight_update(g, mu, infection_time):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def baseline_dog_tracker(g, obs_nodes, infection_times):\n",
    "    query_count = 0\n",
    "    q = min(obs_nodes, key=lambda n: infection_times[n])\n",
    "    while True:\n",
    "        if q not in obs_nodes:\n",
    "            query_count += 1\n",
    "            obs_nodes.add(q)\n",
    "            \n",
    "        found_source = True\n",
    "        for u in g.neighbors(q):\n",
    "            if u not in obs_nodes:  # need to query it\n",
    "                obs_nodes.add(u)\n",
    "                query_count += 1            \n",
    "\n",
    "            if infection_times[q] > infection_times[u]:\n",
    "                # q later than u                    \n",
    "                found_source = False\n",
    "                # found the direction\n",
    "                if infection_times[q] - infection_times[u] == g[q][u]['d']:\n",
    "                    q = u\n",
    "                    break\n",
    "        if found_source:\n",
    "            source = q\n",
    "            expected = min(infection_times, key=lambda n: infection_times[n])\n",
    "            print('**Found** source: {}'.format(source))\n",
    "            assert source == expected            \n",
    "            break\n",
    "    return query_count\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "infp = 0.7\n",
    "g = prepare_standard_infection_graph(g, infp)\n",
    "\n",
    "p_on_mu = abs(infp-0.5) + 0.5\n",
    "fraction = 0.05\n",
    "K = 100\n",
    "M = 5\n",
    "N = 50\n",
    "cnts = []\n",
    "inputs = []\n",
    "for i in range(N):\n",
    "    print('experiment {}'.format(i+1))\n",
    "    obs_nodes, infection_times, tree = make_input(g, infp, fraction, sampling_method='late_nodes')\n",
    "    mu = estimate_source_likelihood(g, K, M, obs_nodes, infection_times)\n",
    "    \n",
    "    inputs.append(\n",
    "        (obs_nodes, infection_times, tree, mu)\n",
    "    )\n",
    "\n",
    "    c = mwu(g, obs_nodes,\n",
    "            infection_times, \n",
    "            p_on_mu,\n",
    "            mu,\n",
    "            max_iter=float('inf'), debug=False)\n",
    "    cnts.append(c)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# manual inspection\n",
    "bad_example1 = 11  # zero division\n",
    "bad_example2 = 2  #  many querys\n",
    "normal_example = 1\n",
    "obs, times, tree, mu = inputs[normal_example]\n",
    "# mu = estimate_source_likelihood(g, K, M, obs_nodes, infection_times)\n",
    "mwu(g, obs, times, \n",
    "    p_on_mu=abs(infp-0.5) + 0.5,\n",
    "    mu=mu,\n",
    "    max_iter=float('inf'), debug=True, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "s1 = pd.Series(cnts)\n",
    "s1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "cnts = []\n",
    "for i in tqdm(range(N)):\n",
    "    obs_nodes, infection_times, tree = make_input(g, infp, fraction, sampling_method='late_nodes')\n",
    "    c = baseline_dog_tracker(g, obs_nodes, infection_times)\n",
    "    cnts.append(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "s2 = pd.Series(cnts)\n",
    "s2.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "print('late nodes sampling')\n",
    "cmp = pd.DataFrame.from_dict({'multiplicative': s1.describe(), 'baseline': s2.describe()})\n",
    "cmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "print('unfiorm sampling')\n",
    "cmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "infection_times = obs_infection_times\n",
    "nodes = obs_nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def reverse_infection():\n",
    "    u = max(nodes, key=lambda n: infection_times[n])\n",
    "    tmax = infection_times[u]\n",
    "    O = copy(nodes)\n",
    "    infection_times = copy(infection_times)\n",
    "    queue = [u]\n",
    "    found = False\n",
    "    time = max(infection_times.values())\n",
    "    m = defaultdict(set)\n",
    "    while not found:\n",
    "        obs_i = filter(lambda n: infection_times[n] == time, O)\n",
    "        for u in obs_i:\n",
    "            for v in g.neighbors(u):\n",
    "                if u not in m[v] or (infection_times[u] < time - g[u][v]['d']):\n",
    "                    m[v]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "nx.ran"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}