class03_degdist_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS446/546 - Class Session 3 - Degree Distribution\n",
    "\n",
    "In this class session we are going to plot the degree distribution of the undirected human\n",
    "protein-protein interaction network (PPI), without using `igraph`. We'll obtain the interaction data from the Pathway Commons SIF file (in the `shared/` folder) and we'll \n",
    "manually compute the degree of each vertex (protein) in the network. We'll then compute\n",
    "the count `N(k)` of vertices that have a given vertex degree `k`, for all `k` values.\n",
    "Finally, we'll plot the degree distribution and discuss whether it is consistent with the \n",
    "results obtained in the Jeong et al. article for the yeast PPI.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by loading all of the Python modules that we will need for this notebook.  Because we'll be calling a bunch of functions from `numpy` and `matplotlib.pyplot`, we'll alias them as `np` and `plt`, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import collections\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import igraph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 1:  load in the SIF file as a `pandas` data frame using `pandas.read_csv`.  Make sure the column names of your data frame are `species1`, `interaction_type`, and `species2`. Save the data frame as the object `sif_data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sif_data = pandas.read_csv(\"shared/pathway_commons.sif\",\n",
    "                            sep=\"\\t\", names=[\"species1\",\"interaction_type\",\"species2\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 2:  restrict the interactions to protein-protein undirected (\"in-complex-with\", \"interacts-with\"). The restricted data frame should be called `interac_ppi`. Then we will make a copy using `copy` so `interac_ppi` is independent of `sif_data` which will be convenient for this exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "interaction_types_ppi = set([\"interacts-with\",\n",
    "                             \"in-complex-with\"])\n",
    "interac_ppi = sif_data[sif_data.interaction_type.isin(interaction_types_ppi)].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 3: for each interaction, reorder `species1` and `species2` (if necessary) so that\n",
    "`species1 < species2` (in terms of the species names, in lexicographic order).  You can make a boolean vector `boolean_vec` containing (for each row of the data frame `interac_ppi`) `True` if `species2 > species1` (by lexicographic order) for that row, or `False` otherwise.  You can then use the `loc` method on the data frame, to select rows based on `boolean_vec` and the two columns that you want (`species1` and `species2`). Thanks to Garrett Bauer for suggesting this approach (which is more elegant than looping over all rows):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "boolean_vec = interac_ppi['species1'] > interac_ppi['species2']\n",
    "interac_ppi.loc[boolean_vec, ['species1', 'species2']] = interac_ppi.loc[boolean_vec, ['species2', 'species1']].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since iterating is reasonably fast in Python, you could also do this using a `for` loop through all of the rows of the data frame, swapping `species1` and `species2` entries as needed (and in-place in the data frame) so that in the resulting data frame `interac_ppi` satisfies `species1 < species2` for all rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rowid, row in interac_ppi.head().iterrows():\n",
    "    if row['species1'] > row['species2']:\n",
    "        interac_ppi['species1'][rowid] = row['species2'] \n",
    "        interac_ppi['species2'][rowid] = row['species1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 4: Restrict the data frame to only the columns `species1` and `species2`. Use the `drop_duplicates` method to subset the rows of the resulting two-column data frame to only unique rows. Assign the resulting data frame object to have the name `interac_ppi_unique`.  This is basically selecting only unique pairs of proteins, regardless of interaction type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "interac_ppi_unique = interac_ppi[[\"species1\",\"species2\"]].drop_duplicates()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 5: compute the degree of each vertex (though we will not associate the vertex degrees with vertex names here, since for this exercise we only need the vector of vertex degree values, not the associated vertex IDs).  You'll want to create an object called `vertex_degrees_ctr` which is of class `collections.Counter`.  You'll want to name the final list of vertex degrees, `vertex_degrees`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "allproteins = interac_ppi_unique[\"species1\"].tolist() + interac_ppi_unique[\"species2\"].tolist()\n",
    "vertex_degrees_ctr = collections.Counter(allproteins)\n",
    "vertex_degrees = list(vertex_degrees_ctr.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print out the vertex degrees of the first 10 vertices, in whatever the key order is. Pythonistas -- anyone know of a less convoluted way to do this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A1BG': 19,\n",
       " 'A1CF': 27,\n",
       " 'A2M': 169,\n",
       " 'A2ML1': 6,\n",
       " 'A4GALT': 3,\n",
       " 'A4GNT': 14,\n",
       " 'AAAS': 68,\n",
       " 'AACS': 4,\n",
       " 'AADAC': 1}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict(list(dict(vertex_degrees_ctr).items())[0:9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print out the first ten entries of the `vertex_degrees` list.  Note that we don't expect it to be in the same order as the output from the previous command above, since `dict` changes the order in the above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[19, 27, 169, 6, 3, 14, 68, 4, 1]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vertex_degrees[0:9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 6: Calculate the histogram of N(k) vs. k, using 30 bins, using `plt.hist`.  You'll probably want to start by making a `numpy.array` from your `vertex_degrees`.  Call the resulting object from `plt.hist`, `hist_res`. Obtain a numpy array of the bin counts as element zero from `hist_res` (name this object `hist_counts`) and obtain a numpy array of the bin centers (which are `k` values) as element one from `hist_res` (name this object `hist_breaks`). Finally, you want the `k` values of the centers of the bins, not the breakpoint values. So you'll have to do some arithmetic to go from the 31 `k` values of the bin breakpoints, to a numpy array of the 30 `k` values of the centers of the bins.  You should call that object `kvals`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAURElEQVR4nO3dfWxd933f8fdnUuw8R3LMeZ4kTEojpFCMLlE5WUWKYIg2SXaKygPcQMYwa6lQAYuypVuHRG6AqktqwN5DvRhLHKixajkL/DA3hYXFqarZLoIBtWw6fpLtOGJsJ5IgW2wkO+2COlX63R/3x/iGJSmSl7ykovcLuLjnfM/v3Ps9R6Q/PA/3OlWFJOn89vfmuwFJ0vwzDCRJhoEkyTCQJGEYSJKAxfPdwExdfPHFtXLlyvluQ5LOKY8++uhfVNXA2Po5GwYrV65kaGhovtuQpHNKku+OV/c0kSTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSOIc/gdyLlbu+NqVxL97w4TnuRJIWBo8MJEmGgSTJMJAkYRhIkjAMJEkYBpIkDANJEoaBJIkphEGSvUlOJjk8zrLfSlJJLm7zSXJzkuEkTyZZ2zV2W5Ij7bGtq/6LSZ5q69ycJLO1cZKkqZnKkcFtwOaxxSQrgI3A97rKVwCr22MHcEsbexGwG7gcWAfsTrK0rXML8Btd6/2d95Ikza2zhkFVfQM4Nc6im4BPAtVV2wLcXh0PAUuSXApsAg5W1amqOg0cBDa3ZW+vqoeqqoDbgat62yRJ0nTN6JpBki3A8ap6YsyiZcDRrvljrTZZ/dg49Yned0eSoSRDIyMjM2ldkjSOaYdBkjcDvw38zuy3M7mq2lNVg1U1ODAw0O+3l6SfWTM5Mvg5YBXwRJIXgeXAN5P8A+A4sKJr7PJWm6y+fJy6JKmPph0GVfVUVf39qlpZVSvpnNpZW1UvAfuBa9tdReuBV6vqBHAA2JhkabtwvBE40Jb9IMn6dhfRtcC9s7RtkqQpmsqtpXcAfw68J8mxJNsnGX4f8DwwDPwB8DGAqjoFfBZ4pD0+02q0MV9q63wH+PrMNkWSNFNn/Z/bVNU1Z1m+smu6gJ0TjNsL7B2nPgRcdrY+JElzx08gS5IMA0mSYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEliCmGQZG+Sk0kOd9X+S5JvJXkyyR8nWdK17Lokw0meS7Kpq7651YaT7Oqqr0pyqNXvSnLBbG6gJOnspnJkcBuweUztIHBZVf0C8G3gOoAka4CtwHvbOl9IsijJIuDzwBXAGuCaNhbgRuCmqno3cBrY3tMWSZKm7axhUFXfAE6Nqf1pVZ1psw8By9v0FuDOqnqtql4AhoF17TFcVc9X1Y+AO4EtSQJ8CLinrb8PuKrHbZIkTdNsXDP4deDrbXoZcLRr2bFWm6j+TuCVrmAZrY8ryY4kQ0mGRkZGZqF1SRL0GAZJPg2cAb4yO+1Mrqr2VNVgVQ0ODAz04y0l6byweKYrJvnXwK8AG6qqWvk4sKJr2PJWY4L694ElSRa3o4Pu8ZKkPpnRkUGSzcAngV+tqh92LdoPbE1yYZJVwGrgYeARYHW7c+gCOheZ97cQeRC4uq2/Dbh3ZpsiSZqpqdxaegfw58B7khxLsh34H8DbgINJHk/yRYCqehq4G3gG+BNgZ1X9uP3V/3HgAPAscHcbC/Ap4D8kGaZzDeHWWd1CSdJZnfU0UVVdM055wv9gV9X1wPXj1O8D7hun/jydu40kSfPETyBLkgwDSZJhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgSWIKYZBkb5KTSQ531S5KcjDJkfa8tNWT5OYkw0meTLK2a51tbfyRJNu66r+Y5Km2zs1JMtsbKUma3FSODG4DNo+p7QLur6rVwP1tHuAKYHV77ABugU54ALuBy4F1wO7RAGljfqNrvbHvJUmaY2cNg6r6BnBqTHkLsK9N7wOu6qrfXh0PAUuSXApsAg5W1amqOg0cBDa3ZW+vqoeqqoDbu15LktQnM71mcElVnWjTLwGXtOllwNGuccdabbL6sXHq40qyI8lQkqGRkZEZti5JGqvnC8jtL/qahV6m8l57qmqwqgYHBgb68ZaSdF6YaRi83E7x0J5PtvpxYEXXuOWtNll9+Th1SVIfzTQM9gOjdwRtA+7tql/b7ipaD7zaTicdADYmWdouHG8EDrRlP0iyvt1FdG3Xa0mS+mTx2QYkuQP4p8DFSY7RuSvoBuDuJNuB7wIfacPvA64EhoEfAh8FqKpTST4LPNLGfaaqRi9Kf4zOHUtvAr7eHpKkPjprGFTVNRMs2jDO2AJ2TvA6e4G949SHgMvO1ockae74CWRJkmEgSTIMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiS6DEMkvz7JE8nOZzkjiRvTLIqyaEkw0nuSnJBG3thmx9uy1d2vc51rf5ckk29bZIkabpmHAZJlgH/DhisqsuARcBW4Ebgpqp6N3Aa2N5W2Q6cbvWb2jiSrGnrvRfYDHwhyaKZ9iVJmr5eTxMtBt6UZDHwZuAE8CHgnrZ8H3BVm97S5mnLNyRJq99ZVa9V1QvAMLCux74kSdMw4zCoquPAfwW+RycEXgUeBV6pqjNt2DFgWZteBhxt655p49/ZXR9nHUlSH/Rymmgpnb/qVwH/EHgLndM8cybJjiRDSYZGRkbm8q0k6bzSy2mifwa8UFUjVfU3wFeBDwBL2mkjgOXA8TZ9HFgB0Ja/A/h+d32cdX5KVe2pqsGqGhwYGOihdUlSt17C4HvA+iRvbuf+NwDPAA8CV7cx24B72/T+Nk9b/kBVVatvbXcbrQJWAw/30JckaZoWn33I+KrqUJJ7gG8CZ4DHgD3A14A7k/xeq93aVrkV+HKSYeAUnTuIqKqnk9xNJ0jOADur6scz7UuSNH0zDgOAqtoN7B5Tfp5x7gaqqr8Gfm2C17keuL6XXiRJM+cnkCVJhoEkyTCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAk0WMYJFmS5J4k30rybJJfSnJRkoNJjrTnpW1sktycZDjJk0nWdr3Otjb+SJJtvW6UJGl6ej0y+BzwJ1X188A/Bp4FdgH3V9Vq4P42D3AFsLo9dgC3ACS5CNgNXA6sA3aPBogkqT9mHAZJ3gF8ELgVoKp+VFWvAFuAfW3YPuCqNr0FuL06HgKWJLkU2AQcrKpTVXUaOAhsnmlfkqTp6+XIYBUwAvxhkseSfCnJW4BLqupEG/MScEmbXgYc7Vr/WKtNVP87kuxIMpRkaGRkpIfWJUndegmDxcBa4Jaqej/w/3j9lBAAVVVA9fAeP6Wq9lTVYFUNDgwMzNbLStJ5r5cwOAYcq6pDbf4eOuHwcjv9Q3s+2ZYfB1Z0rb+81SaqS5L6ZMZhUFUvAUeTvKeVNgDPAPuB0TuCtgH3tun9wLXtrqL1wKvtdNIBYGOSpe3C8cZWkyT1yeIe1/+3wFeSXAA8D3yUTsDcnWQ78F3gI23sfcCVwDDwwzaWqjqV5LPAI23cZ6rqVI99SZKmoacwqKrHgcFxFm0YZ2wBOyd4nb3A3l56kSTNnJ9AliQZBpIkw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJLELIRBkkVJHkvyv9v8qiSHkgwnuSvJBa1+YZsfbstXdr3Gda3+XJJNvfYkSZqe2Tgy+ATwbNf8jcBNVfVu4DSwvdW3A6db/aY2jiRrgK3Ae4HNwBeSLJqFviRJU9RTGCRZDnwY+FKbD/Ah4J42ZB9wVZve0uZpyze08VuAO6vqtap6ARgG1vXSlyRpeno9MvjvwCeBv23z7wReqaozbf4YsKxNLwOOArTlr7bxP6mPs85PSbIjyVCSoZGRkR5blySNmnEYJPkV4GRVPTqL/UyqqvZU1WBVDQ4MDPTrbSXpZ97iHtb9APCrSa4E3gi8HfgcsCTJ4vbX/3LgeBt/HFgBHEuyGHgH8P2u+qjudSRJfTDjI4Oquq6qllfVSjoXgB+oqn8JPAhc3YZtA+5t0/vbPG35A1VVrb613W20ClgNPDzTviRJ09fLkcFEPgXcmeT3gMeAW1v9VuDLSYaBU3QChKp6OsndwDPAGWBnVf14DvqSJE1gVsKgqv4M+LM2/Tzj3A1UVX8N/NoE618PXD8bvUiSps9PIEuSDANJkmEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJoocwSLIiyYNJnknydJJPtPpFSQ4mOdKel7Z6ktycZDjJk0nWdr3Wtjb+SJJtvW+WJGk6ejkyOAP8VlWtAdYDO5OsAXYB91fVauD+Ng9wBbC6PXYAt0AnPIDdwOXAOmD3aIBIkvpjxmFQVSeq6ptt+i+BZ4FlwBZgXxu2D7iqTW8Bbq+Oh4AlSS4FNgEHq+pUVZ0GDgKbZ9qXJGn6ZuWaQZKVwPuBQ8AlVXWiLXoJuKRNLwOOdq12rNUmqo/3PjuSDCUZGhkZmY3WJUnMQhgkeSvwR8BvVtUPupdVVQHV63t0vd6eqhqsqsGBgYHZellJOu/1FAZJ3kAnCL5SVV9t5Zfb6R/a88lWPw6s6Fp9eatNVJck9UkvdxMFuBV4tqp+v2vRfmD0jqBtwL1d9WvbXUXrgVfb6aQDwMYkS9uF442tJknqk8U9rPsB4F8BTyV5vNV+G7gBuDvJduC7wEfasvuAK4Fh4IfARwGq6lSSzwKPtHGfqapTPfQlSZqmGYdBVf1fIBMs3jDO+AJ2TvBae4G9M+1FktQbP4EsSTIMJEmGgSSJ3i4g/8xbuetrUxr34g0fnuNOJGlueWQgSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRfVDcr/EI7Sec6jwwkSYaBJGkBnSZKshn4HLAI+FJV3TDPLc26qZ5OAk8pSeqvBXFkkGQR8HngCmANcE2SNfPblSSdPxbKkcE6YLiqngdIciewBXhmXruaR9M5ipiKqR5peDFcOj8tlDBYBhztmj8GXD52UJIdwI42+1dJnpvh+10M/MUM1+2nWeszN87Gq0z4eufK/oRzp1f7nH3nSq9z3ec/Gq+4UMJgSqpqD7Cn19dJMlRVg7PQ0pyyz9l3rvRqn7PvXOl1vvpcENcMgOPAiq755a0mSeqDhRIGjwCrk6xKcgGwFdg/zz1J0nljQZwmqqozST4OHKBza+neqnp6Dt+y51NNfWKfs+9c6dU+Z9+50uu89Jmqmo/3lSQtIAvlNJEkaR4ZBpKk8ysMkmxO8lyS4SS7FkA/LyZ5KsnjSYZa7aIkB5Mcac9LWz1Jbm69P5lk7Rz3tjfJySSHu2rT7i3Jtjb+SJJtferzd5Mcb/v18SRXdi27rvX5XJJNXfU5/dlIsiLJg0meSfJ0kk+0+oLap5P0uRD36RuTPJzkidbrf2r1VUkOtfe9q92UQpIL2/xwW77ybNswx33eluSFrn36vlafn9+nqjovHnQuTH8HeBdwAfAEsGaee3oRuHhM7T8Du9r0LuDGNn0l8HUgwHrg0Bz39kFgLXB4pr0BFwHPt+elbXppH/r8XeA/jjN2Tft3vxBY1X4eFvXjZwO4FFjbpt8GfLv1s6D26SR9LsR9GuCtbfoNwKG2r+4Gtrb6F4F/06Y/BnyxTW8F7ppsG/rQ523A1eOMn5d/+/PpyOAnX3lRVT8CRr/yYqHZAuxr0/uAq7rqt1fHQ8CSJJfOVRNV9Q3gVI+9bQIOVtWpqjoNHAQ296HPiWwB7qyq16rqBWCYzs/FnP9sVNWJqvpmm/5L4Fk6n7xfUPt0kj4nMp/7tKrqr9rsG9qjgA8B97T62H06uq/vATYkySTbMNd9TmRe/u3PpzAY7ysvJvsh74cC/jTJo+l81QbAJVV1ok2/BFzSphdC/9PtbT57/ng7xN47euplkn762mc7PfF+On8hLth9OqZPWID7NMmiJI8DJ+n8x/E7wCtVdWac9/1JT235q8A7+9Hr2D6ranSfXt/26U1JLhzb55h+5rTP8ykMFqJfrqq1dL6tdWeSD3YvrM6x4YK893ch9wbcAvwc8D7gBPDf5red1yV5K/BHwG9W1Q+6ly2kfTpOnwtyn1bVj6vqfXS+tWAd8PPz3NK4xvaZ5DLgOjr9/hM6p34+NY8tnldhsOC+8qKqjrfnk8Af0/lhfnn09E97PtmGL4T+p9vbvPRcVS+3X76/Bf6A1w/557XPJG+g8x/Yr1TVV1t5we3T8fpcqPt0VFW9AjwI/BKd0yqjH6jtft+f9NSWvwP4fj977epzczslV1X1GvCHzPM+PZ/CYEF95UWStyR52+g0sBE43HoavUtgG3Bvm94PXNvuNFgPvNp1eqFfptvbAWBjkqXttMLGVptTY66l/As6+3W0z63trpJVwGrgYfrws9HOTd8KPFtVv9+1aEHt04n6XKD7dCDJkjb9JuCf07nG8SBwdRs2dp+O7uurgQfa0dhE2zCXfX6r64+A0Lmu0b1P+//7NFtXos+FB52r9N+mc17x0/Pcy7vo3MHwBPD0aD90zmHeDxwB/g9wUb1+R8LnW+9PAYNz3N8ddE4H/A2dc5PbZ9Ib8Ot0LsgNAx/tU59fbn08SecX69Ku8Z9ufT4HXNGvnw3gl+mcAnoSeLw9rlxo+3SSPhfiPv0F4LHW02Hgd7p+tx5u++d/ARe2+hvb/HBb/q6zbcMc9/lA26eHgf/J63cczcu/vV9HIUk6r04TSZImYBhIkgwDSZJhIEnCMJAkYRhIkjAMJEnA/wc2IVWuB7eIWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins=30\n",
    "hist_res = plt.hist(np.array(vertex_degrees), bins=nbins)\n",
    "hist_counts = hist_res[0]\n",
    "hist_breaks = hist_res[1]\n",
    "kvals = 0.5*(hist_breaks[0:(nbins)]+hist_breaks[1:(nbins+1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print the k values of the bin centers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  60.98333333,  180.95      ,  300.91666667,  420.88333333,\n",
       "        540.85      ,  660.81666667,  780.78333333,  900.75      ,\n",
       "       1020.71666667, 1140.68333333, 1260.65      , 1380.61666667,\n",
       "       1500.58333333, 1620.55      , 1740.51666667, 1860.48333333,\n",
       "       1980.45      , 2100.41666667, 2220.38333333, 2340.35      ,\n",
       "       2460.31666667, 2580.28333333, 2700.25      , 2820.21666667,\n",
       "       2940.18333333, 3060.15      , 3180.11666667, 3300.08333333,\n",
       "       3420.05      , 3540.01666667])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kvals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print the histogram bin counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.4898e+04, 1.3460e+03, 4.0800e+02, 1.7200e+02, 7.8000e+01,\n",
       "       3.5000e+01, 2.0000e+01, 2.5000e+01, 5.0000e+00, 9.0000e+00,\n",
       "       5.0000e+00, 1.0000e+00, 2.0000e+00, 2.0000e+00, 0.0000e+00,\n",
       "       4.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00, 3.0000e+00,\n",
       "       1.0000e+00, 1.0000e+00, 2.0000e+00, 0.0000e+00, 0.0000e+00,\n",
       "       0.0000e+00, 0.0000e+00, 1.0000e+00, 0.0000e+00, 1.0000e+00])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hist_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 7: Plot `N(k)` vs. `k`, on log-log scale (using only the first 14 points, which is plenty sufficient to see the approximatey scale-free degree distribution and where it becomes exponentially suppressed at high `k`.  For this you'll use `plt.loglog`.  You'll probably want to adjust the x-axis limits using `plt.gca().set_xlim()`.  To see the plot, you'll have to do `plt.show()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAASSklEQVR4nO3dUYhc133H8e8/iuwuLtG2sShopSA5EkvU6kHu1H5wStM0QRKJYlWBxmoSSCosXKqQhyIiQSilLchB5CGmTlJRGxFIZUQQqlQ7bEIcIwguaBW1lW0hIgSttXqQXFdKExZsq/8+7O71arNazezO3Ttn5/uBfZhzZ2b/kjz+cc//zDmRmUiSBPCepguQJPUOQ0GSVDEUJEkVQ0GSVDEUJEkVQ0GSVHlv0wUsxP33359r165tugxJKsrZs2ffyMyVs10rOhTWrl3L6Oho02VIUlEi4j/vdM3pI0lSpchQiIjtEXH45s2bTZciSUtKkaGQmacyc8+KFSuaLkWSlpQiQ0GSVI+iG80LdeLcGIdGLnL1xjirBgfYt2WYHZuHmi5LkhrTt6Fw4twYB46fZ/ztWwCM3RjnwPHzAAaDpL7Vt9NHh0YuVoEwZfztWxwaudhQRZLUvL4Nhas3xjsal6R+0LehsGpwoKNxSeoHfRsK+7YMM7B82W1jA8uXsW/LcEMVSVLz+rbRPNVMdvWRJL2rb0MBJoLBEJCkd/XU9FFE3BcRoxHxyaZrkaR+VGsoRMSzEXEtIl6ZMb41Ii5GxKWI2D/t0leAY3XWJEm6s7rvFI4AW6cPRMQy4GlgG7AR2BURGyPi48BrwLWaa5Ik3UGtPYXMPB0Ra2cMPwRcyszLABHxHPAo8OvAfUwExXhEvJCZ/1dnfZKk2zXRaB4CXp/2+ArwcGbuBYiILwBv3CkQImIPsAfgAx/4QL2VSlKf6alGM0BmHsnMf5nj+uHMbGVma+XKWU+TkyTNUxOhMAasmfZ49eRY2zxkR5Lq0UQonAE2RMS6iLgHeAw42ckbeMiOJNWj7iWpR4GXgeGIuBIRuzPzHWAvMAJcAI5l5qsdvq93CpJUg8jMpmuYt1arlaOjo02XIUlFiYizmdma7VrPNZolSc0pMhScPpKkehQZCjaaJakeRYaCJKkeRYaC00eSVI8iQ8HpI0mqR5GhIEmqh6EgSaoUGQr2FCSpHkWGgj0FSapHkaEgSaqHoSBJqhQZCvYUJKkeRYaCPQVJqkeRoSBJqoehIEmqGAqSpIqhIEmqFBkKrj6SpHoUGQquPpKkehQZCpKkehgKkqSKoSBJqhgKkqSKoSBJqhgKkqSKoSBJqhQZCn55TZLqUWQo+OU1SapHkaEgSaqHoSBJqhgKkqTKe5suQAI4cW6MQyMXuXpjnFWDA+zbMsyOzUNNlyX1HUNBjTtxbowDx88z/vYtAMZujHPg+HkAg0FaZE4fqXGHRi5WgTBl/O1bHBq52FBFUv8yFNS4qzfGOxqXVB9DQY1bNTjQ0bik+hgKaty+LcMMLF9229jA8mXs2zLcUEVS/+qZRnNEfAj4MnA/8KPM/FbDJWmRTDWTXX0kNS8ys743j3gW+CRwLTN/Z9r4VuAbwDLgHzPzyWnX3gN8JzM/d7f3b7VaOTo62v3CJWkJi4izmdma7Vrd00dHgK0zilkGPA1sAzYCuyJi4+S1TwHPAy/UXJckaRa1hkJmngbenDH8EHApMy9n5lvAc8Cjk88/mZnbgM/e6T0jYk9EjEbE6PXr1+sqXZL6UhM9hSHg9WmPrwAPR8RHgJ3Avcxxp5CZh4HDMDF9VF+ZktR/eqbRnJkvAS+189yI2A5sX79+fZ0lSVLfaWJJ6hiwZtrj1ZNjbfM8BUmqRxOhcAbYEBHrIuIe4DHgZAN1SJJmqDUUIuIo8DIwHBFXImJ3Zr4D7AVGgAvAscx8tcP39ThOSapBrd9TqJvfU5CkzjX5PYVaeKcgSfUoMhRsNEtSPYoMBUlSPQwFSVKlyFCwpyBJ9SgyFOwpSFI9igwFSVI9igwFp48kqR5FhoLTR+oVJ86N8ciTL7Ju//M88uSLnDjX0TZeUs/pmV1SpdKcODfGgePnGX/7FgBjN8Y5cPw8gEeJqlhF3ilIveDQyMUqEKaMv32LQyMXG6pIWjhDQZqnqzfGOxqXSlBkKNhoVi9YNTjQ0bhUgiJDwUazesG+LcMMLF9229jA8mXs2zLcUEXSwtloluZpqpl8aOQiV2+Ms2pwgH1bhm0yq2iGgrQAOzYPGQJaUoqcPpIk1cNQkCRVigwFVx9JUj2KDAVXH0lSPYoMBUlSPQwFSVLFUJAkVQwFSVLFUJAkVQwFSVKlyG0uImI7sH39+vVNlyI14sS5MfdcUi2KvFPwewrqZ1Mnvo3dGCd598Q3jwJVN3QUChHxGxHx2xHxQEQUGShS6TzxTXW66/RRRKwA/gLYBdwDXAd+DfitiPhX4JuZ+eNaq5RU8cQ31amdnsL3gO8Av5+ZN6ZfiIjfBT4fEQ9k5jN1FCjpdqsGBxibJQA88U3dcNdQyMyPz3HtLHC2qxVJmtO+LcMcOH7+tikkT3xTt7TdF4iIv5nxeFlEfLf7JUmay47NQxzcuYmhwQECGBoc4ODOTa4+Uld0siR1TUQcyMyDEXEvcAw4V1NdkubgiW+qSycriP4M2BQRB4BTwI8z869rqUqS1Ih2Vh89OO3hN4B/AH4CnI6IBzPzp3UVJ0laXO1MH319xuP/ATZOjifw0W4XJUlqRjurj/5wMQoBiIgdwCeA9wHPZOYPFut3S5La6ClExOciIua4/sGI+PAc15+NiGsR8cqM8a0RcTEiLkXEfoDMPJGZjwNPAJ9p/48hSeqGdqaP3g/8W0RMfSdh6hvN64E/AN4A9s/x+iPA3zPxBThgYjkr8DTwceAKcCYiTmbma5NP+erkdUnSIrrrnUJmfgN4EDgKrAT+aPLxGPD5zPx0Zv5sjtefBt6cMfwQcCkzL2fmW8BzwKMx4WvA9+/UwI6IPRExGhGj169fb+OPKElqV1vfU8jMW8APJ3+6YQh4fdrjK8DDwJeAjwErImJ9Zn57lloOA4cBWq1WdqkeSRLtLUn9qzkuZ2b+bbeKycyngKe69X6SOudZDf2tnTuFX84ydh+wm4l+w3xCYQxYM+3x6smxtnjIjlSPqbMapvZVmjqrATAY+kQ7PYWvT/0wMW0zAHyRiT7AA/P8vWeADRGxLiLuAR4DTrb7Yg/ZkerhWQ1qa5uLiPjNiPg74D+YuLt4MDO/kpnX2njtUeBlYDgirkTE7sx8B9gLjAAXgGOZ+Wq7RUfE9og4fPPmzXZfIqkNntWgdnoKh4CdTNwlbMrMX3TyCzJz1x3GXwBe6OS9pr32FHCq1Wo9Pp/XS5qdZzWonTuFvwRWMfHdgasR8fPJn/+NiJ/XW56kxbRvyzADy5fdNuZZDf2lnW0ueu4sZhvNUj2mmsmuPupfkVnuUv9Wq5Wjo6NNlyFJRYmIs5nZmu1az90FSJKaU2QouPpIkupRZCj4PQVJqkeRoSBJqoehIEmqFBkK9hQkqR5FhoI9BUmqR5GhIEmqh6EgSaq0dfJar3GbC6n3eVhPmYq8U7CnIPW2qcN6xm6Mk7x7WM+Jc22fpaWGFBkKknqbh/WUy1CQ1HUe1lMuQ0FS193pUB4P6+l9hoKkrvOwnnK5+khS13lYT7k8ZEfSkuAS2PbNdchOkXcKkjTd1BLYqRVPU0tgAYOhQ/YUJBXPJbDdYyhIKp5LYLvH6SNJPafT/sCqwQHGZgkAl8B2zjsFST1lPltkuAS2ewwFST1lPv2BHZuHOLhzE0ODAwQwNDjAwZ2bbDLPQ5HTR35PQVq65tsf2LF5yBDogiLvFNwlVVq63CKjWUWGgqSly/5As4qcPpK0dLlFRrMMBUk9x/5Ac5w+kiRVDAVJUsVQkCRVDAVJUsVQkCRVDAVJUqVnQiEiHoiIZyLie03XIkn9qtZQiIhnI+JaRLwyY3xrRFyMiEsRsR8gMy9n5u4665Ekza3uO4UjwNbpAxGxDHga2AZsBHZFxMaa65AktaHWbzRn5umIWDtj+CHgUmZeBoiI54BHgdfqrEWS5qPTA39K10RPYQh4fdrjK8BQRLw/Ir4NbI6IA3d6cUTsiYjRiBi9fv163bVK6mPzOfCndD3TaM7M/87MJzLzg5l5cI7nHc7MVma2Vq5cuZglSuoz8znwp3RNhMIYsGba49WTY22LiO0RcfjmzZtdLUySppvvgT8layIUzgAbImJdRNwDPAac7OQNPGRH0mLoxwN/6l6SehR4GRiOiCsRsTsz3wH2AiPABeBYZr7a4ft6pyCpdv144E9kZtM1zFur1crR0dGmy5C0hC3F1UcRcTYzW7Nd85AdSZpDvx340zOrjzrh9JEk1aPIULDRLEn1KDIUJEn1KDIUnD6SpHoU2WjOzFPAqVar9XjTtUgqUzdXFS2lFUpFhoIkLcTUnkZTW1hM7WkEdPw/826+Vy8ocvpIkhaim3saLbX9kYoMBXsKkhaim3saLbX9kYoMBZekSlqIbu5ptNT2RyoyFCRpIbq5p1Gv7Y904twYjzz5Iuv2P88jT77Y8dkPNpol9Z2pBnA3Vgx1870WqhtNbzfEk6Ql4pEnX2Rsll7G0OAAP9n/0erxXBviFTl9ZKNZkn5VN5reRYaCjWZJ+lXdaHoXGQqSpF/Vjaa3jWZJWiK60fQ2FCRpCVnooUBOH0mSKkWGgquPJKkeRYaCq48kqR5FhoIkqR6GgiSpYihIkiouSZWkhs11nOdiH/VpKEhSg+ba2RRY9KM+DQVJatDdjvO80zVDYZqI2A5sX79+fdOlSNKCzGdn0zqP+iyy0ez3FCQtFXPtbNrEUZ9FhoIkLRVz7WzaxFGfRU4fSdJS0c7Opou5+sjjOCWpzyy54zglSfUwFCRJFUNBklQxFCRJFUNBklQxFCRJlZ75nkJE3Ad8E3gLeCkzv9twSZLUFYu90+lC1HqnEBHPRsS1iHhlxvjWiLgYEZciYv/k8E7ge5n5OPCpOuuSpMUytQvq2I1xknd3Oj1xbqzp0mZV9/TREWDr9IGIWAY8DWwDNgK7ImIjsBp4ffJpt28LKEmFutsuqL2m1lDIzNPAmzOGHwIuZeblzHwLeA54FLjCRDDMWVdE7ImI0YgYvX79eh1lS1LXzGcX1CY10Wge4t07ApgIgyHgOPDpiPgWcOpOL87Mw5nZyszWypUr661UkhaoiZ1OF6JnVh9l5i8z84uZ+ed3azJHxPaIOHzz5s3FKk+S5qWJnU4XoolQGAPWTHu8enKsbZ6nIKkUOzYPcXDnJoYGBwhgaHCAgzs39ezqoyaWpJ4BNkTEOibC4DHgTxuoQ5IWxY7NQz0bAjPVvST1KPAyMBwRVyJid2a+A+wFRoALwLHMfLXD93X6SJJq4HkKktRnltx5Ct4pSFI9igwFG82SVI8iQ0GSVI+iewoRcRP42QLfZgUwn3moTl/X7vPbfd79wBsd/P6lZL7/ZnVarJq6/Xu68X5+hsqzITNnn2rJzGJ/gMNNvUenr2v3+R08b7Tpv/+S/91Lranbv8fPUPP/7TTxM9ffUenTR3fcDmMR3qPT17X7/G78mZa6Xvw7Wqyauv17/Az1pzv+HRU9fdTPImI077CkTNLd+RmaXel3Cv3scNMFSIXzMzQL7xQkSRXvFCRJFUNBklQxFCRJlSa2zlYNImIH8AngfcAzmfmDhkuSihIRHwK+zMSX2n6Umd9quKRGeKfQwyLi2Yi4FhGvzBjfGhEXI+JSROwHyMwTmfk48ATwmSbqlXpNh5+hC5n5BPAnwCNN1NsLDIXedgTYOn0gIpYBTwPbgI3ArojYOO0pX528LqnDz1BEfAp4HnhhccvsHYZCD8vM08CbM4YfAi5l5uXMfAt4Dng0JnwN+H5m/nSxa5V6USefocnnn8zMbcBnF7fS3mFPoTxDwOvTHl8BHga+BHwMWBER6zPz200UJxVg1s9QRHwE2AncSx/fKRgKS0RmPgU81XQdUqky8yXgpYbLaJzTR+UZA9ZMe7x6ckxSe/wMzcFQKM8ZYENErIuIe4DHgJMN1ySVxM/QHAyFHhYRR4GXgeGIuBIRuzPzHWAvMAJcAI5l5qtN1in1Kj9DnXNDPElSxTsFSVLFUJAkVQwFSVLFUJAkVQwFSVLFUJAkVQwFqYsiYu3MbZqlkhgKkqSKoSDVJCIeiIhzEfF7TdcitctdUqUaRMQwE/v0fyEz/73peqR2GQpS960E/hnYmZmvNV2M1Amnj6Tuuwn8F/DhpguROuWdgtR9bwF/DIxExC8y85+aLkhql6Eg1SAzfxkRnwR+OBkM7tevIrh1tiSpYk9BklQxFCRJFUNBklQxFCRJFUNBklQxFCRJFUNBklQxFCRJlf8HuhDg+2seZDkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "inds = np.where(hist_counts > 0)\n",
    "plt.loglog(kvals[inds],\n",
    "           hist_counts[inds], \"o\")\n",
    "plt.xlabel(\"k\")\n",
    "plt.ylabel(\"N(k)\")\n",
    "plt.gca().set_xlim([50, 2000])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 8:  Do a linear fit to the log10(N(k)) vs. log10(k) data (just over the bins where the count is nonzero).  You'll want to use `scipy.stats.linregress` to do the linear regression.  Don't forget to log10-transform the data using `np.log10`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinregressResult(slope=-2.562269147407539, intercept=8.776071095837441, rvalue=-0.9735828277703321, pvalue=2.6789389286928163e-14, stderr=0.13437166140568355)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.linregress(np.log10(kvals[inds]), np.log10(hist_counts[inds]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Slope is -2.6 with SE 0.13, i.e., gamma = 2.6 with a 95% CI of about +/- 0.26.\n",
    "\n",
    "Now let's compute the slope for the degree distribution Fig. 1b in the Jeong et al. article, for the yeast PPI.  The change in ordinate over the linear range is about -6.5 in units of natural logarithm. The change in abscissa over the linear range is approximately log(45)-log(2), so we can compute the Jeong et al. slope thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.09\n"
     ]
    }
   ],
   "source": [
    "jeong_slope = -6.5/(np.log(45)-np.log(2))\n",
    "print(\"%.2f\" % jeong_slope)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How close was your slope from the human PPI, to the slope for the yeast PPI from the Jeong et al. article?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll do the same thing in just a few lines of igraph code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = igraph.Graph.TupleList(interac_ppi_unique.values.tolist(), directed=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXzU1b3/8deZmexAAgl7AmFfRBEJm7hVsa6oV62KrbdaKtdWvVZve2+tdrm32nrbW3+3VK6WCqLWota6gIogboAghlV2CGFHEsISEsg+5/dHwpCdSTKT7yzv5+Ph4zHf813mk6/DJyef75lzjLUWERGJLC6nAxARkcBTchcRiUBK7iIiEUjJXUQkAim5i4hEICV3EZEI5HE6AIC0tDSbmZnpdBgiImFl9erVBdbaro3tC4nknpmZyapVq5wOQ0QkrBhj9jS1T2UZEZEIpOQuIhKBlNxFRCKQkruISARSchcRiUBK7iIiESisk/uugpMs3XEYr1fTFouI1BYS49xb6+UVe5j9+S4yuiRwx5g+fGt0Ot06xTsdloiI40woLNaRlZVlW/MlptKKKhZuOsTcL/fyRe5R3C7DFUO7MWVsHy4Z3BW3ywQhWhGR0GCMWW2tzWpsX8B77saYYcBDQBrwkbX22UC/x2nxMW5uPL83N57fm9zDxbyWvY83Vu9n0eY8eiXHc9uYDG7LyqBXSkKwQhARCUl+9dyNMbOB64F8a+2IWu1XA38E3MDz1tqnau1zAS9Za79ztuu3tufemPJKL4u35DH3y70s3VGAy8BlQ7pxx5gMLh/aDY+7/R4zbDtUxAuf7yIx1sPPrx+GMfpLQkQCJxA99znAM8BLtS7qBmYAVwL7gWxjzDxr7WZjzA3AD4CX2xJ4a8R6XFx7bk+uPbcn+46e4rXsfby+ah/TtubTrWMct2VlcPuYDDK6JAbl/b1ey2c7DjN72S6W7ijA4zJUei2ZaYn884TMoLyniEh9ftfcjTGZwLune+7GmAnAr6y1V9VsPwpgrf1trXPes9Zed7ZrB7Ln3pjKKi8fb83n1ex9fLotHwtcNDCNKWP7MGlYd2I9be/Nl1ZU8dbaA8xatouc/GK6dYzjuxdmMmVsH/7t9XV8nnOEt+6/kHN6Jbf9BxIRofmee1uS+63A1dba79ds3wWMA94AbgbigK+stTOauN40YBpAnz59Ru/Z0+TkZgF18HgJr6/ax+vZ+zhYWEpqUiy3jk7n9jEZ9O/aocXXyy8q5a8r9vDXlXs5erKc4T078f2L+3H9eb18vzSOFJdx7fSlJMV6mP/gRSTFhfUgJREJEe2a3K21D7Q0wGD33BtT5bUs2X6YuV/u5aOt+VR5LeP7d2HK2D5cdU4P4mPczZ6/9dAJZi3dxTvrDlLh9XLF0G5Mvag/4/t3abS2vmLnEe58/gtuHpXOH24bGawfS0SiSLBGyxwAMmptp9e0tSSwycDkgQMHtiGM1nG7DN8Y2o1vDO1G/olS/r56P69m7+WhV9eRkhjDzaPSmTI2g0HdO/rOOV1Pn7V0F8tyCkiIcXP7mAzumZh51l7/hAGpPHj5IKZ/tIOJA1O5+YL0YP+IIhLF2tJz9wDbgSuoTurZwJ3W2k0tDcKJnntjvF7L8p1HmPvlXhZtPkRFlSWrb2fuGNuH8kovs5blsvPwSbp3qq6n3zm2DymJsX5fv7LKy53Pr2TjgULeffCiVpWBREROa3NZxhgzF7iM6rHrecAvrbWzjDHXAv9L9VDI2dbaJ1sTYKgk99oKisv4x+r9vJq9j10FJwEY0bsT37+oP9ee27PVD2G/Lizh2j8upWdyAm/dfyFxnubLPyIiTQlIzT2YQjG5n2atZdWeY7gMXNCnc0DGqn+0JY+pL67i7gsz+dUN5wQgShGJRs0ld0cnDjPGTDbGzCwsLHQyjGYZYxiT2YXRfRt/UNoaVwzrzvcm9mPO8t0s2nQoINcUEanN0eRurZ1vrZ2WnBx9Y7//45ohnNs7mZ+88RUHjpc4HY6IRJiwnvI3nMV53PxpyiiqvJaH5q6lssrrdEgiEkFUlnFQZloST/7TCFbtOcb/Lt7hdDgiEkFUlnHYjef35rasdGZ8msPnOQVOhyMiEUJlmRDwqxvOYUDXDvzotXUcLipzOhwRiQBK7iEgMdbDM3eO4kRJBVNfzKbwVIXTIYlImFNyDxFDe3Ti2e9cwNavi7hr9koleBFpEz1QDSGXD+2uBC8iAaEHqiHmimFK8CLSdirLhCAleBFpKyX3EKUELyJtoeQewpTgRaS19EA1xCnBi0hraMrfMPHRljx+8Nc1dEmK5eoRPfjmOd0Zm9kFj1t/fIlEK83nHiFW7DzCrGW5LN1RQFmll+SEGK4Y2o1vntOdSwZ3JTFWC2+LRJNgraEq7WzCgFQmDEjlVHklS7YXsGjzIT7ems+baw8Q63Fxz8RM/v2qobhdgZl3XkTCl5J7GEqM9XD1iB5cPaIHlVVesncf4++r9vHnz3LZfqiI6VNG0TE+xukwRcRBKtiGOY/bxYQBqTx9+/k8cdMIluwo4Ob/W87eI6ecDk1EHKTRMhHkO+P78vL3xpJfVMaNM5bxRe4Rp0MSEYfogWoE2l1wkqkvZrPnyCl+dcM5XDakK1VeS6XXUuW1lJRXkVtQzLZDxWzPKyInv5hJw7rz8+uHBWydWBEJPj1QjTKZaUm8+cOJPDh3LY+/vbHJ42LchgFdO9ArJZ7Zn++ic2IMD14xqB0jFZFgUXKPUMkJMcz+bhYLN+VxsrwSj8vgrvkv1u2iX1oSmWlJxLhdWGv5t9fX84cPt9MnNZEbz+/tdPgi0kZK7hHM43Zx3Xk9z3qcMYbf3nIu+4+V8JM3viK9cwKj+3ZphwhFJFg0WkYAiPO4+fNdo+mVHM+9L61mz5GTTockIm2g5C4+nZNimX33GLzW8p1ZKzl4vMTpkESklZTcpY7+XTvw4j1jOX6ygm8/v5L8E6VOhyQiraBx7tLAyIwU5nxvDPknSpnyly84XFTmdEgi0kIa5y5NWpl7hLtfyKZHcjzXntuD4T2TGdKjIymJMcTHuIn3uDQrpYiDNM5dWmVc/1Rm3Z3Fr9/dwnOf5VLlrdsRiPW4+N0t53HTKA2dFAk1Su7SrAsHpLHgoYspq6xiR14x2w4VUVxWSWlFFYs25/GTN9bTKyWBsf00dFIklKgsI61WeKqCf3r2c46dLOetH04kMy3J6ZBEokpzZRkVTKXVkhNjmP3dMVjgnjnZ5OQXAWCt5fmludw1ayUl5VVnvc6RYj2wFQk0JXdpk8y0JJ7/5ywKSyq4bvoynl+aywN/W8sT721h6Y4C/rFmf7Pnf7g5j3G/+Yh9RzVFsUggKblLm2VlduGDH13MxIFpPPHeFhZs/JpHrxnKyPRkZi3bhdfbdOnv4615VHot6/cfb8eIRSKfHqhKQHTrGM+s72bx7ldf0yM5njGZXeiVksCDc9eyeEseVw7vzhe5R0mMdTOkR0fiY9xA9bqwANsPFcF5Tv4EIpFFyV0CxhjD5JG9fNvXjOhB75QEZny6k9ey9/HR1nwAYt0upk8ZxciMZHbXrBi19VAR1loefm0dFw3qyq2j0x35GUQihcoyEjQet4upF/Vj/b7jLMsp4PHrhvHsty+gb2oiTy3YwrIdBQD0S0tiW14R+46W8Pa6g/zHP75ieU6Bw9GLhDdHe+7GmMnA5IEDBzoZhgTRlLF9KCypYPLIngzs1hEAY+C+v67hD4u2k5wQww0jezH94x18tDUPgG4d4/jh39bw6Y8vIyUx1snwRcKWoz13a+18a+205ORkJ8OQIEqIdfPwlYN9iR3gm8N7MKhbBw6dKGVcvy4M69kJa+HlL/aQ1iGOGd++gOOnKli0Oc/ByEXCm8oy0u5cLsP936j+a23CgFSG9qhO/LmHTzK+fxdGZaSQ3jmBBRu+djJMkbCmB6riiMkje1Fe5eW6c3uSEOMmPsZFaYWXCQNSMcZw7bk9eeHzXRSWVJCcEON0uCJhRz13cYTbZbgtK4OkOA8ul2Fw9+re+4T+qUD1SJuKKsuHKs2ItIqSu4SEC/p0JjM1kX4189Ocn5FCZmoij7+9gRmf5BAKcyCJhBMldwkJj147lHceuAhjDFA9Zv6Ve8dzyaCu/H7hNj7bfrjBOX9cvIN56w+2d6giYUHJXUJCnMfdoLbeOyWBZ+68gF7J8fzp47q999KKKp75ZAdPvLuZiipve4crEvKU3CWkxXpc3HfZAFbvOcYXuUd97ZsOFlJRZckvKuODjYccjFAkNCm5S8i7LSuDTvEe3qw1w+SaPdUTjXXrGMcrK/c4FZpIyFJyl5AXH+Pm4sFd+Wz7YV9pZvWeY2R0SeDmC9JZveeYX/PGi0QTJXcJC5cO7kp+UZlvgrHVe48xuk9nxvXvQkWVZc3eY06HKBJSlNwlLFw6uCsAn247zD/WHOBwURnj+qeS1bczLgMrc484HKFIaNE3VCUsdO8Uz7CenZi5ZCcny6q4cEAqt45OJ8btYkTvZBZvyefSIV0Z3VcLdYuAeu4SRv718oGcm57Cdef1ZMadFxDjrv74XjKoK5u/PsGtz63gUGGpw1GKhAb13CVsXHNuT645t2eD9h9NGsSI3snc99fVrNx1hBvP7+1AdCKhJSg9d2PMTcaYvxhjXjPGfDMY7yFymsftYtKwbnSI85C9+yh7j5ziqQVbOVlWyf5j1a+rmlnHVSQS+d1zN8bMBq4H8q21I2q1Xw38EXADz1trn7LWvg28bYzpDPwPsCiwYYvU5XG7uKBvZ/62ci+vfrmPSq9lQNckdh85yXOf7eSGkb0Y3quT02GKtJuW9NznAFfXbjDGuIEZwDXAcGCKMWZ4rUMer9kvEnRjMzvjtfgW3/5wcx4ra77VuiO/yMnQRNqd3z13a+0SY0xmveaxQI61NhfAGPMqcKMxZgvwFLDAWrsmQLGKNOu2MRkUllQw9aL+zPgkh5e/OPPN1Zz8YnLyiyit8DKit1b+ksjX1pp7b2Bfre39NW0PApOAW40x9zV2ojFmmjFmlTFm1eHDDWf8E2mpbh3jeey64fRIjuemUb187cbAjrxiJj29hOv/tMzBCEXaT1BGy1hrpwPTz3LMTGAmQFZWlp52SUCN7tuFt++fyMdb89l0oJANBwp9+yqrvHjcGgUska2tn/ADQEat7fSaNhHHnZ+RwiNXDmZ4r04cOF7ia6/9WiRStTW5ZwODjDH9jDGxwB3APH9PNsZMNsbMLCwsPPvBIq10+5iMOtu5h086FIlI+/E7uRtj5gIrgCHGmP3GmKnW2krgAWAhsAV43Vq7yd9rWmvnW2unJSfrAZcET3rnRObcM4bvTewHwPSPd3CkuAyAkvIq7nt5NZsOqoMhkcXv5G6tnWKt7WmtjbHWpltrZ9W0v2+tHWytHWCtfTJ4oYq03mVDuvHz64cBsHbvcZ58fwsAb609wAebDvGXJblOhicScI4+VVJZRtqTMYYnbqr+/t2n2w5TWeXlg03Vqzit2nOM332wlc9zCvjt+1u0ILeEPRMKH+KsrCy7atUqp8OQKLFgw9f84JU1PH7dMJ54b0ujx/zxjvO58fzePDh3LcdPlfPy1HHtHKXI2RljVltrsxrbp/FgEnXG9queFviJ97bQMd7D0B4dGxyzYEN1j37++oMs3VHQrvGJBIKSu0Sd1A5xdOsYB8C4fqn0S0tqcMzuIydVmpGwppq7RKVeKQkAjOqTgjHVbZcNqV7tqX9aErsKTpJ3osx3/NId+ha1hBdHk7uGQopTKqq8AJzTqxNlFdWvb8/KYOm/f4OpF/ejrNLLR1vzfMffNetLLQQiYUVlGYlKv7h+OEO6dyQrswtDamru6Z0TyeiSyMCuHQB47K2Ndc5ZsPHrRq+Vd6LU98tCJFQouUtUGtc/lYUPX0KHOA8PXzmY1/9lAuemV/8FOSazi69EU9v6fcc5XFTG3iOnfG2FJRWM+81HPNnEqBsRp6jmLlEvxu3yjaABcLkMD3xjYIPj1u47zpgnF/PtWV/42gpqvun66bb84Acq0gKquYs0YljPhqs27anpse87embiseOnygFIjNVyxBJa9IkUaURSnIdfTh7O6L6d6ZeWxIvLd/M/i7b79nu9FpfLsDznSM3xbqdCFWmUau4iTbhnYj/OS0+hY3wMQ3rU7ck//s5Gyiu9/OHD6oSfoJ67hBgldxE/XDG0G49dO4xJw7oD8LeVezlVXunb73GZBl96qqzyUuXVF6HEGXqgKuIHl8tw7yX9uWtCX1/b3qNnRs18vDWfOct31zln4GMLuOXZ5e0VokgdeqAq0gId48+UX2545vM6+95c03ARsnX7jgc9JpHGqCwj0gK9khOa3Odxm3aMRKR5Su4iLdAjOZ6N/3lVo/s8LiV3CR1K7iIt1CHOw8OTBjdoz959jD1HqtdnrdR0BOIwJXeRVnho0qBG25/5OAeAUxVV7RmOSAMaLSPSRiMzUnyvuyTFAtULb4s4SaNlRNrol5OHE+uu/qf0dWEpR0+Wc0rJXRymsoxIG8W6XSx+5FIA5q0/yAW//pDi0sqznCUSXEruIm0UH+OiT2pinbVYd+QX+V7P+CSnwTlzv9zLvPUH2yU+iU5K7iJtFOuunjSsU0KMr+2r/WeeI/1+4bYG5zz65gb+de5aoHoSstO+9dxyfvfB1mCFKlFEyV2kjWI91f+MOsad+fZq/TVXC0sqGj13V8FJ+v/sfRZsqF7lKXv3Mf7v051BilSiiZK7SBudTu6pHWJ9bTsPn2RkejLfrZmL5sCx6jngT5RWsPHAmV59Tn4xAH9fvb+9wpUooeQu0kank/vkkb3qtA/q3pFbRqcD8MrKPQDc9twKrv/TMt8xnROrSzlHT5Y32bsXaQ2Ncxdpo7ia5H7xoK78/b4JvvbySi8ZnRMBeGXlXkrKq9h6qKjOufNrHqoePF7CyP9c1E4RSzTQOHeRNqo9p8zpnjjALaPTSam1veFAw07Miyuqe/T5RWVBjFCikcoyIq302LXDSEmMwZgzyT3Oc2a5vUsHd62z77Y/r2jX+CS6KbmLtNK9l/Rn3S++WactPqbhWqrTp4xqr5BEfJTcRQIoPqbhP6mERhK+SLApuYsEUGM998YSvkiw6VMnEkAxbvXcJTQouYsEQe1RMrV7849cOZhrz+1x1vMzf/peg2+5irSEkrtIgL06bTwfPHSJbzsh9kxyPz8jpc4Imua8s04Ti0nrec5+iIi0xPj+qXW2Y2uVaizg70qrjZV4RPylT49IkKV1iPO9HtazI8dP+TfNQKxbC25L66nnLhJkCbFudj91nW+7R3I8ACPTk1m/v+mpN9Rzl7bQ3DIi7exXN5zDuw9exNxp4xvsG9itg+91cVklP3trA6fKtaqTtJyjPXdr7XxgflZW1r1OxiHSnjrEeRjRu+F8SvddOoDCknLfNMCvZu8DoHdKAndfmIkF3MbUeUAr0hT93ScSIn5y1ZA6D19P+/3CbZz/X4sY8cuFTHjqIwcik3Ck5C4SItwuQ7+0pEb3VVRVL8Xn78NYESV3kRDQq+Yha98mkrtISym5i4SAV+6tfrhaex1WqK631zf8Fx/w1f7j7RKXhC8ld5EQcHq9jw7xdZN7x/iGYx5OlVfxwue7qfJaqry2zr6KKi/eem0SnZTcRUKAq2ZKgqTYusm8Q1zjA9rcLsO43ywm64kP67QPemwBP3ptXXCClLCi5C4SAlw1Xff6ybyxnjvAiZIKCorLOdbIA9Z56yN7TprCUxVUVHmdDiPkKbmLhIDTZZmkesm9sfnhARZtzgt2SCFr5H8t4oevrHE6jJCn5C4SAk6XZWI9Lj798WW+do+mIGjUh1H8y81f+uSIhIDaswBn1gyHTE2KJcbl/+Rh1upBqpyh5C4SAtz15nj/8OFL+PCRS/G0YGZIDZKR2pTcRUKAq15yH9S9I12SYv0uy5RWVDHgZ+/7tm95djkvrdgdwAgl3Ci5i4SA+sn9NH/LMruPnKyzvXrPMX7xzqY2xyXhK+DJ3RjT3xgzyxjzRqCvLRKpTBP/Ev3puX+yLZ9Fm5p/wHiitIKFmw61JrSgyjtRyrIdBU6HEZH8Su7GmNnGmHxjzMZ67VcbY7YZY3KMMT8FsNbmWmunBiNYkUhVv+Z+2g0je5313HteyObpD7c3e8xDc9fyLy+vZt/RU62KL1iu/9MyvjNrpdNhRCR/e+5zgKtrNxhj3MAM4BpgODDFGDM8oNGJRImmyjIjM1LY/dR1LHjo4jZdf8+R6qReVhlaX/45XFTmdAgRy6/kbq1dAhyt1zwWyKnpqZcDrwI3Bjg+kajQRG73ifW0roL647+vp8pryS2orsm/8PmuRo9596vI/lZrNGpLzb03sK/W9n6gtzEm1RjzHDDKGPNoUycbY6YZY1YZY1YdPny4DWGIhL+meu6nxbha90/1jdX7OXqy3Lf9ysq9jR7zwN/Wtur6EroCvsyetfYIcJ8fx80EZgJkZWVphK5EtbMNinHXjHfv0SmeQydKW3RtTwu+CCWRoy099wNARq3t9Jo2EfFT58QY4Ow997iaskxTE4k154qnP6uz/fSibfzb6+tbfB0JL21J7tnAIGNMP2NMLHAHMK8lFzDGTDbGzCwsLGxDGCLh660fTuTJfxrhmxWyKWkd4vjV5OG8cM+YFr9H7bIMwPSPc/jHmv0tvo6EF3+HQs4FVgBDjDH7jTFTrbWVwAPAQmAL8Lq1tkXfmrDWzrfWTktObrgSvEg0yExL4tvj+vp17N0T+5HeOTHIEUmk8He0zBRrbU9rbYy1Nt1aO6um/X1r7WBr7QBr7ZPBDVVEAmnGJzkNevWNsdYyc8lOjp86+7ESOhydfkBlGRHn/H7hNi749YdnPe6L3KP85v2t/OytDe0QlQSKo8ldZRmR0Hd61aOi0kqHI5GW0MRhIuIXTRcfXpTcRcJMp1YMh/TH80tz+cFfVzdoPz1K01I3u5dVVjHp6c9YnqOJv0KRau4iYeadBy7itqx0Jg5M5f5vDAjYdZ94bwsLNjacOdLQ+DDNPUdOkZNfzC/naWrhUBScLoCfrLXzgflZWVn3OhmHSDjpl5bE724d6due8clOB6ORUKWyjIg0y1eWqVdzVw0+tDnacxeR0PPSit0M7NYBlzGM75/qK8oomYcXR5O7MWYyMHngwIFOhiEitdRenm/3U9f5Xtd/oCqhTePcRSLQdyf4N6WBX5ooy/h2a9LJkKSau0gEOttEZC3R1GgZ9eRDm5K7SAQ62xTCraFUHl6U3EUikDuAPfdfzttY/aJedj/do9+eV8wnW/NZv+84mw42/Z2V4rJK5q3Xcn7tRQ9URcJceucE9h8rqdMWqI777oKTbM8rbnRf7bLMPXOyz5xT6yFsbY++uYH56w8yoGsS5/TSc7Zg0wNVkTC37D8ub5BQ3fWy+7dGp7fq2pVer+91W2vsXx+v/gV0qryqTdcR/6gsIxKBAlWWcbdyYW5xnv7PiUQgE6C6TO3FtfUlpvCi5C4SgeqXZVrL466V3ANyRWkvmn5AJAK5A9Rtq/9Loqi0goLicuI8riZ78jn5RXSIi8HlgtJyLz2S44n1qB/Z3jRaRiRCJMa6fQ8rA1WWqX0day1X/b8lHCwsBWDBQxc3es6kp5fU2f7W6HR+/62RjR4rwaPRMiIRIvuxSb7X9b/EFIiSigVfYm+Jz7WYhyP0t5JIhEiKO/OHeP2Ou7eVT0Nbe15tgforQlpGyV0kCrQ2R1d6z5zY2vncldudoeQuEgVsK7N7ZVXtLzG1jpK7M5TcRaJAaxNz7Z57ff4m7aZmlZTgUnIXiQLN5OhmVVbVOrFe79/fPwYCOIeZtICSu0gUaO2D0X+s2e97vX5/3RkfT5VX+nWNQydKKavUfDLtzdHkboyZbIyZWVjY9DShIuK/UX1SGm0f3z+1VdebuSS3yX23PrfCr2uUVni5/5W1rXp/aT2NcxeJIHPvHc+qx8+Md7/7wkwWP3Ipd40/s+ze3+4d1+5xLd6S1+7vGe00/YBIBImPcRMf4/Ztx3lcDOzWwfe6rNJLp/gYp8KTdqSau0iUCcYSfBJ6lNxFooymaI8O+t8sEiVOj5dRzz06KLmLRBmNO48OSu4iUUYTeUUHJXeRKKOyTHRQcheJFjVFd5VlooOSu0iU0URe0UHJXSRK9EyJd/T9V+05BsDO/GIARv3XokaPm/jUx8z4JAeA5TsLyPzpe1z+h0957K0N7RNohNDcMiJR4tVp45k+ZVSdxap/ctWQdo/jvQ1fA3DsVEWj+w8cL+H3C7cB8PzSXQDkHj7JKyv3tk+AEUJzy4hEiZ7JCdwwsledtksHd3UoGgk2lWVERCKQkruISARSchcRiUBK7iIiEUjJXUQkAim5i4hEICV3EZEIpOQuIhKBlNxFRCKQkruISARSchcRiUBK7iIiEUjJXUQkAim5i4hEIE+gL2iMSQL+DygHPrXWvhLo9xARkeb51XM3xsw2xuQbYzbWa7/aGLPNGJNjjPlpTfPNwBvW2nuBGwIcr4iI+MHfnvsc4BngpdMNxhg3MAO4EtgPZBtj5gHpwOn1sKoCFqmIRISlOwoY9Nj7ddoyf/peg+PO1jbosff5+fXD2XigkNdX7ad3SgIHjpdw+dBuHDxegrWwLa+IDnEebh+TwY8mDeIvS3K5//KBLM85gstl6J+WxJzlu+mZHE+/tCTiPG4uGpQGwPdfzGbxlnz+9YpBPHLlYN/7Wmt57rNcrhjWjR//fT078or58VVD+N7ETB5+bR3FZZVcf14vbhrVu0H8H27OIynWzYUD0zhcVMZv39/CneP6kJXZpdX3syl+JXdr7RJjTGa95rFAjrU2F8AY8ypwI9WJPh1YRzN/GRhjpgHTAPr06dPSuEUkwDrEeSguq/Rtj8xIYf2+40F5r4oqG5Br/OKdTb7tA8dLAPh4a36d44rLKpm1bBfz1x8kv6iM5MRYfv3uZgAyuiSw72hJneN3P3UdXq9l8Zbq60z/aEed5L4jv5j//mAr//3BVl/br9/dzKBuHXh73UEAFm/JbzS53/vSKt97nCqv5M21B7hoUFpQkntbHqj2BhHEyI8AAAQhSURBVPbV2t5f0/YmcIsx5llgflMnW2tnWmuzrLVZXbtqqS8RJ/VOSWDjf17F4kcuAaBLUiw/umKQw1EF1onS6jVbyyrPFBSKSyubOrxJVd7GfzFVVHlbF1iQBPyBqrX2JHBPoK8rIiL+a0vP/QCQUWs7vabNb8aYycaYmYWFhW0IQ0SCwdL20ok4py3JPRsYZIzpZ4yJBe4A5rXkAtba+dbaacnJyW0IQ0RE6vN3KORcYAUwxBiz3xgz1VpbCTwALAS2AK9bazc1dx0RCR9WHfew5u9omSlNtL8PvN/YPn8YYyYDkwcOHNjaS4iISCMcnX5AZRmR0KWee3jT3DIiIhFIyV0kgqnzHb2MDYG/vYwxh4E9QDJQf1xk/bb622lAQVADbPx9g3Xu2Y5tbn9T+6L9vvpzXEvuXVPtZ2uLpHvqz7H6rLbu3Jbc177W2sa/BWqtDZn/gJlna2tke5VTsQXj3LMd29z+pvZF+33157iW3Dt/72n9tki6p/4cq89q+9/X2v+FWlmmsekK6rc1OaVBkLXlfVty7tmObW5/U/ui/b76c1xL7l1T7f62BZs+q8ERDvfVJyTKMm1hjFllrc1yOo5Io/saeLqnwaH72rhQ67m3xkynA4hQuq+Bp3saHLqvjQj7nruIiDQUCT13ERGpR8ldRCQCKbmLiESgiEvuxpgkY8yLxpi/GGO+7XQ8kcIY098YM8sY84bTsUQKY8xNNZ/T14wx33Q6nkhhjBlmjHnOGPOGMeYHTsfjlLBI7saY2caYfGPMxnrtVxtjthljcowxP61pvhl4w1p7L3BDuwcbRlpyX621udbaqc5EGj5aeE/frvmc3gfc7kS84aKF93WLtfY+4DZgohPxhoKwSO7AHODq2g3GGDcwA7gGGA5MMcYMp3pFqNNru1YhzZmD//dV/DOHlt/Tx2v2S9Pm0IL7aoy5AXiPNkxJHu7CIrlba5cAR+s1jwVyanqU5cCrwI1UL9SdXnNMWPx8TmnhfRU/tOSemmr/DSyw1q5p71jDSUs/q9baedbaa4CoLc2Gc/LrzZkeOlQn9d7Am8Atxphnce5ryuGs0ftqjEk1xjwHjDLGPOpMaGGrqc/qg8Ak4FZjzH1OBBbmmvqsXmaMmW6M+TNR3HP3ayWmcGKtPQnc43QckcZae4Tq2rAEiLV2OjDd6TgijbX2U+BTh8NwXDj33A8AGbW202vapG10XwNP9zQ4dF+bEc7JPRsYZIzpZ4yJBe4A5jkcUyTQfQ083dPg0H1tRlgkd2PMXGAFMMQYs98YM9VaWwk8ACwEtgCvW2s3ORlnuNF9DTzd0+DQfW05TRwmIhKBwqLnLiIiLaPkLiISgZTcRUQikJK7iEgEUnIXEYlASu4iIhFIyV1EJAIpuYuIRCAldxGRCPT/AUTgdiCMJhpQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs, ys = zip(*[(left, count) for left, _, count in \n",
    "         g.degree_distribution().bins()])\n",
    "plt.loglog(xs, ys)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FittedPowerLaw(continuous=False, alpha=2.980799940666456, xmin=227.0, L=-5488.236776994901, D=0.033363389572659274, p=0.2818062734691993)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "igraph.statistics.power_law_fit(g.degree())"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}