tutorial_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**SCM sampler tutorial**\n",
    "\n",
    "*Author*: J.-G. Young <jean-gabriel.young.1@ulaval.ca> <br/>\n",
    "*Date*: May 2017.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's assume that you have managed to compile the source (check the [README](README.md) file otherwise). We will now go over 2 examples of what can be done with the sampler.\n",
    "\n",
    "First, let's import some required modules: `subprocess` (to call the C++ sampler from `python`), `networkx` and `pyplot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import subprocess as sp\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns  # for aesthetic\n",
    "%matplotlib inline\n",
    "sns.set_style('ticks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rejection sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will begin by sampling from a trivial model: $\\mathrm{SCM}(\\mathbf{d},\\mathbf{s})$ with $\\mathbf{d}=(2,2,2,1,1)$ and $\\mathbf{s}=(3,3,2)$, which has 12 non-isomorphic instances, shown below:\n",
    "\n",
    "![](img/simple_instances.png)\n",
    "\n",
    "The support is not very constrained, meaning that the rejection sampler works well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with open('/tmp/degrees.txt', 'w') as f:\n",
    "    print('2 2 1 2 1', file=f)\n",
    "with open('/tmp/sizes.txt', 'w') as f:\n",
    "    print('3 3 2', file=f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Sample:\n",
      "0 1 4 \n",
      "0 1 3 \n",
      "2 3 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "command = ['bin/rejection_sampler', '--degree_seq_file=/tmp/degrees.txt', '--size_seq_file=/tmp/sizes.txt', \n",
    "           '--seed=0']  # for reproductibility\n",
    "proc = sp.run(command, stdout=sp.PIPE)\n",
    "print(proc.stdout.decode())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the sampler generated instance # 8 (node identifiers are off by one, since the sampler represents nodes as 0 indexed integers).\n",
    "\n",
    "To verify that the sampler indeed generates instance uniformly, we will compute the distribution of the length of the cycles in the bipartite graph representation of the simplicial complexes, using `networkx`. It is enough to verify the validity of the sampler, since there are only two possibilities:\n",
    "\n",
    "![](img/shape.png)\n",
    "\n",
    "i.e., one cycle of length 4 or one cycle of length 6, represented with equal probability (6 out of 12 instances have a cycle of length 4, the rest have a cycle of length 6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def facet_list_to_graph(facet_list):\n",
    "    \"\"\"Convert a facet list to a bipartite graph\"\"\"\n",
    "    g = nx.Graph()\n",
    "    for f, facet in enumerate(facet_list):\n",
    "        for v in facet:\n",
    "            g.add_edge('v' + str(v), 'f' + str(f))  # differentiate node types\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rejection_sampling(degree_seq, size_seq, seed=0):\n",
    "    \"\"\"Wrapper around the rejection sampler.\"\"\"\n",
    "    # Write sequence to files\n",
    "    with open('/tmp/degrees.txt', 'w') as f:\n",
    "        print(\" \".join([str(x) for x in degree_seq]), file=f)\n",
    "    with open('/tmp/sizes.txt', 'w') as f:\n",
    "        print(\" \".join([str(x) for x in size_seq]), file=f)\n",
    "    # Call sampler with subprocess\n",
    "    command = ['bin/rejection_sampler', '--degree_seq_file=/tmp/degrees.txt', '--size_seq_file=/tmp/sizes.txt', '-d', str(seed)]\n",
    "    proc = sp.run(command, stdout=sp.PIPE)\n",
    "    # Read output as a facet list\n",
    "    facet_list = []\n",
    "    for line in proc.stdout.decode().split(\"\\n\"):\n",
    "        if line.find(\"#\") == -1 and line != \"\\n\":\n",
    "            facet_list.append([int(x) for x in line.strip().split()])\n",
    "    return facet_list[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cycle of length 4: 96\n",
      "Number of cycle of length 6: 104\n"
     ]
    }
   ],
   "source": [
    "num_iter = 200\n",
    "\n",
    "lengths = []\n",
    "for i in range(num_iter):\n",
    "    g = facet_list_to_graph(rejection_sampling([2,2,1,2,1],[3,3,2],i))\n",
    "    lengths.append(len(nx.cycle_basis(g)[0]))\n",
    "\n",
    "    \n",
    "print(\"Number of cycle of length 4:\", len([x for x in lengths if x == 4]))\n",
    "print(\"Number of cycle of length 6:\", len([x for x in lengths if x == 6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the sampler is correct (there are of course random fluctuations around the mean!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MCMC sampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MCMC sampler generates multiples instances (efficiently), so we will first need to construct a function that parses its output.\n",
    "\n",
    "For the ensemble used above, the output looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Sample:\n",
      "0 1 \n",
      "2 3 4 \n",
      "1 2 3 \n",
      "# Sample:\n",
      "1 2 \n",
      "0 2 3 \n",
      "1 3 4 \n",
      "# Sample:\n",
      "2 3 \n",
      "1 3 4 \n",
      "0 1 2 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "command = ['bin/mcmc_sampler', 'datasets/simple_facet_list.txt', '-t', '3', '--seed=42']\n",
    "proc = sp.run(command, stdout=sp.PIPE)\n",
    "print(proc.stdout.decode())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function below calls the sampler and separates the samples in a list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rejection_sampling(command, seed=0):\n",
    "    # Call sampler with subprocess\n",
    "    proc = sp.run(command, stdout=sp.PIPE)\n",
    "    # Read output as a facet list \n",
    "    facet_list = []\n",
    "    for line in proc.stdout.decode().split(\"\\n\")[1:-1]:\n",
    "        if line.find(\"#\") == 0:\n",
    "            yield facet_list\n",
    "            facet_list = []\n",
    "        else:\n",
    "            facet_list.append([int(x) for x in line.strip().split()])\n",
    "    yield facet_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 1], [2, 3, 4], [1, 2, 3]]\n",
      "[[1, 2], [0, 2, 3], [1, 3, 4]]\n",
      "[[2, 3], [1, 3, 4], [0, 1, 2]]\n"
     ]
    }
   ],
   "source": [
    "# Same command as in cell [8], but parsed:\n",
    "for x in rejection_sampling(command):\n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now redo the cycle calculation, using the MCMC sampler instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cycle of length 4: 100\n",
      "Number of cycle of length 6: 100\n"
     ]
    }
   ],
   "source": [
    "command = ['bin/mcmc_sampler', 'datasets/simple_facet_list.txt', '-t', str(num_iter), '--seed=0']\n",
    "\n",
    "lengths = []\n",
    "for facet_list in rejection_sampling(command):\n",
    "    g = facet_list_to_graph(facet_list)\n",
    "    lengths.append(len(nx.cycle_basis(g)[0]))\n",
    "\n",
    "    \n",
    "print(\"Number of cycle of length 4:\", len([x for x in lengths if x == 4]))\n",
    "print(\"Number of cycle of length 6:\", len([x for x in lengths if x == 6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we have a 50/50 split between both outcomes, as expected. Only this time, the sampler is much faster than the first time around"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's finish by looking at a more complex dataset, where the MCMC is actually needed: the [crime simplicial complex](datasets/crime_facet_list.txt).\n",
    "\n",
    "Since there are still few simplicial properties that are easy to compute, we will instead look at a graphical property  in the bipartite graph representation of the pollinator simplicial complex.\n",
    "\n",
    "The distribution of the degree assortativity coefficient is as good as any.\n",
    "\n",
    "First, we compute the coefficient of the real dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "crime_facet_list = []\n",
    "with open(\"datasets/crime_facet_list.txt\", 'r') as f:\n",
    "    for l in f:\n",
    "        crime_facet_list.append([int(x) for x in l.strip().split()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real assortativity -0.267639538869\n"
     ]
    }
   ],
   "source": [
    "real = nx.degree_assortativity_coefficient(facet_list_to_graph(crime_facet_list))\n",
    "print(\"Real assortativity\", real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find what the distribution looks like in the randomized ensemble:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "command = ['bin/mcmc_sampler', 'datasets/crime_facet_list.txt', '-t', '500', '--seed=0', '-c', '-b', '500', '-f', '500']\n",
    "# We cheat by giving somewhat short decorrelation times, this is just so that the example runs fast enough\n",
    "coeff = []\n",
    "for facet_list in rejection_sampling(command):\n",
    "    g = facet_list_to_graph(facet_list)\n",
    "    coeff.append(nx.degree_assortativity_coefficient(g))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f171bb046a0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAESCAYAAAD+NaylAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGXexvHvlEx6SALppJEEAoTea4CICKLCIoIF66or\nvth2bViwwGvbdUWwvOhaESmiAoq0CAgaCJ3QMhBI7723yXn/iGSNAoaQmTOZ/D7XlYtkMmfOnQPk\nnueU52gURVEQQgghhM3Rqh1ACCGEEOYhJS+EEELYKCl5IYQQwkZJyQshhBA2SkpeCCGEsFFS8kII\nIYSNkpIXQgghbJSUvBBCCGGj9OZ64aqqKp566ikKCgqorKzkwQcfpH///jzxxBOUlZXh6+vLP//5\nTwwGg7kiCCGEEB2axlwz3n3//fdkZmZy7733kpGRwd13303//v2Jjo5mypQpvPbaa4SFhXHjjTea\nY/VCCCFEh2e2kv+tAwcOsHjxYtLS0vjhhx9wcHBg//79LF++nLfeeuuyXqu6uppjx47h5eWFTqcz\nU2IhhBDCephMJvLy8oiKisLBwaHFy5ltd/15M2fOJD8/n2XLlnHrrbc2hfP09CQ/P/+Syy5ZsoSl\nS5eaO6IQQgjRLnzxxRcMHjy4xc83e8mvWbOG48eP89hjjzUbeSuKgkajueSy8+bNY968ec0eS0lJ\n4eqrr+aLL77A19fXLJmFaBMv/Xoo6vmv1M0hhGj3srOzufXWW/Hy8rqs5cxW8gkJCXTu3Bl/f396\n9+5NQ0MDjo6OVFVV4ejoSH5+Pt7e3pf9uuffKPj6+tK1a9e2ji1E23H69U2t/DsVQrSRyz1MbbaS\nP3ToEBkZGTz99NPk5+dTUVFBTEwMsbGxTJ06la1btxIdHW2u1QshWmFTXHKbvM41I0La5HWEEFfG\nbNfJz549m/z8fG655Rb+9re/sWDBAu6//35WrVrFjBkzKC4uZsqUKeZavRBCCNHhmW0kbzAY+Ne/\n/vWHxz///HNzrVIIIYQQvyEz3gkhhBA2SkpeCCGEsFFS8kIIIYSNkpIXQgghbJSUvBBCCGGjpOSF\nEEIIGyUlL4QQQtgoKXkhhBDCRknJCyGEEDZKSl4IIYSwUVLyQgghhI2SkhdCCCFslJS8EEIIYaOk\n5IUQQggbJSUvhBBC2CgpeSGEEMJGSckLIYQQNkpKXgghhLBRUvJCCCGEjZKSF0IIIWyUlLwQQghh\no/RqBxBCWJ7J1IAxtZhDxlwOJeZyNrMUB4MOjQbs7XQ42uuJCHQnPNAdnVbGAkK0V1LyQnQgZZW1\nrN5mZOveFCqq6wHQajUE+bhiamigoKSa8so6TA0KKdll/JKQRd/wLvTu1hkHg/y6EKK9kf+14o/y\n0uDpSRDSp/Hr+hro2gPu/l/Q6i7vdd5+AF7+ru2yxX8PQ69tu9frIOrqTXy3+xyrthmpqKqjcycH\nxg7syoDu3vQN74Kzox0Am+KSASitqOXomXxOnCtgz7Fs9p/MZVRfP6LCuqj3QwghLpuUvLgwv27w\n7Kr/fv3+Y/DLOhj9F/UyAWx4T0r+Mu07kc373ySQW1iJs6Mdd1/Xm6mjQ7HTX/wNm5uzgdH9/Bna\ny4cT5wo5mJjLzkMZFJbVMLqvP1qtxoI/gRCitaTkRcuED4Dsc42fb/0U4jaA0gDDroVr7oG0U/DR\nfNDbgUYLD7134depr4P3HoGiHKiphL88Asb9jW8qxs1ufM4TMfDMSvjshebPyzoLqSfhrfvgkWWw\n+vXGZU11MOluGH5d42t7+DTmKS2Eq++En1ZDbU3jmxY7e0tsLatgalD4YtNJ1sSeRq/TcMPYMGZN\n7I6rk6HFr2Gw09G/uxfdAjrx/c/nSDiTT0lZDVcPD8be7jL26gghVCFn1Ig/V18Hh7ZBaF/IT4d9\nm+C5NfD8Wtj7HRRkQkk+3PosPLMKegyFX7658GulnYKyAnj+K3jqC6gohTE3wp5fd+lnGME7CAqz\n//i8qX8DR9fGgk/c17jeZ1fD01/CN4uhtrrxcIJOD49/Ci4ekJvSmMnJDZKOWG6bqaykvIYFy35h\nTexpfDs78a+Ho/nrDVGXVfC/5eZsYMb4cIJ9XUnNKWPt9tOUVtS2cWohRFsz60j+zTffZO/evdTV\n1XHvvfeyf/9+Dh06hLOzMwD33HMP48aNM2cE0VpZZ2HhrMbP00/B1Adg0MTGgs86C4t+HXVXVUBe\nOrh6wMpXoa4WirJh5LQLv65fGFSUNI64B10No6Y1lnJlaeMbhQNbG5e90PN+6+wROH3wvxmVBijO\nbfw8tG/jn+5eENz7v59XlbXd9rFixtQiXvkknvySaob19uWRmwfi8usx9ythsNMxZVQovxzJ5MiZ\nfL7bfZYbJ0RgkBG9EFbLbCW/b98+Tp48yapVqyguLub6669n1KhRLFq0iJ49e5prtaKt/PaY/OK/\ngU/If7/Xbxz89bXmz194U+NIu/+ExuPmtdUXfl0Hp8YT8U7Fw48rIH4j/M9SGHkDHNgCx3+Gx/4D\n9o4Xft5vRc+EG+b9cR06/YU/7wCOGPN4+eO91NWZuH1KT2aMj2jT4+dajYbR/QNAA0dO57NlbwpT\nRoWi1cgxeiGskdl21w8YMIC33noLADc3N+rq6igr6xgjKZtz83xY9RrUVEFoHzi5p/FzRWk8bl5b\nDeXFjW8EaqvhUGzjLv4LOZcAezZArxFw+wuNI3KAETfAztWNx9PtHS/+PFPjZV+E9W9cT4OpcZ2f\nLTDvNmgH9p3I5sX/7MFkUnjmrmHMjOluthPkRvbxJ9DHhZTsMvYeyzbLOoQQV85swxy9Xo9e3/jy\na9asITo6mry8PN5++23Kysrw8fHhueeew93d3VwRRFvxDoIhk+HbJTDricYT7RbeBBpN4650gwNM\nugsW3w8evnDtfY3lP+K6C7/W6jdg8yeNX894rPHPTl0aR/kjrr/080L7wAvT4IVvoddIePEvjW82\nYm4z4wawbpvikklKL2bL3lS0WpgyMpTC0uqmy+HMQavVcPWwYL768TQHE3Pp3MmB7kEeZlufEKJ1\nNIqiKOZcwbZt23jvvff4+OOP2bt3L6GhoYSHh7Ns2TKys7N5/vnnL7rskiVLWLp06QW/FxsbS9eu\nXc0VW1haWRG8eiu8vOHyrsW3Zo+MavzzrZ/NuprFqw4RG5+KXq9l6qhQ/L1czLq+3yoqrearH09j\nalCYPi4cH08nAK4ZEWKxDEJ0BOnp6cTExFx295n17Ppdu3bx7rvv8uGHH+Lm5sbEiRMJDw8HICYm\nBqPReMnl582bR2JiYrOP2NhYc0YWajiwtbHgb3nGdgreQvadyCZ2XyoGOx3Xj+lm0YIH8HBzYOKw\nYEwNClvjU6g3NVh0/UKISzNbyZeVlfHqq6+ybNkyPDwad+PNnTuX9PR0AOLj44mIiDDX6kV7Mmgi\nLNoIvUepnaRdMaYW8drn+9FpNUwdHYpvZ2dVcoT4udEvogsl5bXsOyHH54WwJmY7Jr9x40ZKSkp4\n9NFHmx77y1/+wqOPPoq9vT3Ozs688sor5lq9EDYtM7+cFz/cQ12diWtGhKhW8OcN6+3LucxSDhnz\nCO8q59kIYS3MVvKzZs1i1qxZf3h8+vTp5lqlEB1CcVkNLyzbQ2lFLXNv7Ic1XLxmp9cxbmBX1u86\ny48H0rh5UiR6ncy1JYTa5H+hEO1IbZ2Jlz/aQ1ZBBbMmdmeyFZ3gFujjSmSIB/nF1Xyz44zacYQQ\nSMkL0a783zcJGFOLmTA4kFsnRaod5w9G9fXH0V7Pl1sSycgrVzuOEB2elLwQ7cSWvSls2ZtCt4BO\njbvprXCWOQeDnrEDAqirb+Ddr45g5it0hRB/omPN+SmEFbmcyWpyiyr5evsZ7O10jOrrx/b9aWbL\ndaXCAjoxuKcP+0/msO9kDkN7+aodSYgOS0byQli56tp6NsWlYGpQmDg0CDdn675drkaj4a6pvdBq\n4OMNxzHJtfNCqEZKXggrpigKW+NTKausZUhPH4L93NSO1CJBvm5MHBZMem45W/amqB1HiA5LSl4I\nK3b0TD6p2WUE+bgypJeP2nEuy62TInEw6FixOZHK6ovcsEgIYVZS8kJYqcLSauISsnAw6IgZEmiV\nJ9pdioebAzMmRFBcXsPa7XJJnRBqkJIXwgqZGhRi96VialAYPygQJwc7tSO1yrSxYXi6OfDtziTy\ni6vUjiNEhyMlL4QVOnAyh9yiKnoEe9AtoJPacVrNwV7PnMmR1NaZ+PyHk2rHEaLDkZIXwsrkFFay\n/1QOLo52jOkfoHacKzZ+cBAhfm7sOJBGWk6Z2nGE6FCk5IWwIvWmBrbtS0VRIGZIIPZ27f/Wuzqt\nhlsm9aBBgZVbE9WOI0SHIiUvhBU5cDKH4rIa+oZ3oau3q9px2szwKD+6BXRi1+EMUrJL1Y4jRIch\nJS+ElSgqreZgYh4ujnYMi7KtWeI0Gg23TopEUWDlFhnNC2EpUvJCWAFFUdhxMJ0GRWHMgAAM+va/\nm/73hvTyITzQnd1HMknOktG8EJYgJS+EFUhMLSIzv4IQPze6+bffs+kv5fxoHuDLLadUTiNExyAl\nL4TKqmvr+flIJnqd1ibOpr+UQZHe9Ajy4JejWZzNKFE7jhA2T0peCJXFJWRRXWtiaC8f3JwNascx\nK41Gwy2/juZXbJbRvBDmJiUvhIqyCyo4ca6Qzp0c6BvhpXYcixjQw4ueIZ7sPZ4to3khzExKXgiV\nKIrC7iOZAIztH4BO277mpm8tjUbDrIndAVgda1Q5jRC2TUpeCJWcSS8hp7CSsIBO+Hu5qB3Hogb2\n8CY80J1fjmbKLHhCmJGUvBAqqK0zEZeQhVajYUQfP7XjWJxGo2HWVd1RFBnNC2FOUvJCqGDDrrOU\nVdbSJ7wLnVzs1Y6jiqG9fAnxc+Ong+lk5perHUcImyQlL4SFlZTXsDrWiL1Bx+Ce3mrHUY1Wq+Gm\nmO40KPBV7Gm14whhk6TkhbCwFZtPUVldz5CePjgY9GrHUdXIfv4EeDnz4/40cosq1Y4jhM2RkhfC\ngtJyyti0J4UAL2eiwrqoHUd1Oq2GmTHdMTUorP1RRvNCtLWOPYwQwsK+2HSKhgaFO67tTUl5jdpx\nzGZTXHKLn2tqUHBzNrB5TwreHk44O9o1fe+aESFtnk2IjkRG8kJYSFJ6MT8fzaR7kDvDbewuc1dC\np9UwsIc3pgaFQ8Y8teMIYVPMOpJ/88032bt3L3V1ddx7770MHTqUJ554grKyMnx9ffnnP/+JwWDb\n03gKcd7yTY3TuN52TU80mo4x8U1LRQZ7sO9kDsfPFjAo0htHe9nJKERbMNtIft++fZw8eZJVq1bx\n0Ucf8corr/D6668zY8YMVq9eTUBAAOvXrzfX6oWwKqeSC9l/Mofe3TrTv3vHmL72cuh0WgZ096Le\n1MCR0zKaF6KtmK3kBwwYwFtvvQWAm5sbdXV17NmzhwkTJgAQExPD7t27zbV6IazK5z+cBGDOZBnF\nX0zvbp1xtNdz9Ew+1bX1ascRwiaYreT1ej3Ozs4ArFmzhujoaKqqqnBwcADA09OT/Pz8S77GkiVL\n6NGjR7OPmJgYc0UWwiyOnM7j6Jl8Bvbwpne3zmrHsVp6nZb+3b2oq28g4UyB2nGEsAlmP/Fu27Zt\nrF69mvnz52Nn99+zZhVF+dMRzbx580hMTGz2ERsba+7IQrQZBYUvzh+LnxypchrrF9WtM/YGHUfO\n5FFbZ1I7jhDtnllLfteuXbz77rt8+OGHuLm54ezsTFVVFQD5+fl4e3fc2b5Ex1Bb18DJ5EKGR/kS\nEeihdhyrZ7DT0S/ci5paE8fOymheiCtltpIvKyvj1VdfZdmyZXh4NP5yGzNmTNNIfOvWrURHR5tr\n9UKoTkGhoqoOgFsmySi+pfqGd8Gg13LYmCfH5oW4Qma7TmXjxo2UlJTw6KOPNj326quv8tRTT/Hx\nxx8TGhrKlClTzLV6IcympRO9jKqso66+gVB/NxJTikhMKTJrLlthb9DRJ7wLB07lsikuhWnRYWpH\nEqLdMlvJz5o1i1mzZv3h8c8//9xcqxTCqpw/pjy4p4/KSdqf/t29OHomn7XbT3PNiOAOP8e/EK0l\nM94JYQaZ+eXUmxT0Og3eHk5qx2l3HAx6+oV3obis5rKmyBVCNCclL4QZHDiZCzSeSCZap193L5wc\n9Kz98QzVNXJsXojWkJIXoo3lFFaSmlOGXqdBr5P/Yq3lYNBz3ZhuFJfXsPGXZLXjCNEuyW8gIdrY\ngVM5gIzi28K0sWE4O+j5esdpGc0L0QpS8kK0ofziKs5lluLr6YROJ9PXXikXJwPXjw2jpLyW738+\np3YcIdodKXkh2tDBxMZj8YN6+qBBSr4tXN80mj9DlYzmhbgsUvJCtJHSihrOpBXTuZMDwb6uasex\nGS6OdtwQHU5pRS3f7kxSO44Q7YqUvBBt5LAxDwUY0MNb7jTXxm4Y2w13F3u+2XGaorJqteMI0W60\naIYJo9HIunXrKCsrQ1GUpsdffvllswUToj2pqqnnZHIhrk52hHd1VzuOzXFysGP21T14/+ujrNyS\nyAMz+qkdSYh2oUUj+YceegiDwUCvXr3o3bt304cQolHCmXzqTQr9IrzQaWUUbw6Thgfj38WZzXtS\nyMgrVzuOEO1Ci0by/v7+PPzww+bOIkS7VFdv4mhSPvYGHb1CPdWOY7P0Oi23X9uLVz/dx2cbT/D0\nHUMtnqEtZ9+7ZkRIm72WEBfTopK/4YYbeO+99+jfvz96/X8XGTJkiNmCCdFenEwupKbWxOCePtjp\n5dp4cxrZx48ewR78cjSLU8mFRIbImyohLqVFJb9hwwaSk5PZtWtX02MajYYvvvjCbMGEaA9MDQqH\njXnodRr6hndRO47N02g03DW1N0+9s5uPvzvOqw+OlpMchbiEFpV8WVkZ27ZtM3cWIdqdpPRiyirr\n6BPWGUd7uVOaJfTu1plhvX3ZezybXxKyGNXXX+1IQlitFp14N2TIENLS0sydRYh2RVEUDhlz0dB4\na1RhOXdd1xu9TssH3yZQWV2ndhwhrFaLhh47duzgk08+wd3dHb1ej6IoaDQaduzYYeZ4QlivjLxy\n8ourCe/aCTdne7XjdCgBXi7cOCGClVsTWbE5kb/eEKV2JCGsUotK/r333kOrbT7oLy0tNUsgIdqL\nQ8Y8APp391Y5Scc0MyaCnYfS2bArifGDuhIm8xMI8Qct2l0fGBhIVVUVmZmZZGZmkpyczBNPPGHu\nbEJYrYKSalKzy/Dr7IyPp5PacTokg52OuTP60qDAO18dwdSg/PlCQnQwLRrJL1y4kN27d1NQUEBA\nQAAZGRn89a9/NXc2IazWkdO/juJ7yLF4NfXv7k30gK7sPJTOprhkrh0VqnYkIaxKi0byCQkJbNq0\nicjISL799ls+/PBD2V0vOqyK6joSU4vo5GIg1M9N7Tgd3j3X98bZQc9nG09QWCrz2gvxWy0q+fPH\n400mEyaTiX79+nHkyBGzBhPCWh07k09Dg0L/CC+5RtsKeLg5cMe1vaisrmfxykM0yG57IZq0aHd9\nZGQkn3zyCVFRUdx5550EBgZSUVFh7mxCWJ3qmnoSzhbgYNDRI1hmWzO3y5lGNsjXlYOJubzyaTyD\nIn2afU+mkBUdVYtKfsGCBZSWluLk5MSGDRsoKSmRuexFhxS7P+03U9jKnZqthUaj4aohQazamsje\n49n4dXHGv4uL2rGEUN2f/pZSFIWFCxfi5uaGXq9n+vTp3Hnnnfj4+PzZokLYFFODwrqfktBpNfQJ\n66x2HPE7jvZ6rh4WDAps2ZtKVU292pGEUN2flrxGo8HBwYG4uDiqq6tpaGho+hCiI9l3Ipus/Aq6\nB3ng5GCndhxxAf5eLgzt7UtFVR2x+1JRFDk+Lzq2Fu2uX7lyJf/5z3+a/YfRaDScPHnSbMGEsDbr\nfkoCoF+EXDZnzQZFepOZV05KdhnxJ3IY1ttX7UhCqKZFJb9//35z5xDCqp1JK+ZYUgEDe3jTuZOD\n2nHEJWg0Gq4aGsRXP55m/8kc7O10cuKd6LBaVPJLliz5w6VCer2ekJAQJk2adNHLiIxGI3PnzuXO\nO+/ktttu4+WXX+bQoUM4OzsDcM899zBu3Lgr+wmEsIDzo/gbosPILaxUOY34M04OdtwwNoxvdpzh\n56OZ/PDLOSaPlIlyRMfTopIvLS1l+/btjBw5Eq1WS1xcHIMHD+bw4cPs3LmTV1555Q/LVFZW8vLL\nLzNixIhmjy1atIiePXu23U8ghJnlF1ex63AGQb6uDOjuxeY9KWpHEi3QycW+seh3JvHu2qPYG3RM\nGBykdiwhLKpFJZ+Zmcn69etxcmqco7uyspJnnnmG999/n9mzZ19wGYPBwAcffMAHH3zQ9JhcWy/a\no+92n8XUoDBtbJhMftPOeLg5cP2Ybnz/8zkWrzwEaJgwOFDtWEJYTItKPi0tDb3+v0+1t7fn7Nmz\nmEwmampqLvzCen2zZaCx5N9++23Kysrw8fHhueeew9394neOWrJkCUuXLm1JRCHMoqqmnk17UnB3\nsSd6YFe144hW6OLuyIv3jeC5//uFf395kNOpRdx9fZTMcyA6hBaV/Lhx45g0aRJ9+vQB4NSpU4wZ\nM4ZvvvmG6OjoFq9s9uzZhIaGEh4ezrJly3j77bd5/vnnL/r8efPmMW/evGaPpaenExMT0+J1CnEl\nYvelUlFVxy2TIjHY6dSOI1qpe5AHbz4SzSufxPPdz+c4nVbMk7cPwcvDUe1oQphVi0r+scceY/r0\n6Zw+fZqGhgYefPBBevTogclkQqdr+S++iRMnNn0eExPDggULLj+xEBZialBY/9NZ7PRaJsvZ2e1e\ngJcL/3xoLO+sPcKOA+k8/OYOHp7Vn2FRfmpHE8JsWrS/ymQycfjwYfbu3cv+/fs5ffo0wGUVPMDc\nuXNJT08HID4+noiIiMuMK4TlxB/PIquggvGDAnF3tVc7jmgDDvZ6Hrt5IHNn9KWqpp6FH8ez4IM4\nUrPlrprCNrVoJP/CCy9QUlLCsGHDAPjuu+84fPgwzz777EWXOXbsGK+99hoZGRno9Xo2b97Mbbfd\nxqOPPoq9vT3Ozs4XPCtfCGvx7c5fL5sb203lJKItaTQaJo8MpVe3zny47hgHT+Vy2JjHlBEh3Dwp\nEjdng9oRhWgzLSr5pKQkVqxY0fT1Lbfcwi233HLJZaKiovj888//8PikSZMuM6IQlmdMLeLEuUIG\nRXoT5Cv3jLdFwb5uvHTfCPadyOE/64/x3c/n+PFAGteN7sYN0WG4OknZi/avRSVfX1/f7Pi7zF0v\nbN26X0fx06PDVU4izEmj0TC0ty8Deniz8ZdzfBV7mlXbjKzfdZapo0OZFh0uI3vRrrX47PoZM2Yw\nfPhwAPbu3SsjcmGzcosq2X00kxA/N/pGdFE7jrAAO72WG8aGMWl4MJviklm7/QxrYk/z/c/nmD2x\nB1NHd5NL7kS71KKSnzt3LiNGjODo0aMoisKLL75I3759zZ1NCFVs2HWWhgaFadEy+U1H42DQMy06\nnGtGhLApLplVW418tOE4P8Qlc891vVEURf5NiHalRSX/1FNP8eqrrzJgwICmx+69995ms9kJYQsq\nq+vYsjcFTzd7xg6QyW86qvNlP2FwEF9uOcXGX5JZ+HE8QT6uxAwJlFsNi3bjkiW/fv16Vq5cyenT\np7n11lubHq+qqqK4uNjs4YSwtC17U6msrmfG+AjZPStwczZw//S+TBkZygffJnDImMfq2NNcMzwY\n387OascT4k9dsuSvv/56hg0bxj/+8Y9mM89ptVrCw+WEJGFbTKYGNuxKwt6gY/LIELXjCCsS6OPK\ni/eN4JVP9rHnWBbf7ExibP8AenfrrHY0IS7pT4cqPj4+LFu2DD8/P4YOHYqdnR0nTpyQ41LC5vxy\nNIvcoipiBgfK5VPiDzQaDQMjvbluTONJeDsOprP9QBoNDYra0YS4qBbtj3zyyScxGo1kZGTwyCOP\nYDQaefLJJ82dTQiLURSFr3ecRqOBaXLZnLiEQB9XborpThd3R06cK2TnoXQURYpeWKcWlXx2djYx\nMTFNs9YtXLhQbhsrbMqxpALOpJcwoo8ffl3kWKu4NDdnA9Ojw5qKft+JHLUjCXFBLT6zSFEUtm7d\n2nTXubq6OrOFEsLSvt5xBoDp42QUL1rGYKdj6uhQ3JwN7DuZw7GzBWpHEuIPWlTyQ4cOZeDAgbi6\nutK9e3eWLVtGt24yn7ewDSnZpew/mUOvUE8igz3VjiPaEWcHO64b3Q0Hg46fDqZzNrNE7UhCNNOi\n6+T/8Y9/cN999+Hm1jiH95QpU/Dzk9szCtvQNIWtjOJFK7i72jN1dDe+3ZnElj0pzJgQgZe73Kde\nWIdLlvz//d//cf/99/P4449f8Gz6119/3WzBhLCEwtJqth9IJ8DLmaG9fNWOI9opH08nJg0L5vtf\nzhG7L5WZEyLQ6WSeBaG+S5Z8ZGQkDQ0NDB06VC6ZEzbpu91nqTc1MC06HK1W/o2L1gvxd6N3t84c\nP1tA/IkcRvSRvZ1CfZcs+fvvvx+NRtN0ecjvP7/xxhvNn1AIM6msrmPjL8l0cjEwfnCg2nGEDRjZ\nx4+0nDIOJeYS6u8ms+IJ1V2y5E+dOtXs6zlz5lzwHvFCtAeb4pKbfX3YmEdFVR1De/myfX+aKpmE\nbTHY6ZgwOJBvdyYRuy+Nm67qLtMjC1Vd1r8+2WUvbIXJ1MBhYy56nZY+4TI1qWg7AV4u9IvoQnF5\nDXuOZakdR3Rwl1XyMquTsBWJqUVUVNfTu1tnHAwtushEiBYbHuWHu6s9R8/kk5lfrnYc0YHJSF50\nOA2KwiFjHlqNhv4RXdSOI2yQXqcl5tfzPHYfyZQBklDNJYcw0dHRzYq9oKCAcePGoSgKGo2GHTt2\nmDufEG3uXGYJxWU19AzxxEVuRCPMxLezMxGB7pxOK8aYWkyPYA+1I4kO6JIlv2LFCkvlEMIiFEXh\n4KlcAAZ091I5jbB1w6P8OJtRwp5jWYR17YRerp0XFnbJkg8ICLBUDiEsIiOvnNyiKroFdMLDzUHt\nOMJCfn/A08iyAAAgAElEQVRlhaW4ORvoG96FQ8Y8DhvzGNzTR5UcouOSt5WiQzk/ih/Yw1vlJKKj\nGNTTBweDjoOJuVRWy429hGVJyYsOI6ewkrTccgK8XPDxdFI7jugg7O10DO3lS119A/FyS1phYVLy\nosPYf7LxF+zgnjKKF5bVq1tn3F3tOXG2gMLSarXjiA5ESl50CEnpxSRnleLb2YkALxe144gORqfV\nMLKPHwoQfzxb7TiiA5GSFx3Cqm1GAIb09JH5HoQqQvzc8PZwJCmjhIKSKrXjiA7CrCVvNBq56qqr\nWL58OdB4nf0999zDTTfdxEMPPURtba05Vy8EAClZpcQlZOHt4Uigj6vacUQHpdFoGPLr7Yz3nZRj\n88IyzFbylZWVvPzyy4wYMaLpsddff50ZM2awevVqAgICWL9+vblWL0ST/47ifWUUL1QV7OvaOJpP\nLyElq1TtOKIDMFvJGwwGPvjgA7y9/3uSU3x8PBMmTAAgJiaG3bt3m2v1QgCQllPG7iMZdAvoRLCf\njOKFujQaDUN6No7mV25NVDmN6AjMdmcOvV6PXt/85SsqKnBwaJyAxNPTk/z8/Eu+xpIlS1i6dKm5\nIooOYHWsEUWB2RO7U1Iuh4eE+oL9XPHycOTno5mkZpcS5OumdiRhwyx64p2dnV3T5+fnv7+UefPm\nkZiY2OwjNjbW3DGFjcjMK+eng+kE+7oyrLef2nGEAM6P5n1QFFi11ah2HGHjLFryzs7OVFU1nlWa\nn5/fbFe+EG1txeZEGhSYfXUPtFo5Fi+sR4ifG2FdO7HrSAap2XJsXpiPRUt+zJgxTSPxrVu3Eh0d\nbcnViw4kOauUnw6n0y2gEyP7+KsdR4hmNBoNsyf2QFFg9bbTascRNsxsJX/s2DHmzJnDN998w2ef\nfcacOXO47777WLVqFTNmzKC4uJgpU6aYa/Wig1v+w0kUBeZM7imjeGGVhvX2JcTPjV2H08nKr1A7\njrBRZjvxLioqis8///wPj1/oMSHaUmJKIXuPZ9MzxJNBkXJISFgnjUbDzJgI3lh+gLXbT/M/M/ur\nHUnYIJnxTtic5T+cAmDOlJ5yXbywaqP6BeDXxZnYfWkyC54wCyl5YVOOnsnj8Ok8BnT3ok9YF7Xj\nCHFJOq2GGydEUG9q4JsdSWrHETZISl7YDEVR+HzjSQBum9xT5TRCtMz4QYF06eTApj3JlJTXqB1H\n2BgpeWEz9p3I4VRKESP6+NE9yEPtOEK0iJ1ey/Tx4dTUmtiw66zacYSNkZIXNqHe1MBHG46j1Wq4\n7ZpIteMIcVmuHhZMJxcD3+0+S2V1ndpxhA2Rkhc2YVNcMhl55UwaHizThIp2x8Gg54axYVRU1/P9\nz+fUjiNsiJS8aPfKK2tZsfkUTg56bp0ko3jRPk0ZGYqzg55vdyZRXVOvdhxhI6TkRbu3apuRsso6\nZl3VnU4u9mrHEaJVnB3tuG5MGKUVtfwQl6x2HGEjpORFu5aZX853u8/i4+nE1NHd1I4jxBW5fmw3\nHO31fL3jDDV1JrXjCBsgJS/atU++O0G9SeHOqb0w2OnUjiPEFXF1MjB1dCjFZTVs3pOsdhxhA6Tk\nRbuVkJRPXEIWPUM8GdVXbkIjbMMNY8NwMOhY++MZamU0L66QlLxol+rqG3hv7VE0GvjrDVEyfa2w\nGZ1c7JkyMpTC0mq2xqeqHUe0c1Lyol36ducZ0nLKuGZ4iEx8I2zOtHFhGOx0fBVrpK5eRvOi9aTk\nRbuTXVDByq1G3F3suX2KTF8rbI+HqwOTR4SQX1JN7L40teOIdsxst5oVoi1siktu9rWiKHz/8zlq\n60yMHRDA7iOZquQSwtz+Mj6cjb+cY9U2IxMGB8qJpaJVZCQv2pWzGSWkZJfR1duF7oHuascRwmw8\n3RyYOrob+cVVMgueaDUpedFu1NaZ2HU4A61WQ/SArnKynbB5M2MicHbQsybWSHmVzGkvLp+UvGg3\n9hzLoqK6nkE9vHF3lZnthO1zdTJwY0x3yirr+Hr7abXjiHZISl60C2k5ZSQkFeDhZs/ASG+14whh\nMVNHh+Lp5sC6n85SUFKldhzRzkjJC6tXXVtP7P40tBqYOCQIvU7+2YqOw8Gg55ZJkdTWmfhyS6La\ncUQ7I78thdXbdTiDiqo6hvTyxcvDSe04QljcVUMC6ertwtb4VNJzy9SOI9oRKXlh1c6kF2NMLcbH\n04mBPWQ3veiYdDott0/pSUODwmcbT6odR7QjUvLCahWWVrPjYDp6nYarhgSh1crZ9KLjGh7lR88Q\nT+ISsth/MkftOKKdkJIXVsnUoPDvLw9SU2tiZB9/OZtedHgajYa5N/ZDp9Xw/tdHqa6tVzuSaAek\n5IVVWrH5FIeNeYT4uREV1lntOEJYhRA/N6ZFh5FTWMnqbUa144h2QEpeWJ3449ms3mbEt7MTMUMC\nZdIbIX5j9sQeeHs48vX2M6Rkl6odR1g5mbteWJWs/AreXHEAg17L03cMxZhapHYkIczi9/dluBxD\nevry/S/nWPRRPO8/FSPnq4iLsuhI/tixY4wdO5Y5c+YwZ84cXn75ZUuuXli56tp6Xvk0norqeh6Y\n0Y9uAZ3UjiSEVQrxd6NbQCeyCirYtk/uOS8uzqIj+crKSiZNmsQzzzxjydWKdkBRFN5Zc4RzmaVM\nGh7MVUOD1I4khFUb08+ftJwyPtpwnH4RXvh4yhwS4o8sOpKvqKiw5OpEO7J80yl2HEynR5AH903r\no3YcIayei5OB0f38qaiq4/XP91FX36B2JGGFLFrylZWVHDhwgLvuuovbbruNuLg4S65eWKlNccms\n3mbEr7Mzz949TO6bLUQL9QzxZNygrhhTi/n0+xNqxxFWyKK76yMjI7n//vuZNGkSKSkp3HnnnWze\nvBmDwXDB5y9ZsoSlS5daMqKwsPgT2by39ghuzgZeuG+4XA8vxGXQaDTMndGPM2nFrPspid7dOjOi\nj5/asYQVsehIPiwsjEmTJgEQHBxMly5dyMm5+MxN8+bNIzExsdlHbGyspeIKMzOmFvH65/vR63U8\nf88w/Lu4qB1JiHbH0V7PU7cPwWCnY/HKg2QXyGFR8V8WLflvvvmGTz75BICCggIKCgrw8fGxZARh\nJc5llvDih3uoqzPx+G2D6BHsqXYkIdqtYD83HvhLHyqq63nts31U18hseKKRRXfXX3XVVTz++ONs\n2bKF+vp6FixYcNFd9aJ9u9Q1wPnFVaz7KYnqWhPjB3WluKzmiq4ZFkLAVUODOXGukK3xqSz6JJ7n\n7xmGnV7Ob+noLFryrq6uvP/++5ZcpbAyvy/4XqEyZa0QbWXujf0oKa8l/kQ2r3++n6duH4JOJxOb\ndmTyty8sRgpeCPPS67Q8eftg+oZ3Yc+xbN5adYiGBkXtWEJFUvLCIrLyK/hWCl4IszPY6Xj27mH0\nCPZgx4F03vv6qBR9ByYlL8zubEYJ635KorbOxITBgVLwQpiZo72eF/46nBA/NzbFJfOqnIzXYUnJ\nC7M6lpTPprhkNBoN144MpWeInEUvhCW4OBlY9MAo+oR1IS4hiyeX7ia3qFLtWMLCpOSFWSiKQlxC\nFjsPZeBgr2d6dBjBfm5qxxKiQ3FzNvDifSOYNDyYs5kl/H3xT5xKLlQ7lrAgKXnR5iqr69j4SzIH\nE3Pp5GJgxvhwvOXmGUKowk6v5cEb+3H/9D6UVtTy9Ls/89WPpzGZZK77jkBKXrSp9Nwy/r74J5Kz\nSunq7cKN4yPo5CJT1QqhJo1Gw9TR3Xjhr8NxcbLj0+9P8PiSXaRklaodTZiZRa+TF7Yt/kQ2//ri\nAJXV9fSP8GJEHz+0Wo3asYQQvxrQw5t3n5jAsm8T2HEgnUf+vYPZE3vwl/ERxLbRfemvGRHSJq8j\n2oaUvLhi9aYGlv9wkrXbz2DQa/n7LQOprjWpHUsIcQGuTgb+fssgxvQP4J01R1i+6RSx+9LoG9GF\nsIBOaDTyxtyWSMmLK5JbWMkby/dzKqUIv87OPHH7YMK7uss0tUJYuaG9fOn1RGdWbD7Fxp/PsXlP\nCr6dnRjV1x/fzs5qxxNtREpetFpcQiaLVx2moqqO6AFdmXtjX5wc7NSOJYRoIRdHO+6b1oepo0J5\n7fP9nM0oYe32M4T6uzGkpw9eHnLCbHsnJS8uW2V1HR+uO8bW+FQMdjoeuqk/Vw0Nkt18QrRT/l4u\nTB4RQmZ+OXEJWZzLLOVcZimh/m4M7eVLF3dHtSOKVpKSF5cl4Uw+b608SG5RFd38O/H3WwcS5CvX\nvwthC/y7uPCXceGk5ZQTfyK7WdkPjvSRS2HbISl50czFjqXX1Tew93gWR07no9HA4EhvBvfy4cS5\nQk6ck8k1hLAVGo2GIF9XAn1cSMspI/5ETlPZB/m4MrinD35d5Jh9eyElL/5Uem4Z2w+kU1pRi7uL\nPTFDAuXEHCFsXGPZuxHo40p6bjkHTuWSmlNGak4Z/l2cGRTpQ6CPixyms3JS8uKiqmvr+floJqeS\ni9AA/SO8GNrbFzu9zKEkREeh0WgI9HEl0MeVrPwK9p/KITW7jMzdZ/HycGRQpDfd/OXSO2slJS/+\nQFEUjKnF/Hw0k6qaerq4OzB+UCDecqatEB2aXxdnrhvdjbyiSg4k5pKUXsKmuBQ8XO0Z0MOb7kHu\nakcUvyMlL5rJK6rkp8MZZBdUotNqGBHlR7/uXuhk5johxK+8PJy4ZngIRWXVHDyVizG1iB/3pxF/\nPJu6+gYmDQ/B0V7qxRrI34IAoKS8huWbTjWdeBcW0ImRff1xczaomksIYb08XB2IGRLE0N6+HDmd\nx/Gzhfxn/XFWbjVyzfBgpo7uJpffqUxKvoOrrK5j3c4kvtl5hqoaEx6u9ozpH0Cgj6va0YQQLWAN\ns0u6OhkY3S+AwT19qK418f3uc6zdfoZvdyYxZkAAN4wJIzxQduWrQUq+g6qtM7EpLplV24yUVtTS\nycXAnMm90Go1smteCNEqDgY906LD+cu4cHYcTOfbnUnsOJDOjgPp9Aj2YMrIUEb388dgp1M7aoch\nJd/BVFTV8UNcMut/SqKorAYnBz23XRPJ9WPDcLTXW8WoQAjRvhnsdFw9LJiJQ4M4mJjLd7vPceBU\nDokpRfxn/TEmDg0iZkiQ7DG0ACn5DiK/uIoNu87yQ1wyVTX1ONrrmTE+nOnjwuV+70KINnOhgcKw\n3r5EBntw/GwBJ5ILWbv9DGu3n8HH04meIZ6EB7pjf4HRvdy29spJyduwelMD+07ksGVvCgdP5dCg\ngIerPTdd1Z3JI0JwdpSbyQghLKOTiz0j+/oztLcv5zJLOJVcRGpOGTmFlew6nEGQryvhXd0J8XOT\n3fltSErexjQ0KBhTi/j5aCY7D6ZTVFYDQPcgdyYND2HcwK7yH0gIoRq9TktEoAcRgR6UV9aSmFqE\nMbW4aepcnVZDsK8bIf5uDIvyxcPVQe3I7ZqUvA2orq3nxLlC9h7LYs+xLApLG4vd2dGOqaNDuXpY\nMKH+nVROKYQQzbk4GRgU6cOgSB8KSqpJSi/mTHoxZzNLOJtZwvYDaUQEujO4py/9IroQEeghM25e\nJin5dqiyuo6k9BKOJeVz5Ew+iSmF1JsUAFyd7LhqSBAj+vrRP8JLRu1CiHahcycHOnfyZUgvH4rL\nakjOLqWsoo4T5wowphazYjPYG3T0CvGkT3gXIkM8iejqjoNMunNJsnWsmKIo5BVXkZ5TTlpuWdO7\n3PTccpTGTkejaZy4pm+4FwN7eBMV1hmdTt7pCiHaJ41Gg4ebAx5uDlwzIoTyqjqOnM4j4Uw+R8/k\nc8iYxyFjHgBaDQT5utEj2INQ/06E+LkR7OuKi5NM4nWexUt+8eLFxMXFUVtby4svvkifPn0sHcFq\n1NWbKCqrobishsLSanKLKsktrGr8s6iSjNxyqmtNzZZxctDTJ6wL4V3diQzxpE9YZ1ycDGyKSyar\noIKsggp1fhghhDADF0c7RvX1Z1RffwCKSqs5drYAY2oRiSlFJKUXk5xV2myZzp0cCPBywa+LM76d\nnfHr4oyXuyOdOzng7mLfoQZCFi35PXv2kJCQwMqVKzEajbz44ot88cUXloxwQTV1Jiqq6lB+HR4r\nyvkPhYbzHw2NH6bf/Flvavj1Q6G+voGaOhO1v37U1JmoqjFRVVNPVU09ldV1lFfVUV5ZS1llHaUV\ntVRU1V00k8FOR4CXM129XQn0dqGrtyuhAW74d3FBK5PVCCE6KA83B8b0D2BM/wCg8Sqi1OwykrNK\nSMkqIzm7lJSsUo7+OvL/Pa2m8Ux/d1d7XJ0MjR/OBpwd9DjaN3442OtxMOiw0+uw02sx2Gmx0+nQ\n6RonC9PptOi0GrRaDVpN458aDWjQoP31/YO9nc4q9ihYtOT37t1LTEwMAN27dyc3N5eqqiocHVs+\nt7HJ1Diyzc7ObpNMNbUmHl+yi8rqixduW9LptLg46nF2NODnbaCTiz1uzgY6OdvT2d3h1+NSjrg6\n2f3u1o0K1JaQmVlywdctyGub7SHaTkZV479V+bsRonXS01tWUQagu6+W7r6dgMaTjGvqTOQXVZFb\nXEleURVFpTUUl//6UVZMamod1bX1Zsuu0cDDswYQFdalTV7vfOed78CWsmjJ5+XlERkZ2fS1p6cn\n+fn5BAYGXvD5S5YsYenSpRf83q233mqWjEK0ubg71E4ghFDBw7Ft/5p5eXkEBwe3+PkWLXk7u+aT\nryiK8rvRanPz5s1j3rx5zR6rrq7m2LFjeHl5odO1/zPHY2JiiI01w78EGyPbqWVkO7WMbKeWke3U\nMpbYTiaTiby8PKKioi5rOYuWvJeXFwUFBU1fFxYW0qXL5e3KcHBwYPDgwW0dTVVdu3ZVO0K7INup\nZWQ7tYxsp5aR7dQylthOlzOCP8+ipxiOHTu26d3O8ePHCQwMxMFBZjMSQgghzMGiI/moqCgiIyOZ\nPn06Op2ORYsWWXL1QgghRIdi8evkH3/8cUuvUgghhOiQdC+88MILaofo6IYNG6Z2hHZBtlPLyHZq\nGdlOLSPbqWWsdTtplPMzwAghhBDCpnScuf2EEEKIDkZKXgghhLBRUvJCCCGEjZKSF0IIIWyUlLwQ\nQghho6TkhRBCCBtl8clwOiKTycSLL76I0WgE4I033vjDnfeWLl3Krl27qK+vZ/z48fzP//xPi5az\nJS35eX/44Qc++ugjoPG61H/84x/k5OTw9NNPU1dXR319PfPnz6dPnz4Wz28prd1OAB999BHfffcd\nAC+88AJ9+/a1YHLLupLtBJCfn8/kyZNZunSp1V4DfaVau41MJhPPPvssKSkp1NbW8sQTTzB06FCL\n57eUK9lOqv8OV4TZffXVV8rzzz+vKIqi/Pjjj8rjjz/e7PupqanKvHnzFEVRlPr6emXixIlKdnb2\nny5na/7s562urlbGjx+vlJeXKw0NDcrMmTOVxMRE5dVXX1VWrFihKIqi7N+/X7n77rstnt2SWrud\nUlNTlTlz5igNDQ1KQkKC8tZbb6kR32Jau53Oe/zxx5Xp06cre/bssWhuS2rtNvrmm2+UBQsWKIqi\nKKdPn1amT59u6egW1drtZA2/w2V3vQXs3buXmJgYAMaMGUN8fHyz7wcGBvL2228DUFJSgk6nw9nZ\n+U+XszV/9vPa29uzbt06nJ2d0Wg0dOrUifLycjw8PJrublhSUoKHh4fFs1tSa7fT1q1buf7669Fo\nNERFRfHwww+rEd9iWrudAOLi4nBxcaF79+4Wz21Jrd1GU6ZM4cknnwTAw8ODiooKi2e3pNZuJ2v4\nHS4lbwF5eXl4enoCoNfrMZlMmEymPzxv4cKFXHvttTz44IO4uLi0eDlb0ZKf19XVFQCj0Uh2djZR\nUVHcfvvtbN68mWuvvZYFCxbw0EMPWTy7JbV2O2VlZZGZmcnf/vY37rjjDk6dOmXx7JbU2u1UW1vL\nu+++y6OPPmrxzJbW2m1kMBhwdHQE4NNPP2Xq1KmWDW5hrd1O1vA7XI7Jt7E1a9awZs2aZo9d6Jep\nRqP5w2PPPvssDz74IHfccQd9+/bFzs6uRcu1R1eynZKTk3nsscd4/fXXMRgMLF26lMmTJzN37ly2\nb9/Oa6+9xjvvvGO27JbUltuptraW0tJS3n//ffbv38/8+fP5+uuvzZbdktr639PNN9/c9EvbVrTl\nNjrviy++4NixY7z//vttH1glbbmdrOJ3uMUPEHRA8+fPV3bs2KEoiqLU1NQo48ePb/b9zMxM5ciR\nI01fv/DCC8r333//p8vZmpb8vFlZWco111zTbHvdc889SkJCQtNyY8eOtUxglbR2Oy1ZskTZuHFj\n09fDhg0zf1gVtXY7zZo1S5k5c6Yyc+ZMZdiwYcrkyZMVo9FosdyW1NptpCiKsnr1auXOO+9Uqqqq\nLJJVTa3dTtbwO1x211vA2LFjiY2NBWDHjh2MHDmy2fdLSkp46aWXmnblHDt2jNDQ0D9dzta05Oed\nP38+CxYsaHZWeFBQEMeOHQPg5MmTBAcHWyawSlq7nUaPHs3u3buBxl2Kfn5+lgmsktZup5UrV7J6\n9WpWr17NuHHjWLBgARERERbLbUmt3UZpaWmsWLGCd955BwcHB4vlVUtrt5M1/A6Xu9BZgMlkYv78\n+Zw+fRpHR0f+9a9/4evry7JlyxgyZAgDBgxg2bJlbN26FUVRGD9+PA8++OBFl7NVf7ad3N3dmTZt\nWrP/RHfeeSd9+vTh6aefpra2loaGBp577jkiIyNV/EnMq7XbKSYmhiVLlhAfH09VVRXPPvss/fv3\nV/EnMa8r2U7nPfXUU0yfPt2mL6FrzTY6cuQI33//Pf7+/k2P/+c//2m2K9+WtHY7jRs3TvXf4VLy\nQgghhI2S3fVCCCGEjZKSF0IIIWyUlLwQQghho6TkhRBCCBslJS+EEELYKCl5ISwoPT2dqKgo5syZ\nw80338y0adN44403qK6uVjuaxW3YsIGGhoZLPufgwYOkpaUBsGjRoqb5EC7kt99ft25d2wUVoh2T\nS+iEsKD09HRuueUWfvrpJwBqa2t55ZVXyM7O5r333rui125oaECrtf737edzXn311WzcuBG9/uKz\naz/77LNMmTLlsiYRMZlMTJkyhc2bN7dFXCHaNZm7XggVGQwG5s+fz1VXXcWZM2cICwvj9ddfJyEh\nAZPJRO/evXnmmWdQFIX58+djNBoJCQlBo9EwYsQIhg8fzgMPPEDv3r0JDAxk7ty5F1xeo9Hw8ccf\ns3XrVrRaLT4+PixcuLDpJiPnLV26lJ9++gmtVouvry9vvPEGGo2G5557jnPnzlFbW0uvXr146aWX\niI+P55///CcGg4Gamhrmz5/PoEGD2LFjB++88w729vbY29uzaNEifH19mTBhAtdddx3JycmEhYWR\nkpLCnXfeydKlS9myZQtr1qxpWubf//43e/fuZdOmTRw9epSnn36ad999lwceeIB//etfPPvsswwY\nMABonHTkrrvu4sMPP+SBBx5g3bp1ZGRkcPfdd+Ph4cHo0aOZPn06AM8//zyRkZHccsstFv+7FkIN\n1v+2XwgbZ2dnR8+ePTEajWzatIni4mKWL1/Ol19+SW5uLlu2bOHnn3/mzJkzfPXVVzz55JNNRQyQ\nlJTEvffey4MPPnjR5Y8ePcrOnTtZvnw5y5cvx9vbmy+//LJZjvr6euzs7Pjss89YuXIlpaWl7N69\nG6PRyKFDh1i5ciVff/01gwYNori4mE8//ZS77rqLFStWsHjxYvLy8qiqquKZZ55h8eLFLF++nOjo\naP797383rSMgIIDFixc33Snwk08+wd3dnerqat59912WL1+Ov78/69evZ+LEifTs2ZOnnnqKESNG\nNL3Gddddx6ZNmwAoKCggKSmJUaNGNX1/3rx5eHp68tFHHzF79uymm/CYTCZ2797NddddZ56/SCGs\nkIzkhbACGo0GOzs74uPjOXjwIHPmzAGgvLycjIwMAAYNGtQ0Cv/t9Jlubm6EhYUBjcewL7R8VlYW\nycnJ3HHHHQBUVVUxcODAZhn0ej319fXcdddd6PV6kpKSKCoqYvjw4Tg5OXHPPfcwYcIEJk2ahKen\nJ5MmTeKNN97g6NGjjBs3jmuuuYaTJ0/i5eXVNN3piBEjWLVqVdM6zo++f8/e3p6HH34YrVZLRkYG\n3t7eF91W1157LbNnz+bpp59m06ZNTJ48+aK7/IcMGUJhYSHp6emkpqYyaNAgm7u7nBCXIiUvhMpq\na2s5cuQIzzzzDPv27eOmm27innvuafacDz/8sNnx9t/ervK384VrNJoLLv/pp58yYcIEnn/++Yvm\n2Lt3L+vXr2ft2rW4uLgwd+5cABwdHVm7di1Hjx5lx44dTJs2jRUrVjBt2jTGjh3L7t27Wbp0Kdu2\nbWPGjBnNXlNRlGa5LzS3eVpaGm+++SYbNmzA29ubRYsWXWpz4eXlRVBQEEePHuWHH37gqaeeuuTz\nZ86cybp168jOzmbmzJmXfK4QtkZ21wuhorq6Ov73f/+XUaNG0bVrVwYNGsS2bduor68H4N133yUp\nKYnQ0FCOHj2Koijk5uaSkJBwwde72PIDBw5k165dVFRUAPDll1+yf//+ZsuWlJTg7e2Ni4sLqamp\nJCQkUFtbS0JCAitXrqRfv348/PDDhIeHYzQaefvtt6mrq+P666/n0UcfZf/+/YSGhpKfn09WVhYA\nu3btol+/fhfMqtFoqK6uprS0FGdnZ7y9vSksLOSXX36htra22XN+77rrrmPNmjWUlJQQFRXV7Hta\nrZaampqmr6dNm8bWrVtJTExkyJAhf/p3IoQtkZG8EBZWWFjInDlzUBSFiooKxowZw/z58wG4+uqr\nOXLkCLNnz0aj0RAZGUlQUBBBQUGsW7eOGTNmEBERwcCBA9HpdH947Ystb2dnx6233sqcOXMwGAx4\neXn94dj06NGj+eSTT5g9ezahoaE89NBDvP/++7zzzjvs2rWL9evXA+Dv78+YMWOoqalh7ty5ODs7\nU8YOZk0AAAECSURBVF1dzWOPPYajoyOLFi3ioYcewmAw4OLictGR+ZgxY7j55pt55513iIyMZMaM\nGfj7+/PII4/w0ksvER0dzahRo1i4cGHTm5bf/pwvv/wy999//x9e19vbGx8fH2688UY+++wz3N3d\nCQ4Opnfv3s32gAjREcgldEK0A2VlZcTGxjJt2jRMJhPTp09n4cKFzY7NiwsrLS3lpptuYsWKFXh6\neqodRwiLkpG8EO2As7Mz+/bt47PPPkOr1TJ27Fgp+Bb4+uuv+fTTT/n73/8uBS86JBnJCyGEEDZK\nTrwTQgghbJSUvBBCCGGjpOSFEEIIGyUlL4QQQtgoKXkhhBDCRv0/T2cQx4+3ttAAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f171ba7f208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "sns.distplot(coeff, norm_hist=True)\n",
    "plt.plot([real, real], [0,30], color=\"#ff5b1e\")\n",
    "plt.text(-0.29,25, 'Real system', color=\"#ff5b1e\")\n",
    "plt.ylim(0,30)\n",
    "plt.xlabel('Degree assortativity')\n",
    "plt.ylabel('Histogram')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The assortativity of the real system therefore seems to be typical. \n",
    "\n",
    "**Unrelated note**: The fact that we are looking at a bipartite representation of a simplicial ensemble more or less explain this behavior, so there is nothing worthwhile here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further information: \n",
    "\n",
    "\"*Construction of and efficient sampling from the simplicial configuration model*\"<br/>\n",
    "[J.-G. Young](http://jgyoung.ca), [G. Petri](https://github.com/lordgrilo), F. Vaccarino, and [A. Patania](https://github.com/alpatania) (2017)<br/>\n",
    "Phys. Rev. E, **96**, 032312 (2017)\n\n",
    "Read it on: [arXiv](https://arxiv.org/abs/1705.10298) | [Phys. Rev. E](https://doi.org/10.1103/PhysRevE.96.032312)\n"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}