notebooks/choicerank-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ChoiceRank to understand network traffic\n",
    "\n",
    "This notebook provides a quick example on how to use ChoiceRank to estimate transitions along the edges of a network based only on the marginal traffic at the nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import choix\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Generating sample data\n",
    "\n",
    "First, we will generate sample data.\n",
    "This includes\n",
    "\n",
    "1. generating a network,\n",
    "2. generating a parameter for each node of the network,\n",
    "3. generating samples of choices in the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_items = 8\n",
    "p_edge = 0.3\n",
    "n_samples = 3000\n",
    "\n",
    "# 1. Generate a network.\n",
    "graph = nx.erdos_renyi_graph(n_items, p_edge, directed=True)\n",
    "\n",
    "# 2. Generate a parameter for each node.\n",
    "params = choix.generate_params(n_items, interval=2.0)\n",
    "\n",
    "# 3. Generate samples of choices in the network.\n",
    "transitions = np.zeros((n_items, n_items))\n",
    "for _ in range(n_samples):\n",
    "    src = np.random.choice(n_items)\n",
    "    neighbors = graph.successors(src)\n",
    "    if len(neighbors) == 0:\n",
    "        continue\n",
    "    dst = choix.compare(neighbors, params)\n",
    "    transitions[src, dst] += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The network looks as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAFBCAYAAACrYazjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVGX/P/D3AZcBQUCUcl9we9xT0MwNc8ctt9TMhRa/\n+a0nK3ctl/wpj0aauaVPGpqVigu5oLmG5hK4oGhuZIo7CIgigyzz+f2B8nVgBplhZg7L+3VdXjVn\nuc891sWb8zn3uW9FRARERESkCju1O0BERFScMYiJiIhUxCAmIiJSEYOYiIhIRQxiIiIiFTGIiYiI\nVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIRg5iIiEhFDGIiIiIVMYiJiIhUxCAmIiJS\nEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUhGDmIiISEUMYiIiIhUxiImIiFTEICYiIlIRg5iIiEhF\nDGIiIiIVMYiJiIhUxCAmIiJSEYOYiIhIRQxiIiIiFTGIiYiIVMQgJiIiUhGDmIiISEUMYiIiIhUx\niImIiFTEICYiIlJRCbU7QERE1hUTE4O1gYG4fPYskhIT4eTigrpNmmCknx8qVKigdveKPUVERO1O\nEBGR5YWHh2ORvz927tqF/gC8U1LgDOARgDAHB2wVQc8ePTB2yhR4e3ur3Nvii0FMRFQErVy+HDPG\nj8ckrRYjReBm4JgEAIGKgvkODpgVEIDRY8bYupsEBjERUZGzcvlyzBs/Hr8lJ6N2Ho6PAtDN0RGT\nGMaqYBATERUh4eHh6OPjg8PZQng4gH0AkgGUB/AOgGnP7Y8C0M7REdtDQ+Hl5WW7DhNHTRMRFSWL\n/P0xSavNcSc8BcA/ABIB7AKwGMBvz+2vDWCiVotF/v626Shl4R0xEVERERMTg3rVq+NqSorBZ8LP\nXALQGcCvAJo/tz0egKdGg8vR0RxNbUO8IyYiKiLWBgaiH2A0hD8EUAZAI2SWpZtn218OQD9FwdrA\nQKv1kXJiEBMRFRGXz55Fy5QUo/uXAkhC5rPizwGEGzjGW6vF5chI63SQDGIQExEVEUmJiXB+wTEK\ngA4ABgH4xcB+ZwCPEhIs3TXKBYOYiKiIcHJxwaM8HpsOwNHA9kcAnN1ye8JMlsYgJiIqIuo2aYIw\njSbH9lgAGwA8BqBD5mjpIAB9DbQR7uCAuo0bW7OblA1HTRMRFRFRUVFoVq8ebuh0egO27gMYCOAs\nAAFQB8AXAHpnO5+jptXBO2IioiLgjz/+QOfOnSE6HVZn21cewO/IDNoEAGHIGcIAsEZR0MvXlyFs\nY7wjJiIqxNLS0jBz5kz85z//gU6nAwCUBXASyNP0ls9wZi31cBlEIqJC6vLlyxg2bBhOnDiht/0R\ngPaKgkMiJs01PSsggCGsApamiYgKGRHBypUr8corr+QIYQAY+8knmPbNN2jn6IiFigJjLyPFA1ig\nKGjHBR9UxdI0EVEhEhsbi/feew/btm3Lsa9ixYpYs2YNunTpAgA4ceIEFvn7Y0dICPopCry12qz1\niMOfrkfcy9cXY6dM4Z2wihjERESFxK5du+Dn54d79+7l2Ne/f3+sXLkS7u7uOfbFxsZibWAgLkdG\n4lFCApzd3FC3cWOMGDWKA7MKAAYxEVEBp9VqMXHiRCxZsiTHvjJlyuDbb7+Fn58fFEVRoXeUXxys\nRURUgJ0+fRrDhg3DhQsXcux79dVXsW7dOnh6eqrQM7IU3hETERVQIoLmzZsjIiJCb7u9vT2++OIL\nTJs2DSVK8H6qsOOoaSKiAkpRFDRq1EhvW61atXD48GHMmDGDIVxEMIiJiAqgmzdvok6dOtiwYQMG\nDRoEAPDz80NERARat26tcu/IkhjEREQFzIoVK1CzZk2UKFEC0dHRWL9+Pfbv34/Vq1fD2flFCx1S\nYcNnxERWFhMTk/nqyNmzSEpMhJOLC+o2aYKRfn58dYT0JCcno3v37jhy5AgmTpwIf39/tbtENsAg\nJrKS8PBwLPL3x85du9AfgHdKStZkCmFPJ1Po2aMHxk6ZAm9vb5V7S2rbv38/+vbtC41Gg3379qFZ\ns2Zqd4lshEFMZAUrly/HjPHjMUmrxUgRGFpmPQFAoKJgvoMDZnF6wWJLp9PBz88PP/74I/r374/1\n69dzEFYxwyAmsrCVy5dj3vjx+C052aQJ9znXb/Fz6dIldOzYEQkJCVi/fj369u2rdpdIBRysRWRB\n4eHhmGEghJcC8AagAfBOtnNqA/gtORkzxo83OIE/FU1z5sxBgwYNULlyZdy7d48hXIwxiIksaJG/\nPyZptTnuhCsD+ALAu0bOqw1golaLRRycU+TFx8ejWbNmmDFjBubPn4/w8HCULVtW7W6RiliaJrKQ\nmJgY1KteHVdTUgw+EwYyw/gWgNUG9sUD8NRocDk6mqOpi6iNGzdixIgR8PDwwMGDBzk1JQHgHTGR\nxawNDEQ/wGgIv0g5AP0UBWsDAy3XKSoQ0tPT0bt3bwwZMgQjRozAtWvXGMKUhUPziCzk8tmzaJmS\nkq82vLVaRERGWqhHVBCEh4ejW7duSE9Px++//4727dur3SUqYHhHTGQhSYmJyO+cR84AHiUkWKI7\nVAB89tlnaNWqFZo3b46YmBiGMBnEICayECcXFzzKZxuPAFy4cgVHjx6FTqezRLdIBTdv3kTt2rWx\nZMkSrFq1Cvv27YNGo1G7W1RAMYiJLKRukyb408gP2wwAKU//mQ7gydN/zy4UwJlLl9CmTRvUqFED\n48ePR3h4ODimsvBYvnw5atasiZIlSyI6Ohp+fn5qd4kKOI6aJrKQK1eu4JX69XFDp8sxYGvW0z/K\nc9tmAJj+3Od4ZL7mZOgpc61atTB48GAMHjwYTZo0gaIoBo4iNXGeaDIXg5jIAo4cOYK33noL96Oj\n8SWAcWa0EYDMYNa+4Lj69etnhfK//vUvM65ElrZ3717069cPDg4O2Lt3L+eJJpMwiInyISMjA3Pm\nzMGsWbOynumWBXASyNP0ls9EAWjn6IhPZ8zAyZMnsX37dmi1L4pkoHHjxnjvvffw8ccfm9N9yqfn\n54keMGAANmzYADs7PvEj0/D1JSIz3bhxA2+//TYOHTqkt/0RgHYADiNvYfxsrunnF35ISkrCjh07\nsGHDBoSEhCA1NdXguZGRkTh16lR+vgaZ6fl5ordu3copKsls/NWNyAxbtmxB06ZNc4QwAAiAewBa\nAPgamassGRIPYIGioJ2BBR+cnJwwZMgQbN26FTExMVizZg18fX0Nrsqzc+dOjBs3joO6bIjzRJMl\nsTRNZILk5GR89tlnWLFixQuPrV69Olo1boyQXbvwhqLgtfT0rPWIw5+uR9zL1xdjp0yBl5dXnq4f\nFxeHrVu3YsOGDThw4ADc3Nywb98+bNmyBRs2bEB6ejrefPNNDB48GE2bNuWgLguLj4/H66+/jnPn\nzmHevHkYN86c0QBE+hjERCaIi4tD06ZNcevWLb3t9vb2yMj4vxeS3Nzc8Oeff6JOnTqoV68eenTr\nBu2DB3iUkABnNzfUbdwYI0aNytec0vfu3cOFCxfg4+MDABARREREYMOGDdiwYQNKly6dNairQYMG\nZl+HMj0/T3RoaChq1qypdpeoiGAQE5kgIiIC7du3x6NHmVN3lChRAmXLlkV8fHzWMfb29ti9ezc6\nd+6Mu3fv4l//+hfu378Pe3t7m/VTRBAWFoYNGzZg48aNKFeuXFYo165tyjAySk9PxxtvvIGQkBC8\n9957+O677zggiyyK/zcR5dGsWbPQokULNG/eHJMnT4anpyeaNm2qF8IAsHDhQnTu3BkAcPDgQfj4\n+Ng0hAFAURS0atUKCxYsQHR0NJYuXYo7d+6gbdu28PLywldffYXr16/btE+FUXh4ODw8PHDo0CH8\n/vvvWLlyJUOYLE+IKFeJiYnSrFkzsbe3l8WLF4uISGpqqrzzzjuCzLFZWX9Gjx4tOp0u69z33ntP\nvv32W7W6nkNaWprs27dP3n//fXF3d5fWrVvLN998I7du3VK7awXOp59+KoqiSOfOnUWr1ardHSrC\nGMREudi5c6c4ODhIpUqVJCoqKmv7kiVLcoRw+/bt5cmTJ3rn16xZU86dO2frbudJamqqhISEyMiR\nI8XV1VU6dOggy5Ytk5iYGLW7pqobN25IrVq1pGTJkhIYGKh2d6gYYBATGZCRkSEjR44URVHkrbfe\nkoyMjKx9+/btE3t7e70QrlGjhsTGxuq1cfXqVfHw8NC7Qy6otFqtBAcHy9ChQ8XFxUW6dOkiq1at\nkvj4eLW7ZnUPHjyQDz74QG7evCnLli2TEiVKSP369eXOnTtqd42KCQYxUTbXrl2TqlWrikajkeDg\nYL19ly9fFjc3N70QdnJyksjIyBztrFq1SoYMGWKrblvM48ePZePGjTJgwAApW7as9OzZU9auXSuJ\niYlqd83iQkNDpVq1agJAXF1dRVEUmTJlitrdomKGow6InrNixQrUrl0brq6uuHXrlt5EDYmJiejT\npw8SnlsvWFEU/PTTT2jUqFGOtg4cOIDXX3/dJv22JEdHRwwaNAibNm3CzZs3MXToUAQFBaFq1aro\n378/Nm7ciMePH6vdzXxJTU3F1KlT4ePjg+joaADAgwcPMH78eMydO1fl3lGxo/ZvAkQFwZMnT6RT\np05iZ2dn8I4oPT1dunfvnuO5sL+/v8H2dDqdVKxYUe+5cmEXHx8vq1evlm7duomLi4sMHjxYtm7d\nWugGMl28eFFatGiR478lAPH29pb09HS1u0jFDIOYir2wsDBxcXERV1dXCQ8PN3jMuHHjcvzQHjZs\nmNHnvxcuXJDq1asXiufD5oiJiZHvvvtOOnbsKK6urjJ8+HDZuXNnjsFqBYlOp5Ply5eLg4ODwRD+\n9NNPC90vFVQ0MIipWJs8ebIoiiJdunQxGiI//PBDjh/aLVu2lOTkZKPtLl26VPz8/KzV7QLl9u3b\n8u2330qbNm3E3d1d3n33Xdm7d6+kpaWp3bUs9+7dk169ehkM4EqVKsnevXvV7iIVY5xZi4ql+Ph4\ndOjQARcvXsTSpUsxevRog8cdPXoUHTt21Fv9qFKlSjhx4gQqVqxotP0BAwagX79+ePvtty3e94Is\nOjoaQUFB2LBhA65fv44BAwZgyJAhaNu2rUUnwoiJicHawEBcPnsWSYmJcHJxQd0mTTDSzy/HtKEh\nISHw8/NDTExMjnYGDBiAFStWwN3d3WJ9IzKZ2r8JENnali1bpHTp0lK1alW5fv260eOuX78uHh4e\nendPGo3GaPn6mYyMDClXrpzcvHnT0l0vVKKiomTu3LnStGlTqVSpkowdO1aOHTuWr3J9WFiYDOvX\nT1w1GnlHo5HlgKwDZDkgfg4O4qrRyLB+/SQsLEweP34s//u//2vwLtjJyUl++OGHIvvogAoXBjEV\nGxkZGTJ06FBRFEX8/Pz03g3O7smTJ9KsWbMcP8B/+eWXF17n1KlTUq9ePUt2vdC7cOGCzJw5U+rX\nry/Vq1eXCRMmyMmTJ00KwhXLlsnLjo6yUFEkHhAx8CcekAWKIh4ajVR86SWDIdy6desiNYiOCj8G\nMRULUVFRUrFiRXFwcJCQkJA8nbN06VK9H+DTpk3L03kBAQEyZsyY/HS3yNLpdHLmzBmZOnWqeHp6\nSu3atWXatGkG38N+3oply6SWo6NcMRLA2f9cAeRlQJTn/vvZ29vLrFmzCtSzayIRPiOmYmDRokUY\nN24cmjZtioMHD6Js2bJ5Oq9Fixa4cuUKSpcujbZt22Lz5s15es7Zs2dPvPPOOxgwYEB+u16kiQhO\nnjyZtWyjs7MzhgwZgsGDB6Nu3bpZx4WHh6OPjw8OJyfj2bpRzgCerbQsALQAPgSw6Ln2owC0APAQ\ngKenJ9atW4dXX33V6t+LyFQMYiqyUlJS0LVrVxw5cgQzZszA9OnT83yun58ffvrpJ1y6dAmKoqB8\n+fJwcnJ64XlpaWkoX748rl69ygFAJtDpdDh+/HjWso0vv/wyhgwZgjfffBNfjBsHr+BgfGLkR9Vj\nABUB7ALQJtu+AAA/VK+O45GRcHZ2tu6XIDITg5iKpCNHjqBHjx4oWbIkDh48iCZNmuT53GXLluGj\njz7Cjh074Ovra9J1jx49ig8//BCnT582tcv0VEZGBg4fPowNGzYgKCgI2vh43BSBm5Hj1wCYjcw7\n4OziAXhqNLgcHZ1jNDVRQVFC7Q4UR6a8ekGm++yzz/DNN9/A19cXwcHBKFEi7/+bHzlyBP/+978x\ne/Zsk0MYKLzTWhYk9vb28PHxgY+PD2pWr47zM2bA7bnXx7JbC2CEkX3lAPRTFKwNDMS4CROs0V2i\nfGMQ21B4eDgW+ftj565d6A/AOyUFzgAeAQjbsgV1Z8xAzx49MHbKFHh7e6vc28InJiYG7du3x9Wr\nVxEYGIgRI4z9eDbs7t276Ny5M/r06YNp06aZ1YcDBw5g3LhxZp1LOUWdO4fWuYTwdQCHAKzOpQ1v\nrRYRkZGW7hqRxXDRBxtZuXw5+vj4wCs4GFdTUrAqJQUfABgG4AMAq7VaXE1JQYvgYPTx8cHK5ctV\n7nHB9+TJE/z9998AgI0bN6Jq1apIS0vDtWvXTA7h9PR0NG/eHFWrVsXmzZvN6o9Wq0VYWBjatWtn\n1vmUU1JiInJ7svsjgLYAqudyjDOAR88t1EFU0PCO2AZWLl+OeePH6436NMQNwKci6J2cjG7jxwMA\nRo8ZY5M+FjYXL17E0KFDkZiYiEaNGmHHjh344IMPsGzZMrPa69ixI5KSknDx4kWzZ4A6duwYmjRp\nkudR2fRijmXL4lEu+38EMPUFbTwC4Oxm7AkzkfoYxFYWHh6OGdlCOBXA/wLYByABgCeAuQC6P91f\nG8BvycloN348mnt7w8vLy+b9LqhEBKtXr8bHH3+M5ORkAJnTKu7duxedOnUyq82xY8fi2LFjOHPm\nTL5ClM+HLSc1NRVr165F8K5deIzMqlF2RwHcBjDwBW2FOzigYePGFu8jkaWwNG1li/z9MUmr1bsT\nTgdQDcBhAInIHPH5JoDo546pDWCiVotF/v4262tB9+DBAwwZMgTvvfdeVggDmaNsY2NjzWpz7dq1\nWLx4MX7++Wc0bNgwX/3bv38/gzifUlJSsGzZMtSpUwfvv/8+4hISEIzMX1izWwtgAIAyubQXD2Cr\nCEaMGmWF3hJZBl9fsqKYmBjUq14dV1NSjL568UxTADMB9HtuW3F49SKvI8iPHj2Kt956C9evX8/R\nxptvvokVK1bA1dXVpGtHRETAy8sL48aNw7x58/L1PR4+fIhKlSohNjYWDg4O+WqrOEpOTsZ///tf\nzJ8/H7dv39bb5wjgSwDmDIFbqCg41a8ffjTzuT+RTag1pVdx8NW8eeKn0bxwOr67gDgAcsnAPj8H\nBwmYP1/tr2JxeZ28/9ixYzJ79myxt7fPMWewo6OjrFq1yqyJ++Pi4qRMmTLSqVMni3yfHTt2SMeO\nHS3SVnHy6NEjmT9/fo7FNZBtaspyJUrkeXpLvWkuHR1fuEgHkdoYxFb0/rBhsvwFPyzSAOkMyBgj\n+5cBMnr4cLW/ikXldfL+rxVFXO3s9OYLfvanWbNmcuHCBbOun5GRITVq1JCqVavmuvCDKT777DOZ\nPXu2RdoqDhITE2XOnDni7u5uNIBLly4tH330kURHR5s113QtR0dZsWyZ2l+V6IUYxFY0tFcvWZfL\nDwsdIIMB6QlIupFjfgRkaK9ean8VizH0AzUekDcAKQNIDUB+fsHk/Z988omkpKSY3Yfu3buLo6Oj\nxMbGWux7NWvWTI4cOWKx9oqq+Ph4mTFjhri6uhoNYAcHB/nss8/k9u3beuc++wVuQS6/wMU9/QXu\nZYYwFSIMYit60R2xHyCdAHmSyzHLAGlWr55cuXJF7a+Tb2FhYfKygbuaIU//JAPyByAugPyVLYzL\nAuLq6io7d+7MVx8+//xzsbOzs2i58v79+1K2bFlJTU21WJtFTWxsrEydOlWcnZ2NBrCTk5NMmjRJ\n7t27Z7Sd8PBwebt/f3HVaMTPwUGWPf1ldRn+75HG2/37sxxNhQqD2Ipye0b8P4C0BuTxC0psg5/e\nDSqKIr1795YDBw4U2sXMh/XrJwsVRe/7PQakFCBRz20bAciUbH8PAYAM7NEjX9ffsmWLKIoi33//\nvYW+UaagoCDx9fW1aJtFxd27d2X8+PFSpkwZowHs4uIiX3zxhdy/fz/P7cbExEjA/PkyevhwGdqr\nl4wePlwC5s+XmJgYK34bIutgEFvRvXv3xFWjyVFGu/40XB0AcXr6xzlbSfZZmU1j4AdXs2bNJDAw\nMF/lWVsz9ndx+mlJ+vltXwPSx8DfhatGY/YP2osXL0qJEiWssk7wmDFjJCAgwOLtFmY3b96UsWPH\nikajMRrAbm5u8uWXX0pCQoLa3SVSFd8jtiIPDw/07NEDaxRFb3s1ADoAycic9ecRMtdMHZrt/NX4\nvzVXnxcREYFRo0ahRo0amD17ttnv0NrS2sBA9ANyvMaVBCD7FBplgRyzKT0/eb+pkpKS0KpVK3h5\neZk981ZuOJGHvuPHj6NWrVpYtGgRUlJScuyvUKEC/vOf/+D69ev44osvTH7tjKioYRBb2dgpUzDP\nwcHgEm25iQLwtaMjJs+ahebNmxs85u7du5g+fTqqVauG999/H+fPn893f63l8tmzaGngh7ITMn8J\neV4iYHB+YW+tFpdNnLxfp9PB29sbGo0Ghw8fNuncvLh16xZiY2PRtGlTi7ddWDVq1Mjg2r8vv/wy\nFixYgH/++QeTJk3i+sBETzGIrczb2xuzAgLQzdExz2EcBaCLRoNZAQGYPn06Tpw4gdDQUPTr1w+K\nkvMeOSUlBd9//z0aNWqErl27YteuXdDpdBb9HvllbPL+usicaezv57adAWBojitzJu8fPHgw/vnn\nH5w6dcqk5RDz6uDBg/Dx8TF7fuqiRKvV4ttvv0X9+vVRqVKlrO1VqlTB4sWLcfXqVXz66acoUya3\nubCIih/+9LCB0WPGYFJAANo5OmKhohicrg/InElrgaLgtdKl8VCjQb+BmbPoKoqC9u3bY8uWLYiK\nisLYsWPh5ORksI29e/fC19cXDRs2xIoVK/SmglSTk4uLwcn7HQH0BzAdmaX6PwBsBzDcwLGmTt4/\nf/58bNmyBbt379YLBks6cOCA2XNcFxWPHz/G119/DU9PTxw4cADBwcE4ceIE2rdvj++++w5RUVH4\n6KOPOOMYkTFqP6QuTkx59WLixInSvXt3oxNOPHjwQL7++mupXr260cEwAMTd3V2mTp0qt27dsvG3\n1ZfbCPLn3yOuDsh6IyPITZllbM+ePaIoinzzzTdW+046nU6qVatm9sQihV1iYqLMnTtXPDw8ZNCg\nQRIREaF2l4gKJc41rYLY2NjM+ZUjI/EoIQHObm6o27gxRowalTW/clpaGnx8fNC7d29MnjzZaFvp\n6ekIDg7GwoULcfToUaPHdezYEQcOHLD4d8krU+bdNsSUebejo6NRp04dDBo0COvWrTOrv3nx999/\no127drh165bBRwZF1YMHD/Dtt99i8eLF6Nq1K6ZNm4YGDRqo3S2iwkvt3wTIuOjoaHnppZfk0KFD\neTr+zz//lKFDhxqcl3nUqFGSmJho5R7nztB7xHn987WiyNv9+7/wGk+ePJHy5ctL48aNrf59Vq5c\nKcOGDbP6dQqK+/fvy7Rp08Td3V1GjRolly5dUrtLREUCnxEXYFWrVsXq1avx1ltv5ekVpZYtW+Ln\nn3/GtWvXMGnSpKzXQipVqoTHjx+jZs2a+PTTT/HPP/9Yu+sG5WcE+SwRVKlb94XHvvbaa9DpdAgL\nCzOrj6YoLq8txcTEYNKkSahbty5iYmIQFhaGH374AXXz8N+DiPJA7d8E6MUmTZqU6/NiY5KSkmTp\n0qWybt06ERG5fv26TJw4Udzd3aV///5y+PBhm8/SZc7k/c/mmlYURYKCgoy27efnJyVLlpSoqCir\nfw+dTiceHh7yzz//WP1aarl165Z88skn4ubmJh9++KFcv35d7S4RFUkM4kIgNTVV2rRpI3PnzrVI\ne48ePZIlS5ZInTp1xMvLS9atW2fTeZK/nDlTygLy1dOBWsYm7w9QFCkL/QUfSpUqJfv378/R5rJl\ny0RRlHzPRZ1X586dk1q1atnkWrYWHR0tH374obi5ucknn3yi+kA/oqKOpelCoGTJkvjll1+waNEi\ni0xK4eTkhA8//BAXL17E9OnTsWrVKtSsWRP+/v6Ij4+3QI9z9zApCQ+R+cpSJQBDACwHsO7pP/00\nGnhqNIjo1w8zFyyAnb191rn29vY5BkYdO3YMH330EWbOnAlfX1+r9x8ommXpf/75B6NHj0azZs1Q\npkwZXLhwAQsXLrTaq19E9JTavwlQ3u3cuVOqVKlilYntT58+LaNGjRI3Nzf54IMPrPZKjlarzbEG\nbemSJUXzdM5tDSCzZs7U+45r1qwRIHN5PEVRJDQ0NGvfnTt3RKPRSJ8+fazSX2P69u0rP//8s02v\naS2XL1+WUaNGSbly5WTatGkWXR6SiF6MQVzITJo0Sbp162axBe2zu3PnjkyfPl08PDykR48esmfP\nHos+R/7pp59yTPyffWk8Q6O7V65cKZcvX5aBAwdKyZIl5dy5c5KWliaVKlWSOnXqWO3vw5D09HRx\ndXWVO3fu2Oya1nD+/Hl56623pHz58jJr1iwuvkCkEgZxIWPp58XGaLVa+f7776Vhw4bSqFEj+f77\n70Wr1ea73Xbt2umF7ieffJJjhZ7Hjx+/sA0HBwdp1aqVODs72/y1rPDwcGnQoIFNr2lJZ86ckUGD\nBomHh4fMnTtX9dfaiIo7PiMuZEqWLIn169db7HmxMRqNBu+++y4iIyOxYMECbNmyBdWrV8eMGTNw\n7949s9rU6XTo0KEDSpYsmbVt9OjRSE9P1zvuRXNC//7779BoNPjzzz/x22+/oWzZ7Os3WVdhfT58\n8uRJvPHGG+jWrRtatWqFq1evYsqUKTb/+yMifQziQqhKlSr44Ycf8vx+cX4oioIuXbpg586dCA0N\nRUxMDOqY9UCIAAAgAElEQVTXr49Ro0bhzJkzJrVlZ2eHCRMmIC0tDdOnT8f48eNRv359k4P4559/\nxoMHD+Dm5oZBgwYhNTXV5O+VH4Vtfunjx4+jZ8+e6Nu3Lzp16oSrV69i3LhxXHyBqKBQ+5aczDd5\n8mSrPi825v79+zJ37lypXLmydOzYUbZt25bnPnz00UdSrly5rM9paWl6ZWk7O7tczz99+rTY29vL\nhAkTJDExUVxdXaVBgwY2+zt48uSJODs7S3x8vE2ulx+hoaHSuXNnqV69uixfvlxSUlLU7hIRGcAg\nLsTS0tKkTZs2MmfOHFWun5qaKj/99JN4eXlJnTp1ZMmSJfLo0aNcz3FxcZHx48dnfdZqtXpBXKpU\nKaPnJiQkSJkyZeT111/P2nbjxg1xcHCQDh065Pv75MWhQ4ekRYsWNrmWOXQ6nezbt0/at28vnp6e\n8v3338uTJ0/U7hYR5YJBXMjduHFDXnrpJb1XemxNp9PJ4cOHpX///uLu7i4TJ06U6OjoHMdt375d\n7Ozs9AZjPXr0SC+Iy5QpY/AaGRkZUqtWLalatWqOu9/IyEgpWbKkDBo0yLJfzICZM2fKhAkTrH4d\nU+l0OgkJCZHWrVtLvXr1ZO3atZKWlqZ2t4goDxjERUBISIhUrlxZ7t27p3ZX5OrVq1nTIg4ePFiO\nHz+eta9Zs2bSpk0bvePj4+P1gtjFxcVguz169BAHBwej77iGhoaKnZ2dfPzxx5b7Mga0b99edu3a\nZdVrmEKn00lwcLB4eXlJw4YNZf369ZKenq52t4jIBAziImLy5MnStWtXmz8vNiYxMVEWLFggNWrU\nkNatW8vq1atFURT5448/9I6LiYnJsX5ydl988YXY2dlJWFhYrtcMCgoSRVHE39/fot9FJPM95tmz\nZ4tGo5EHDx5YvH1TZWRkSFBQkDRt2lSaNWsmmzdvLjD/7YnINAziIiItLU3atm2r2vNiY9LT02Xz\n5s3y0ksviZ2dnQQEBOhNHHH79m29IH755Zf1zg8ODhZFUeT777/P0/UWL14siqJIYGCgRb9HkyZN\n9O7ajx07ZtH28yo9PV1+/vlnadCggXh7e8v27dttvnAHEVkWg7gIKQjPi41xcnISPz8/eeutt8TN\nzU3+/e9/S1RUlFy/fl0viKtUqZJ1zsWLF6VEiRLywQcfmHStqVOnip2dnYSEhFik79nv2hVFsfmo\n6dTUVAkMDJS6detKmzZtZPfu3QxgoiKCQVzE7Nq1q8A8L34mKChI7O3ts2bmunnzpkyZMkXKly8v\nnTt31gu5GjVqiIjI48ePxdXVVVq2bGnWNf38/MTe3l7Cw8Pz3f8NGzbo9dGWo6afPHkiK1eulJo1\na4qPj48cOHCAAUxUxHBCjyKme/fuGDFiBIYPHw6dTqd2dwAAs2bNQrt27aDRaAAAlStXxty5c3H9\n+nW0bdtW79hnk3l4e3ujVKlSZs8etnr1anTt2hVt27bF33//na/+HzhwQO+zLWbVSklJwbJly1Cn\nTh1s2rQJa9euxcGDB9GxY8ccq08RUeHGIC6CvvzySyQnJ8Pf31/truDu3bs4f/485s+fn2Ofo6Mj\nBg4cqLft2rVraNy4MaKionDy5EmUKlXK7Gvv2LEDjRo1QrNmzXD//n2z27FlECcnJ+Obb76Bp6cn\ndu3ahY0bN+K3337L8QsLERUdDOIiqESJEvjll1+wePFihIaGqtqXSZMm4aWXXoK3t7fB/dmnt3R1\ndcW5c+dQqlQpfPnll/jrr7/MvradnR2OHz8ODw8PNGjQAMnJySa3cePGDVy5ciXrc4kSJdCuXTuz\n+2RMUlISvvrqK9SqVQuHDh3C9u3bsX37drRq1cri1yKigoVBXERVqVIFgYGBGDZsGGJiYlTrx6ZN\nmzBmzBij+7MH8f3797Fw4UL8/fffqFKlCjp16oRu3bph9+7dEBGTr1+iRAlERkYCABo3bpzjei9y\n8OBBvc+vvvqqRedoTkxMxJw5c1CrVi2cPHkSe/fuxZYtW9C8eXOLXYOICjYGcRH27Hnx22+/rcrz\n4nXr1iE1NRWTJ082ekxaWpreZ3d3d3zyySfw8PDA9OnTce3aNQwdOhSTJk1Cw4YNsXLlSmi1WpP6\n4ejoiL/++guxsbF49dVXTfq72L9/v95nS5WlExISMHPmTNSuXRsXL15EaGgo1q9fj8aNG1ukfSIq\nPBjERdyXX36JlJQUVZ4Xz5kzBx07dsz1OW/2O9T69evrfS5dujRGjRqFiIgILFmyBDt27ED16tXx\n+eef486dO3nuS/ny5XH69GmcO3cOvXr1AgDExcXleo6IWPz58P379zF16lTUrl0bN27cwLFjx/Dj\njz/iX//6V77aJaLCi0FcxD17XrxkyRKbPi++efMmLl26hICAAKPH7NmzB0+ePNHbZmwJREVR8Prr\nr2Pbtm04cuQIHjx4gIYNG2L48OE4depUnvrk6emJP/74A3v27EGvXr1Qr149fPfdd0aPj4qKws2b\nN7M+azQavPrqq3m6VnZ3797FhAkTULduXcTHx+PkyZNYtWoVateubVZ7RFSEqP3+FNnG7t27pXLl\nynL37l2bXG/IkCFSuXJlo/s3bdokAKR8+fJ67+h26tQpz9eIj4+XefPmSZUqVaR9+/aydevWPM2z\nPHv2bL3JOYKCggwe99133+n1rXPnznnu2zM3b96Ujz/+OGsSkxs3bpjcBhEVbQziYmTq1KnSuXNn\nqy8KkJGRIRqNRubNm2dw/6lTp8TR0VEv5J796d69u8nXS01NlV9++UVatmwpnp6esmjRInn48KHB\nYxMSEsTFxSXH0ov79u0TEZF79+7JV/PmyfvDhkndihXFARDl6XGmzGF97do1GTNmjLi5uclnn30m\nt2/fNvl7EVHxwCAuRtLS0qRdu3Yye/Zsq17nv//9r5QsWdLgMnx3796VqlWrGgxhANKrV698Xfvo\n0aMyaNAgKVeunIwbN06uXbuW45gdO3aIvb293nUdHBzE18dHXDUaeUejkeWArANkOSCDAdEA0qND\nhxcuPBEVFSXvvvuulCtXTiZPniwxMTH5+j5EVPQxiIuZmzdvyssvvywHDx602jVq164tvr6+Oban\npKRI69atjYYwAHnjjTcs0odr167JuHHjpFy5cjJw4EA5cuSI3tSQa9as+b/yNCBlAQkAJB4QMfAn\nHpCvFUVednSUFcuW5bjexYsXZcSIEeLu7i5ffPGFxMXFWeR7EFHRxyAuhqz5vPjq1auiKIr89ddf\nett1Op2MHDnSYCn6+c+DBg2yaH8ePnwoixYtklq1aknLli3ll19+kdTUVBERCQgIEAWQlwG5YiSA\ns/+5Akit58I4MjJShgwZIhUqVJDZs2frrSxFRJQXDOJiatq0aVZ5Xty/f3+pXr16ju0BAQE5QrhL\nly56d6YAZOjQoRbtzzPp6ekSHBwsHTp0kCpVqsh//vMf2bdvn7iXKGE0hC8/LUkPNxDGHhqNdOzY\nUV566SWZN2+e0WfSREQvwteXiqmZM2ciNTXVou8X63Q67NixA59++qne9pCQEEyYMEFvW506dbBh\nw4Ycs2UZe30pv+zt7dG3b1/8/vvv+PXXX3H+/HkM6NEDU9PTYewFoo8AtDSwvTaACSkpSL5/H3//\n/TcmTpwIZ2dnq/SbiIo+BnEx9ez94qVLl+aYxtFcy5YtAwD8+9//ztp24cIFDB06VC9wXVxcsH37\ndri5ueWY0KNkyZIW6UtumjdvjoCAACh2dvAzcsx6AG4AOhnZ/w6AS1eumDV/NRHR8xjExVilSpWw\nZs0avP3227h79y5++OGHfAVLQEAAfH19YWeX+b9VXFwcevfujYcPH2YdY2dnhw0bNqBevXoAck5x\naa074uzWBgain6LAzcC+hwBmAFiAzHq5IeUA9FMUrA0MtFIPiai4YBAXc127dsWQIUPQpEkTvPPO\nO/j444/NaufSpUuIjo7OmkkrLS0Nb775Zo61gL/++mt069Yt63OrVq0wd+5cODg4YNy4cfD19TX/\ny5jg8tmzaJmSYnDfdADvA6j0gja8tVpcfrqgBBGRuWxz+0EF1sWLF7F582bExsYCAFatWoUOHTpg\n+PDheTr/3r17WLp0KY4ePYqaNWvC09MTAPDpp5/mmKf5nXfewdixY/W2vfLKK3jllVcwb948fP75\n53B1dbXAt3qxpMREGHqqGwFg39N/vogzgEcJCRbtFxEVPwziYq5KlSrQaDR62z744AN4eXnlaSGC\n1atXY/bs2QAylxkMDQ3FhQsXsHTpUr3j2rZti2XLlkFRFIPtpKWl2awsDQBOLi54ZGB7KIDrAKoh\nsyydBCADwF8ATmQ79hEAZzdDxW0iorxjabqYc3JywsaNG/XCODk5GYMGDXrh82KdToeVK1dmfY6M\njMTGjRv1BmsBQLVq1bB582aULl3aaFvp6ek2Gaj1TN0mTRCW7RcQAPgfAH8j8474DIAPAPQCsMdA\nG+EODqjLZQuJKJ8YxIQmTZpg8eLFetvOnz+fI1Cz27NnD65du5b1uVSpUvjll1/0RkKXKVMG27Zt\ng4eHR65tpaen2/SOeMSoUdgKIHthWQPA47k/Tk+3lct2XDyArSIYMWqUlXtKREUdg5gAAO+++y6G\nDRumt2316tVYu3at0XNWrFih99nBwQEJ2Z6Z/vjjj2jatGmu19bpdNDpdFmjrW3Bw8MDPXv0wBoj\npfJnZgAw9DewRlHQy9cXFSpUsEr/iKj4UCT7jApUbCUlJcHLywuXLl3K2ubo6IgTJ07keF588+ZN\n1KhRAxkZGUbbmz17Nj7//PMXXjc1NRVOTk5ITU01v/NmCA8PRx8fHxxOTjY6qYchUQDaOTpie2go\nvLy8rNU9IiomGMSkJzIyEi1btkTKc6/2NGzYEDt27MCmjRtx+exZJCUm4ur16wiLjDT6nu2QIUPw\n888/Gx2c9bzk5GSUL19elckxVi5fjnnjx+O3PIZxFIBujo6YFBCA0WPGWLt7RFQMMIgph1WrVuG9\n997L+uwIQLG3x+CSJeGdkpL52g6A3wH8CkABoH3u/BYtWuDQoUNwdHTM0/USExNRtWpVvYk/bGnl\n8uWY9sknmJSaincBg5N8xAMIVBR85eCAWQxhIrIgBjHlICIYMWIEflq3Ds7InODiHRgOqAQAqwDM\nRmY4v1yxIsLDw1G5cuU8Xy8uLg5169ZFXFycBXpvnlatWiEyLAwCoC+ADkDWLxzhDg7YKoJevr4Y\nO2UKy9FEZFF8j5hyUBQF3i1aYP9PP+GQSFbJ1gfAnwBKIvMd2yoALgAYD+ANAO0AjBw+3KQQBmw/\nYjq7+/fv4+TJk3j2tHsjAPTuDTsROLu5oWHjxpg3ahQHZhGRVTCIKYfw8HD4T5uGw8+FMJBZgl4G\nGFwooTaAwwDaLVmCAYMGmXTXaOvJPLILDg7WG3TWoGFDrN+2TbX+EFHxwteXKIdF/v6YpNUaHLyU\n23OM2gAmarVYZOLSiraezCO7TZs26X0eNGiQSj0houKIQUx6YmJisHPXLow0MnRgCjInumiHzOkg\nsxspgh0hIVlzV+eFmqXp+Ph47N+/X28bg5iIbIlBTHrWBgaiHwwPzJoP4CqAW8hcnag3gH+yHWPO\n8oBqlqaDg4P1ZgJr0KABGjRooEpfiKh4YhCTntyWB/QGUAaZg7VGAGgDIMTQcSYuD6hmaTooKEjv\n88CBA1XpBxEVXwxi0mNseUBDFBh+Zmzq8oBqlabj4+Oxb98+vW0sSxORrTGISY+x5QETkbkC0RNk\nLgv4EzJHSXc3cKypywOqVZr+9ddf9crS9evXR8OGDW3eDyIq3hjEpMfY8oBpAD5H5kCtCgCWInNW\nLUMjq01dHlCt0nT2svSgQYPyNCUnEZElMYhJj7HlAcsDCEPmnXE8gKMAXjdwvjnLA6pxR5yQkMCy\nNBEVCAxi0pPX5QGNMWd5QDWeEf/6669IS0vL+lyvXj00atTIpn0gIgIYxGTA2ClTMM/BAVEmnhcF\nYL6DA8ZOmWLSeWqUplmWJqKCgkFMOXh7e2NWQAC6OTrmOYyfLQ84KyDA5EURbF2afvDgAfbu3au3\njWVpIlILg5gMGj1mDCYFBKCdoyMWKkqOZ8bPxAP4WlHQLh9r9Nq6NL1t2za9snTdunXR2ITBZURE\nlsQgJqNGjxmD7aGhONWvH2ppNHjHwQHLAawDsBzAEACVARzp2hXbQ0PNXqPX1qVplqWJqCDh6kuU\nKy8vL/y4eTNiY2OxNjAQEZGR+G3nTtyLj8cTZE7oMfx//idfa/TasjSdmJiIPXv26G1jWZqI1MQ7\nYsqTChUqYNyECVixdi36vv02UvB/s2qdPn06X23b8o5427ZtSE1Nzfpcu3ZtNGnSxCbXJiIyhEFM\nJnvllVf0PkdEROSrPVs+I2ZZmogKGgYxmaxZs2Z6n/MbxLYqTScmJuK3337T28ayNBGpjUFMJmvQ\noIFeKfnGjRuIi4szuz1blabj4+PRvXv3rDtgT0/PHL9UEBHZGoOYTFaqVKkca/aeOXPG7PZsVZqu\nWbMmlixZAldXVwQGBmLWrFksSxOR6hjEZJbsz4nzM2DLlqOmN23ahDfeeAMjR47EsGHDbHJNIqLc\nMIjJLJZ8TmzLUdNBQUF8LkxEBQqDmMxi6SC2xR3xjRs3cOnSJXTq1Mnq1yIiyisGMZmladOmep8v\nXLiAlJQUs9qyVWl68+bN6NOnD0qVKmX1axER5RWDmMzi6uqKmjVrZn3OyMjAuXPnzGrLVqVplqWJ\nqCBiEJPZLFWetkVp+ubNm7hw4QI6d+5s1esQEZmKQUxms1QQ26I0vWXLFpaliahAYhCT2Sx5R2zt\n0jTL0kRUUDGIyWzZ3yU+c+YMdDqdye1YuzR9+/ZtnD9/Hl26dLHaNYiIzMUgJrNVqVIF5cqVAwAo\nioKmTZsiISHB5HasXZrevHkzevfuzbI0ERVIDGIym6Io2Lx5M+bMmYPSpUvjjz/+gLu7u8ntWLs0\nHRQUhIEDB1qtfSKi/GAQU774+Phg9OjRSElJQVJSklltWLM0fefOHURGRqJr165WaZ+IKL8YxJRv\n5cuXR6lSpRASEmLW+dYsTW/evBm9evVC6dKlrdI+EVF+MYjJIipWrIi9e/eada41S9McLU1EBR2D\nmCyiQYMGOHHihFnnWqs0fefOHZw9e5ZlaSIq0BjEZBHt27fHP//8Y9a51ipNb9myBT179oRGo7F4\n20RElsIgJovo3bs3EhMTzX6P2BqlaZaliagwYBCTRTRs2BB2dnY4evSoyedaozR99+5dnDlzBt26\ndbNou0RElsYgJotxd3fH9u3bTT7PGqXpLVu2wNfXl2VpIirwGMRkMXXr1sXx48dNPs8apelNmzax\nLE1EhQKDmCymVatWuHz5ssnnWfqOOCYmBqdOnWJZmogKBQYxWYyvry9iY2NNPs/Sz4iflaUdHBws\n1iYRkbUwiMliOnTogIyMDFy5csWk8yxdmuZoaSIqTBjEZDElSpRA2bJlsW3bNpPOs2RpOiYmBidP\nnkT37t0t0h4RkbUxiMmiqlevjkOHDpl0jiVL01u3bkX37t1ZliaiQoNBTBbVvHlznDt3zqRzLFma\nZlmaiAobBjFZVOfOnXHr1i2TzrFUaTo2Nhbh4eHo0aNHvtsiIrIVBjFZVK9evfDkyRM8ePAgz+dY\n6o74WVna0dEx320REdkKg5gsytXVFRqNBjt27MjzOZZ6RsyyNBEVRgxisrhKlSph3759eT7eEqXp\n+/fvIywsDL6+vvlqh4jI1hjEZHGNGjXCqVOn8ny8JUrTW7duRbdu3ViWJqJCh0FMFte+fXtcu3Yt\nz8dbojTNsjQRFVaKiIjanaCi5e+//0bt2rXzXHK2s7NDWloa7O3tzbpeXFwcatWqhdu3b6NMmTJm\ntUFEpBbeEZPFeXp6wt7eHqGhoS88VqfTAYDZIQwAwcHB6Nq1K0OYiAolBjFZRYUKFRASEvLC41iW\nJqLijkFMVpHXtYnzO2I6Pj4ex44dQ8+ePc1ug4hITQxisorXXnstT6sw5XfEdHBwMLp06cKyNBEV\nWgxisoqePXsiLi7uhcfltzQdFBSEgQMHmn0+EZHaGMRkFa+99hp0Oh3Onz+f63H5KU3Hx8fj6NGj\n6NWrl1nnExEVBAxisgo7Ozu4uLi8cG3i/JSmf/31V3Tq1AlOTk5mnU9EVBAwiMlqatWqhcOHD+d6\nTH5K0xwtTURFAYOYrKZFixb466+/cj3G3NJ0QkIC/vjjD5aliajQYxCT1XTp0gV37tzJ9RhzS9PP\nytLOzs7mdo+IqEBgEJPV+Pr6IjU1Fffv3zd6jLmlaZaliaioYBCT1Tg5OcHBwQG//vqr0WPMKU0/\nePAAhw8fRu/evfPbRSIi1TGIyaqqVKmCAwcOGN1vTml627ZteP3111mWJqIigUFMVtW4cWNEREQY\n3W9OaZplaSIqShjEZFU+Pj64fv260f2mlqYTExNx6NAhlqWJqMhgEJNV9e3bF48fP0ZqaqrB/aaW\nprdt2wYfHx+ULVvWUl0kIlIVg5isqlq1aihRogT27dtncL+ppWmWpYmoqGEQk9V5eHhg9+7dBveZ\nUppOTEzE77//zrI0ERUpDGKyuvr16yMsLMzgPlNK09u3b0eHDh3g4uJiye4REamKQUxW99prryEq\nKsrgPlPuiFmWJqKiiEFMVtenTx/Ex8dDp9Pl2JfXZ8QPHz7EwYMH0adPH2t0kYhINQxisroWLVoA\nAM6ePZtjX15L09u3b0f79u3h6upq8f4REamJQUxWZ2dnB1dXV4NrE+e1NM2yNBEVVQxisolatWrh\njz/+yLE9L6Xphw8f4sCBA+jbt6+1ukdEpBoGMdlEy5YtceHChRzb81Ka3rlzJ9q1a8eyNBEVSQxi\nsolu3brh3r17ObbnpTTNsjQRFWUMYrKJbt26IS0tDbdv39bb/qI74qSkJOzfv59laSIqshjEZBMa\njQYODg5YtGiR3vrEL3pGvGPHDrRp0wZubm626CYRkc2Ztv4ckRni4uLQvn17aLVazJ8/H2XLlsWD\nBw+gKMoLS9MsSxNRUcc7YrK6cuXK4datW1mfHz58iGvXrgHIvTSdlJSEffv2sSxNREUag5isTlEU\nNGvWTG/b6dOnAeRemt65cydat26NcuXKWb2PRERqYRCTTWQP4oiICAC5j5pmWZqIigMGMdmEsSA2\nVpp+/Pgx9u7dizfeeMMm/SMiUguDmGwityA2dEe8c+dOvPrqq3B3d7dJ/4iI1MIgJpto0KCB3p3v\njRs3EBcXZ7Q0zbI0ERUXDGKyiVKlSqFhw4Z62yIiIgyWph8/fow9e/awLE1ExQKDmGzGUHnaUGk6\nJCQErVq1Qvny5W3ZPSIiVXBCD7IZY8+Jswfxpk2bWJYmomKDd8RkM6+88ore59OnT+coTScnJ+O3\n335Dv379bN09IiJVMIjJZpo2bar3+eLFi3jy5IneHfGuXbvg7e3NsjQRFRsMYrIZFxcX1KxZM+tz\nRkYG4uLi9IKYo6WJqLhhEJNNZX9OnJCQkFWa1mq12L17N8vSRFSsMIjJprI/J37w4EHWHfGuXbvg\n5eWFChUqqNE1IiJVMIjJprLfET98+DAriIOCgjBw4EA1ukVEpBpFRETtTlDxcePGDVSrVg0AoADQ\nAGj76quoUqMG1m/dihOnTqFBgwaq9pGIyJZ4R0w2defOHbiULAkNgCEAFgAYefw4Wq5fj35paWjT\nogXe7t8f4eHhKveUiMg2eEdMNrNy+XLMGD8eE5KT4QfAzcAxCQACFQXzHRwwKyAAo8eMsXEviYhs\ni0FMNrFy+XLMGz8evyUno3Yejo8C0M3REZMYxkRUxDGIyerCw8PRx8cHhw2E8HoAXwKIBlARQCCA\nNk/3RQFo5+iI7aGh8PLysll/iYhsic+IyeoW+ftjklabI4T3ApgCYA2AJACHANR6bn9tABO1Wizy\n97dNR4mIVMA7YrKqmJgY1KteHVdTUnI8E24D4D0AfrmcHw/AU6PB5ehovl9MREUS74jJqtYGBqIf\ncg7M0gE4ASAGQB0A1QD8G8CTbMeVA9BPUbA2MNDKPSUiUgeDmKzq8tmzaJmSkmP7PQBpADYDOAIg\nAsBpAP/PQBveWi0uR0Zas5tERKphEJNVJSUmwtnAdoen//wYgAcy73w/AxBi4FhnAI8SEqzTQSIi\nlTGIyaqcXFzwyMB2VwBVsm1TjLTxCICzm6G3jomICj8GMVlV3SZNEKbRGNznB2AxgFhkTuSxEEBv\nA8eFOzigbuPGVusjEZGaOGqarCq3UdPpAMYC+BmZperBAOYBKPXcMRw1TURFHe+Iyao8PDzQs0cP\nrFFyFp5LAFiKzLvh28i8Iy6V7Zg1ioJevr4MYSIqsnhHTFaX28xaueHMWkRUHPCOmKzO29sbswIC\n0M3REVF5POfZXNOzAgIYwkRUpDGIySZGjxmDSQEBaOfoiIWKAmMvI8UDWKAoaMcFH4iomGBpmmzq\nxIkTWOTvjx0hIeinKPDWajPfE0bm6OitIujl64uxU6bwTpiIigUGMakiNjYWawMDcTkyEo8SEuDs\n5oa6jRtjxKhRHJhFRMUKg5iIiEhFfEZMRESkIgYxERGRihjEREREKmIQExERqYhBTEREpCIGMRER\nkYoYxERERCpiEBMREamIQUxERKQiBjEREZGKGMREREQqYhATERGpiEFMRESkIgYxERGRihjERERE\nKmIQExERqYhBTEREpCIGMRERkYoYxERERCpiEBMREamIQUxERKQiBjEREZGKGMREREQqYhATERGp\niEFMRESkIgYxERGRihjEREREKmIQExERqYhBTEREpCIGMRERkYoYxERERCpiEBMREamIQUxERKQi\nBjEREZGKGMREREQqYhATERGpiEFMRESkIgYxERGRihjEREREKmIQExERqYhBTEREpCIGMRERkYoY\nxAkdAcsAAAAKSURBVERERCr6/y2v8L9Dh1+LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x102302a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx.draw(graph, with_labels=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we aggregate the all the transitions into *incoming traffic* and *outgoing traffic*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "incoming traffic: [ 465.  682.  155.   83.  369.  593.  148.  505.]\n",
      "outgoing traffic: [ 397.  367.  365.  359.  384.  399.  365.  364.]\n"
     ]
    }
   ],
   "source": [
    "traffic_in = transitions.sum(axis=0)\n",
    "traffic_out = transitions.sum(axis=1)\n",
    "\n",
    "print(\"incoming traffic:\", traffic_in)\n",
    "print(\"outgoing traffic:\", traffic_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Estimating transitions using ChoiceRank\n",
    "\n",
    "ChoiceRank can be used to recover the transitions on the network based only on:\n",
    "\n",
    "1. information about the structure of the network, and\n",
    "2. the (marginal) incoming and outgoing traffic at each node.\n",
    "\n",
    "ChoiceRank works under the assumption that each node has a latent \"preference\" score, and that transitions follow Luce's choice model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params = choix.choicerank(graph, traffic_in, traffic_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can attempt to reconstruct the transition matrix using the marginal traffic data and the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "est = np.zeros((n_items, n_items))\n",
    "for src in range(n_items):\n",
    "    neighbors = graph.successors(src)\n",
    "    if len(neighbors) == 0:\n",
    "        continue\n",
    "    probs = choix.probabilities(neighbors, params)\n",
    "    est[src,neighbors] = traffic_out[src] * probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True transition matrix:\n",
      "[[   0.  233.    0.    0.   94.   70.    0.    0.]\n",
      " [ 148.    0.    0.    0.  143.   76.    0.    0.]\n",
      " [   0.    0.    0.    0.    0.  365.    0.    0.]\n",
      " [   0.    0.    0.    0.    0.    0.    0.  359.]\n",
      " [  83.    0.  155.    0.    0.    0.    0.  146.]\n",
      " [  83.  168.    0.    0.    0.    0.  148.    0.]\n",
      " [ 151.    0.    0.    0.  132.   82.    0.    0.]\n",
      " [   0.  281.    0.   83.    0.    0.    0.    0.]]\n",
      "\n",
      "Estimated transition matrix:\n",
      "[[   0.  232.    0.    0.  102.   63.    0.    0.]\n",
      " [ 150.    0.    0.    0.  134.   83.    0.    0.]\n",
      " [   0.    0.    0.    0.    0.  365.    0.    0.]\n",
      " [   0.    0.    0.    0.    0.    0.    0.  359.]\n",
      " [  84.    0.  155.    0.    0.    0.    0.  146.]\n",
      " [  83.  169.    0.    0.    0.    0.  148.    0.]\n",
      " [ 149.    0.    0.    0.  133.   83.    0.    0.]\n",
      " [   0.  280.    0.   84.    0.    0.    0.    0.]]\n",
      "\n",
      "Difference:\n",
      "[[ 0.  1.  0.  0. -8.  7.  0.  0.]\n",
      " [-2.  0.  0.  0.  9. -7.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [-1.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [ 0. -1.  0.  0.  0.  0.  0.  0.]\n",
      " [ 2.  0.  0.  0. -1. -1.  0.  0.]\n",
      " [ 0.  1.  0. -1.  0.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(\"True transition matrix:\")\n",
    "print(transitions)\n",
    "\n",
    "print(\"\\nEstimated transition matrix:\")\n",
    "print(np.round_(est))\n",
    "\n",
    "print(\"\\nDifference:\")\n",
    "print(np.round_(transitions - est))"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}