sociopatterns.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SocioPatterns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import igraph as ig\n",
    "import dionysus as d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use(\"fivethirtyeight\")\n",
    "plt.rcParams[\"figure.figsize\"] = 10, 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2.0.7.dev0'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = ig.read(\"data/sociopatterns/infectious/infectious.graphml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH UN-- 10972 415912 -- \\n+ attr: id (v), name (v), id (e), time (e)'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(max(g.es[\"time\"]) - min(g.es[\"time\"])) // (3600 * 24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temporal partitioning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "from zigzag import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "wins = sliding_windows(g, 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(wins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zigzag persistence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqcAAAF9CAYAAAA5qMHQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt4lPWd9/HPMIEQQg4QcuIQUIicCkXQBGVFBQRiapEConVXlwpBqFvhKUEotG5buuG0HFppPIB2H9d21YAKBcH6EE8BEy1iEBGj4CykISeceCckZAjz/JEyOhUyd2AmmcP7dV1eV5n7O8NvOj8mn/xOt8VutzsFAAAA+IEO7d0AAAAA4ALCKQAAAPwG4RQAAAB+g3AKAAAAv0E4BQAAgN8gnAIAAMBvEE4BAADgNwinl1BSUtLeTYAfoB9Aoh+gGf0AF9AXfItwCgAAAL9BOAUAAIDfIJwCAADAbxBOAQAA4DcIpwAAAPAbhFMAAAD4DcIpAAAA/AbhFAAAAH6DcAoAAAC/QTgFAACA3whr7wYAAID2ZTMcWnHAUNmZJiV3sWr5yCj1jerY3s1CiCKcAgAQwmyGQ3fuqdZxo8n12PuVjXp5UhwBFe2CaX0AAELY0sIat2AqSceNJq04YLRTixDqCKcAAISogrJ67Tpx9qLXvjDOtXFrgGaEUwAAQlBBWb2+t/v0Ja+X159vw9YAXzMVTgsKCnT33Xdr8ODBio2N1XPPPee65nA49Oijj+rGG29Uz549NXDgQM2ePVsnTpxwe42zZ88qOztbV199tXr27Km7775bpaWlbjUnTpzQzJkz1bNnT1199dVavHixGhsbvfA2AQDABQVl9bpj92k5W6hJ6Gxps/YA32QqnNbV1WnIkCFauXKlIiIi3K6dOXNGH374oRYtWqQ333xTf/zjH1VaWqrp06fr3LmvpwSWLl2qHTt2aMuWLdq1a5cMw9DMmTPV1NS8zqWpqUkzZ85UbW2tdu3apS1btmj79u1atmyZF98uAAChzWY4NO210/I0LnpVNJuh0D5M7dafOHGiJk6cKEmaP3++27WYmBi9/PLLbo+tX79eo0eP1tGjRzV06FDV1NTo2Wef1aZNm3TrrbdKkp544gkNGzZMb7zxhsaPH6+9e/fqyJEjOnTokHr37i1J+uUvf6mf/OQn+vnPf67o6OgrfrMAAISygrJ6Tf/LaTV4SKYRVmn5yKi2aRTwD3yy5tQwmnf4xcbGSpIOHjwoh8OhcePGuWp69+6tgQMHqrCwUJJUVFSkgQMHuoKpJI0fP15nz57VwYMHfdFMAABCxtbPDWXuPq36ppbrOneQ8m7rzjFSaDdeP+e0sbFRy5cv1+TJk9WrVy9JUkVFhaxWq+Li4txq4+PjVVFR4aqJj493ux4XFyer1eqquZiSkhIvv4O2eW0EDvoBJPoBmgVqP/jrlx304OFwSS2vI+0e1qSnv9uohNqTCtC32mYCtS/4g9TU1BavezWcnjt3TllZWaqpqdGf/vQnb770JXl6g5erpKTEZ6+NwEE/gEQ/QLNA7QcFZfWa/86ld+VfkNhZeu17yYyYmhCofSFQeG1a/9y5c3rggQd0+PBhvfLKK+revbvrWkJCgpqamlRdXe32nMrKSiUkJLhqKisr3a5XV1erqanJVQMAAMy7sCu/pSWmHSTd3idcr30vgWAKv+CVcOpwODRr1iwdPnxYO3bsUGJiotv1ESNGqGPHjsrPz3c9VlpaqqNHjyo9PV2SlJaWpqNHj7odL5Wfn6/w8HCNGDHCG80EACBkPP6RXZkegqlF0o7J3fXHCT0IpvAbpqb1a2trdezYMUnS+fPndfLkSRUXF6tbt25KTk7W/fffrw8++EB/+tOfZLFYVF5eLkmKjo5WRESEYmJi9C//8i969NFHFR8fr27dumnZsmUaOnSobrnlFknSuHHjNHjwYD344INasWKFvvzyS/3iF7/Qfffdx059AABa4TfvV2vNoYYWayySNo+N1pjkiBbrgLZmKpx+8MEHuuOOO1x/zsnJUU5Oju655x4tWbJEu3btkiRX0Lxg06ZNuvfee13PsVqtmjVrlhoaGjR27Fg9/vjjslqtkiSr1arnn39eixYt0uTJk9W5c2fNmDFDv/71r73xPgEACAmPf2T3GEytFmn7pO4EU/gli91ub+kGESGLxc6Q6AdoRj+AFBj9wMyIKcH0ygVCXwhkPjnnFAAAtC0zwVQimML/EU4BAAhwZqbyJSnaKoIp/B7hFACAAFZQVq8l79WZqv3Tbd09FwHtjHAKAEAAm/u23WNNt04W7ZzMdD4Cg9dvXwoAAHzLZji04oChE7XndLKupZNMpexhnbXsurgWawB/QjgFACCA2AyH7txTreNGk8dagikCEeEUAIAAsrSwxlQwXXl9pB78TmwbtAjwLsIpAAABYuvnhnadOHvRaxZJUWFSbGercv8phvWlCFiEUwAA/JjNcGjJuzUqOHVWX527dN2Ufp31h1uZwkfgI5wCAOCnbIZDmbsqdfJMyzdzDO8g/fK66DZqFeBbHCUFAICfWnHA8BhMJWlcz3D1jerYBi0CfI9wCgCAnzpqd3isCbNIK0fHtEFrgLZBOAUAwA9t/dxQ8ekWFpmqeRPUEzdFM2qKoMKaUwAA/IzNcGjOW1+1WNPZKm29jbs+IfgQTgEA8CM2w6GJf65US/d96t3Fop23xzNiiqBEOAUAwE8UlNVrxl9O60wLZ+x3CSOYIrgRTgEA8AMFZfW6Y/fpFkdMJenFCd0IpghqbIgCAKCdbf3cUKaJYJoWZ2WNKYIe4RQAgHa09XNDD3jY/CQ1rzN9ijtAIQQQTgEAaCdmduV3kHR7n3DWmSJksOYUAIB2YDMcyny1qsWpfKtF2j6J46IQWhg5BQCgjdkMh76/u0on6y4dTS0imCI0MXIKAEAbujBi2lIwlaTNY6MJpghJhFMAANqImXNMrRbpyZuiNa1/VNs1DPAjhFMAAHzIZji04oChwlP1+t8zLdemdLVqx+Q4Nj4hpBFOAQDwEZvh0J17qnXcaGGo9O+uirLq5UkEU4BwCgCAl10YLX3jb2dV2eDpaP3mEVOCKdCMcAoAgBe1ZrRUkiLDxFQ+8A0cJQUAgBetOGCYDqaS9MKE7gRT4BsIpwAAeFFZS1vx/8HK6yM5Lgr4B0zrAwBwmS6sLS0706TkLlYtHxml5C5WU89deX2kHvxOrI9bCAQewikAAJfhYmtL95efVUInz8+d2S+MYApcAuEUAIDLcLG1pSfrzutkXcvPuypCeuLWRB+2DAhsrDkFAOAytGZt6QWRYdLLmQk+aA0QPBg5BQDgMphdW3pBZJhFL0zoxs58wANGTgEAuAzLR0YpKcLcj9GUrlbtuzOenfmACYRTAAAuw/sVDTpV7/nuT70jO3DIPtAKhFMAAFqpoKxeD7z1VYs1HSTdktxJOzN6EEyBVmDNKQAArXDva2XaWdryiGlKVyujpcBlIpwCAGDS1J2lyq/wXEcwBS4f0/oAAJiQ+3kHU8G0o0QwBa6AqXBaUFCgu+++W4MHD1ZsbKyee+45t+tOp1M5OTkaNGiQkpKSlJmZqSNHjrjV2O12ZWVlKSUlRSkpKcrKypLdbnerOXz4sG6//XYlJSVp8ODBWrVqlZxO5xW+RQAArsxv3q/W02XhpmofHxvt49YAwc1UOK2rq9OQIUO0cuVKRUR8+xiMjRs3atOmTVq1apX27t2r+Ph4TZ06VYZhuGpmz56t4uJi5eXlKS8vT8XFxZo7d67r+ldffaWpU6cqISFBe/fu1cqVK/W73/1Ojz32mBfeJgAAl8dmOLTmUIMki8fa7GGdNa1/lO8bBQQxU2tOJ06cqIkTJ0qS5s+f73bN6XQqNzdXCxYs0JQpUyRJubm5Sk1NVV5enmbNmqWjR4/q9ddf1+7du5WWliZJWr9+vTIyMlRSUqLU1FS9+OKLqq+vV25uriIiIjRkyBB9+umn+v3vf6+HHnpIFovnLwUAALzJZjiU+WqVqdqdk7tzjingBVe85tRms6m8vFzjxo1zPRYREaEbb7xRhYWFkqSioiJ17dpV6enprprRo0crMjLSreaGG25wG5kdP368ysrKZLPZrrSZAACYZjMcuucvVUrbVqGTdZ7PMiWYAt5zxbv1y8vLJUnx8fFuj8fHx6usrEySVFFRobi4OLfRT4vFoh49eqiiosJV07Nnz2+9xoVr/fr1u+jfX1JScqVv4ZJ8+doIHPQDSPSDUFJab9FDh8N1ssHT+I1TXTs4tXZwoxJqT4ouElr4Trh8qampLV4P+KOkPL3By3VhuQFCG/0AEv0g1Kx+87RONtR7rMseFqFl18W1QYvgb/hO8K0rDqeJiYmSpMrKSvXp08f1eGVlpRISEiRJCQkJqq6ultPpdI2eOp1OVVVVudVUVla6vfaFP1+oAQDAl7ILKvXiscYWazhgH/CtK15z2rdvXyUmJio/P9/1WENDg/bv3+9aY5qWlqba2loVFRW5aoqKilRXV+dWs3//fjU0NLhq8vPzlZycrL59+15pMwEAaNGw50r11KctB9OrogimgK+ZCqe1tbUqLi5WcXGxzp8/r5MnT6q4uFgnTpyQxWLRvHnztHHjRm3fvl0ff/yx5s+fr8jISE2fPl2SNHDgQE2YMEELFy5UUVGRioqKtHDhQk2aNMk1LD59+nRFRERo/vz5+vjjj7V9+3Zt2LBB8+fPZ6c+AMBnbIZD8c+U6kTLuVQxYef18iSCKeBrpqb1P/jgA91xxx2uP+fk5CgnJ0f33HOPcnNz9fDDD6u+vl7Z2dmy2+0aNWqUtm3bpqior89627x5sxYvXqxp06ZJkjIyMrR69WrX9ZiYGL300ktatGiRbr31VsXGxurHP/6xHnroIW+9VwAA3BSU1Stz92lTtTd0ayKYAm3AYrfbuQXTRbDYGRL9AM3oB8HJZjh0bV6FPB8UJYVZpLyR9bpl+ACftwv+j+8E37riNacAAAQam+HQ6G3mgqkkPXFTtHpFMJYDtAXCKQAgpNgMh67Lq1C9yWS6c3J3bkkKtCHCKQAgpNzycoUcJmu58xPQ9ginAICQkV1QqS/PmavdMjaaYAq0g4C/QxQAAJ7YDIcefudLvXHK85hpB0k7GDEF2g3hFAAQ1GyGQ5mvVulknedFpimRFu3IiOfIKKAdEU4BAEHFZji04oChsjNNigqzqKiyUdVnze20L76rp49bB8ATwikAIGjYDIfu3FOt40ZTq5+7c3J3H7QIQGuxIQoAEDRWHDBaHUy7hrErH/AnjJwCAILGccPkVvy/uzmxo165PcFHrQFwORg5BQAEjQqzJ+tL6t3Fot/e1M2HrQFwOQinAICgcf68uSn9W5I7aeft7MoH/BHT+gCAoJBdUKmTZzzXpcVZ9fLkeN83CMBlYeQUABDwHv/Irqc+bfRYl9hZeurWuDZoEYDLxcgpACBgfPMM0+QuVi0fGaW+UR31s/fqPD739j7hykmPYSof8HOEUwBAQLjYGab7y89qUEyYPG2Dujmxo/44oYdvGwjAKwinAICAcLEzTE/WndfJupan82M7il35QABhzSkAICCUnWn9XZ8irdKbUxKYygcCCOEUABAQkrtYW1UfY5X2TSWYAoGGcAoACAjLR0YpKcLzj634zhbNuDpCbxFMgYDEmlMAQEA4WXtO5R7uAHVVlFUvT4ojlAIBjJFTAEBAmPu2Xc4WrlstIpgCQYBwCgDwezWN51V2puVR0ydviiaYAkGAcAoA8Gs1jec17bUqNV1i2NQiacvYaE3rH9Wm7QLgG6w5BQD4rY9On9X3d1fr9NmLJ9Mwi/TKpO4akxzRxi0D4CuEUwCAX7lwi9JD1Q36pObbobRjBymigxTb2arcf4ohmAJBhnAKAPAbNsOhzFerdLLu0utLM1M66w+3xrVhqwC0JdacAgD8xv17q1sMppJU3dDydQCBjZFTAIBfuGVrqQ5+5bkuqZV3igIQWAinAIB2N+C/SlVlYkA0Mqz5TlEAghfT+gCAdlNQVq/YZ8wGU4temNCds0yBIMfIKQCgXWz93NADb5mYx5d0e59w5aTHEEyBEEA4BQC0uYKyetPBdHiMRX+c0MPHLQLgL5jWBwC0KZvh0B27T5uqtUp69rZ43zYIgF9h5BQA0CZshkNLC2u068RZU/XhkoqmJzCVD4QYwikAwOfMHK7/TeGSymf18m2jAPglpvUBAD634oBhOphKzSOmAEITI6cAAJ+5MJX/qsmpfEnaOZnjooBQRjgFAPiEzXBo0s4qnao3P2K6ZWy0xiRH+LBVAPwd0/oAAJ9YWljTqmCa2auDpvXn7k9AqCOcAgB84o2/mZ/Kz+zVQc9NTPZhawAECq+E06amJq1YsULDhw9XYmKihg8frhUrVujcuXOuGqfTqZycHA0aNEhJSUnKzMzUkSNH3F7HbrcrKytLKSkpSklJUVZWlux2uzeaCABoIxduSXqmyVx99rDOBFMALl4Jpxs2bNDmzZu1atUqFRUVaeXKlXrqqae0bt06V83GjRu1adMmrVq1Snv37lV8fLymTp0qwzBcNbNnz1ZxcbHy8vKUl5en4uJizZ071xtNBAC0gYKyemWaPGBfag6my66L82GLAAQar2yIKioq0uTJk5WRkSFJ6tu3rzIyMvTXv/5VUvOoaW5urhYsWKApU6ZIknJzc5Wamqq8vDzNmjVLR48e1euvv67du3crLS1NkrR+/XplZGSopKREqamp3mgqAMCHpr1mLpjGhVv0f2/txuYnAN/ilZHT0aNH65133tGnn34qSfrkk0/09ttv67bbbpMk2Ww2lZeXa9y4ca7nRERE6MYbb1RhYaGk5oDbtWtXpaenu71uZGSkqwYA4L8e/8iuBhP7n7aMjdbnP+xJMAVwUV4ZOV2wYIFqa2uVnp4uq9Wqc+fOadGiRZo9e7Ykqby8XJIUH+9+f+T4+HiVlZVJkioqKhQXFyeLxeK6brFY1KNHD1VUVFzy7y4pKfHGW2jz10bgoB9Aoh94Ulpv0ZK/hsvzmMd5DT9/SiUlp9qiWV5HP8AF9IXL52k23CvhdNu2bfqf//kfbd68WYMGDdKhQ4e0ZMkSpaSk6L777vPGX3FJvpruZykBJPoBmtEPWrb1c0Nz/vqVqdotY2OVGqDHRdEPcAF9wbe8Ek5/8Ytf6KGHHtK0adMkSUOHDtWJEye0fv163XfffUpMTJQkVVZWqk+fPq7nVVZWKiGh+RZ1CQkJqq6ultPpdI2eOp1OVVVVuWoAAP7BZji04oChd0rrVWbyxKgtY6M5xxSAR15Zc3rmzBlZrVa3x6xWq86fb1581LdvXyUmJio/P991vaGhQfv373etMU1LS1Ntba2KiopcNUVFRaqrq3NbhwoAaF82w6E791TrxWPmgmmYpA+nJxBMAZjilZHTyZMna8OGDerbt68GDRqk4uJibdq0SXfffbek5rWj8+bN07p165SamqoBAwZo7dq1ioyM1PTp0yVJAwcO1IQJE7Rw4UJt2LBBkrRw4UJNmjSJoXMA8CMrDhg6bpg7xHREd6v+a1yc+kZ19HGrAAQLr4TT1atX6ze/+Y1++tOfqqqqSomJibr//vu1ePFiV83DDz+s+vp6ZWdny263a9SoUdq2bZuior7+TXrz5s1avHixa3lARkaGVq9e7Y0mAgC8wGY49Mrxeo91Vot0YFoCoRRAq1nsdruzvRvhj1jsDIl+gGb0g2Y2w6FReRU657k0KNeX0g9wAX3Bt7yy5hQAENxshkM3vmQumA7sqqALpgDajlem9QEAwefCjvxjXzl0qPqcGk3MsyWGSy9M5oQVAJePcAoA+JYLO/LNbnySpPSEMD05tjvrTAFcEcIpAOBbWrMjX5Jm9gvTE7cm+rBFAEIFa04BAG5shkNv/K3BdH2MRQRTAF5DOAUAuNgMhzJfrVJlg7mDXCIkvTWNNaYAvIdwCgBwWVpYo5N1503Vjuhu1bvTOcsUgHex5hQA4PLG30zcj1RS9rDOWnZdnI9bAyAUEU4BAJKke18r0xkTe6BWXh+pB78T6/sGAQhJhFMAgG54oVRH6jzXZQ/rTDAF4FOEUwAIYTbDoTF5Fao1UctUPoC2wIYoAAhRNsOhm18xF0wtEsEUQJtg5BQAQpDNcGjMSxWqNXnO/uax0b5tEAD8HSOnABBiCsrqNSLPXDDtIGnL2GhN6x/l83YBgMTIKQCEFJvh0A9eOy0zR+xHSCqb1cvXTQIAN4ycAkCIsBkOjd5WobPmztjXu9O58xOAtkc4BYAQYDMcGp1XoXoTwXREd6s+5M5PANoJ0/oAEKRshkMrDhgqO9Okz2vOqd7Ec+Zc00lrxsT7vG0AcCmEUwAIQjbDoTv3VOu4YXI7vqSrIkQwBdDumNYHgCC04oDRqmAqSS9nssYUQPsjnAJAEDpunGtV/Zax0awxBeAXCKcAEGRshkMHqxymai3iHFMA/oU1pwAQRGyGQ//0UoXOmTjI9PY+4cpJj2HEFIBfIZwCQJCwGQ6N214pM0tNV14fqQe/E+v7RgFAKxFOASAI2AyHJv65QtWNnms7SwRTAH6LNacAEASy3jqt8gZztVsnd/dtYwDgChBOASDAZRdUqrDC8+786I7SzsndNSY5og1aBQCXh2l9AAhQNsOhcS9VqNrEGtOBXaXCGb183ygAuEKEUwAIQDbDofS8CpmZyQ+T9MJkDtgHEBiY1geAAPSDPVWmgmm4RXplcneOiwIQMBg5BYAAYjMcytxRoZNnPdemJ4TpybEEUwCBhXAKAAHAZji05N0a7Tl5VudN1F8VIe3JTPR5uwDA2winAODnbIZDd+6p1nEzp+v/3cuZrDEFEJgIpwDgpy6Mlu7921mdNTNc+nc7WWMKIIARTgHAD9kMhzJ3VerkGafp51gt0vZJnGMKILARTgHAD604YLQqmKZEWrQjI54RUwABj3AKAH7oqN1hunbL2GhN6x/lw9YAQNvhnFMA8EMff+n5dqSSdHNiR4IpgKDCyCkAtDOb4dCKA4bKzjQpuYtVk3t31DkTM/rRYdJvb+rm+wYCQBsinAJAO7rYMVEvHqv3+DwO2AcQrAinANCOVhwwWnV+qSStvD5SD34n1kctAoD2xZpTAGhHZWdaF0wze3UgmAIIal4Lp6dOndKDDz6o/v37KzExUenp6XrnnXdc151Op3JycjRo0CAlJSUpMzNTR44ccXsNu92urKwspaSkKCUlRVlZWbLb7d5qIgD4neQu1lbVPzcx2UctAQD/4JVwarfbNWnSJDmdTr3wwgsqLCzU6tWrFR8f76rZuHGjNm3apFWrVmnv3r2Kj4/X1KlTZRiGq2b27NkqLi5WXl6e8vLyVFxcrLlz53qjiQDgl5aPjFLfrua+ireMjfZxawCg/Xllzelvf/tbJSUl6YknnnA91q9fP9f/djqdys3N1YIFCzRlyhRJUm5urlJTU5WXl6dZs2bp6NGjev3117V7926lpaVJktavX6+MjAyVlJQoNTXVG00FAL9hMxx69P2vVNXQ8r1JI6zSY2M4yxRAaLDY7XbztyC5hPT0dI0fP15lZWV6++23lZSUpPvuu09z5syRxWLRF198oREjRmjv3r0aOXKk63l33XWXunfvrscff1zPPvusli5dqhMnTshisUhqDrW9e/fWqlWr9M///M8X/btLSkqutPkA0OZK6y3KKg5XhaOlUVOnfj3grCYntRxeASCQeBpw9MrI6RdffKEtW7Zo/vz5WrBggQ4dOqRHHnlEkpSVlaXy8nJJcpvmv/DnsrIySVJFRYXi4uJcwVSSLBaLevTooYqKikv+3b4aUWW0FhL9AM180Q9+9lqlKhyN33q8g6SuYVJsZ6ty/ylGY5IjvPr34vLxfYAL6Au+5ZVwev78eV177bV69NFHJUnf/e53dezYMW3evFlZWVne+CsAIGjUNJ5X/t++HUwlqVu49PkPe7VxiwDAf3hlQ1RiYqIGDhzo9tg111yjkydPuq5LUmVlpVtNZWWlEhISJEkJCQmqrq6W0/n1KgOn06mqqipXDQAEuprG85r2WtUl7wDVwcIJfwBCm1e+BUePHq3PPvvM7bHPPvtMffr0kST17dtXiYmJys/Pd11vaGjQ/v37lZ6eLklKS0tTbW2tioqKXDVFRUWqq6tz1QBAILsQTN+vdFyy5roe3PEJQGjzyrT+/PnzNXHiRK1du1Y/+MEPVFxcrCeffFI///nPJTWvHZ03b57WrVun1NRUDRgwQGvXrlVkZKSmT58uSRo4cKAmTJighQsXasOGDZKkhQsXatKkSazrABDQLuzK/3+lZ2U4Lr0HtXcXi1aOjmnDlgGA//FKOB05cqSee+45/epXv9KaNWvUu3dv/exnP9Ps2bNdNQ8//LDq6+uVnZ0tu92uUaNGadu2bYqK+vpolM2bN2vx4sWaNm2aJCkjI0OrV6/2RhMBoF3YDIe+v7tKttpv77gfEddR/aKsqm44r6Qu1uYzT6MYOQUQ2rxylFQwYiceJPoBml1uP7AZDmW+WqWTdd8OpnHhHfTB9ERFd2KNaaDg+wAX0Bd8yysjpwAAdzbDoYxdlfrbmYv//n9NjJVgCgAXQTgFAC97/CO7lrxX12JN7658/QLAxfDtCABeYjMcunNnhY7Xt1x3VVTz+lIAwLcRTgHAC2yGQze/UiH7pU+JkiSFd5BenhTHxicAuAQWPAGAF0zc4TmYStK4nuEEUwBoASOnAHAFbIZD382rMFXLOaYA4BkjpwBwmVoTTGM7Sjtvj2fUFAA8IJwCwGUaYTKYdrFIb05JIJgCgAmEUwC4DLHPlMrMHUxGdLdq/zSCKQCYxZpTAGiFgrJ6Ze4+bbr+jSlJPmwNAAQfRk4BwKTWBtPsYZ192BoACE6EUwAwaUorgunMfmFadl2cD1sDAMGJcAoAJkzdWapzJmuzh3XWE7cm+rQ9ABCsWHMKAB78c1EHHW00Vzvnmk6MmALAFSCcAkALej5TqjMKN1U7s1+Y1oyJ93GLACC4Ma0PABdRUFav2GdKdUaSZPFYP7NfGFP5AOAFhFMA+Aet3ZWf2asDwRQAvIRwCgDsGuP7AAAchElEQVT/4PutPC7quYnJPmwNAIQW1pwCwDfc8EKpmkzWrrw+Ug9+J9an7QGAUMPIKQD83S1bS3WkzlztzH5hBFMA8AFGTgFAzSOmZoPpnGs6sSsfAHyEcAog5A16tlSnTJ6wv2VstKb1j/JtgwAghBFOAYQsm+HQd/MqTNevvD6SYAoAPsaaUwAhqXXB9Lyyh3VmjSkAtAFGTgGEpBEmg2mYpLxRZ3XL8D6+bRAAQBIjpwBCUL9nSuU0WfvX6QnqFWG2GgBwpQinAEKCzXBozpunlfhMqewmn7Nzcnf1jero03YBANwxrQ8g6NkMhzJfrdLJuvOmn/Ph9ASCKQC0A0ZOAQS9pYU1rQqmW8ZGE0wBoJ0QTgEENZvh0K4TZ03Xc1wUALQvwimAoFVQVt/qc0w5LgoA2hfhFEBQ2vq5oczdp03XZ/bqQDAFAD9AOAUQdArK6vXAW1+Zrp/ZL0zPTUz2YYsAAGaxWx9AULEZjlaNmO6c3F1jkiN82CIAQGswcgogaDz+kb1Va0wJpgDgfxg5BRAUfvN+tdYcajBdTzAFAP/EyCmAgLf1c6NVwXTL2GiCKQD4KcIpgIC29XOjVZufOMcUAPwb4RRAwGrtrnzOMQUA/0c4BRCwWrMrf841nQimABAA2BAFIODYDEerduVvGRvNVD4ABAifjJyuW7dOsbGxys7Odj3mdDqVk5OjQYMGKSkpSZmZmTpy5Ijb8+x2u7KyspSSkqKUlBRlZWXJbrf7ookAApDNcGjSn8sJpgAQxLweTt977z394Q9/0NChQ90e37hxozZt2qRVq1Zp7969io+P19SpU2UYhqtm9uzZKi4uVl5envLy8lRcXKy5c+d6u4kAApDNcGjCnytVWHnO9HNGdxfBFAACjFfDaU1NjebMmaPHHntMsbFfr+1yOp3Kzc3VggULNGXKFA0ZMkS5ubmqra1VXl6eJOno0aN6/fXXtWHDBqWlpSktLU3r16/Xnj17VFJS4s1mAggwNsOhyTurVNngNP2czF4dtHtKLx+2CgDgC15dc3ohfI4dO1arVq1yPW6z2VReXq5x48a5HouIiNCNN96owsJCzZo1S0VFReratavS09NdNaNHj1ZkZKQKCwuVmpp60b/Tl8GVUAyJftDeSusteuhwuMoazP4u7dSPks9q3lXnvfrZ0Q8g0Q/wNfrC5btUprvAa+H0v/7rv3Ts2DE9+eST37pWXl4uSYqPj3d7PD4+XmVlZZKkiooKxcXFyWKxuK5bLBb16NFDFRWXXl/m6Q1erpKSEp+9NgIH/aD9rX7ztE421Juun3NNuNaM6e3VNtAPINEP8DX6gm95JZyWlJToV7/6lXbv3q2OHTt64yUBQL95v1ovHuOWpAAQSrwSTouKilRdXa3Ro0e7HmtqatK+ffv09NNP691335UkVVZWqk+fPq6ayspKJSQkSJISEhJUXV0tp9PpGj11Op2qqqpy1QAIHXPzy/X8F+Y3P2UP60wwBYAg4JUNUZmZmdq3b5/efvtt13/XXnutpk2bprffflsDBgxQYmKi8vPzXc9paGjQ/v37XWtM09LSVFtbq6KiIldNUVGR6urq3NahAgh+Wz83WhVMR3eXll0X58MWAQDaildGTmNjY91250tSly5d1K1bNw0ZMkSSNG/ePK1bt06pqakaMGCA1q5dq8jISE2fPl2SNHDgQE2YMEELFy7Uhg0bJEkLFy7UpEmTWNcBBDmb4dCKA4bKzjQpKsyiV0+eNf3czF4d9NzEZB+2DgDQltrsDlEPP/yw6uvrlZ2dLbvdrlGjRmnbtm2Kivr6DMLNmzdr8eLFmjZtmiQpIyNDq1evbqsmAmgHNsOhO/dU67jR1Orn3poggikABBmL3W43f3BgCGEnHiT6QVuY8+ZpvXjM/G78C0Z3V5udY0o/gEQ/wNfoC77VZiOnAHAxxw3za0svWHl9pB78TqznQgBAwCGcAmhXFfXnW1X/4fQE9Y3iyDoACFZevX0pAJhlMxya8+ZpldeZX2u6ZWw0wRQAghwjpwDaXGs3QVkkbR4brWn9ozzWAgACG+EUQJtbccAwHUzTE8L05NjujJgCQIggnAJoc++Wm7slaVqcVXsyE33cGgCAP2HNKYA2lV1QqRN1nk+w693Foqdu5a5PABBqGDkF0CZshkMzX6vQJ195rr0luZM2jollKh8AQhDhFIDP2QyHbt1RqdMm7ko6PMailyfH+75RAAC/RDgF4FM2w6GbX6mQ3eG5Nraj9OxtBFMACGWsOQXgMzbDoYk7zAXTwTEd9OYUDtgHgFDHyCkAn7n7L1UqNzGVP+eaTlozhhFTAADhFICP/Ob9ah2p8Xxr0uxhnbXsOnblAwCaEU4BeJXNcGjFAUMvHvN8lunK6yP14Hdi26BVAIBAQTgF4DWtuS3p8BgLwRQA8C1siALgNWZvS9pJ7MoHAFwcI6cArpjNcGhpYY12n/C8+8ki6aXJ3dmVDwC4KMIpgCtiMxyatLNKp+o9b34Kk/TK5O4akxzh+4YBAAIS4RTAFbn7L56DadcwKSMlQstHRjFiCgBoEeEUwGUze1xURkqEnrq5exu0CAAQ6NgQBeCyPP6RXWsOeT4u6qooq5aPjGqDFgEAggEjpwBaxWY4dN/rlfrQ7vRYe3ufcOWkxzCVDwAwjXAKwLSvNz95DqY3J3bUHyf0aINWAQCCCdP6AEx7uMBuale+JP32pm4+bg0AIBgRTgF4ZDMcuucvVXqjrNFU/Zax0UzlAwAuC9P6AFpkMxzKfLVKJ+vMjZiuvD5S0/qzAQoAcHkIpwBa9HCB3VQw7dbJov8e140D9gEAV4RwCuCibIZDD7/zpd445WixziLpz9z1CQDgJYRTAN9idio/MsyiFyYwWgoA8B7CKYBvWVpY4zGYWiTtuzOejU8AAK9itz4ANwVl9dp14qzHuqQIC8EUAOB1jJwCkNQ8lb/k3RrtPuk5mFokbb6Zc0wBAN5HOAWggrJ63fX6adWd81yb0LmDnrkllnWmAACfIJwCIc5mODTtL6fV0NRyXUdJ709PYCofAOBTrDkFQtz9e6s9BlNJmtA7nGAKAPA5wikQwh7/yK6Dpz0n08TO0srRMW3QIgBAqGNaHwgyNsOhFQcMlZ1pUnIXq5aPjLroiOfWzw0tea+uxdfqHm7R6IROykmPYdQUANAmCKdAELEZDt25p1rHja9HQ9+vbNTLk+LcwmVBWb0eeOurFl8rPlwq+WFPn7UVAICLYVofCCIrDhhuwVSSjhtNWnHAcP25oKxed+w+7fG1/nBrd6+3DwAATxg5BYJI2ZmLrx89dabJ7RxTp4fXGR5j4agoAEC7IJwCQSS5i/Wij0d1tCjz1SqPtySVpNiO0rO3xXu7aQAAmOKVaf1169bp1ltvVZ8+fdS/f3/NnDlTH3/8sVuN0+lUTk6OBg0apKSkJGVmZurIkSNuNXa7XVlZWUpJSVFKSoqysrJkt9u90UQgJCwfGaWrotwD6lVRVh2vaTQVTNMTwvTmFM4yBQC0H6+E03feeUcPPPCA9uzZo+3btyssLEx33nmnvvzyS1fNxo0btWnTJq1atUp79+5VfHy8pk6dKsP4ei3c7NmzVVxcrLy8POXl5am4uFhz5871RhOBkNA3qqMeGxOjlK5WxXS0KKWrVROSrTrylaeJfOnmxI7ak5lIMAUAtCuvTOtv27bN7c9PPPGEUlJS9O677yojI0NOp1O5ublasGCBpkyZIknKzc1Vamqq8vLyNGvWLB09elSvv/66du/erbS0NEnS+vXrlZGRoZKSEqWmpnqjqUBQsxkOPVRQo/+tbV57WuNo0lOfmjvH9Lc3dfN18wAA8Mgna05ra2t1/vx5xcbGSpJsNpvKy8s1btw4V01ERIRuvPFGFRYWatasWSoqKlLXrl2Vnp7uqhk9erQiIyNVWFh4yXBaUlLii7fg89dG4AikfvDzox113DA/8tnR4tQN3Zr0f65yqPHUFyo55cPGBbhA6gfwHfoBLqAvXD5PA44+CadLlizRsGHDXCOg5eXlkqT4ePdNFvHx8SorK5MkVVRUKC4uThaLxXXdYrGoR48eqqiouOTf5asRVUZrIQVeP6j9rFJSo+n696cxjW9GoPUD+Ab9ABfQF3zL6+H0Zz/7md59913t3r1bVuvFdw4D8I3ojhbPRX83PMZCMAUA+B2vHsK/dOlSbd26Vdu3b1e/fv1cjycmJkqSKisr3eorKyuVkJAgSUpISFB1dbWczq83bjidTlVVVblqALQs/8RZU3UxYRwXBQDwT14Lp4888ogrmF5zzTVu1/r27avExETl5+e7HmtoaND+/ftda0zT0tJUW1uroqIiV01RUZHq6urc1qECuLhbtpaq3kTd4JgOeutOjosCAPgnr0zrL1q0SM8//7z++7//W7Gxsa41ppGRkeratassFovmzZundevWKTU1VQMGDNDatWsVGRmp6dOnS5IGDhyoCRMmaOHChdqwYYMkaeHChZo0aRLrOgAP5uaX6+BXnuuyh3XWsuvifN8gAAAuk1fC6ebNmyXJdUzUBY888oiWLl0qSXr44YdVX1+v7Oxs2e12jRo1Stu2bVNUVJTb6yxevFjTpk2TJGVkZGj16tXeaCIQtKbuLFX+pfcMunw4ndFSAID/80o4NXMXJ4vFoqVLl7rC6sXExsbqySef9EaTgKBlMxxaccBQ2ZkmffFlo06aWGaaPawzwRQAEBB8cpQUAN+wGQ5lvlpl6lak38RUPgAgUHh1tz4A31paWNPqYLry+kgftQYAAO8jnAIB5N1y8wfsW9QcTB/8TqzvGgQAgJcxrQ8EiIKyep1udHoulHRrgvRSZi8ftwgAAO9j5BQIAFs/N5S5+7Sp2pn9wgimAICAxcgp4Oe2fm7ogbdMHGIqaefk7hqTHOHjFgEA4DuMnAJ+7PGP7KaDaZxVBFMAQMAjnAJ+auvnhpa8V2e6fu/UBB+2BgCAtkE4BfyQzXCYHjGVmnflc8g+ACAYEE4BP1NQVq/v5pm4H+nfZfbqwHFRAICgQTgF/MjjH9lN78qXmm9L+tzEZB+2CACAtsVufcBPPP6R3fQa04FdpcIZHBcFAAg+hFOgndkMh7LePK3CynOmn/PCZDY/AQCCE+EUaEc2w6EJf65UZYO5Oz9JbH4CAAQ31pwC7WhpYU2rgumcazqx+QkAENQYOQXa0a4TZ03XZg/rrGXXxfmwNQAAtD/CKdAObIajVcdFrbw+khFTAEBIYFofaGMEUwAALo2RU6CNXduKYLpzcneNSY7wYWsAAPAvjJwCbcRmOBT7TKnOm6wnmAIAQhHhFGgDrb0lafawzgRTAEBIYlof8DGb4WjVLUlZYwoACGWMnAI+xOYnAABah3AK+MjlTOUTTAEAoY5pfcDLbIZDWW9Uq7CqyfRz2PwEAEAzwingRTbDoYl/rlB5g/nnfDg9QX2jOvquUQAABBCm9QEvWnHAIJgCAHAFCKeAF9gMh+a8eVpbj9Wbfs6WsdEEUwAA/gHT+sAVshkO3bmnWscN82tMV14fqWn9o3zYKgAAAhPhFLhCKw4YrQqmbH4CAODSCKfAFXrZ5FT+4JgO+p/bejCVDwBACwinwBUY9adSOUzW7v9Bsk/bAgBAMCCcApfBZjg0dmuFapzm6reMjfZtgwAACBKEU6CVCsrqlbn7tKnaMIv0xE3RbH4CAMAkwinQCjbDoe+ZDKaSVPWvvXzYGgAAgg/nnAKtkLGzQiZn8pnKBwDgMhBOAZOyCyr1N5Nn7G8Zy1Q+AACXg2l9wITfvF+tpz5t9FjXQdIH3JIUAIDLxsgp4MHWzw2tOdTgsa6ThWAKAMCVYuQUaMFfv+ygBw9/Zar2vWkEUwAArhQjp8Al2AyHHv64k8e6Dmq+JSnBFACAK+eX4XTz5s0aPny4EhMTdfPNN2vfvn3t3SSEoEff/0pnnS3/E4m0Nk/lj0mOaKNWAQAQ3PwunG7btk1LlizRT3/6U7311ltKS0vTjBkzdOLEifZuGkJITeN5/b/Ssy3WRFmlfVOZygcAwJv8Lpxu2rRJP/zhD3X//fdr4MCBWrNmjRITE/X000+3d9MQQqwWKdxqueT1uE4WvUMwBQDA6yx2u93smeI+19jYqOTkZG3ZskV33nmn6/FFixbp448/1q5du771nJKSkrZsIkLIZ7UW3fdhZzmc7iG1R9h5bf7uWfWK8Jt/OgAABIzU1NQWr/vVbv3q6mo1NTUpPj7e7fH4+HhVVFRc9Dme3uDlKikp8dlrIzCkSvq/+kxLP+sqw3FeHSzS9fGdlJMew4hpiOH7ABL9AF+jL/iWX4VTwN8M6OrUwRlJ7d0MAABChl+tOY2Li5PValVlZaXb45WVlUpISGinVgEAAKCt+FU47dSpk0aMGKH8/Hy3x/Pz85Went5OrQIAAEBb8btp/R//+MeaO3euRo0apfT0dD399NM6deqUZs2a1d5NAwAAgI/5XTj9wQ9+oNOnT2vNmjUqLy/X4MGD9cILLyglJaW9mwYAAAAf87twKkmzZ8/W7Nmz27sZAAAAaGN+teYUAAAAoY1wCgAAAL9BOAUAAIDfIJwCAADAbxBOAQAA4DcIpwAAAPAbhFMAAAD4DYvdbne2dyMAAAAAiZFTAAAA+BHCKQAAAPwG4RQAAAB+g3AKAAAAv0E4BQAAgN8gnAIAAMBvhHQ4zczMVGxsrNt/P/rRj9xq7Ha7srKylJKSopSUFGVlZclut7vVHD58WLfffruSkpI0ePBgrVq1Sk4nJ3QFus2bN2v48OFKTEzUzTffrH379rV3k+AlOTk53/q3f80117iuO51O5eTkaNCgQUpKSlJmZqaOHDni9hpmvhvgXwoKCnT33Xdr8ODBio2N1XPPPed23VufOz8T/J+nvjBv3rxvfUdMmDDBrebs2bPKzs7W1VdfrZ49e+ruu+9WaWmpW82JEyc0c+ZM9ezZU1dffbUWL16sxsZGn7+/QBfS4VSS7r33Xh09etT13/r1692uz549W8XFxcrLy1NeXp6Ki4s1d+5c1/WvvvpKU6dOVUJCgvbu3auVK1fqd7/7nR577LG2fivwom3btmnJkiX66U9/qrfeektpaWmaMWOGTpw40d5Ng5ekpqa6/dv/5i8fGzdu1KZNm7Rq1Srt3btX8fHxmjp1qgzDcNV4+m6A/6mrq9OQIUO0cuVKRUREfOu6Nz53fiYEBk99QZJuueUWt++IF1980e360qVLtWPHDm3ZskW7du2SYRiaOXOmmpqaJElNTU2aOXOmamtrtWvXLm3ZskXbt2/XsmXLfP7+Al1IH8KfmZmpIUOGaM2aNRe9fvToUaWnp2v37t0aPXq0JGn//v3KyMjQe++9p9TUVG3ZskX//u//rk8//dTVwdesWaOnn35aH3/8sSwWS5u9H3jP+PHjNXToUP32t791PTZy5EhNmTJFjz76aDu2DN6Qk5Oj7du3a//+/d+65nQ6NWjQIM2ZM0eLFi2SJNXX1ys1NVW//vWvNWvWLFPfDfBvvXr10urVq3XvvfdK8t7nzs+EwPOPfUFqHjk9ffq0nn/++Ys+p6amRgMGDNCmTZt01113SZJOnjypYcOGKS8vT+PHj9df/vIX3XXXXTp06JB69+4tSXr++ef1k5/8RCUlJYqOjvb9mwtQIT9yunXrVl199dUaPXq0li9f7vYbclFRkbp27ar09HTXY6NHj1ZkZKQKCwtdNTfccIPbb17jx49XWVmZbDZb270ReE1jY6MOHjyocePGuT0+btw41+eOwPfFF19o0KBBGj58uH70ox/piy++kCTZbDaVl5e7ff4RERG68cYb3f7de/puQGDx1ufOz4TgsX//fg0YMECjRo3ST37yE1VWVrquHTx4UA6Hw62/9O7dWwMHDnTrCwMHDnQFU6m5L5w9e1YHDx5suzcSgMLauwHtacaMGerTp4+SkpL0ySef6Je//KUOHz6sl156SZJUUVGhuLg4t990LRaLevTooYqKCldNz5493V43Pj7eda1fv35t82bgNdXV1WpqanJ9jhfEx8e7PncEtuuuu06///3vlZqaqqqqKq1Zs0YTJ07Uu+++q/Lyckm66OdfVlYmydx3AwKLtz53fiYEhwkTJuiOO+5Q37599b//+79asWKFvv/97+uNN95QeHi4KioqZLVaFRcX5/a8b/6cqKio+FZ/iouLk9Vq5XvCg6ALpytWrNDatWtbrNmxY4duuukm/eu//qvrsaFDh6pfv34aP368Dh48qBEjRvi4pQDay2233eb25+uuu04jRozQH//4R11//fXt1CoA/mLatGmu/z106FCNGDFCw4YN0549e/T973+/HVsWGoIunM6bN8+1/uNSvjnE/k3XXnutrFarjh07phEjRighIUHV1dVyOp2u35SdTqeqqqqUkJAgSUpISHAb6pfk+vOFGgSWC7/ZXuxz5TMNTl27dtWgQYN07Ngxfe9735PU/Hn36dPHVfPNz9/MdwMCS2JioqQr/9z5mRCckpOT1bNnTx07dkxS82fZ1NSk6upq9ejRw1VXWVmpG264wVXzj8t8LszM0RdaFnRrTuPi4nTNNde0+F+XLl0u+tzDhw+rqanJ9SWVlpam2tpaFRUVuWqKiopUV1fnWnOUlpam/fv3q6GhwVWTn5+v5ORk9e3b14fvFL7SqVMnjRgxQvn5+W6P5+fnu601Q/BoaGhQSUmJEhMT1bdvXyUmJrp9/g0NDdq/f7/bv3tP3w0ILN763PmZEJyqq6tVVlbmygcjRoxQx44d3fpLaWmpa9Oc1NwXjh496na8VH5+vsLDw5md9cC6ZMmSf2/vRrSH48eP68knn1RkZKQaGxtVVFSkBQsWqFevXlq+fLk6dOigHj166P3331deXp6GDRum0tJSLVy4UCNHjnQdHdK/f38988wzOnTokFJTU7V//3794he/0IIFC/ghFcCioqKUk5OjpKQkde7cWWvWrNG+ffv02GOPKSYmpr2bhyu0fPlyderUSefPn9dnn32m7OxsHTt2TOvXr1dsbKyampq0YcMG9e/fX01NTVq2bJnKy8u1YcMGhYeHm/pugP+pra3VJ598ovLycj377LMaMmSIoqOj1djYqJiYGK987vxMCAwt9QWr1apf/epX6tq1q86dO6dDhw7p3/7t39TU1KQ1a9YoPDxcnTt31qlTp7R582YNHTpUNTU1WrhwoaKjo/XLX/5SHTp0UL9+/bRjxw7t3btXQ4cO1SeffKJFixZpxowZuuOOO9r7/wK/FrJHSZ08eVJZWVk6cuSI6urq1KtXL02cOFFLlixRt27dXHV2u12LFy/Wq6++KknKyMjQ6tWrFRsb66o5fPiwFi1apAMHDig2NlazZs3SI488wpEhAW7z5s3auHGjysvLNXjwYP3Hf/yHxowZ097Nghf86Ec/0r59+1xTctddd52WLVumQYMGSWqeql25cqX+8Ic/yG63a9SoUVq7dq2GDBnieg0z3w3wL2+//fZFQ8E999yj3Nxcr33u/Ezwfy31hXXr1unee+9VcXGxampqlJiYqJtuuknLli1zWxZ49uxZLV++XHl5eWpoaNDYsWP1n//5n241J06c0KJFi/TWW2+pc+fOmjFjhn79618rPDy8Td5noArZcAoAAAD/E3RrTgEAABC4CKcAAADwG4RTAAAA+A3CKQAAAPwG4RQAAAB+g3AKAAAAv0E4BQAAgN8gnAIAAMBv/H+j1QWgdnJ9+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for w in wins[0:1]:\n",
    "    (f, t) = presence_times(w)\n",
    "    zz, dgms, cells = d.zigzag_homology_persistence(f, t)\n",
    "    d.plot.plot_diagram(dgms[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.plot.plot_diagram(dgms[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pool = multiprocessing.Pool(processes=multiprocessing.cpu_count())\n",
    "dgms = pool.map(zigzag_network, wins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pool.terminate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sliced Wasserstein kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgms0 = [dgm[0] for dgm in dgms if dgm != []]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgms1 = [dgm[1] for dgm in dgms if len(dgm) > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sliced_wasserstein import diagram_array, SW_approx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[      0.        ,  168498.26575157,  142213.96446884, ...,\n",
       "         181004.87138341,  980779.00902599,  604002.56924413],\n",
       "       [ 168498.26575157,       0.        ,  108974.83210031, ...,\n",
       "         292878.72367602,  864985.91166401,  466984.82235283],\n",
       "       [ 142213.96446884,  108974.83210031,       0.        , ...,\n",
       "         223518.93608983,  951003.34397253,  553228.56514372],\n",
       "       ...,\n",
       "       [ 181004.87138341,  292878.72367602,  223518.93608983, ...,\n",
       "              0.        , 1153158.29742231,  755045.92126984],\n",
       "       [ 980779.00902599,  864985.91166401,  951003.34397253, ...,\n",
       "        1153158.29742231,       0.        ,  607164.48834291],\n",
       "       [ 604002.56924413,  466984.82235283,  553228.56514372, ...,\n",
       "         755045.92126984,  607164.48834291,       0.        ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gram0 = np.array([[SW_approx(dgms0[i], dgms0[j], 10) for i in range(len(dgms0))] for j in range(len(dgms0))])\n",
    "gram0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0.        , 1782.59876486, 1976.144841  , ..., 1674.27089231,\n",
       "        4008.08473571, 3074.53612959],\n",
       "       [1782.59876486,    0.        ,  814.44479983, ..., 1389.37805531,\n",
       "        3609.13022428, 1694.74694713],\n",
       "       [1976.144841  ,  814.44479983,    0.        , ..., 1665.50397679,\n",
       "        4034.59792213, 2097.51910373],\n",
       "       ...,\n",
       "       [1674.27089231, 1389.37805531, 1665.50397679, ...,    0.        ,\n",
       "        2953.12379337, 1952.36875879],\n",
       "       [4008.08473571, 3609.13022428, 4034.59792213, ..., 2953.12379337,\n",
       "           0.        , 3326.13065883],\n",
       "       [3074.53612959, 1694.74694713, 2097.51910373, ..., 1952.36875879,\n",
       "        3326.13065883,    0.        ]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gram1 = np.array([[SW_approx(dgms1[i], dgms1[j], 10) for i in range(len(dgms1))] for j in range(len(dgms1))])\n",
    "gram1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dgm2array(dgm):\n",
    "    return np.array([[p.birth, p.death] for p in dgm])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn_tda import SlicedWasserstein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bottleneck distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client = Client()\n",
    "distmat = []\n",
    "for dgm in dgms0:\n",
    "    distmat.append(client.map(lambda x: d.bottleneck_distance(x, dgm), dgms0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        , 37.67822266, 71.04492188, ...,         inf,\n",
       "                inf, 19.10028076],\n",
       "       [37.69311523,  0.        , 70.44287109, ...,         inf,\n",
       "                inf, 37.14916992],\n",
       "       [70.62719727, 70.30273438,  0.        , ...,         inf,\n",
       "                inf, 69.04138184],\n",
       "       ...,\n",
       "       [        inf,         inf,         inf, ...,  0.        ,\n",
       "                inf,         inf],\n",
       "       [        inf,         inf,         inf, ...,         inf,\n",
       "         0.        ,         inf],\n",
       "       [19.10028076, 37.14916992, 69.04138184, ...,         inf,\n",
       "                inf,  0.        ]])"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distmat = np.array(client.gather(distmat))\n",
    "distmat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lozeve/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/lozeve/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/lozeve/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/lozeve/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.svm import OneClassSVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = OneClassSVM(kernel='precomputed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0,\n",
       "       1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0,\n",
       "       1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,\n",
       "       0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(clf.predict(gram)+1)//2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import AgglomerativeClustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lozeve/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/site-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "gram1 = 1/gram1\n",
    "gram1[gram1 == np.inf] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = AgglomerativeClustering(n_clusters=2, affinity=\"precomputed\", linkage=\"complete\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AgglomerativeClustering(affinity='precomputed', compute_full_tree='auto',\n",
       "            connectivity=None, linkage='complete', memory=None,\n",
       "            n_clusters=2, pooling_func=<function mean at 0x7f3dd41c29d8>)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(gram1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0,\n",
       "       1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0,\n",
       "       0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(gram1)\n",
    "clf.labels_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "distmat[distmat==np.inf] = 1e100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1,\n",
       "       1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0,\n",
       "       1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0,\n",
       "       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1])"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(distmat)\n",
    "clf.labels_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}