Generalized degree distribution of the Random Simplicial Complex (RSC) model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import numpy as np\n",
    "from scipy import random\n",
    "import pandas as pd\n",
    "import copy\n",
    "import random\n",
    "from collections import OrderedDict, Counter\n",
    "from multiprocessing import Pool\n",
    "import itertools\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_my_simplicial_complex_d2(N,p1,p2):\n",
    "    \n",
    "    \"\"\"Our model\"\"\"\n",
    "    \n",
    "    #I first generate a standard ER graph with edges connected with probability p1\n",
    "    G = nx.fast_gnp_random_graph(N, p1, seed=None)\n",
    "    \n",
    "    if not nx.is_connected(G):\n",
    "        giant = list(nx.connected_components(G))[0]\n",
    "        G = nx.subgraph(G, giant)\n",
    "        print('not connected, but GC has order %i ans size %i'%(len(giant), G.size())) \n",
    "\n",
    "    triangles_list = []\n",
    "    G_copy = G.copy()\n",
    "    \n",
    "    #Now I run over all the possible combinations of three elements:\n",
    "    for tri in combinations(list(G.nodes()),3):\n",
    "        #And I create the triangle with probability p2\n",
    "        if random.random() <= p2:\n",
    "            #I close the triangle.\n",
    "            triangles_list.append(tri)\n",
    "            \n",
    "            #Now I also need to add the new links to the graph created by the triangle\n",
    "            G_copy.add_edge(tri[0], tri[1])\n",
    "            G_copy.add_edge(tri[1], tri[2])\n",
    "            G_copy.add_edge(tri[0], tri[2])\n",
    "            \n",
    "    G = G_copy\n",
    "             \n",
    "    #Creating a dictionary of neighbors\n",
    "    node_neighbors_dict = {}\n",
    "    for n in list(G.nodes()):\n",
    "        node_neighbors_dict[n] = G[n].keys()           \n",
    "                \n",
    "    #print len(triangles_list), 'triangles created. Size now is', G.size()\n",
    "        \n",
    "    #avg_n_triangles = 3.*len(triangles_list)/G.order()\n",
    "    \n",
    "    #return node_neighbors_dict, node_triangles_dict, avg_n_triangles\n",
    "    #return node_neighbors_dict, triangles_list, avg_n_triangles\n",
    "    return node_neighbors_dict, triangles_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_p1_and_p2_correction(k1,k2,N):\n",
    "    p2 = (2.*k2)/((N-1.)*(N-2.))\n",
    "    p1 = (k1 - 2.*k2)/((N-1.)- 2.*k2)\n",
    "    if (p1>=0) and (p2>=0):\n",
    "        return p1, p2\n",
    "    else:\n",
    "        raise ValueError('Negative probability!')\n",
    "        \n",
    "def one_realization(ii):\n",
    "    print('ITERATION %i'%ii)\n",
    "    \n",
    "    #Generating a simplicial complex with the given parameters\n",
    "    node_neighbors_dict, triangles_list = generate_my_simplicial_complex_d2(N,p1,p2)\n",
    "\n",
    "    return node_neighbors_dict, triangles_list\n",
    "    \n",
    "def parse_results(results):\n",
    "    \n",
    "    degree_list = []\n",
    "    tri_degree_list = []\n",
    "    \n",
    "    for res in results:\n",
    "        node_neighbors_dict, triangles_list = res\n",
    "        \n",
    "        #Adding the degree sequence of this realization to the global list\n",
    "        degree_list += [len(neighs) for neighs in node_neighbors_dict.itervalues()]\n",
    "    \n",
    "        #creating a triangles degree sequence and adding it to the global list\n",
    "        vertex_in_a_triangle = []\n",
    "        for tri in triangles_list:\n",
    "            vertex_in_a_triangle += tri\n",
    "            \n",
    "        tri_degree_list += Counter(vertex_in_a_triangle).values()\n",
    "\n",
    "    degree_counter = Counter(degree_list)\n",
    "    triangle_degree_counter = Counter(tri_degree_list)\n",
    "\n",
    "    return degree_counter, triangle_degree_counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running in parallel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Simplicial Complex\n",
    "N = 2000\n",
    "k1 = 20\n",
    "k2 = 6\n",
    "p1, p2 = get_p1_and_p2_correction(k1,k2,N)\n",
    "\n",
    "n_simulations = 80\n",
    "n_processes = 8\n",
    "\n",
    "#################################################\n",
    "\n",
    "iteration_numbers = range(n_simulations)\n",
    "\n",
    "pool = Pool(n_processes)                         \n",
    "results = pool.map(one_realization, iteration_numbers)\n",
    "degree_counter, triangle_degree_counter = parse_results(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAADQCAYAAACa9N1LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+UHNV15z9XILCwhZrhpyUsoRnAxo6cMCMl/MghiTVD\nNnEcMJ5B6032JD4bTccxSTYLmUEntsEsZ+WRYe14Q8gMTrJx1t6IbntNHNuxp8E2OQQcNOPYAmRj\ndwtkJCQBmgZhBELS2z9e1UxNddVMV3V1V/+4n3Pq9Myreq/e666+/X58371ijEFRFKUTWJJ2BRRF\nURqFGjxFUToGNXiKonQMavAURekY1OApitIxqMFTFKVjUIOnKErHkJrBE5FBEZkSESMikyLSXUWe\nEU+eoogMJ1GuoiidQSoGT0QGgRxQAIac5KKIZBbIMwaMOXkGgDww7qTHLldRlM5B0thpISJFYNoY\nM+RJmwEmjDGjIXkMsM173jF2I8YYiVuuoiidw8mNvqHT2+oGsr5T9wKDQIVhcoalJWC771SxxnKX\nA8ux78NbgKeBE9W3RlGUBrMEOBt4zBjzetTMDTd4WKME1oB5KQLXB2UwxpSAnoBTWWA6brnAjcAt\nC1VWUZSmpBf4btRMaRi8rpD0MlDVXJuI9AL3YI1cXw3l3glMAG8GdkxNTXHeeedVU4W6ceTv7Gh8\n2e/kUq1Hq/OBzzwEwN/+3pUp10RJkv3799PX1wfwXJz8aRi8QyHpixo7Z9h6D3aImgc2GmPKccs1\nxhwGDosIAOeddx4rV65crBp1pbTp7wFSr0erc/dv299BfR/bllhTT2ms0roGyi8X6aFyODqLM4+3\nG9uV7TPGDHmMXexyFUXpHBpu8Jz5uBJzshGXfmyvLYxJYIcxpscYM+0/WUO5TcWKxz/Lisc/m3Y1\nWp7t33uB7d97Ie1qKE1GGkNasCumOUdGUgC2YHtm4+4Fjqh4wBgz5MzZdQN5R2s3D2OMa9AWLbfZ\nOe2ZB9OuQlvwyJ6X066C0oSkYvCMMXkRyWJXWV0xcZ/TS3Ppw87VAax3XkdCipQI5SqKkhClUome\nnh5axXN6alvLjDETxpg+Y4wYYwb8w1RjTNYVFDvXStgRpVxFUZKjUCjQ29ubdjWqRp0HKIoSm6mp\nKfr7+2f/L5fLlErNO6BKaw5PUZQGceXHH2Bv+UhF+qrMMh66+V01lV0oFBgbs9vZ8/k8mzdvJpfL\n0d3dnD471OA1Gc/+6meASm2NEo1PvmdN2lVoGvaWj/DUx99dkX7BzV+puexSqURvby+jo6Pk83mm\npqaa1tiBGjxFaSs+9uXHeWLfS1Vfv2n84dm/377ydG55zzuqzjs9bafHs9kshUKBYrE4z9jl83ky\nmcy8IW/a6Bxek6E6vGRQHV79KRQKAAwMDNDb20s+Pyd3zWazDA0NUS6Xw7KngvbwmgzV4SVDp+rw\ngnpoCw1dt2cvj32vRx99lJGREUZGRshkMoyOjjIyYpVjY2NjZDLN54ZSe3iKosSiUCgwMDAAwPDw\nMOVyebaXl8lkOPPMM9OsXiDaw1OUNmdVZllgL29VZlnsMsvlMuVymfXr18+mDQ8PMz4+zuBgxWao\npkENnqK0ObVKT4LIZDIVuyvGx5t/B6cOaRVF6Ri0h9dkqA4vGVSHlz7uAkYzoT08RVE6BjV4TYbq\n8JJBdXhKEDqkbTJUh5cMnarDUxZGe3iKonQMavAURekY1OApitIxqMFrJnbmOfWl3bzhwDTcdRns\nbJnYQ4rSEuiiRRrszMODd8DzP4Sz3gpX3WTT77+NJdf/Hay+HPY8DPfdYNPXNe9WnWZFdXiNodVi\nWqjBazQ783D/bXDNX8wZti/9AZjj8N5xWHuVvW7tVfaar46owVOaFo1poSzMg3dYQ7b2KjhpqX29\n9i/hpWetAfSy+nLbC1Qiozq8xtBqMS3U4DWa538YbNhOeaPt7XnZ87Ad8iqReWTPy6rF87Izb+eF\nP3ZGovPDhUKBDRs2ANbD8dq1a5va4OmQttGc9VZryNyhK9j/T11h5+y8Q937boCNH02vrkp7EDSN\nktD8sMa0aAAikgEwxjSX/+hquOomO2d37V/Of/iuvg2AE/f+DnJkBjnlNHjXR3T+TonG126G/Tvn\npz23C4b+d+X8cO53Ycffzl133jr4tY9XfauFYlpks1l27NhBuVwml8s1zTxfakNaERkUkSkRMSIy\nKSJRfhZywHBImSbgGEuu5jWybhDOuhj+73+C28+xixIbP2rT1w3y2ulree2sd8DRn8Lxo2nXVmkH\nXjkUPI3yyqGaig2LaeGmT01Nkcvl2Lx5c033SZJUengiMog1WtuArUAWKIrIGQv12hyjmAX6gcmA\nS7qBMjDqS9+RRL0T4cQJOPgEdP8S/MfPBV5iTjoVVl4Kj38JrvzjBldQaWmCemh3XRY8jXL22+AD\n8UM1hsW06OrqYnTUfgWbbnhrjGn4ARSBnC9tBhhbIM8wYDzHSMA148BkjPqsBMzevXtN3dnzb8bc\ncrox39seeLpYLJpisWjMv3zSXnfoqfrXqQ2ZfR8VY76fM+aT64wpfduYY0ft6yfX2fQayGQyZnJy\ncvZ/wORyc2UWi0XT29s775pa2bt3r/v9X2li2J6GD2md+bduxzh5uRcInbAyxkwYY8QYI9heXBDr\ngWkRyYhIc0wa+Nl1HyxZChddvfB177jWvj5xX/3rpLQ36wbttMlXRyqnUWKyUEwLgImJCbLZLLlc\nrqni0qYxpHX7uP616yJwfQJlDwIjACICkDfGDAVdLCLLgeXAuTXetzqMgV1fhu5fhmXBIexmfeG9\n51Z488/BE1+CK/+oIdVrJ1wN3pZmG1KlhTNHnBQLxbQoFApMTk4yORk065QuaSxadIWkl4FaA1lm\nsIa0x+kJDgGDIhIWXeRGYC8wXeN9q2P/Tph5Ci55T+glpz3z4JxPvHdcC3unoLynIdVrJ1SHlx6T\nk5MUCgV6enro6emZDeXYDKRh8MKWhmqO2usMeQeMMSXn/zx2AaNiRdfhTmAV0Jjh764vgyyBt727\nuuvffo19feIf61cnRUmYsbExZmZmKBaLFIvFpurppWHw3Pk3/1ijh8phbhKUYE6758UYc9gYsw84\nUIf7VrLrH2HNlfDGs6q7vqsbznunHdYqilIzDTd4Tu+rhB1ueukHYu93EZFeEZkJWKzYAJRN2iLl\n556E534Al/xmtHzvuBaeeRRefKY+9VKUDiIt4fEoMCwiI46hyuFbuRWRYSe9Kowx01hDmnPK7XcE\nxyNA+srHXc6w9JLfiJbv7bpaqyhJkYrw2BiTF5EsVkQ8BhSAPnfuzaGPBWQqIeX2OQsUWawBnQYG\njDGFZGpeA7u+DOdvgNNXLnhZRVzaM3vslp/HvwSXf6i+dWwj1B+eEkRqe2mNMRPAxALnXYMYdO6M\nRfI1FzNPw7P/DgO3xcvfdSH8+BvW04XrMFT32CpKZFrSeUBLsTMPk7cAAtN/D6evWtBYzerwum+d\ny79vCt7/D+oJOQKqw1OCUINXT1y3PO+9u2pjVRGX9sE74Jq71BNyRFSD1zjK5TKZTM2qsoagDkDr\nSZB342v+wqZXS5jDUPWErDQBpVKJe++9d97/zg6nRSmXy2zbtq1eVQtEDV49ScJYuQ5DvagnZKVJ\nGB8f5/rr53aERolxkclkKBaL9apaIGrw6kkSxuqqm+wwePeDcPx1+3rfDXORzhQlRfzD2agxLoaG\nhmb95zUCNXj15Kqb4IvDtRkrr6eL/3423Ps7NXu6UJQkyOfzDA3N3z8QNcZFf38/uVzVctua0UWL\nerJuEO6/3brSPjJje3aLGKsKHZ5bzrpB+MaH4Tvj8LaI4uUORHV49WdycnLWQ4pLnBgXPT09DVv4\n0B5ePTlxAn56AN65CW6ZgQ89UlvPbPUV1u37vsY4d1GUMNwA3F68MS62bdvG5ORkVR6PBwcHmZgI\nleQmivbw6kn5aXj9FetKu0oqdHheVl9mX5/+V1hzRe31a2M6WYe3afzhirSNl5zD8FU9VZ2vhkKh\nwODgYEUa2BgXhw4dIp/PMzIyMns+n8+TyWQqHIJ2d3c3bPFCe3j15OAu+3rO26vOMs8fXsXJLjj7\nksqFEKUC9YdXX/r7+2eD9rh4Y1xks1m2bt06ey6bzTI0NES5XOnDI6i3WC+q6uGJyAVYn3HdwJnA\nC9iN+gVjzEv1qlzLc/AJ+3p2ghKSNZdbQfOJ47DkpOTKVdqG7dnLazpfDUG9skKhMLsAMTw8TDab\nJZ/PMzg4yNjYWOgcXT6fZ3g4zGVlsoT28ETkdBHZKiIvYN2v57FRxrLOax6YEZEficiNInJ6Q2rc\nShzcBStWwxsSfGtWXwGvvQQHHkuuTEWJwcDAwOy83WIxLjKZDGeeeWZgOcVisWE7NSoMnmPo/grr\nqLMPuBlYb4xZ4hxdxpglwBnA1VgHAB/EGr8bG1LrVuG5H8A5lyRb5hrn1/lpHdYq6TI4ODjPoBlj\n5hmu8fHxRb0dFwqFCmlLPQnq4U1jQyaeYYy52hhzjzHmu/6LjDEvGmPuN8Z8whhzIXAR8Asisr3O\ndW4Njr8Ozz+ZvMFbcb7tNe7512TLVZQY1Noza3RUs6A5vD5jzItRC3J82V0vIitqr1YbcKhkJSQR\nDV6gDs/Pmsuh+E0bBa3KfYudhurwGkM2m2ViYiLWHFy5XG7YYoVLhcGLY+ySzN82uAsWSffwwO7H\n/f52a1TPbOwDoyheuru75+2lXQivRAVs79CfVm9iy1JE5AIReZcuVoRwcJeNUHbWxZGyrXj8s3Na\nvDBchwRP67A2jO3fe2FWi6fUl1ZxDQW1CY/7sDFmSyKyAVhrjPlMMtVqAw7uslHHli6LlC1Ug+fl\n7LfCsi6rx+v9zzEr2N6oBk8JIrbBM8Z8wfPvUyKyVkS2G2M2JVCv1ufgrvoMZ8HO262+XHt4ihKR\nxHZaGGN2Y3t9yuuvwqFipB0WkVlzOczshsP763cPRWkzajZ4IrJRRDaLyNeBxsilm53nnwRzItIe\n2sisdvbSai9PUaqmZoPnaPHuAT4BDNRepTYgxh7ayLz5nbD0NN1XqygRWHAOT0ROD9srKyKXAsYY\n8+/YPwrODo0tyVezxTj4BCxZGksyUpUOD2yMjPM36I6LEFSH1zjaKYjPgIj8XtAJZ/fFiyJynTOs\n/VNgPOjajuO5H1g5yklL63ufNVfYPbVHKj1QKEojaLcgPi8C4yLysyHnLwWyni1mn0i2esGISEZE\nmvcn5eATsVdoq9LhuRw/CqecBtvWwl2XWS8qCqA6vEbRVkF8jDEF4A7gARFZDhVeVPLA+oXKCENE\nBkVkSkSMiEyKSBRPjTlCFkhqLLd2XjsM5T2xDd6C/vC87Mzb4/3/AB8+CL++zcbAVaMHqD+8RtF2\nQXyMMaNYhwIPOI4BZoBRYAoYMsYE+3xZABEZxBqtAuC6Sigu1msTkW4RGQMCdxvHLTdRnnNCMNZL\ng+eSRMxbRamBtgjiIyIXGGOe8iUPYQ1cL3Y1dtzR3cVlDMg7xhQgLyIz2AWP0aAMIjLM4nOEkctN\nnHruofWiAbqVlGmXID4lETkuIl93hq7XAcexvSoBflSLsXN6W91UGq97gdAIN8aYCWOMGGME66uv\n5nJFZLmIrATOrb4Fi3BwF5y8DDIXJFZkIBqgW0mRdgricz12Xq6XuV6RwRqZMjAhIoeA+2O6d3ff\nAX8/t+jcOy5xyr0RuKWGe1ZycBec8zZYUudwIW6A7mv+wvbs9jxs/9/40freV2l+/vbdlWkX/ypc\n+UfVna+CqEF8stksO3bsoFwuk8vl5i1sNDKIT5B7qDx2MQIAEVmLFRT3MmcIvwAYESkBk8aYP4hw\nz66Q9DJQS582Trl3Yj02n4udp6ydg7vgwo2xs1etw3PDPX51xMpgTjsTfm1MA3Q7qA6vvrhBfLzu\nnbxBfDKZDKOjo4yMjMwawqmpKaanp9m8eTNTU1Oz+ZoqiI8zfJ3X33REx/3AzxN9d8WhkPRaB/CR\nyzXGHAYOV6sbWpRXDsHL++u7pcyLG6D7vhtg15fhHdc15r5Kc/OBr9R2vgqiBPHp7u5mdHR0Np+f\npgjisxDGmO86urshY8xFEbO782/+lvdQORxthnKrJ4EtZZF0eC6rL4dXy7pg4UF1ePWn2iA+vb29\ndHd3UyqV2LhxI2NjY/PKaWQQn4YH4jbGlJyh8BBWPuLSj2co3SzlRiKBFdqqNHh+3ADdex6u/+pw\ni6AavPozODhINptlfHx8NoiPF+8K7sTEBLlcjlwuN6+Xl3oQHxH5sYj8lziFicifOl5TFmMUGBaR\nERHpFZEcvhVWERl20qOwaLl1Y2ceHvwEIPB/3tdYAXBXN7zxHNjzSOPuqShU5+24UCgwOTkZuGrb\n6CA+QUPaq4FNIvKCiNwtIr8SltnZdXGdiGwXkeNY47LoSquzMJIFNmH1fRls8CDv0LOPBWQqNZSb\nPDvzdpfD+z4DH3mu8bseRGwvTz2nKA3GDeKzEJOTkxQKBXp6eujp6WFgwE77pxHER/zd0NkTIv3A\nCHZI6MpSXMORwa6KZrDavAlgrEYxcmo4Wry9e/fuZeXKldELuOsya+TWXjWXtvtBu4L6oWi9riN3\nvwuAZR98IFodHv5L+PoW+JMnYMWqaHnbkGv+/JsA3PfHob/XSkI00lvKvn37WLVqFcAqY8y+qPlD\n5/CcfbQFmDV+/cwZukNYfdu0Meb+GPVuL5ph14M7j/eTR2DF+xp3X6XjaRXXUFDlooXX+CkBuLse\nvD28mLseqtbh+TnvnbD0jXYe72fU4KkOTwkiaNFihbOt7Lgzj/fPIqJPz0JcdRP8v9+3w9jjr9vX\n+26w6Y3ipJPh/PU6j6coCxC0aHEPVky8GzvxfzV2f61OhoSxbhBOXwXbfxtuP8fO3W38aKxdD7F0\neC5rroADj8OrGgtddXhKEEFD2n5gyhizAWY35U8BORHpjrl/tv15tQxrfhHe//maiomlw3NZfZkN\nHvTMo3Bh45b6mxHV4SlBBPXwMnjm64wxZazUo4sQP3Qdz+tH4IUfw3k/k249Vq0HOUn1eIoSQtjW\nsnljAWfRAmLMpXcEB5+wPatzUzZ4p77JRjNTg6cogdTZh1GHsP8x+5p2Dw+sHOaZHXDsaNo1UZSm\nQw1eEhx4DE55U/2dflbD6svg2BHY//20a6IoTUeYDi8rIl3YubwdulCxCPsfsx5SEnD6GVuH5/IW\njyOB82PFV2oLVIenBBFk8O4HNgI343g8dryQGGyc2kljzPcaV8UmxxgrBVnXJGLf5edaZwJ7HoEr\n/jDt2tSVKz/+AHvLR+alnXLyEo4eO+FJsS67VmWW8dDN72pg7ZRmJMjj8QDMejruxerw+rB7ZgeA\nacdhZom5HuBfN6rCTceLP4HXXkxswWJWg9d9a/xCVl8OT/6zNcZJOTdNkSDD5vLUx+e7K7/g5q9U\npLnpirLQXtrdWPHxF9w0x9PxeqwRvBQrVxkGOtfgzS5YrEukuJp0eC5LToZjr8JtXXZ721U3tbTr\n973lI2rEIrBQWMSFqCbgjv8+rifjTZs2VcS4mJ6erjood6OI5ADUGPNd4LvY3RiIyApiBuJuGw48\nBkhNXo4TZWceSt+0Abq9wX2gpY1erZxy8pJAA6lD3fiMjo4yNjZGV1cXW7duZXx8nPHxcbq7uymX\ny2zfvr21DZ4fY8yL2Dm/zmX/TuhaazVwzcCDd8A1d805MnADdH91pOkN3kJD11o5euxEx/QSo/bU\n4lAqlRgbG5u919jY2GyPz+1hNjLAdrU03MV723HgsfQFx16awVVVTHTo2jq4hi6ft05u+/v76e7u\nnmfkyuWK8NGpowavFl57GQ7thp99f9o1mSNBV1XNTtAwNWzoGqUM0KFuNQwNDbFhwwYymQybN29m\nw4YN88I2bt26tSJgT9qowauFg08AJtEeXs06vBYI0B116BrHKLnDKrcnEmYEO2momySlUoktW7bM\nztENDw8zPT3N0NDQ7Ht/zz33pFnFQNTg1cL+nfa1GbaUufgDdC87A379E001fxd16BpmlJT0CJon\n7O3tbcp5Oy9q8GrhwGPwhhWw4i2JFZmIDs8N0P3A7fAvd8KFGxOpWyvh+sLb4nwxV2WWaa9NUYNX\nE/udBYsExb2J6PBcLhywoSOLD7SE2/eFhq5R8fvDCxv6hhlBndtrT9TgxeXECbul7NLfSrsm4Zy/\n3g5pf1RoGoMXFiUPmmvoqnN77YkavLjM7IbXf9pckhQ/S06CnnfBjyetgU7AuUEU6qmrSwod6nYW\navDicqCJfOAtxEVXw2NfgP3fg5WXNvTWYYsTF3/4a4kNXWsl6lBXaW1SM3giMghswTooKABZY8yC\nmwAXy+OcD1om2maMGU2q7oCdv5MlzbOlLIweZ8HiR5MNN3hhNNPQVeksUjF4HsO0DdiKdUJQFJEz\nnBgacfN0A2Uct1YediTeiAOPw5kXwtJkeyU16/D8vOlsWNlrDd4vjSx+fZug/vCUINLq4Y0BeU+v\nKy8iM9jeW1hPrJo8PVh3VRN1qvccB3baoDmtwEUD8O1t8MohOK0r7dq0BGFze6fq6m1L03CD54R9\n7Mb20LzcCwwSYPAi5FkPFNzrjTHTCVZ9jldfhPIe6PvdxItORIfn56Kr4dtjVp7SIAHyjw+mGybR\nr8OLSpDxev34CS76s6/p6m0Lk0YPz30C/fN1ReD6GvN0Yw3gCIDjqDRvjBkKKlRElgPLgXOrqfgs\nBx63r+cm4wPPS6I6PJeVl8KyLvjRN+pi8MJWYxfqDdWbesSlXXqShoBpddIweGFjqjI2Jm4teTLY\n+boBY0zJnfcTkXFjjL93CHAjcEsVdZ5jZx4KtwAC3/gzeO2lptq2FciSk2xg7h/fXxd5ykJbxXRx\nQmkm0vjJOhSSHmbsqs5jjBFjzIC7cmuMyWOHu8Mh+e8EVmFXfRdnZx7uvw2uvRs+8hy8+077/858\nVdlT5aIBeOV5ePa7addEUVIjjR6ed0XVO0TtoXLIWkselxLYeUD/CrAx5jBwWKrdGvbgHdYLSQs6\n17TyFLGrtav60q5N26Fb0VqDhhs8Z6hZAoawWjqXfiCwq1RNHhHpxYm45lus2ACUw+QukWhh55q8\n8Uw4Yy18Z9wuYLRBrIs0WGhnhi5mND9pyVJGsXNrRawB24LtvY27F4jIMHYubqiaPMaYacco5kRk\nHJjGRlkbwRrK2mmAc83EdXguO/Nw/DW4/u86ItZFvXR4ujOjtUnF4Blj8iKSxcpMxrAGrM+306IP\nu+JadR5jTJ9j7LJYmzGNNZreXmF8rroJvpiF68ab1rlmKA/eAe/9q8SH45//zp5EvZwoSj1JbWuZ\nIw4OFQg7q6pZX9qCeTz56sO6QZi8BXK/C0dmbM9u40cT7SHVRYcHiQzHw+Qnbzr1ZHbd9h84aUnz\nxMCtVYentCfqPCAKhw/AS89A/63wi39Sl1vURYcHiQzHF5KfNJOxg/ro8OKgixnNhRq8KJS+ZV+7\nfyXVasSiBWJdtDK6mNEaqMGLQulb1qHmee9MuybR8ca6eP4HsPSNcPXtbblgkQa6mNEa6F6ZajHG\nGry1v9RwR5qJsW4QPvQIfPBhOPoyvPJC2jVSlIbSot/cFHj+STi8D3pacDjr55xLrCfkf7sHjh2t\nKsux4yfqXClFqT86pK2W2fm7X67rbeqmw/Nz2R/A5wbhiS/BOyt9NgStyLaS/KSZ/eFd+fEHgMrh\nri5k1B81eNVS+haccYE92oGejXDWxfDwXbBuqCLymjoESIawxQxdyEgHNXjVcPx12P0vDZngr5sO\nz8+SJfALvw9f+W+w5xFYc/nieVqIZtHhBfXY1LClh87hVcPeaTh6uO7DWbA6vLpp8fz87PvtqvMj\ndzXmfg3kkT0vN40WT2ketIdXDaVvATJftNsOnHIa9H0AHvoUzDw1O1zXBYp0UJFy/VGDVw2lb8LK\nn2vPeBA/vxke+nP4m1+Flw9izrqYvz9pEGjyaGxtiAb/rj9q8BbjtcPwzKNwxR+mXZO6cOv/uocP\nn3Y2J193D6y+HNnzMJvyH+SxU97LBTdXXt+MK7Kthgb/Tg81eIvx9L/CiWMNmb9Lg9/ji5w8eM88\nLyqnDd7Nn3z+D7nzttvTrVybEnVXhg51k0MN3mIUvwknvwHecllDbtcwHZ7Dytf3BHpRWfn6ngbV\noD40sw4vKjrUTQ5dpV2M0resQVj6hrRrUhf2LV1tHQl42fOwTVeUNkN7eGHszFtX6M//yEYm25lv\nLx0e8JXvP8vXX/kNPu33ovLFLNteuY5P170G9aNZdHhR0Lm9+qMGLwg3OpnflRLU3ejVS4MX5rwT\nruDTGy91vKj8EN50DhzezyGzvC71aBStqMHTub36owYviFaOThZC2Faxiz/8NS743DLgIwCceuQo\n/3TKn/E/T5mAI1krTFaaEp3bi44avCBaOTpZRAK/NPsugM/0w+c22eH88z/UKGcpEnWoqz2/cNTg\nBdGA6GT1InzoGoGVl8LbfgOefggG/6Yjopw1M1GHutrzC0cNXhAt7A59IS8nkXjuB9bYtdGwvtPR\nnp8avGDmuUP/YV2ik4XRaB1eKM8/2dLD+nbS4SWF9vzU4IWzbrCpezJRh66RnXeGDevfdC7cdZnO\n6zUBSclYOqnnpwavyahWhxd16Br26x5K0LD+H37L7jr59W1NP6/Xijq8qCQVOCjs2bj4w19rO0Oo\nBq/J8OvwElmEiEPQsP7U5fDev2qJeb1W1OElRVI9vyBD6D6PreqePjWDJyKDwBagFygAWWNMqdY8\nccptBuYM2402wfNARenJJRp3wj+s/9gZ4fN6//55eOjTlUPdnXmra9QhcMOoZ8jIBfWcAeUvW3oS\nR14/XnV6vQ1nKgbPMUo5YBuwFcgCRRE5wxhTjpsnTrlRCOttRflQTzl5CUePBTvYTGJCOfLQNQph\n83pLlkLhFnjfX88f6v7kO/Dk14N3rECwIQwzkFHTlQrquXVtoQWRKOmLGc5jh2sMLWqMafgBFIGc\nL20GGKt8mRBBAAAKUUlEQVQlT5xynWtWAmbv3r1mIdaM/lPN6UmUESc9Eb6fM+aT64wpfduYY0ft\n6yfXGXPHW+3fXkrfNuZ/nB+cfsfbgsv5yk3JpH8/Z37zUw+Y3/zUA/V7L9qIVnh+3fS9e/cawAAr\nTQzb0/AenohksKqLrO/UvcAgMBonT8xylwPLgTcD7N+/f8G6Hzv8Avv27aspPYkyAJYcmeH8D322\nIv3NK94QeH0inHkFrLsBPv9f4VARunrgshvgq38KJ68B731PXgPPvxic/twhuPZuOPVCOPCcfb3s\n1rlFklrT//F2Xj2+BaB+70UbcfaSnwY+S1D5/iX1/MZN93xH43l6imMlazmwc2sG6PaljwAzcfPE\nLPdWJ48eeujRWseVLdHDA8ICQ5SBTA154pR7JzABrAL+DVgPPBtybbtwLjCN/YE4kHJd6kmntBM6\nq61vBnYAP4mTOQ2DdygkPcwoVZsncrnGmMPAYZkLQv2sMaatx0Ceth5o57Z2SjuhY9t6LE7+NDwe\nu6ulfkVoDxAmH6kmT5xyFUXpIBpu8IzVxJWAId+pfiAfN0+ccj0cBj7mvLY7ndLWTmknaFurRpzJ\n+4bi0cuNYsXBW7ArqT2O4UJEhoEBY8xQhDyLXqMoSueSShAfY0weKx/ZBExh59n6fEapD2usqs5T\nZbmKonQoqfTwFEVR0kDDNLYxIpJxBNmK0rKISLeIJOL2puMNnogMisiUiBgRmUzqjW0ScsBw0Il2\nareIjHjaUnTmf/3XtHR7nR+vcRGZcdow5cxZ+69r6Xb6cX6wpwjeKRW5rR1t8DyLHAXmVneLrd4r\ncn4Rx7Ar1EHn26bdTjvHsG0ZwK7Ijzvp7jXt0N4ccD32iz+AFd/mRGT2M26TdvrJEaCljd3WRm8t\na6aDmM4GmvnA9ui8W3BG2rndThvHfGlj9tFuj/Y6X3gDDAZ8juPt0s6Ado84bZrXzlra2rE9PI+z\ngXHfKdfZQEtijJkwxogxRpgTY8/STu12hjAlYLvvVNFzTbu015Va+XFdo7VLOwEQkV7sD9cQvue4\nlrZ2rMFjbkeGX7JSJHxfbjvQNu02xpSMMT3GmGnfqSx2bym0QXuNMWVjzDZjTNmZrhgRkUnntPul\nb/l2+sgBowGfLdTQ1k42eHGcDbQDbdtuEekVkSnsF8Kd12m39vZjdab92C+8u4e8bdopIjmgZIzZ\nFnJJ7LZ2ssGL48SgHWi7djsrmDnsal4JWGvmxOZt1V5nyqLPmbLoAu53TrVFO53FiH4qt4h6id3W\nTjZ4nepsoK3a7czj7ca6RuozxgyZ+e78W769Ts91JODUdmy7oQ3a6bABa7hc+Y3BtnHY+b+XGtra\nsQbP1OZsoGVpw3ZPAjtC5vLaqb1jATqz2S94G7VzHCu78R4l5mRHpVra2ulhGkexWqYic84GglZ/\n2o22aLfza98N5INEuMburYYWb68xZlpESsCkiIxiezgDWAmSN6RBS7cT5hnuWUSkjDV03lXqeG1N\nW2uT9oF9aKawOqdJoDftOiXYthkCdHjt0m4qNYfzjnZqr/NlzmFXIo3TlsGA61q6nSFtn8Knw4vb\nVnUeoChKx9Cxc3iKonQeavAURekY1OApitIxqMFLGY9roxnnmBKRsRb3cDEPp02Tnv+LjlC4UffP\nOTswGpq3lXFcL5lWdy/lRw1eSjhi0hnsBmmwG58nsCryEWB3uz1sHrxbopTmZAPMykTahk7X4aWC\n03tztwQNmPn6Ildfdj92qb2nwdWrO8aYgbTroCxKL8HeWVoa7eGlwxh2+8yQ39iBFZoCW4Fur4NH\nmHXuOekZAk86BtJ7zZTjHbfb4xF2Jmh7UsTyegOGp67XWe+QPNDxqCdPv3sPZx+sCTlmotTTuW7Y\n0+ZF6xInb4S6+D3y9jpTGF7npAu9t9Xep6rrItLPnMcZ9z79Tptad7olbVFhJx5YQfBUjHy9WJFl\nETvsdR0kGjyiS6wYc9I5N4b1ETbpXDccs7ycp7x+j/DTFX0OO+fc/N3++vjL8/zf7ztmy41YzxHm\nRLkjWNW9qeb9rjZvhLp43xtvecWg9yLgva32PlVdF/E568bncNSp27znpxWP1CvQaQdz3msrlONV\n5J0CigHp/i+lqz7v911nfIYmanl+j7tTfkMS8mVZ0OCFtHMGyFRbT8/7Oum7xjU8oQYvSt4I79lM\nQHkj/vss8t5W+9ksel3E58xtd8Y5Jp3yWn7Xhg5pG4/ry6voP+EZ/niPcedcBvtrPubPh13s8A9h\nyqZyuDw7AR2zvHkbs411U9Tnu85daInldNIZ7vVih/vlCPV0h5/zgr0YYyYI8Pzso6q81dbFGQpn\n/NeZcP9u897bCPeJ+hlWSx+23V1Yg9qFdbkV5IyzpdBFiwZjjCmJCAQvRmxlvsubLZ6/3fRx1wj6\nEZGMmXONFLS65v3iJ1Geu8DSj13Vczfzx8IxFCPANo+xrqqehHvBBRvwZiEDXG3eqHXZEXBJkPH1\n3zfqfar9DKtlvfPq/ijnY5TRlKjBS4cycw/VLP4elHdy20NYbAN8D2W1so/Y5Tn1G8F+YfPYeaoS\nAb3XxXC+wDmgYIypCMm3WD2dH5FGUUtdgj6XsM+q2vtU+xlWi9szzGJ/mEdiGs6mQ4e06TAB9EqA\nSyOXgBVCtxdQNsZMew8nPWrPqqbyHAPl9sZ6jDGjTq8srr7Olen4fZxVW8+S738vFT8uIfdYLG/U\nugTdt5rPKep9knom3B472CmFCebcLW0JydJSqMFLAacHUwbuCZI+OA9dzpenjJUJjHplATKn6dsU\nsQ61lud+mR71pV8fpR7OPcfxzNvFrKfbw5nXK3Z+VBaTUVSVN0JddmA/33k9VQn2WlxBtfdJ+plw\ncJ/HgnMP1/nmbHDzBSQ4TS+U1yFtemzEEReLSIE5zZM7JzZB5S/0ZpyVUc+cTRb7pQwaBi5G7PKM\ndUrpGm03XKIbXAZgSEQKZhGlvmPwh3G+YP4fAKfXuGg9nSHeKNYz8BTW/XmPU/aCk+0R80apy7ij\nq8thFwL6sYawmqFhtZ9NVdeJyLDz/6h/6sTHBqyzTW8dx4F+5wcg7L3swj6vzb0zI+1l4k4+mFvJ\nc6UJM1gJgKvFGsSne8I+VK5MwL2+23fNPBmIL90fvLiW8nqZ0/cVsV8Mt00Gx/moP7+3HsxJNcKO\n7mrr6VzndQo55byH/YQ4Qo2TN0JdBpmT2Ew6+fzBswPf24j3qeYzdN/nBXV0/vp50t1yvXKjfl8d\nMml/pxY71AGootQBp0c1r4frDDVnsL2sMIlKveuEsXNzccvoN84KuogMu2V505sZncNTlPowhp2u\n6IfZ6GruvGxsgxMXx9gOESyViYJXazjk/N1NDXKkRqJzeIpSH/qwBm7SIx8pYUNJpiHvWI8dqtYs\nHhYbOCePnfMsAtPGmIXiyDYNOqRVlDriEQj7FwKUFFCDpyhKx6BzeIqidAxq8BRF6RjU4CmK0jGo\nwVMUpWNQg6coSsegBk9RlI5BDZ6iKB2DGjxFUTqG/w9Ru28vW1FqQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f20d39f7250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(4.5,3));\n",
    "plt.rcParams['xtick.major.width'] = 1.2\n",
    "plt.rcParams['ytick.major.width'] = 1.2\n",
    "plt.rcParams['axes.linewidth'] = 1.2\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "norm_k1 = 1.*np.array(degree_counter.values())/sum(degree_counter.values())\n",
    "norm_k2 = 1.*np.array(triangle_degree_counter.values())/sum(triangle_degree_counter.values())\n",
    "\n",
    "#################### plot pdf\n",
    "ax.plot(degree_counter.keys(), norm_k1, 's-', label='$k_1$',\n",
    "        clip_on=True, mfc='white');\n",
    "ax.plot(triangle_degree_counter.keys(), norm_k2, 'o-', label='$k_2$',\n",
    "        clip_on=True, mfc='white');\n",
    "\n",
    "\n",
    "#################### real mean\n",
    "sum_of_numbers = sum(number*count for number, count in degree_counter.iteritems())\n",
    "count = sum(count for n, count in degree_counter.iteritems())\n",
    "mean_k1 = float(sum_of_numbers) / count\n",
    "\n",
    "sum_of_numbers = sum(number*count for number, count in triangle_degree_counter.iteritems())\n",
    "count = sum(count for n, count in triangle_degree_counter.iteritems())\n",
    "mean_k2 = float(sum_of_numbers) / count\n",
    "\n",
    "ax.vlines(mean_k1,0,0.2, linewidth=1.5, linestyles='--',\n",
    "          color=u'#1f77b4', label=r'$\\langle k_1 \\rangle$');\n",
    "ax.vlines(mean_k2,0,0.2, linewidth=1.5, linestyles='--',\n",
    "          color=u'#ff7f0e', label=r'$\\langle k_2 \\rangle$');\n",
    "\n",
    "#################### expected mean\n",
    "ax.vlines(k2,0,0.2, linewidth=1.5, linestyles='-', color='lightgray', zorder=1, label=r'$\\approx$');\n",
    "ax.vlines(k1,0,0.2, linewidth=1.5, linestyles='-', color='lightgray', zorder=1);\n",
    "\n",
    "ax.tick_params(axis='both', which='major', labelsize=16)\n",
    "ax.set_xlabel('Generalized degree, $k_\\omega$', size=18)\n",
    "ax.set_ylabel('$P(k_\\omega)$', size=18)\n",
    "\n",
    "ax.set_ylim(0,0.2);\n",
    "ax.set_xlim(0,40);\n",
    "\n",
    "plt.legend(fontsize=14, loc='upper right', handlelength=1, frameon=False, borderpad=0)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}