sociopatterns_daily.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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "datadir = \"data/sociopatterns/infectious/INFECTIOUS_cumulative_daily_networks/\"\n",
    "files = [os.path.join(datadir, f) for f in os.listdir(datadir)]\n",
    "gs = []\n",
    "for f in files:\n",
    "    gs.append(ig.read(f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH U-W- 242 1092 -- \\n+ attr: id (v), weight (e)'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gs[0].summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(gs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Persistent homology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from filtration import wrcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_diagram(graph, weight=\"weight\"):\n",
    "    filt = wrcf(graph, weight=weight)\n",
    "    pers = d.homology_persistence(filt)\n",
    "    return(d.init_diagrams(pers, filt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 36s, sys: 128 ms, total: 1min 36s\n",
      "Wall time: 1min 36s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "dgms = [get_diagram(g) for g in gs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "with Client() as client:\n",
    "    dgms = client.map(get_diagram, gs)\n",
    "    dgms = client.gather(dgms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApYAAAF9CAYAAAC+mSjWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X90VPWd//HXkIQwhCkxMSS0QORHUoQvVIMrHPDHFo9Sv7QbUal6usfFalAsPdJdQMKmaml2EWXZYqXxB+2hW+1pBTnitit8dw8U+RlF6tJdUWMjsWgSQ2KSCZmQZDLfPzBjJjOTuZO5M3Nn5vk4x3PKZ+7M/dz7IeTVufd937bW1laPAAAAgAiNiPcEAAAAkBwIlgAAADAFwRIAAACmIFgCAADAFARLAAAAmIJgCQAAAFMQLAEAAGAKgmUKq6mpifcUMAhrYj2sifWwJtbDmlhPvNaEYAkAAABTECwBAABgCoIlAAAATEGwBAAAgCkIlgAAADAFwRIAAACmIFgCAADAFARLAAAAmIJgCQAAAFMQLAEAAGCK9HhPAACSTZ2zR5UnnarvdGv86DRVlDhU6MiI97SGJZJjieS9L//Zqe8fdepCr0eZ6Tb9dL5DBaPTteJwm1pcduW83aCqa8ZqwXi7qcf2Wt15/eOJ83J7pDSb9E9XZemB/5Nt6L2SDH/eM//T6jcuyW/s0Ccu/f7jPu9+F39lhM+f+7Xe85VhnwfATLbW1lZPvCeB+KipqVFRUVG8p4EBWBPrCXdN6pw9umVfsz50ur1jkx1pemVRbsKFy0iOJZL3vvxnp+59vd1vfISkgZEq3SbtWZQzrHAZaH5j06Q2t/+2j/+Vb7gM9N4Jo22Szaaz57+YYbDPCxYOIxXPcMm/XdYTrzXhUjgAmKjypNMncEjSh063Kk864zSj4YvkWCJ57/ePBt5mcBTr9UgrDreF/Dyj8wsUAqWL3yKGeu/ZTo9PqBzq86IRKiXpswvR+VwgHARLADBRfWfgNNEQZNzKIjmWSN57odf4hbS2YYapYPMLxD1oOuG8N5Zu2XeOcIm4I1gCgInGj04LOF4QZNzKIjmWSN6bmW4LuU2/sZnD+zUWbH6BpA2aTjjvjaX/bu7Rbf/vnHr7uMMN8UOwBCygztmjsoMteuBUpsoOtqjO2TPkdt98rWnI7RA/FSUOTXb4Bo/Jji+KOxJJJMdSUeLQhCzfXzETskYYeu9P5wfeZvAvrDSbVHXN2JCfJ0lH6l2avbNBk174RLN3NugbEzL85jcmyG/E/sKafoHOy4TRNuWO9H1fsDs/F38lOr96bZKWXz5G6SOMB3PAbARLIM76CwF21rr0Vnuadta6dMu+Zr/QOHC7ww3dQbdDfBU6MvTKolwtnWLXtQUjtXSKPSELdyQTjsXjGfrPQVw1bpRyB+1ibJqUk+k7dmmmNGFM6IebHKl3qXRfiz7qcKu9x6OPOtwqe71dXb2+l42zR9m0ZtYo7zeUaTb/wh0p8Hn50VUOtQ76Uey2SWXFI/0+78FZ2X6/fEfIf9s7LjP+4BabpKprL9Gd00Ybfg8QDTxuCIizoYocnr8+J+ztEH+FjoykWZPhHkvlSafOdvoGybOdHkN/XytPOtU8KKS1uSUNurWxsUuGPm/F4TYNvm2zT9K5C75jZzs9OnPepuZloaurB5+X2Tsb/O7FdHukfZ+4/T5v9s4Gv0KkPvlv++VffRJyHhKhEtbCN5ZAnBktckimohAkv2gU/gz381rDKGgZ7s9TsH0EKi4yuq2RIiZCJayGYAnEmdEih2QqCkHyi0bhz3A/LzuMAp/h/jwF20eg4iKj24YqYiJUwooIlkCcGS2QSKaikERF8ZRxkRb+5A+6nzI3Q37FNsE+b3ChzoOXj9LgjDZC0qDpKTdDWlZsH9YaV10z1m8f6Tbp0ZIsv88Ltu3gQqRgRUwSoRLWlbZu3brH4j0JxEdLS4tyc3PjPY2Ul52ZppsnZqq5q0/2vgu65itZ+tm12X4FEgO3y80cobnjRgbcDubq/znpL5460titjzrceuezXu37S5dunpip7Ey+NR4skr+v/9vSrZ+/79LAC8HdHmnLPIcy00YM+XPSX6jzWbdHF/qktm6P/lDfo8qrsvR+e5/kkfJHp6nsq5k60NDr815Xn/SfZ7v0ZlNP2Gs8yZGhBfkZOtzQ7d3H5nkOVf7xvN/fmeUzxmjxpEyfbV9YeIlfB6EZOZlqu9CrE+d852nFUMnvE+uJ15rQ0jGF0YLLelgT6+lfk7KDLdpZ6/J7fekUe9IU6ljF7J0N+qjD/17HSWPSdGppwZA/J6He2+/Lv/pEnQYfxD7cNY7078yeMy599w8tPkVBVgyVEv92WREtHQHAwiieip1wCmGG+95wuvsMd40j+TuTSKESGIhgCQAGUDwVO+EUwgz3veF09xnuGg/37wyhEomM51gCgAEVJQ6daOr2eZZosOKROmePKk86Vd/p1vjRF7eJxb2w0dhvLI5l8D4eLcnS/YfafZ49mW6THrx8lGbvbFCLy66ctxtUdc1YTRiT7vfestfbfZ4TOUKBC2Pufb3dby55o2xq6vpix5Mdad6CnlDnYPBxLCu2G/47049QiURHsAQAA/q7rVSedKqh062CIAGjv8hnYJg40dQd9e470dhvLI4l2D6evfZL+tHJ82q70KexmSP04OWjVHHi/Odh06aODre+tbdF4+wj1OD6Ikbu/4sr4MPHGzp9C2Bum3ox3H3/qFMXej3KTLfpp/MdumrcKJ81XlZs18ojbSHPQbDjeHrBWO143zXk35l+hEokA4IlABhkpAtNvDokRWO/sTiWYPvYe7bHp9hm9s6GgN1zBoZKSX4de/p9/6jTGyb73TbV4TcmyefYyg62RNQZa8f7rqQr1AGGwj2WAGCieBX5RGO/sTgWo/sIp3tOIOEU6wwUi85YhEokE4IlAJgoXkU+0dhvLI7F6D7C6Z4TSDjFOgNFuzMWoRLJJuRP6qxZs5Sdne3337e//W3vNtu3b9fs2bOVn5+v66+/XkePHo3qpAHAquLVISka+43FsRjdR6BuNSMkFdh9f43lBrn1c6guNmbMbzjnilCJZBTyHssDBw7I7f7iq/yGhgb99V//tW655RZJ0u7du7Vu3Tr9y7/8i+bNm6ft27dr6dKlOn78uCZOnBi9mQOABRkt8kmE/cbiWIzuY8F4u5699kv6/lGnunr7NCp9hH4636GC0elacbjNW+RTdc1YNXT2+hXlBLqX8ki9SysOt6n1Qp+yP3/v4O43hY4MVVyZ5fN5FVdm+c0v3HNFqESyCrvzzubNm/XUU0/pvffek91u1w033KCZM2fqqaee8m5TUlKi0tJSPfroo6ZPGOahU4L1sCbWw5pYQ6Cq6wmjbZLNprPnv7j/crIjzVDVen/rx8GPNNqzKMcnXBrdLhzJGCr5ObGehOi84/F49Ktf/Up33HGH7Ha7uru79fbbb2vhwoU+2y1cuFDV1dWmThQAkLoCVV2f7fT4hErpi4rtUFYcbvOrMu/1XBwfznZGJWOoBAYK63FDBw4cUF1dne6++25JUnNzs9xut/Ly8ny2y8vL06effjrkZ9XU1IQ5VUQD62A9rIn1sCbxV3suU5KxoqHa5g7V1DQPuU2Ly66LkW7weK/Pehvdzoj959K0/t2Rcg/4PJs8eqy4W3M8HyvR/5rxc2I90ViTUN+ChhUsf/nLX6qkpESzZs2KaFJS6Ikh+rh0YT2sifWwJtYw5ZMWvdXuMrZt7hgVFU0acpuctxvU0eH/KKAce7rPehvdLpQ9Z1xa/16LBn7SxW8qc5Lim0p+TqwnXmtiOFg2NTXpP/7jP7R582bvWG5urtLS0tTU1OS37bhx48ybJQDEWH97vtpzmZrySUvM2jImgkjaPAZ6r6SQY8uK7Xr9E5cau774rNyRkj1jhN89lkaq1quuGRvw3snBrR+NbidJL//ZGbBoiMvfSCWGg+Wvf/1rZWZm6rbbbvOOjRw5UldccYUOHDjgrRKXLl4y/5u/+RtzZwoAMeJbKJKmt9pdMWnLmAgiafMY6L3HGrr8CnACjb3+sUvnLvh+XmuP9MS8Mdp7tifsqvUF4+3asyjHr6J8cEGO0e1e/rPTp/d4Z69H977erhNNPXruXRehEinDULD0eDz6t3/7N916660aM2aMz2vf+973dP/992vOnDmaO3eufvGLX6ihoUH33HNPVCYMANEWr7aMiSCScxOsAEfyhBxrHBQqJcntkX508rxP68dwLBhv16mloSu7jWz3/aOBC4aqTvtevidUItkZCpaHDh3Sn//8Zz333HN+r916661qaWnRk08+qcbGRl1++eV66aWXNGnS0Pe3AIBVxastYyKI5NwEe28k2iJs9WgWIy0jCZVIBYaC5XXXXafW1tagr99333267777TJsUAMRTvNoyJoJIzk2w90ZibIStHs2SmW5T5xDhklCJVBFWVTgAJKPBHVgeLcnSiaZun8u2sWjLGC1mFtssK7brWOMFn/sfJ2SNMHRuKkocOtbQ9fml7ovyM6WMdN8CnEAPPs/PlM51y+dexWBFNEaP12ghkZFz9dP5Dp97LAciVCKVECwBpLTBnVXae9y6/1C7nr32S9p7tke1zR2akjsmYavCo1Fs09M36Js5TxgN3Gw2Dbx/MiN9hJ69Nls73nf5FOBI8muPeLajVysOt6nF1asce3rAIhqjx2u0kMjoubptqkMnmnq4pxIpj2AJIKUF66zSXxRSU9Mc8pmIVhadYhv5jRn9vMGdcs6e79OO910B3zt4rNCRoVNL7UM+n8/o8RotJDJ6rvaccem5dwmVgDVuTgGAOGkNUvxhlaKQSMWq2CaSzzOzKMroPsw8Np5TCXyBYAkgpWUHKf6wSlFIpGJVbBPJ55lZFGV0H2YdG6ES8JUc/3ICwDBVXTNW6YNaQQcrCklEFSUOTXb4BiOjhUiB3jthtE0Tsnx/dUTyeZMdFwuCyg626JuvNansYIvqnD2qc/b4jRlh9HiDHVv+KPmNBTs2QiXgj3ssAaQ0o51VElWhI0OvLMr1K4QxUogU7L2Sf2HNcD9vWbFdK4+0mVZEY/R4g83l/kOtkgbcBmEb9P86PkeoBAKztba2hlHOh2QSrwb1CI41sR7WJLrKDrZoZ60r9IaSlk6x6/nrc6K2JsHm0r/ffoRKf/ycWE+81oRL4QCAuDG7QCgSRgp/CJXA0AiWAIC4MbtAKBKhCn8IlUBo3GMJAIibihKHX5ejQJ13wul8NNzOO8uK7UE7LhEqAWMIlgCAuDG7QCiSzjsnmrr19IKxfl2A3m7uJVQCBhEsAQBxVejIMNR5x4hIOu986HT7dQHim0ogPNxjCQBIGpF23qFQB4gMwRIAkDQi7bxDoQ4QGS6FAwDiKlCxjSRDBTiDBSoGCtZ551jjBZ8CoQlZIyjUASJEsAQAxE2gIppYdN6RJHk8fn/+r7NdWlvdTqgEholgCQCIm0BFNGc7PZJ8Q1+gApxgghUDDd7vxf347nf18XafPRMqgfBwjyUAIG7i1Xkn2H4JlUBkCJYAgLiJV+edUPslVALDw6VwACnvSL1LKw63qfVCn7IzR6jqmrFaMN4e72klFKPdbgYLp/POsmK7yg62qPZcpqZ80qKKEofOdvQOa+0C7bcfoRIYPoIlgJR2pN6l0n0t6v38Gmh7j1ul+1q0Z1EO4dIgo91uAjHaeWdZsV0rj7R9vo80vdXu0usfu3SuW95Cm3DWrtCRoZUzs7inEjAZwRJASltxuM0bKvv1ei6On1pKsDTCaLebYIx03ik72OK3j8YL/p9ldO32nHFpbTWhEjAb91gCSGmtF/oCjrcFGYc/o91uorGPQEKtHc+pBKKHYAkgpWVnBv5ncGyQcfgz2u0mGvsIZKi1I1QC0cW/nABSWtU1Y5Vu8x1Lt10cx8X7J8sOtuibrzWp7GCL6pw9fttUlDg02eEb/AJ1u4lEoH3kZ0qD42aaTXq0JCvgnAmVQPRxjyWAlLZgvF17FuVoxeE2tV3o01iqwr2MFuWE1e1mmAbuo7a5Q1Nyx2hZsV3fPdDic69ldrr06Fsdfl17Vs7MoqMOEAMESwApb8F4O4U6AYRTlGOk202k+vdRU9OsoqJJKjvY4lfA09wjqcf3HssPnW6qv4EY4VI4ACCgWBTlRCKcgh5CJRAbBEsAQECxKMqJRDgFPf0IlUB0cSkcSHDD7XiSbKx0Hqw0l0jmE6g7jdlFOeHoP47+zjvLiu2Guvb0I1QC0UewBBJYJB1PkomVzoOV5hLpfGJRlGOU73Fc7LxzoqlbTy8Yqx3vu3zm919nu7inEogTLoUDCWyo4opUYqXzYKW5mDGf/oKZf785T89fnxO3/8MS7Dh2vO/ymd/bzb101AHiiGAJJDCrF1fEipXOg5XmIllvPsNl5Dh4TiUQfwRLIIFZvbgiVqx0Hqw0F8l68xmuUMdBqASsgXssgQRmteKKWBlcjBKoiCNe5yGcNYlFkU9FiUPHGrp0tvOLxDVhtM30cxPoWE582qXvH3XqQq9Hmek2/XS+Q7dNHd5+hzqvhErAOgiWQAKzUnFFrAQrRglUxBGP82B0TWJa5GOzyedJjjZb0E2HI9Cx7P+L6+LDyj/X2evRva+3S9KwwmWgzjsVJQ693dxLqAQshGAJJLhYdDyxklBFHFZgZE3C6WoTicqTTr9H75w932fqfgIdS7N/S3FJ0vePOof9reXgzjt8UwlYD/dYAkgoqVSMkij7CacDzoVeT+iNDCBUAtZkKFg2NDTogQce0NSpU5Wfn6+5c+fq8OHD3tc9Ho82btyo6dOnq6CgQIsXL9bp06ejNmkAqStVilESaT/hdMDJTI/8Mvz+c2mESsCiQgbL1tZWLVq0SB6PRy+99JKqq6v1xBNPKC8vz7vN1q1btW3bNm3atEn79+9XXl6elixZIqcztZ6lByD6KkocmuzwDTKTHRcLeMoOtuibrzWp7GCL6pxBrsVaRLDjMLuoJhb7CbSP3CC3iT5y5eiI1mnPGZfWvzuSUAlYVMh7LJ966ikVFBTo2Wef9Y5ddtll3v/t8XhUVVWlVatWqbS0VJJUVVWloqIi7dq1S/fcc4/5swaQsgIVxywrtmvlkTbLdLsxIlaFV7HYT7B9DK4Kf+TK0Xr23a5hr5P38re++NaTUAlYi621tXXIG17mzp2rG264QfX19Tp06JAKCgp09913q6ysTDabTWfOnNEVV1yh/fv3q6SkxPu+b3/728rJydEzzzwT9YPA8NTU1KioqCje08AArMnwlB1s0c5al9/40in2iAtUWBPzRLJO3FNpbfycWE+81iTkN5ZnzpzRz3/+cz344INatWqV/vSnP+nhhx+WJC1fvlyNjY2S5HNpvP/P9fX1QT+3pqYmknnDJKyD9bAm4as9lynJ/z6/2uYO1dQ0R/z5rIk5hrtO+8+lXbz87fNNpUePFXdrjudjsTzWwM+J9URjTUKF1ZDBsq+vT1deeaUeffRRSdLXvvY11dbWavv27Vq+fHnUJobo4/9hWg9rMjxTPmnRW+3+34RNyR2joqJJEX02a2Ke4azTnjMurX+vRQPrzm3yqOraHL6ptBB+TqwnXmsSsngnPz9fX/3qV33GiouLdfbsWe/rktTU1OSzTVNTk8aNG2fWPAEgqFgVwiSTOmePoSIao9sZEWydvjEhQ7N3NmjSC59o9s4GHam/GD6DXf5+rLibUAlYVMhvLOfNm6cPPvjAZ+yDDz7QxIkTJUmFhYXKz8/XgQMHvPdYdnV16dixY9qwYUMUpgwAvlKxA1EkjHb9Mbs7UKB1+saEDN1/qF39j7ds73GrdF+L1n5tjDb9d0fAeyrneD4Oe98AYiNksHzwwQd10003afPmzbr11lt16tQpPffcc/rhD38oSbLZbFqxYoW2bNmioqIiTZs2TZs3b1ZWVpZuv/32qB8AAEip14EoEka7/kSjO9DgdZq9s0GDn5ne65H++e0On7GBhTrcygdYV8hgWVJSohdffFEbNmzQk08+qQkTJmj9+vW67777vNs89NBDcrlcWrNmjVpbWzVnzhzt3r1bDgeXoQDAaox244lF157WC30ht6H6G0gchnqFL1q0SIsWLQr6us1mU3l5ucrLy02bGAAgOox244lF157szBFq7wkeVAmVQGIxFCwBANZX5+xR5Umn6jvdGj/EfaYVJQ6daOr2ucwdqNjJ6HbhzEeSz9ijJVk+91gORKgEEg/BEgCSQDiFNkaLnSIpigo0n2MNXZLNprPnv7j8faKpW2XT7ao67fsYIkIlkJgIlgCQBMIttDFa7DTcoqhA8znb6ZHk+9Xkh063niFUAkkj5HMsAQDWF4tCm3AEm08gA6MmoRJIbARLAEgCsSi0CUew+QyFUAkkPi6FA0AA/YUntecyNeWTloCFJ1Z6CHukhTaS8eIfo/M51njB537K/FFSRtoIn7F+hEogORAsAWAQ38KTNL3V7gpaeDLcLjRmi7T7kNlddiRJHt/7KTNG2HR30WhtfLuDy99AkiJYAsAg4RSeRNKFxmyRdB8yu8tO5Unn5+fsC2c7PYRKIMlxjyUADBJO4Um8imPMZnbxT7DPI1QCyY1gCQCDhFN4Eq/iGLOZXfwT6hwSKoHkxKVwAJZhZvFIJAIVwkwYbfO7xzLc4pjBjHSmiVXRUDjFP0bWKdDn9SNUAsmLYAnAEqJSPDJMAwthaps7NCV3jE/AG05xzGBGO9PEqmjIaPGP0XUqdGRo5cwsrT7ezuVvIIUQLAFYgtnFI5HqL4SpqWlWUdEk77hZczFaIBTLoiEjxT9G12nPGZfWVhMqgVTDPZYALMFqnWOiLZwCoUCs1lFn4Hz2nHHpu39okXtAqiRUAqmBYAnAEqzWOSbahtOZZiCrddTpnw+hEkhtBEsAllBR4tBkh29oibQ4xsoCHe+E0TZNyBoRcixa56XO2aOygy365mtNKjvYojpnT8B5+80xa4QqShyESgDcYwnAGiLtHJNogh2v5F8gFGjM7PMSVvHUoI468nj0X2e7tLa6nVAJpDiCJQDLiKRzTCIKdrxGx8xktCgnWEcdqr8BSFwKBwDIePEUHXUADIVgCQAwXDxFRx0AQ+FSOIC4MNplxyrdeGLFaDeecM6Bkc9cVmzXscYLPg9i7y/KGYiOOgCGQrAEEHNGC0Ws1I0nFox24wnnHITT4aenz78oZzA66gAYCpfCAcTcUIUiw9kuWQTrxjMwAErhnQOjn3m206PGLvmNDd4PHXUADIVgCSDmIi0UoRuP8XNgZocfnlMJIBSCJYCYi7RQhG48xs+BWR1+CJUAjCBYAgnESGeURGC0y0443XiS4dxUlDg0YbTNZyw/UwE77ywrths6XjM6/BAqARhF8Q6QIJKpkMVolx2j2yXTuZHNpoFPhcxIH6Fnr83Wjvdd3nOwrNiulUfaDB1vpB1+3m7uJVQCMIxgCSQIo51REoXRLjtGtkuWc1N50ulfVHO+Tzved/kcR9nBlrCOd7gdfvimEkC4uBQOJIhUK2QJR7KcGysVNREqAQwHwRJIEKlWyBKOZDk3VilqIlQCGC6CJZAgwilkSTXJcm4CFe9MGG2LqKgpXIRKAJHgHksgQRgtZElFSXVuBhXvXPyzr2gdL6ESQKQIlkACMVrwkoqS4dwEK94JVJRj9vESKgGYgUvhAGAR8SpCIlQCMAvBEgAsIh5FSIRKAGbiUjgAxECds0eVJ52q73Rr/KAHkvePLSu260RTt88zKs0oygm070JHBqESgOkIlgAQZYE6Ax1r6JJsNp97Kk80devpBWN9uuxEWpQTrCvRyplZWlvdTqgEYCqCJQBEWaDOQGc7PfKp/tbF7jmDu+xEY98fOt1afbzdZ++ESgBm4B5LAIiyYEU5gZhdqBNs34RKANFAsASAKAtWlBOI2YU6ofZNqARgppDBcuPGjcrOzvb5r7i42Pu6x+PRxo0bNX36dBUUFGjx4sU6ffp0VCcNAFZW5+xR2cEWffO1JpUdbNGyYrtfp5wJo22akOX7T3A4hTqD91Hn7Am4XaAuPf0IlQDMZugey6KiIv3ud7/z/jkt7Yt/pLZu3apt27Zp27ZtKioq0hNPPKElS5bozTfflMORWO3UACBSwYplAhXlSBpW95xg+3hlUa7f+wsdGVo5M4t7KgHEhKFgmZ6ervz8fL9xj8ejqqoqrVq1SqWlpZKkqqoqFRUVadeuXbrnnnvMnS0AWFywYplgRTnDKdQJto9AHXr2nHFpbTWhEkBsGAqWZ86c0fTp0zVy5EhdddVVeuSRR3TZZZeprq5OjY2NWrhwoXdbu92u+fPnq7q6eshgWVNTE/nsETHWwXpYE+sJZ01qz2VK8r/0XNvcoZqaZlPmY3Qf+8+laf27I+XWF/3GbfLoseJuzfF8rET+q8bPifWwJtYTjTUpKioa8vWQwfKqq67Sz372MxUVFencuXN68sknddNNN+n48eNqbGyUJOXl5fm8Jy8vT/X19RFNDNFXU1PDOlgMa2I94a7JlE9a9Fa7y388d4yKiiaZMicj+9hzxqX177Vo4PeaF7+pzEn4byr5ObEe1sR64rUmIYPljTfe6PPnq666SldccYV+/etf66/+6q+iNjEASEQVJY6odM8JZx901AEQL2E/bmjMmDGaPn26amtrvfddNjU1+WzT1NSkcePGmTNDAAnPaAVzMih0ZOjpBWM1aUyaxmbYNGlMmp5eMDai7jlH6l2avbNBk174RLN3NuhsR69eWZSrpVPsurZgpJZOsXsLdwiVAOIp7GDZ1dWlmpoa5efnq7CwUPn5+Tpw4IDP68eOHdPcuXNNnSiAxNRfwbyz1qXDDd3aWevSLfuakzZc1jl7tPJImz7qcKutx6OPOtxaeaRt2Md7pN6l0n0t+qjDrfbPP690X4vOdvTq+etz9O835+n563MIlQAsIWSwrKio0OHDh3XmzBmdOHFCf/d3f6fOzk7dddddstlsWrFihbZu3apXX31V77zzjh588EFlZWXp9ttvj8X8AVjcUBXMycjs411xuE29vp0f1eu5OD4QoRKAFYS8x/KTTz6r8+RWAAAeNklEQVTRfffdp+bmZl166aW66qqr9J//+Z+aNOniDeIPPfSQXC6X1qxZo9bWVs2ZM0e7d+/mGZYAJAVvKWh260KrMPt4Wy/0BRxvGzBOqARgFSGD5S9+8YshX7fZbCovL1d5eblpkwKQPIK1FDS7daFVmH282Zkj1N7jH0rHZl684ESoBGAl9AoHLCpZCl4CtRQ0u0raSoId7zcmZPgU4Byp939cUCBV14xVus13LN12cZxQCcBqDD0gHUBshdOyz+oKHRl6ZVHusFoXJqJAx/uNCRm6/1C7917J9p6LBTh7FuVowXj7kJ+3YLxdexblaMXhNrVd6NPYzBGqumaszl0QoRKA5RAsAQsKp2VfIih0ZCTkvIdr8PHO3tkQtADn1NKhg6V0MVwO3I5vKgFYFZfCAQtKtYKXZGekAMcoQiUAKyNYAhaUagUvyS47M/A/tWODjAdDqARgdQRLwIISpeAlWQqMoi1YAc6jJVmGzx+hEkAi4B5LwIISoeAlmQqMoi1QAc6jJVmq/ON5Q+ePUAkgURAsAYuyesFLshUYRdvgApyygy2Gzh+hEkAi4VI4gGGhwCgyRs4foRJAoiFYAhgWCowiE+r8ESoBJCIuhQMYlooSh040dftczg2nwKjO2aPKk07Vd7o13oL3kJpt8PEuK7YHPX+ESgCJimAJYFgiKTBKtcKfYMf79IKx2vG+y+f8vd3cS6gEkLAIlgCGbbgFRqlW+BPseHe876JQB0BS4R5LADGXaoU/FOoASBUESwAxl2qFPxTqAEgVBEsApjLSjSdROguZZajjJVQCSCbcYwnANEaLchKhs5CZgh0vhToAkg3BEoBpwinKsXpnIbMNPl6+qQSQjLgUDsA0qVaUM1yESgDJimAJwDSpVpQzHIRKAMmMS+EATBNpNx6ri7RbEKESQLIjWAIwTTIX5UTaLYhQCSAVECwBmCpZi3Ii6RZEqASQKrjHEgAMGG5hEqESQCohWAKAAcMpTCJUAkg1XAoHYizSAhDER7DCpGXFdpUdbPFbT0IlgFREsARiKNICEMRPoMKkZcV2rTzS5reeK2dmaW11O6ESQMohWAIxFEkBCOJvcGFS2cGWgOu5+ni7BmRKQiWAlME9lkAM0ZkmuQRbT0IlgFRFsARiiM40ySXYevYjVAJINQRLIIYqShya7PANI8nUmSbZHal3afbOBk164RPN3tmgb0zI8FvPfoRKAKmIYAnEUH8ByNIpdl1bMFJLp9gp3EkQR+pdKt3Xoo863Grv8eijDrfuP9Sub0wYKdugbQmVAFIVxTtAjCVrZ5pkt+Jwm3o9vmO9HqnqtMtnjFAJIJXxjSUAGNB6oS/kNoRKAKmOYAkABmRnDv3PJaESAAiWABBQf6HO14/ZNXtngx68fJTSB99M+TlCJQBcRLAEgEEGFup0uG36qMOtihPndceUTL9tCZUA8AWCJQAMEqxQ58U/X/AZI1QCgC+CJQAMQqEOAAxP2MFyy5Ytys7O1po1a7xjHo9HGzdu1PTp01VQUKDFixfr9OnTpk4UAGKFQh0AGJ6wguWbb76pHTt2aObMmT7jW7du1bZt27Rp0ybt379feXl5WrJkiZxOp6mTBYBw1Dl7VHawRd98rUllB1tU5+zx655zpN7l976qa8ZSqAMAw2A4WLa1tamsrExPP/20srOzveMej0dVVVVatWqVSktLNWPGDFVVVamjo0O7du2KyqQBIJQ6Z49u2desnbUuHW7o1s5al27690/1N4O655Tua/ELlwvG27X2a2P8PpNQCQBDMxws+4Pjdddd5zNeV1enxsZGLVy40Dtmt9s1f/58VVdXmzdTAAhD5UmnPnS6fcYaL0juAEU5Kw63+YztOePSpv/u8BkjVAJAaIZaOv7yl79UbW2tnnvuOb/XGhsbJUl5eXk+43l5eaqvrw/6mTU1NeHME1HCOlgPa2KO2nOZktIMbdvi6vWe9/3n0rT+3ZFyD+gAbpNHjxV3a47nY7E81sDPifWwJtYTjTUpKioa8vWQwbKmpkYbNmzQ3r17lZGREbOJIfpqampYB4thTcwz5ZMWvdXuf/9kIDn2dBUVFWnPGZfWv9eigd9z2uRR1bU5fFNpIfycWA9rYj3xWpOQl8LfeOMNNTc3a968ecrNzVVubq6OHDmi7du3Kzc3Vzk5OZKkpqYmn/c1NTVp3Lhx0Zk1kCACFY8gNipKHJqQ5ftPXO5IKW1QUU667WKxzp4zLn33Dy0+l8ptkh4r7iZUAoBBIb+xXLx4sa688kqfse9973uaOnWq/v7v/17Tpk1Tfn6+Dhw4oJKSEklSV1eXjh07pg0bNkRn1kAC6C8eGXif34mmbr2yKFeFDvO+/ccQPL43VNrTbXpinkM/OnlebRf6NDZzhKquGatzFxQwVFZde4nmeD6O7ZwBIIGFDJbZ2dk+VeCSNHr0aF1yySWaMWOGJGnFihXasmWLioqKNG3aNG3evFlZWVm6/fbbozNrIAEEKh750OlW5Umnnr8+J06zSh2VJ5062+kbLM92erT3bI9OLS3wjgX7prK/UIfbxgDAOEPFO6E89NBDcrlcWrNmjVpbWzVnzhzt3r1bDofDjI8HElJ9pzvgeEOQcZjLyPkPFSoBAOEZVrD8/e9/7/Nnm82m8vJylZeXmzIpIBmMHx24IrkgyDjMFer8EyoBwHz0CgeipKLEockO33Az2ZGmihK+yY+Foc4/oRIAosOUS+EA/BU6MvTKolxVnnSqodOtgtEXQw2FO7ER7Py/3dxLqASAKCFYAlFU6MigUCeOBp9/vqkEgOjiUjiAlECoBIDoI1gCSHqESgCIDS6FA0g4dc4eVZ50qr7TrfEh7l0lVAJA7BAsASSUcDoaESoBILa4FA4goQzV0WggQiUAxB7BEkBCoaMOAFgXwRJAQqGjDgBYF8ESQEKhow4AWBfFOwASCh11AMC6CJYAEg4ddQDAmrgUDiChESoBwDoIlgASFqESAKyFYAkgIREqAcB6CJYAEg6hEgCsiWAJIKEQKgHAugiWABIGoRIArI1gCSAhECoBwPoIlgAsj1AJAImBYAnA0giVAJA4CJYALItQCQCJhWAJwJIIlQCQeAiWACyHUAkAiYlgCcBSCJUAkLgIlgAsg1AJAImNYAnAEgiVAJD4CJYA4o5QCQDJgWAJIK4IlQCQPNLjPQEAya/O2aPKk07Vd7o1fnSaKkocKnRkECoBIMkQLAFEVZ2zR7fsa9aHTrd37ERTt1bOzNLa6nZCJQAkEYIlgKiqPOn0CZWS9KHTrdXH2zUgUxIqASAJcI8lgKiq73QHHCdUAkDyIVgCiKrxo9OGfJ1QCQDJg2AJIKoqShya7AgcLgmVAJBcCJYAoqrQkaGVM7NkGzROqASA5EOwBBBVe864tLaaQh0ASAUESwBRw3MqASC1ECwBRAWhEgBST8hg+fzzz2v+/PmaOHGiJk6cqBtvvFH79u3zvu7xeLRx40ZNnz5dBQUFWrx4sU6fPh3VSQOwNkIlAKSmkMHyy1/+sn70ox/p4MGDOnDggK677jp95zvf0f/8z/9IkrZu3apt27Zp06ZN2r9/v/Ly8rRkyRI5nc6oTx6A9RAqASB1hQyWixcv1o033qgpU6Zo2rRp+uEPf6gxY8bozTfflMfjUVVVlVatWqXS0lLNmDFDVVVV6ujo0K5du2IxfwAWQqgEgNRma21t9YTe7CK3261XXnlFDzzwgP7whz8oKytLV1xxhfbv36+SkhLvdt/+9reVk5OjZ555Juhn1dTURDZzAJay/1ya1r87Uu4BDxayyaPHirv1f8cF7r4DAEgsRUVFQ75uqFf4//7v/+qmm25SV1eXsrKy9MILL2jmzJmqrq6WJOXl5flsn5eXp/r6+ogmhuirqalhHSwmUddkzxmX1r/XooHx8eI3lTkJ/01loq5JMmNNrIc1sZ54rYmhYFlUVKRDhw6pvb1de/bs0YoVK/S73/0u2nMDkAC4/A0A6GfocUMjR47UlClTdMUVV+jRRx/VrFmz9LOf/Uz5+fmSpKamJp/tm5qaNG7cOPNnC8BSCJUAgIGG9RzLvr4+dXd3q7CwUPn5+Tpw4ID3ta6uLh07dkxz5841bZIArIdQCQAYLOSl8Mcee0w33XSTvvKVr3irvQ8fPqyXXnpJNptNK1as0JYtW1RUVKRp06Zp8+bNysrK0u233x6L+QOIA0IlACCQkMGysbFRy5cv16effqovfelLmjlzpnbt2qUbbrhBkvTQQw/J5XJpzZo1am1t1Zw5c7R79245HI6oTx5A7BEqAQDBhAyWVVVVQ75us9lUXl6u8vJy0yYFwJoIlQCAoRiqCgeQOuqcPao86VR9p1vjR6eposShQkcGoRIAEBLBEoBXnbNHt+xr1ofOL55IeaKpWytnZmltdTuhEgAwJIIlAK/Kk06fUClJHzrdWn28XQNbdBEqAQCBDOtxQwCSU31n4NaLhEoAgBEESwBe40enDfk6oRIAMBSCJQCvihKHJjsCh0tCJQAgFIIlAK9CR4ZWzsySbdA4oRIAYATBEoDXnjMura2mUAcAMDwESwCSePg5ACByBEsAhEoAgCkIlkCKI1QCAMxCsARSGKESAGAmgiWQogiVAACzESyBFESoBABEA8ESSDGESgBAtKTHewIAzFfn7FHlSafqO90aPzpNFSUOFToyCJUAgKgiWAJJps7Zo1v2NetDp9s7dqKpWytnZmltdTuhEgAQNQRLIMlUnnT6hEpJ+tDp1urjdNQBAEQX91gCSaa+0x1wnFAJAIg2giWQZMaPThvydUIlACBauBQOJLjBhTrLiu060dTtdzlcIlQCAKKLYAkksGCFOndNtWvj2x1c/gYAxBSXwoEEFqxQh1AJAIgHgiWQwCjUAQBYCcESSGAU6gAArIRgCSSwihKHJjsCh0tCJQAg1giWQAIrdGRo5cws2QaNEyoBAPFAsAQS2J4zLq2tpqMOAMAaCJZAgtpzxqXv/qGF3t8AAMsgWAIJiFAJALAiHpCOpDO4E01FiUOFjox4T8s0hEoAgFURLJFUgnWieWVRblKES0IlAMDKuBSOpBKsE03lSWecZmQeQiUAwOoIlkgqwTrRNAQZTxSESgBAIiBYIqkE60RTEKJDjZURKgEAiYJgiaQSqBPNZMfFAh4rq3P2qOxgix44lamygy2qc/ZIIlQCABILxTtIKoWODL2yKFeVJ51q6HSrIAGqwn0LjtL0VrtLJ5q6tXJmltZWtxMqAQAJg2CJpFPoyNDz1+fEexqGBSs4Wn2cjjoAgMTCpXAgzoIVHBEqAQCJJmSw3LJli77+9a9r4sSJmjp1qu644w698847Ptt4PB5t3LhR06dPV0FBgRYvXqzTp09HbdJAMglWcNSPUAkASBQhg+Xhw4d17733at++fXr11VeVnp6uW265RZ999pl3m61bt2rbtm3atGmT9u/fr7y8PC1ZskROZ+I/OxCItkAFR/0IlQCARBIyWO7evVt/+7d/qxkzZmjmzJl69tlnde7cOR0/flzSxW8rq6qqtGrVKpWWlmrGjBmqqqpSR0eHdu3aFfUDABJdoSNDK2dmyTZonFAJAEg0Yd9j2dHRob6+PmVnZ0uS6urq1NjYqIULF3q3sdvtmj9/vqqrq82bKZCk9pxxaW01hToAgMQXdlX4unXrNGvWLF199dWSpMbGRklSXl6ez3Z5eXmqr68P+jk1NTXh7hpRwDrE1/5zaVr/7ki5B3xfaZNHjxV3a47nY7E81sDPifWwJtbDmlhPNNakqKhoyNfDCpbr16/X8ePHtXfvXqWlRdbJJNTEEH01NTWsQxztOePS+vdaNLAm3CaPqq7N4ZtKC+HnxHpYE+thTawnXmti+FJ4eXm5Xn75Zb366qu67LLLvOP5+fmSpKamJp/tm5qaNG7cOHNmCSSZYB11HivuJlQCABKWoWD58MMPe0NlcXGxz2uFhYXKz8/XgQMHvGNdXV06duyY5s6da+5sgSQwVJvG/zsu8DMtAQBIBCEvha9evVq//e1v9cILLyg7O9t7T2VWVpbGjBkjm82mFStWaMuWLSoqKtK0adO0efNmZWVl6fbbb4/6AQCJJFTvb25RAgAkspDBcvv27ZKk0tJSn/GHH35Y5eXlkqSHHnpILpdLa9asUWtrq+bMmaPdu3fL4XBEYcpAYgoVKgEASHQhg2Vra2vID7HZbCovL/cGTQC+CJUAgFRAr3AgygiVAIBUQbAEoohQCQBIJQRLIEoIlQCAVEOwBKKAUAkASEUES8BkhEoAQKoKu1c4AKnO2aPKk07Vd7o1fnSaKkocKnRkECoBACmNYAmEqc7Zo1v2NetD5xddck40dWvlzCytrW4nVAIAUhbBEghT5UmnT6iUpA+dbq0+3q4BmZJQCQBIOdxjCYSpvjNwP29CJQAg1REsgTCNH5025OuESgBAqiJYAmGqKHFosiNwuCRUAgBSGcESCFOhI0MrZ2bJNmicUAkASHUESyBMe864tLaaQh0AAAYjWAJh4DmVAAAER7AEDCJUAgAwNIIlYAChEgCA0AiWQAiESgAAjCFYAkMgVAIAYBzBEgiCUAkAQHgIlkAAhEoAAMKXHu8JALFQ5+xR5Umn6jvdGj86TRUlDhU6MgJuS6gEAGB4CJZIenXOHt2yr1kfOt3esRNN3XplUa5fuCRUAgAwfFwKR9KrPOn0CZWS9KHTrcqTTp8xQiUAAJEhWCLp1Xe6A443DBgnVAIAEDmCJZLe+NFpAccLPh8nVAIAYA6CJZJeRYlDkx2+4XKy42IBD6ESAADzULyDpFfoyNAri3JVedKphk63Cj6vCn+7uZdQCQCAiQiWSAmFjgw9f32O9898UwkAgPm4FI6UQ6gEACA6CJZIKYRKAACih2CJlEGoBAAgugiWSAmESgAAoo9giaRHqAQAIDYIlkhqhEoAAGKHYImkRagEACC2CJZISoRKAABij2CJpEOoBAAgPgiWSDqZaVKa7Ys/EyoBAIgNgiWSzjcm2vVvC3M0cgShEgCAWKJXOJJSf7hsveAhVAIAECOGvrE8cuSI7rzzTl1++eXKzs7Wiy++6PO6x+PRxo0bNX36dBUUFGjx4sU6ffp0VCYMGPWNiXZCJQAAMWQoWJ4/f14zZszQ448/Lrvd7vf61q1btW3bNm3atEn79+9XXl6elixZIqfTafqEAQAAYE2GguVNN92kRx55RKWlpRoxwvctHo9HVVVVWrVqlUpLSzVjxgxVVVWpo6NDu3btisqkAQAAYD0R32NZV1enxsZGLVy40Dtmt9s1f/58VVdX65577gn4vpqamkh3DROwDtbDmlgPa2I9rIn1sCbWE401KSoqGvL1iINlY2OjJCkvL89nPC8vT/X19cOeGKKvpqaGdbAY1sR6WBPrYU2shzWxnnitCY8bAgAAgCkiDpb5+fmSpKamJp/xpqYmjRs3LtKPBwAAQIKIOFgWFhYqPz9fBw4c8I51dXXp2LFjmjt3bqQfDwAAgARh6B7Ljo4O1dbWSpL6+vp09uxZnTp1SpdccokmTpyoFStWaMuWLSoqKtK0adO0efNmZWVl6fbbb4/q5AEAAGAdhoLlH//4R33rW9/y/nnjxo3auHGj7rrrLlVVVemhhx6Sy+XSmjVr1Nraqjlz5mj37t1yOBxRmzgAAACsxVCwvPbaa9Xa2hr0dZvNpvLycpWXl5s2MQAAACQWqsIBAABgCoIlAAAATGFrbW31xHsSAAAASHx8YwkAAABTECwBAABgCoIlAAAATEGwBAAAgCkIlgAAADAFwRIAAACmIFimqO3bt2v27NnKz8/X9ddfr6NHj8Z7Silhy5Yt+vrXv66JEydq6tSpuuOOO/TOO+/4bOPxeLRx40ZNnz5dBQUFWrx4sU6fPh2nGaeeLVu2KDs7W2vWrPGOsSax19DQoAceeEBTp05Vfn6+5s6dq8OHD3tfZ01iy+12q7Ky0vt7Y/bs2aqsrFRvb693G9Ykuo4cOaI777xTl19+ubKzs/Xiiy/6vG7k/Le2tmr58uWaNGmSJk2apOXLlw/ZWXE4CJYpaPfu3Vq3bp3+4R/+Qa+//rquvvpqLV26VH/5y1/iPbWkd/jwYd17773at2+fXn31VaWnp+uWW27RZ5995t1m69at2rZtmzZt2qT9+/crLy9PS5YskdPpjOPMU8Obb76pHTt2aObMmT7jrElstba2atGiRfJ4PHrppZdUXV2tJ554Qnl5ed5tWJPY+slPfqLt27dr06ZNeuONN/T444/r+eef15YtW7zbsCbRdf78ec2YMUOPP/647Ha73+tGzv99992nU6dOadeuXdq1a5dOnTql+++/39R58oD0FHTDDTdo5syZeuqpp7xjJSUlKi0t1aOPPhrHmaWejo4OTZo0SS+++KJuvvlmeTweTZ8+XWVlZVq9erUkyeVyqaioSD/+8Y91zz33xHnGyautrU3XX3+9nnrqKW3atEkzZszQk08+yZrEwYYNG3TkyBHt27cv4OusSezdcccduuSSS/TMM894xx544AF99tln+u1vf8uaxNhXvvIVPfHEE/rOd74jydjPxHvvvae5c+dq7969mjdvniTp2LFjuvnmm/Xmm2+qqKjIlLnxjWWK6e7u1ttvv62FCxf6jC9cuFDV1dVxmlXq6ujoUF9fn7KzsyVJdXV1amxs9Fkfu92u+fPnsz5RtmrVKpWWluq6667zGWdNYu/3v/+95syZo3vuuUfTpk3TNddco+eee04ez8XvQViT2Js3b54OHz6s999/X5L07rvv6tChQ7rxxhslsSbxZuT8v/HGGxozZozmzp3r3WbevHnKysoydY3STfskJITm5ma53W6fS0qSlJeXp08//TROs0pd69at06xZs3T11VdLkhobGyUp4PrU19fHfH6p4pe//KVqa2v13HPP+b3GmsTemTNn9POf/1wPPvigVq1apT/96U96+OGHJUnLly9nTeJg1apV6ujo0Ny5c5WWlqbe3l6tXr1a9913nyR+TuLNyPn/9NNPlZubK5vN5n3dZrPp0ksvNfX3P8ESiJP169fr+PHj2rt3r9LS0uI9nZRVU1OjDRs2aO/evcrIyIj3dCCpr69PV155pffWnK997Wuqra3V9u3btXz58jjPLjXt3r1bv/nNb7R9+3ZNnz5df/rTn7Ru3TpNmjRJd999d7ynBwvhUniKyc3NVVpampqamnzGm5qaNG7cuDjNKvWUl5fr5Zdf1quvvqrLLrvMO56fny9JrE8MvfHGG2pubta8efOUm5ur3NxcHTlyRNu3b1dubq5ycnIksSaxlJ+fr69+9as+Y8XFxTp79qz3dYk1iaVHHnlEK1eu1G233aaZM2fqzjvv1Pe+9z3967/+qyTWJN6MnP9x48apubnZe0uJdPHezHPnzpm6RgTLFDNy5EhdccUVOnDggM/4gQMHfO67QPQ8/PDD3lBZXFzs81phYaHy8/N91qerq0vHjh1jfaJk8eLFOnr0qA4dOuT978orr9Rtt92mQ4cOadq0aaxJjM2bN08ffPCBz9gHH3ygiRMnSuLnJB46Ozv9rqykpaWpr69PEmsSb0bO/9VXX62Ojg698cYb3m3eeOMNnT9/3tQ1Slu3bt1jpn0aEoLD4dDGjRtVUFCgUaNG6cknn9TRo0f19NNPa+zYsfGeXlJbvXq1fvOb32jHjh2aMGGCzp8/r/Pnz0u6GPptNpvcbrd+8pOfaOrUqXK73frHf/xHNTY26ic/+YkyMzPjfATJZ9SoUcrLy/P5b+fOnZo0aZK+853vsCZxMGHCBG3atEkjRoxQQUGBDh48qMrKSv3gBz/QnDlzWJM4eO+99/Tb3/5W06ZNU0ZGhg4dOqQf//jHuvXWW3XDDTewJjHQ0dGhd999V42NjfrVr36lGTNm6Etf+pK6u7s1duzYkOf/0ksv1YkTJ7Rr1y7NmjVLH3/8sX7wgx+opKTE1EcO8bihFLV9+3Zt3bpVjY2Nuvzyy/XP//zPWrBgQbynlfT6q78He/jhh1VeXi7p4qWJxx9/XDt27FBra6vmzJmjzZs3a8aMGbGcakpbvHix93FDEmsSD/v27dOGDRv0wQcfaMKECSorK9P999/vLTxgTWLL6XTqn/7pn/S73/1O586dU35+vm677TatXbtWo0aNksSaRNuhQ4f0rW99y2/8rrvuUlVVlaHz39raqrVr1+q1116TJN1888164okngv5uGg6CJQAAAEzBPZYAAAAwBcESAAAApiBYAgAAwBQESwAAAJiCYAkAAABTECwBAABgCoIlAAAATEGwBAAAgCn+Pxk5frUEXZYbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "d.plot.plot_diagram(dgms[0][1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sliced Wasserstein Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sliced_wasserstein import SW_approx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5108.143167337765"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SW_approx(dgms[0][1], dgms[1][1], 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[    0.        ,  5108.14316734,  9472.53497695, ...,\n",
       "         1776.37863582,  5490.20107947,  3987.46566694],\n",
       "       [ 5108.14316734,     0.        , 14124.45002927, ...,\n",
       "         3872.70359252,  1866.27573268,  1759.4397636 ],\n",
       "       [ 9472.53497695, 14124.45002927,     0.        , ...,\n",
       "        10464.74527685, 14731.10894512, 12627.65440869],\n",
       "       ...,\n",
       "       [ 1776.37863582,  3872.70359252, 10464.74527685, ...,\n",
       "            0.        ,  4448.29923611,  2525.56709237],\n",
       "       [ 5490.20107947,  1866.27573268, 14731.10894512, ...,\n",
       "         4448.29923611,     0.        ,  2883.85881363],\n",
       "       [ 3987.46566694,  1759.4397636 , 12627.65440869, ...,\n",
       "         2525.56709237,  2883.85881363,     0.        ]])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gram = np.array([[SW_approx(dgms[i][1], dgms[j][1], 10) for i in range(len(dgms))] for j in range(len(dgms))])\n",
    "gram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn_tda import SlicedWasserstein\n",
    "from sliced_wasserstein import diagram_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "too many indices for array",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-e7ab8731d940>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0msw\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSlicedWasserstein\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0msw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdiagram_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdgms\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtolist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/site-packages/sklearn/base.py\u001b[0m in \u001b[0;36mfit_transform\u001b[0;34m(self, X, y, **fit_params)\u001b[0m\n\u001b[1;32m    515\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0my\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    516\u001b[0m             \u001b[0;31m# fit method of arity 1 (unsupervised transformation)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 517\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mfit_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtransform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    518\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    519\u001b[0m             \u001b[0;31m# fit method of arity 2 (supervised transformation)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/share/virtualenvs/tda-networks--KypeAmE/lib/python3.5/site-packages/sklearn_tda/code.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X, y)\u001b[0m\n\u001b[1;32m    429\u001b[0m                 \u001b[0mdiag_thetas\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist_proj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtensordot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdiagram\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mthetas_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\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[0m\n\u001b[1;32m    430\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_directions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 431\u001b[0;31m                     \u001b[0mlist_proj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdiag_thetas\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mj\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[0m\n\u001b[0m\u001b[1;32m    432\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mproj_\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist_proj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    433\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: too many indices for array"
     ]
    }
   ],
   "source": [
    "sw = SlicedWasserstein()\n",
    "sw.fit_transform(diagram_array(dgms[0][1]).tolist())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import OneClassSVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = OneClassSVM(kernel='precomputed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OneClassSVM(cache_size=200, coef0=0.0, degree=3, gamma='auto',\n",
       "      kernel='precomputed', max_iter=-1, nu=0.5, random_state=None,\n",
       "      shrinking=True, tol=0.001, verbose=False)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(gram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1,\n",
       "       1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0,\n",
       "       1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0,\n",
       "       1, 0, 1])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(clf.predict(gram)+1)//2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import KernelPCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = KernelPCA(n_components=2, kernel='precomputed', remove_zero_eig=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([], shape=(69, 0), dtype=float64)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca.fit_transform(gram)"
   ]
  },
  {
   "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
}