class20_partialcorr_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS446/546 - Class Session 20 - Partial correlation network\n",
    "\n",
    "In this class session we will continue with analyzing the tumor gene expression dataset from the NIH human bladder cancer cohort (*M*=414 tumors), building on what we learned in Class Session 19 (Correlation Network).  In order to keep the analysis simple&ast;, in this notebook will restrict our analysis to a set of *N*=164 genes that are very highly expressed in bladder cancer. Using the 164 x 414 matrix of transcript abundance measurements, we will construct a network based on gene-gene partial correlation coefficients. We will also compare the distribution of partial correlation coefficients to the distribution of Pearson correlation coefficients. Do you think they will be different?  In what way would you expect them to be different?\n",
    "\n",
    "&ast; Here, \"simple\" means that the covariance matrix will be nonsingular, so that we can obtain the partial correlation matrix by inversion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll import all of the python modules that we will need for this exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas  ## data file loading\n",
    "import numpy\n",
    "import sklearn.covariance  ## for covariance matrix calculation\n",
    "import matplotlib.pyplot\n",
    "import matplotlib\n",
    "import pylab\n",
    "import scipy.stats  ## for calculating the CDF of normal distribution\n",
    "import igraph  ## for network visualization and finding components\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the tab-deliminted text file of gene expression measurements (rows correspond to genes, columns correspond to bladder tumor samples). (use Pandas, `pandas.read.csv`, and `as_matrix`). As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using `shape`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_for_network_df = pandas.read_csv(\"shared/bladder_cancer_genes_tcga.txt\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert your data frame to a numpy matrix, using the `pandas.DataFrame.as_matrix` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_for_network = gene_matrix_for_network_df.as_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using `shape`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4473, 414)\n"
     ]
    }
   ],
   "source": [
    "print(gene_matrix_for_network.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the median expression level for each row of your matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_median_expression = numpy.median(gene_matrix_for_network, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter the matrix to include only rows for which the gene's median expression is > 12; matrix should now be 164 x 414; this will enable us to easily compute the partial correlation matrix using the inverse of the covariance matrix. Print the size of the filtered matrix, as a sanity check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_np = numpy.array(gene_matrix_for_network)\n",
    "genes_keep = numpy.where(genes_median_expression > 12)\n",
    "matrix_filt = gene_matrix_np[genes_keep, ][0]\n",
    "matrix_filt.shape\n",
    "N = matrix_filt.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the shape of your filtered matrix, as a sanity check.  It should be 164x414."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(164, 414)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matrix_filt.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the 164 x 164 matrix of gene-gene Pearson correlation coefficients, using `numpy.corrcoef` (this function treats each row as a random variable, so you don't have to do any transposing of the matrix, unlike the situation in R). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_cor = numpy.corrcoef(matrix_filt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the covariance matrix using `sklearn.covariance.empirical_covariance` (from the `sklearn.covariance` package, . Make sure you take the transpose of the `matrix_filt` matrix before passing it to the `empirical_covariance` function!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_cov = sklearn.covariance.empirical_covariance(numpy.matrix.transpose(matrix_filt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `numpy.linalg.inv` to get the inverse matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_cov_inv = numpy.linalg.inv(matrix_cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a double for loop to \"scale\" the negative of the precision matrix, which will give you the partial correlation. Print the dimensions of the matrix you get back, as a sanity check. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(164, 164)\n"
     ]
    }
   ],
   "source": [
    "matrix_pcor = -matrix_cov_inv\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        matrix_pcor[i,j] /= numpy.sqrt(matrix_cov_inv[i,i]*matrix_cov_inv[j,j])\n",
    "print(matrix_pcor.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the correlation coefficients and the partial correlation coefficients of the lower triangle of the matrix (not including the diagonal), as two vectors `cor_values` and `pcor_values`; your resulting vectors should each have length 13,366.  You will want to use `numpy.tri` and `numpy.where` (see class session 19 exercise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13366\n",
      "13366\n"
     ]
    }
   ],
   "source": [
    "cor_values = matrix_cor[numpy.where(numpy.tri(*matrix_cor.shape, k=-1))]\n",
    "pcor_values = matrix_pcor[numpy.where(numpy.tri(*matrix_pcor.shape, k=-1))]\n",
    "print(len(cor_values))\n",
    "print(len(pcor_values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot the histograms of the correlation coefficients (upper triangle only) and the partial correlation coefficients, on the same plot using alpha blending (refer to class session 17 exercise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGWNJREFUeJzt3X+UFeWd5/H3J90oTIIaoDNh+NU4EZX4A5IGUQ+jgSHr\nuFnMbnSDESLqhlUjJpnJydG4x2SdnJ3jeiZuIp5hGDWo0cjoZNxOhhmDgnGIQGgNQfEHMsbFNmwk\njSKuYmzz3T+qqNy03XR103Xr9u3P65x7qB/PvfWhafrbT9VTTykiMDMzA3hP2QHMzKx2uCiYmVnG\nRcHMzDIuCmZmlnFRMDOzjIuCmZllXBTMzCzjomBmZhkXBTMzyzSWHaCvxowZE83NzWXHMDMbVB57\n7LFfR0RTb+0GXVFobm6mra2t7BhmZoOKpP+Tp51PH5mZWcZFwczMMi4KZmaWGXTXFLrz9ttv097e\nzv79+8uOUpjhw4czfvx4hg0bVnYUM6tjdVEU2tvbGTlyJM3NzUgqO86Aiwg6Ojpob29n8uTJZccx\nszpWF6eP9u/fz+jRo+uyIABIYvTo0XXdEzKz2lAXRQGo24JwQL3//cysNtRNUTAzs0NXF9cUurpx\nzfYB/bwvzZsyoJ9nZlar6rIoDCadnZ00NvqfIZd1f1W9Y33s6uody6yG+PTRALrjjjs46aSTOPnk\nk1m0aBEvvPACc+bM4aSTTmLu3Lns3LkTgMWLF3PppZdyyimn8JWvfKXk1GZmv+NfUQfItm3b+MY3\nvsGjjz7KmDFj2LNnDxdeeGH2uu2227jyyiu5//77gWQY7aOPPkpDQ0PJyc3MfqfwnoKkBkk/k/TD\nbvYdLmmVpB2SNklqLjpPUdauXct5553HmDFjABg1ahQbNmzgM5/5DACLFi1i/fr1WfvzzjvPBcHM\nak41Th99AXi6h32XAK9ExIeAG4Hrq5CnJrz3ve8tO4KZ2bsUWhQkjQf+PXBLD03OAW5Pl+8D5mqQ\nDsifM2cO9957Lx0dHQDs2bOH0047jXvuuQeAu+66i9mzZ5cZ0cysV0VfU/hfwFeAkT3sHwe8CBAR\nnZL2AqOBX1c2krQEWAIwceLEXg9axhDSD3/4w1xzzTWcccYZNDQ0MH36dG666SYuuugibrjhBpqa\nmvjOd75T9VxmZn1RWFGQ9Ang5Yh4TNKZh/JZEbECWAHQ0tISAxCvEAcuKldau3btu9qtXLmySonM\nzPqmyNNHpwPzJb0A3APMkfTdLm1eAiYASGoEjgQ6CsxkZmYHUVhRiIirI2J8RDQDC4C1EbGwS7NW\n4MCv1uembWq2J2BmVu+qfp+CpOuAtohoBW4F7pS0A9hDUjzMzKwkVSkKEfEw8HC6fG3F9v3AedXI\nYGZmvfM0F2ZmlnFRMDOzTH3OfTTQs2l6xkwzGyLcU6hRnZ2dZUcwsyHIRWGAvPDCCxx33HFccMEF\nHH/88Zx77rm88cYbbN68mdNOO42TTz6ZmTNnsm/fPvbv389FF13EiSeeyPTp01m3bh2Q3NQ2f/58\n5syZw9y5c0v+G5nZUFSfp49K8uyzz3Lrrbdy+umnc/HFF7Ns2TKWL1/OqlWrmDFjBq+99hojRozg\nW9/6FpJ44okneOaZZ/j4xz/O9u3J0+Ief/xxtm7dyqhRo0r+25jZUOSewgCaMGECp59+OgALFy7k\ngQceYOzYscyYMQOAI444gsbGRtavX8/Chcl9fMcddxyTJk3KisK8efNcEMysNC4KA6jrBK9HHHFE\nnz/DU2qbWZlcFAbQzp072bBhAwB33303s2bNYteuXWzevBmAffv20dnZyezZs7nrrrsA2L59Ozt3\n7uTYY48tLbeZ2QH1eU2hpCGkxx57LDfffDMXX3wxU6dOZenSpcyZM4elS5fy5ptvMmLECB588EEu\nv/xyLrvsMk488UQaGxtZuXIlhx9+eCmZzcwq1WdRKEljYyPf/e7vTwQ7Y8YMNm7c+K623T1bYfHi\nxSxevLioeGZmvfLpIzMzy7goDJDm5maefPLJsmOYmR2SuikK9f4Yhnr/+5lZbaiLojB8+HA6Ojrq\n9gdnRNDR0cHw4cPLjmJmda4uLjSPHz+e9vZ2du/eXXaUwgwfPpzx48eXHcPM6lxhRUHScOAR4PD0\nOPdFxNe6tFkM3EDyrGaAZRFxS1+PNWzYMCZPnnxogc3MrNCewlvAnIh4XdIwYL2kf46IruMzV0XE\nFQXmMDOznAorCpGc4H89XR2WvurzpL+ZWZ0o9EKzpAZJW4CXgTURsambZp+StFXSfZImFJnHzMwO\nrtCiEBHvRMQ0YDwwU9IJXZr8AGiOiJOANcDt3X2OpCWS2iS11fPFZDOzslVlSGpEvAqsA87qsr0j\nIt5KV28BPtrD+1dEREtEtDQ1NRUb1sxsCCusKEhqknRUujwCmAc806XN2IrV+cDTReUxM7PeFTn6\naCxwu6QGkuLz9xHxQ0nXAW0R0QpcKWk+0AnsARYXmMfMzHpR5OijrcD0brZfW7F8NVDOPNdmZvYu\ndTHNhZmZDQwXBTMzy7gomJlZxkXBzMwyLgpmZpZxUTAzs4yLgpmZZVwUzMws46JgZmYZFwUzM8u4\nKJiZWcZFwczMMi4KZmaWcVEwM7OMi4KZmWVcFMzMLFPk4ziHS/qppJ9L2ibpv3fT5nBJqyTtkLRJ\nUnNReczMrHdF9hTeAuZExMnANOAsSbO6tLkEeCUiPgTcCFxfYB4zM+tFYUUhEq+nq8PSV3Rpdg5w\ne7p8HzBXkorKZGZmB1foNQVJDZK2AC8DayJiU5cm44AXASKiE9gLjC4yk5mZ9azQohAR70TENGA8\nMFPSCf35HElLJLVJatu9e/fAhjQzs0xVRh9FxKvAOuCsLrteAiYASGoEjgQ6unn/iohoiYiWpqam\nouOamQ1ZRY4+apJ0VLo8ApgHPNOlWStwYbp8LrA2IrpedzAzsyppLPCzxwK3S2ogKT5/HxE/lHQd\n0BYRrcCtwJ2SdgB7gAUF5jEzs14UVhQiYiswvZvt11Ys7wfOKyqDmZn1je9oNjOzjIuCmZllXBTM\nzCzjomBmZhkXBTMzy7gomJlZxkXBzMwyLgpmZpZxUTAzs4yLgpmZZVwUzMws46JgZmYZFwUzM8v0\nWhQk+fGYZmZDRJ6ewkZJ90o6W5IKT2RmZqXJUxSmACuARcBzkv6HpCnFxjIzszL0WhQisSYizgc+\nR/L4zJ9K+rGkU3t6n6QJktZJekrSNklf6KbNmZL2StqSvq7t7rPMzKw6en3yWnpNYSFJT+FXwFKS\nZytPA+4FJvfw1k7gLyLicUkjgcckrYmIp7q0+9eI+ER//wJmZjZw8jyOcwNwJ/DJiGiv2N4maXlP\nb4qIXcCudHmfpKeBcUDXomBmZjUiT1E4NiKiux0RcX2eg0hqJnle86Zudp8q6efAL4EvR8S2PJ9p\nZmYDL8+F5h9JOurAiqT3S3og7wEkvQ/4B+CLEfFal92PA5Mi4mTgJuD+Hj5jiaQ2SW27d+/Oe2gz\nM+ujPEWhKSJePbASEa8AH8jz4ZKGkRSEuyLi+133R8RrEfF6urwaGCZpTDftVkRES0S0NDU15Tm0\nmZn1Q56i8I6kiQdWJE0Cuj2dVCm9p+FW4OmI+GYPbT544N4HSTPTPB15gpuZ2cDLc03hGmC9pB8D\nAmYDS3K873SSEUtPSNqSbvsqMBEgIpYD5wKXSeoE3gQW9HT9wszMitdrUYiIf5H0EWBWuumLEfHr\nHO9bT1JEDtZmGbAsT1AzMytenp4CwOHAnrT9VElExCPFxTIzszLkuXnteuDTwDbgt+nmAFwUzMzq\nTJ6ewidJ7lV4q+gwZmZWrjyjj54HhhUdxMzMypenp/AGsEXSQ0DWW4iIKwtLZWZmpchTFFrTl5mZ\n1bk8Q1JvlzQCmBgRz1Yhk5mZlSTP4zj/A7AF+Jd0fZok9xzMzOpQngvNXwdmAq8CRMQW4OgCM5mZ\nWUnyFIW3I2Jvl22/7balmZkNankuNG+T9BmgQdIxwJXAo8XGMjOzMuTpKSwFPkwyHPV7wGvAF4sM\nZWZm5cgz+ugNkplSryk+jpmZlSnP3Efr6Ob5CRExp5BEZmZWmjzXFL5csTwc+BTQWUwcMzMrU57T\nR4912fQTST8tKI+ZmZUoz81roypeYyT9O+DIHO+bIGmdpKckbZP0hW7aSNK3Je2QtDV9mI+ZmZUk\nz+mjx0iuKYjktNEvgEtyvK8T+IuIeFzSSOAxSWsi4qmKNn8GHJO+TgH+Jv3TzMxKkOf00eT+fHBE\n7AJ2pcv7JD0NjAMqi8I5wB3pc5k3SjpK0tj0vWZmVmV5Rh/9p4Ptj4jv5/iMZmA6sKnLrnHAixXr\n7ek2FwUzsxLkOX10CXAasDZd/xjJHc27SU4rHbQoSHof8A/AFyPitf6ElLQEWAIwceLE/nyEmZnl\nkKcoDAOmHjilI2kssDIiLurtjZKGkRSEu3roUbwETKhYH59u+z0RsQJYAdDS0vKueybMzGxg5Jnm\nYkKXc/y/Anr9dV2SgFuBpyPimz00awU+m45CmgXs9fUEM7Py5OkpPCTpAZJ5jwA+DTyY432nA4uA\nJyRtSbd9lbSgRMRyYDVwNrCD5LGfvfY+zMysOHlGH10h6T8Cf5JuWhER/5jjfetJhrEerE0An88T\n1MzMipenpwDwOLAvIh6U9AeSRkbEviKDmZlZ9eUZkvo5kpE/o4A/JhkyuhyYW2w0q3nr/qrsBGY2\nwPJcaP48yfWB1wAi4jngA0WGMjOzcuQpCm9FxG8OrEhqpJuptM3MbPDLUxR+LOmrwAhJ84B7gR8U\nG8vMzMqQpyhcRXL38hPAfyUZRvrfigxlZmblOOiFZkkNJBPWXQD8XXUimZlZWQ7aU4iId4BJkg6r\nUh4zMytRnvsUnid52lor8P8ObDzI1BVmZjZI9dhTkHRnujgf+GHadmTFy8zM6szBegoflfRHwE7g\npirlMTOzEh2sKCwHHgImA20V20Vyn8LRBeYyM7MS9Hj6KCK+HRHHA9+JiKMrXpMjwgXBzKwO9Xqf\nQkRcVo0gZmZWvjw3r5mZ2RDhomBmZpnCioKk2yS9LOnJHvafKWmvpC3p69qispiZWT55H7LTHyuB\nZcAdB2nzrxHxiQIzmJlZHxTWU4iIR4A9RX2+mZkNvCJ7CnmcKunnwC+BL0fEtpLzmNWUG9dsr+rx\nvjRvSlWPZ7WnzKLwODApIl6XdDZwP3BMdw0lLSF5JCgTJ06sXkIzsyGmtNFHEfFaRLyeLq8Ghkka\n00PbFRHREhEtTU1NVc1pZjaUlFYUJH1QktLlmWmWjrLymJlZgaePJH0POBMYI6kd+BowDCAilgPn\nApdJ6gTeBBZEhJ/9bGZWosKKQkSc38v+ZSRDVs0GhWpf9DUrg+9oNjOzjIuCmZllXBTMzCzjomBm\nZhkXBTMzy7gomJlZxkXBzMwyLgpmZpZxUTAzs4yLgpmZZVwUzMws46JgZmYZFwUzM8u4KJiZWcZF\nwczMMi4KZmaWKawoSLpN0suSnuxhvyR9W9IOSVslfaSoLGZmlk9hT14DVpI8We2OHvb/GXBM+joF\n+Jv0TzMrSRlPl/vSvClVP6b1rLCeQkQ8Auw5SJNzgDsisRE4StLYovKYmVnviuwp9GYc8GLFenu6\nbVc5cWww2fB8R6Gfv7HTz2O2oWlQXGiWtERSm6S23bt3lx3HzKxulVkUXgImVKyPT7e9S0SsiIiW\niGhpamqqSjgzs6GozKLQCnw2HYU0C9gbET51ZGZWosKuKUj6HnAmMEZSO/A1YBhARCwHVgNnAzuA\nN4CLispiZmb5FFYUIuL8XvYH8Pmijm9mZn03KC40m5lZdbgomJlZxkXBzMwyLgpmZpZxUTAzs4yL\ngpmZZcqc+8isZs3auaKqx9s4cUlVj2fWE/cUzMws46JgZmYZFwUzM8u4KJiZWcZFwczMMi4KZmaW\ncVEwM7OMi4KZmWV885oNiA3Pd5QdwcwGQKE9BUlnSXpW0g5JV3Wzf7Gk3ZK2pK//UmQeMzM7uCIf\nx9kA3AzMA9qBzZJaI+KpLk1XRcQVReUwM7P8iuwpzAR2RMTzEfEb4B7gnAKPZ2Zmh6jIojAOeLFi\nvT3d1tWnJG2VdJ+kCQXmMTOzXpQ9+ugHQHNEnASsAW7vrpGkJZLaJLXt3r27qgHNzIaSIovCS0Dl\nb/7j022ZiOiIiLfS1VuAj3b3QRGxIiJaIqKlqampkLBmZlZsUdgMHCNpsqTDgAVAa2UDSWMrVucD\nTxeYx8zMelHY6KOI6JR0BfAA0ADcFhHbJF0HtEVEK3ClpPlAJ7AHWFxUHjOrTTeu2V71Y35p3pSq\nH3OwKPTmtYhYDazusu3aiuWrgauLzGBmZvmVfaHZzMxqiIuCmZllPPdRHarWOdpZOz3fkVm9cU/B\nzMwyLgpmZpZxUTAzs4yLgpmZZVwUzMws46JgZmYZFwUzM8u4KJiZWcZFwczMMr6juWBlzABpZtZf\n7imYmVnGRcHMzDIuCmZmlim0KEg6S9KzknZIuqqb/YdLWpXu3ySpucg8ZmZ2cIVdaJbUANwMzAPa\ngc2SWiPiqYpmlwCvRMSHJC0Argc+XVSmejdr54qyI5gNCtUeADKYHv9Z5OijmcCOiHgeQNI9wDlA\nZVE4B/h6unwfsEySIiIKzGVWc6pZ0DdOXFK1Y9ngU2RRGAe8WLHeDpzSU5uI6JS0FxgN/LqIQB4e\nalb9HqWL0OAyKO5TkLQEOPCd9bqkZwf4EGMoqBAdIufqG+fqmyrl+uu+vqHuvl5/PsBBusiba1Ke\nDyuyKLwETKhYH59u665Nu6RG4EjgXc94jIgVQGG/3khqi4iWoj6/v5yrb5yrb5yrb4ZKriJHH20G\njpE0WdJhwAKgtUubVuDCdPlcYK2vJ5iZlaewnkJ6jeAK4AGgAbgtIrZJug5oi4hW4FbgTkk7gD0k\nhcPMzEpS6DWFiFgNrO6y7dqK5f3AeUVmyKlWx3I6V984V984V98MiVzy2RozMzvA01yYmVlmSBYF\nSaMkrZH0XPrn+3toN1HSjyQ9LempoqfhyJsrbXuEpHZJy4rMlDeXpGmSNkjaJmmrpMLuTK/V6VNy\n5Prz9Ptoq6SHJOUaIlh0rop2n5IUkqoywiZPLkn/Of2abZN0dy3kSn8urJP0s/Tf8uwq5bpN0suS\nnuxhvyR9O829VdJH+nWgiBhyL+B/Alely1cB1/fQ7mFgXrr8PuAPaiFXuv9bwN3Aslr4egFTgGPS\n5T8CdgFHFZClAfg34GjgMODnwNQubS4HlqfLC4BVVfga5cn1sQPfQ8BltZIrbTcSeATYCLTUQi7g\nGOBnwPvT9Q/USK4VwGXp8lTghaJzpcf6E+AjwJM97D8b+GdAwCxgU3+OMyR7CiTTa9yeLt8OfLJr\nA0lTgcaIWAMQEa9HxBtl50qzfRT4Q+BHBefJnSsitkfEc+nyL4GXgaYCsmTTp0TEb4AD06f0lPc+\nYK4kFZClT7kiYl3F99BGknt3ipbn6wXwlyRzj+2vQqa8uT4H3BwRrwBExMs1kiuAI9LlI4FfViEX\nEfEIySjNnpwD3BGJjcBRksb29ThDtSj8YUTsSpf/L8kP2K6mAK9K+n7aTbwhneSv1FyS3kNyi+iX\nC87Sp1yVJM0k+S3r3wrI0t30KeN6ahMRncCB6VOKlCdXpUtIfqsrWq+50tMMEyLin6qQJ3cukv+D\nUyT9RNJGSWfVSK6vAwsltZOMrlxahVx59PV7sFuDYpqL/pD0IPDBbnZdU7kSESGpuyFYjcBsYDqw\nE1gFLCa5t6LMXJcDqyOifSB/+R2AXAc+ZyxwJ3BhRPx2wALWEUkLgRbgjBrI8h7gmyTf27WmkeQU\n0pkkvapHJJ0YEa+WmgrOB1ZGxF9LOpXkXqsT6uX7vW6LQkT8aU/7JP1K0tiI2JX+EOuuW9oObInf\nzfJ6P8l5ukMqCgOQ61RgtqTLSa5zHCbp9Yjo8QJilXIh6Qjgn4Br0u5rEQZs+pQSciHpT0kK7RkR\n8VbBmfLkGgmcADyc/pLxQaBV0vyIaCsxFyT/BzdFxNvALyRtJykSm0vOdQlwFkBEbJA0nGT+oWqc\n3jqYXN+DvRmqp48qp9e4EPjf3bTZTHJO7sB58Tn8/rTfpeSKiAsiYmJENJOcQrrjUAvCQORSMpXJ\nP6Z57iswS61On9JrLknTgb8F5lfp/HivuSJib0SMiYjm9HtqY5qvyILQa67U/SS9BCSNITmd9HwN\n5NoJzE1zHQ8MB3YXnCuPVuCz6SikWcDeitO++VXjqnmtvUjOLz8EPAc8CIxKt7cAt1S0mwdsBZ4A\nVgKH1UKuivaLqc7oo15zAQuBt4EtFa9pBeU5G9hOcs3imnTbdSQ/zCD5T3ovsAP4KXB0lb6vesv1\nIPCriq9Pay3k6tL2Yaow+ijn10skp7aeSv8PLqiRXFOBn5CMTNoCfLxKub5HMqrvbZJe1CXApcCl\nFV+vm9PcT/T339F3NJuZWWaonj4yM7NuuCiYmVnGRcHMzDIuCmZmlnFRMDOzTN3evGZWLZLeIRkC\n2Aj8AlgU5d91a9Yv7imYHbo3I2JaRJxAMmHZ58sOZNZfLgpmA2sD/ZiEzKxWuCiYDZB0Ft25vHta\nBLNBw0XB7NCNkLSF300rvqbkPGb95qJgdujejIhpwCSS+Wd8TcEGLc99ZHaI0qnL35cuTyeZ3fOP\nI3nAj9mg4p6C2QCKiJ+RzKx7ftlZzPrDPQUzM8u4p2BmZhkXBTMzy7gomJlZxkXBzMwyLgpmZpZx\nUTAzs4yLgpmZZVwUzMws8/8BgorORnOCQqcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4a0bea8e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "matplotlib.pyplot.hist(cor_values, normed=1, alpha=0.5, label=\"cor\")\n",
    "matplotlib.pyplot.hist(pcor_values, normed=1, alpha=0.5, label=\"pcor\")\n",
    "matplotlib.pyplot.legend(loc=\"upper left\")\n",
    "matplotlib.pyplot.xlabel(\"R\")\n",
    "matplotlib.pyplot.ylabel(\"frequency\")\n",
    "matplotlib.pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fisher transform the partial correlation values, using `numpy.log`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_scores = 0.5*numpy.log((1+pcor_values)/\n",
    "                         (1-pcor_values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute a p-value for each gene pair (upper triangle only), using the fact that sqrt(M-N-5) times the fisher Z sore should be approximately univariate normal (with zero mean) under the null hypothesis that a given gene pair's measurements (conditioned on the measurements for all the other 162 genes) are independent. You will want to use `scipy.stats.norm.cdf`, `numpy.abs`, and `math.sqrt` function (see class session 19 exercise)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = gene_matrix_for_network_df.shape[1]\n",
    "\n",
    "P_values = 2*scipy.stats.norm.cdf(-numpy.abs(z_scores)*(math.sqrt((M-N-5))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many gene pairs have a P value less than 0.01? (use `which` and `length`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "706"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(numpy.where(P_values < 0.01)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the sizes of the components in the undirected graph whose edges have `P < 0.05` in the statistical test that you did?  You will need to use `zip`, `tolist`, `list`, and `igraph.Graph.TupleList` (see class session 19 exercise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH UN-- 164 706 -- \\n+ attr: name (v)'"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inds_tri = numpy.where(numpy.tri(*matrix_pcor.shape, k=-1))\n",
    "inds_sig = numpy.where(P_values < 0.01)\n",
    "graph_edge_list = list(zip(inds_tri[1][inds_sig].tolist(),\n",
    "                                                inds_tri[0][inds_sig].tolist()))\n",
    "final_network = igraph.Graph.TupleList(graph_edge_list, directed=False)\n",
    "final_network.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the graph degree distribution on log-log scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4k2WePvD7SZo2aUvTMzRtoJRCkfOhoCAnj6AOioh4\nHnUUFRx3Vl125Nr97bg7s6O7jDrODCgHFR3PIqI4uIiKFAGBQkEUrNAWeqInem7SU/L8/miDUNI2\naZO8eZP7c139oy/Jmy/Q5O77Ps/zfYSUEkRERF1plC6AiIj8EwOCiIicYkAQEZFTDAgiInKKAUFE\nRE4xIIiIyCkGBBEROcWAICIipxgQRETkFAOCiIicClG6gP6Ij4+XqampSpdBRKQqBw8erJJSJvT2\nOFUHRGpqKrKzs5Uug4hIVYQQp115HG8xERGRUwwIIiJyigFBREROMSCIiMgpVQ9SE5H3bc4pwcpt\nuSittcIUbcDyuRlYMDFZ6bLIBxgQRNStzTklWLHpKKxtNgBASa0VKzYdBQCGRBDgLSYi6tbKbbnn\nwsHB2mbDym25ClVEvsSAIKJuldZa3TpOgYUBQUTdiokIdXrcFG3wcSWkBI5BENFFpJRYtysf1U2t\nEAKQ8uc/M+i0WD43Q7niyGd4BUFEF2htt2PFpqP449YfccPYJPzvLeOQfN4Vw+PXDOcAdZBQ5RWE\nEGI+gPnp6elKl0IUUGotrVj65iHszT+Lx65Mx+NXj4BGI3BrphkV9c2Y9uxXqG5qU7pMAJx+6wuq\nvIKQUm6RUj5kNBqVLoUoYBRUNWHh6j04eLoGzy8ejyevzYBGI879eWKUHldkJODDQ8Vot9kVrPTn\n6bcltVZI/Dz9dnNOiaJ1BRpVBgQRedbevLNYsGo3aq1teGvJpVg4KcXp4xZnmlHZ0IIduZU+rvBC\nnH7rGwwIoiD3/oEi3PPKPiQMCMPmZZdjSmpst4+9YmQi4iPD8N6BIh9WeDFOv/UNBgRRkLLbJZ7Z\nehz/+uF3mDYsDh8unY7BceE9Pken1eCWycnYkVuBivpmH1V6MVO0vpvjnH7rSQwIoiBkaW3HI28e\nxJqsfNx92WC8dt8UGA06l567ONMMm13iw0PK3e+/aYLzwehfzUj1bSEBjgFBFGTO1Flx68t78cXx\ncjw9fxR+f9MYhGhd/ygYlhCJqamx+CC7CPL8BRI+9ENpPSLDtDAZ9RAAEgeEwaDT4K1vC1Hd1KpI\nTYFIldNcich1508HjR8Qhpa2dtjswCv3TsEVIxP7dM7FU8z4lw+O4MCpGkwd2v2YhTf8VN6AnT9V\n4slrRuCxq4afO559qhp3rd+H+zccwDtLLkV4KD/e+otXEEQBrOt00MqGFjQ02/Dolel9DgcAuH7s\nIESGhSgyWL1+Vz70Og3uvmzIBcczU2Px1zsm4mhxLR596xDaFJ6KGwgYEEQBzNl0UAngrW8L+3Xe\n8NAQzB9vwj+OlqK+2XcL5yoamrE5pxS3TjY77RN17ehB+MOCsdiRW4kVm44qdgssUDAgiALUkaJa\nlHhxOuhtU8xobrNjy5HSfp/LVW/sOY02ux0PzBja7WPuvHQw/vnq4dh4sBh/+pzrIvqDN+mIAkhL\nuw1bj57B63tO43BRLQQ6rhi68sR00PEpRmQMHID3DxThrkuH9P6EfrK0tuPNfadxzSUDkRof0eNj\nf3PVcJTXt2DVjjycqbViX0ENW3L0AQOCKACcqbPi7X2FeGd/IaoaW5EWH4Gn549CmE6D/9py/ILb\nTJ7qxiqEwOIpZvz+02P4saweIwdF9fucPfnwYDFqLW1YMivNpdr+sGAMjhbXYlPOz1c43BHPPQwI\nIpWSUmJ/QTVe33sK234oh11KXDUyEb+clooZ6fHn+igZdCFea2p388RkPPvZcbx3oAi/mz/aI+d0\nxmaXeOWbAkwwRyNzSIxLz9FqhNMpr46WHAyI3jEgiFTG0tqOzTmleGPvKfxY1gCjQYcHZgzF3ZcO\ncboSesHEZK99GMZGhOLaUYPwUU4JnrpuJMJCtF55ne3HynHqrAWr5o6EEKL3J3Q6U+d8tTdbcriG\nAUHkh5y1sp44OBp/33sa72cXob65HZckReHZhWNx04RkGEK988HsisVTzPjH0TPYfqwcvxhn8spr\nrN+Vj5QYA+aOHujW80zRBqcD9WzJ4RoGBJGfcaxdcIwblNRa8cT7h2GXQIhGYN6YQbh3eioyh8S4\n9du0t8xIj4fJqMd7B4q8EhA5hTXIPl2D//jFKLdWfAPA8rkZF/xbAtwRzx0MCCI/42ztgl0CA/Qh\n2P74bAwyOm9UpxStRmBRphl//eoEimssSInpueGfu9bvKkCUPgSLp5jdfq7j1trKbbkoqbUiRCPw\nx5vHcPzBRVwHQeRnurs/3tjc7nfh4HDr5I79IzYeLPboeYuqLfjs+zO489IhiAzr2++zCyYmY/dT\nV+KZhWPRbpdIS4j0aI2BjAFB5EfONrZA181tFH++b26ODceM9Hh8kF0Mu91zq5df+aYAGiFw3/TU\nfp/rhnFJ0Os0+OCgsntZqAkDgshP5FU24ubVe2Cz26HTXji2oIb75oszzSiptWJ3XpVHzldnacP7\n2UW4cbzJI1dOUXod5o0ehE8Ol6K5yy08co4BQeQH9uadxcLVe9DU0o4Plk7HykXjkRxtgACQHG3A\nMwvH+v1982tHD0R0uA7veqiB39v7C2FpteHBmb0vjHPVrZlm1De3Y/uxco+dM5BxkJpIYR8eLMZT\nm77DkLgIvHbfFJhjwzFpcIzfB0JXYSFaLJiQjLf3FaKmqdVpMz1XtbbbsWFPAWakx2OUyXMrtKel\nxcFk1GPjwWLMH++dKbmBhFcQRAqRUuK5z3Px5AdHMHVoLD5cOh3mWM/OAPK126aY0Wqz46Oc/u02\nt+VIKcrrW/DgzO6b8vWFRiNwy+QU7DpRibJuFtHRzxgQRApobrPhN+8exl+/OonbMs3YcP9Ul7f8\n9GeXJEVhXIoR7/djtzkpJdbtyseIgZGYPSLBwxUCiyanwC6BTTmenXEViBgQRD52trEFd6/fh0+O\nlOJf52Xg2VvGdjtzSY0WZ5rxY1kDviuuc/u5m3NKkPmHL/BjWQMqGlrw8WHPtxIfEheBqamx2Jhd\nzP0iehE4P5VEKuCYqXS0pA6r7pyEZXPS/WI1tCfdOMEEvU6D97LdG6x2rCA/29lgr9bShhWbjmJz\nP29XObMoMwX5VU04VFjr8XMHEgYEkY84ZipZWtvxzkOX4YZxSUqX5BVReh2uH5OELYdLYW11bTpp\na7sd/7nlh4tWkDs6r3ra9WOTYNBpPb6wL9AwIIh8YOPBYvzy1X1IGBCGj5ZdjkmDXWtZrVaLp5jR\n0NKOrUfP9Pi4srpmPP95LqY/+xVqLM63LvVG59XIsBBcPzYJnx5xPcSCEQOCyIscM5X+JYBmKrni\n0qGxSI0Ld3qbSUqJPXlVWPrmQVz+P1/hrztOYnyKEXHdTIv11gryRZNT0NDSjs+PlXnl/IGA6yCI\nvKS5zYblG7/DliOluC3TjD/cPCagBqN7IoTAmOQofPpdGYY+9Q+Yog147Mp0tNnseGPvaZyoaER0\nuA4PzhiKuy8bAnNs+EVdbAHvriC/dGgsUmIM+CC7GDdNUNeaE19hQBB5wdnGFjz094M4eLoGv503\nEo/MTgu4weiebM4pwRfHKgB07IldUmvFU51bfY5NNmLlonGYP94Eve7nfSzO77zqi/2jNRqBRZNT\n8OKXJ1BSa0WyH/e6UgoDgsjD8iobcf9rB1Be34xVd04K2MHonqzclovmdvtFxxMiw/DJry/vNiy9\nufudM7dMSsGfvziBTQeL8dhVw332umoRHNe7RD4SLDOVetPdwHJVY4tfXUmZY8MxLS0OGw9xTYQz\nDAgiDwm2mUo96W5g2R9bli+anILTZy3IPl2jdCl+h7eYiPrgwj2j9RhtisLnxypweXocVt81OSDa\nZvSHmrb6vG7sIKzY9B3ufXU/rK02r499qAkDgshNF+8Z3YyS2mZcNjQWG+6fGjQzlXri6wHn/vj8\nh3LYJNDa+vMe4Cs6B9T9sV5fYkAQucnZntEAUFRjYTicx9cDzn21clsubF12wXOs4FZD/d7En2Yi\nN3U3AFtay/bRatT9/6fnV3CrDQOCyE1qGoCl3vH/s3sMCCI3LZ+bAa1GfXtGk3PL52bAcN6CPYD/\nnw5+ExBCiDQhxCtCiI1K10LUkysyEqERHR8iatozmpxbMDEZzywciyh9x5BsklHP/89OXh2kFkK8\nCuAXACqklGPOOz4PwIsAtADWSymflVLmA3iAAUH+7u/fnkKbTWLLY9MxcpDn9ksm5SyYmAy9ToNH\n3jyEdb/MxJhko1vPv3Das//O2HKXt68gNgCYd/4BIYQWwCoA1wEYBeAOIcQoL9dB5BHNbTZs2HMK\nczISGA4BJiWmo8tuUbXFrec5pj2X1FrP9Z3y1kZHvubVgJBSZgGo7nJ4KoCTUsp8KWUrgHcB3OTN\nOog85cNDxahqbMXDs4YpXQp5mKMNe1GNewHhbNqztzY68jUlxiCSAZzfJL4YQLIQIk4I8TKAiUKI\nFd09WQjxkBAiWwiRXVlZ6e1aic6x2SXWZeVjfIoRl6XFKl0OeZjRoMMAfQiKa9yb3hrI02T9ZpBa\nSnlWSvmIlHKYlPKZHh63VkqZKaXMTEhI8GWJFOQ+/6EMp85a8PDsYX7VcI48xxwT7vYtpkCeJqtE\nQJQAMJ/3fUrnMSK/JaXEyzvzMCQuHHNHD1K6HPISc6wBRW5eQXRMk73wozRQpskqERAHAAwXQgwV\nQoQCuB3AJwrUQeSyfQXVOFJchyUz0y5aA0GBwxwTjuIai1utvxdMTL4gDCLDQgJmmqxXA0II8Q6A\nvQAyhBDFQogHpJTtAH4NYBuA4wDel1L+4M06iPprzc48xEWEYtHkFKVLIS8yx4ajuc2OysYWt543\nMunnGW2XpcUGRDgAXl4HIaW8o5vjWwFs9eZrE3nKj2X12JFbiSevGXHBFpkUeMyxHeMGRdVWJA7Q\nu/w8x7jF2GQjTlY0eqU2JfjNILU7hBDzhRBr6+rqlC6FgsDarHwYdFrcM22I0qWQl5k710IUuznV\ntajaCq1GYNaIeBRWW9DspNuvGqkyIKSUW6SUDxmN7q12JHJXaa0Vnxwuxe1TzYgOD1W6HPKy5JiO\nKwh3p7oWVltgitYjY1AU7BIoqGryRnk+p8qAIPKVV78pgATwwIyhSpdCPhAeGoL4yFC3p7oW1Vgw\nODYc6QmRABAwt5kYEETdqLO04Z39hbhxvOlcGwYKfCkx4W6vpi6qtsAcE460hAhoBAOCKOC9ue80\nmlpteGhWmtKlkA+ZY8NRVO36LSZLazuqGlthjg2HXqeFOTacAUEUyJrbbHhtdwFmj0jAJUlsyhdM\nzDEGlNZaL9qGtDuOMHH0ckpPiGRAEAWyTYdKOpryzebVQ7Axx4aj3S5xps61q4jCzvGKwY6AGBiJ\ngqomtNvsXqvRVxgQRF3Y7BLrduVjXIoR09LilC6HfOznqa6uBURR14BIiESrzX4uONRMlQHBdRDk\nTduPlaGgqgkPz2JTvmCUEuNYLOfaB3xhtQURoVrEhOsAAOmJgTOTSZUBwXUQ5C1SSry0Mx+DY8Mx\nbwyb8gUjU7QBQsDlpn3FNRaYY8PP/TIxzBEQlQwIooCyv6AaR4pqsWQWm/IFq9AQDZKi9Ch24wrC\nMUANAFF6HQZF6XGynAFBFFDWZOUjLiIUt7IpX1BLiXVtLYSUEkXV1nPjDw7piZG8giAKJLllDfjq\nxwrcOz2VTfmCXMfGQb3fYqpqbIW1zQZzzIWbA6Undkx1dadtuD9iQBB1OteU7zI25Qt25lgDyhua\n0dLec9O9c1Nc4y6+grC02lBa1+y1Gn2BAUEE4EydFR8fLsFtU8yIiWBTvmBnjgmHlEBpbc8f8I6u\nr85uMQHqn8nEgCACm/LRhVyd6lp41tL5eAYEUUCqs7Th7X2F+MW4pAtmo1Dwcvwc9DZQXVRjQeKA\nsIvGrOIiQhETrsPJigav1egLqgwILpQjT3I05Xt41jClSyE/MTBKD51W9DpQ3XWKq4MQ4txAtZqp\nMiC4UI48paMp3ynMGpGAUSY25aMOWo1AcrSh9ysIJ1NcHdITI3FC5TOZVBkQRJ7yUU4Jqhpb8Ahb\nelMX5tjwHhfLtbbbcabOetEUV4f0xAGotbThbFOrt0r0OgYEBS2bXWJdVj7GJhsxbRib8tGFUmLC\ne2zYV1prhV2i23GrQBioZkBQ0Np+rBz5VU14eHYam/LRRcyxBpxtakVTS7vTPy/qZoqrAwOCSKWk\nlHh5Z15HU77RbMpHF0vppe23Y5Fcd1cQJqMeEaFaBgSR2hw4VYPDRbVYMnMoQrR8G9DFzL2shSiq\ntiJUq8HAKL3TPxdCYJjKZzLxnUFBac3OPMRGhGLRZLPSpZCf6m0tRFG1Bckxhh67/qp9+1EGBAWd\nn8ob8OWPFbh3WioMoWzKR87FRYTCoNN2uxaiuzUQ50sfGImy+mbUN7d5o0SvY0BQ0HE05fvlNDbl\no+4JIWCONZzrt9RVUY0Fg2OdT3F1SE/oGKjOU+lVhCoDgiupqa/YlI/cYY4Jd7qzXH1zG2otbef2\nr+6O2mcyqTIguJKa+uq13adgl2zKR65JiTGguNpy0Wpox8B1d1NcHQbHhiNUq2FAEPm7OmtHU74b\nxrIpH7nGHBuOhpZ21FkvHEMo6mWKq0OIVoOh8REMCCJ/99a+02hsacfDs9lWg1zjWAvRdaDa8b0r\nv2ioeftRBgQFBUdTvpnD4zHaxFuT5Bpz5yB016muhdUWROlDYDToej1HemIkCqstaG7reXc6f8SA\noKCwOacElQ0teGQ2W3qT686thai+OCC6bjPanfTESEgJ5Fc2ebw+b2NAUMCz2yXWZuVjTHIUprMp\nH7khSq+D0aC7qN1GxxRX1wMCgCpvMzEgKOBtP97ZlG/WMDblI7eZYy/cF8JulyiutvY6xdVhaHwE\nNAI4Wa6+3eUYEBTQHE35zLEGXDeGTfnIfSnR4RfcYqpoaEGrze7yTDi9TovBseGBfwUhhIgRQowW\nQqQJIRgu5PeyT9cgp7AWS2amsSkf9UnHamrrubUQvXVx7WpzTgnK6pqx9WgZLn/2K2zOKfFarZ4W\n0tsDhBBGAI8CuANAKIBKAHoAA4UQ3wJYLaXc4dUqifpozc48xITrcCub8lEfmWPD0dJuR2VDCxKj\n9C4vkgM6wmHFpqNobrcDAEpqrVix6SgAYMHEZO8V7SGu/Eq1EUARgJlSygwp5QwpZaaU0gzgWQA3\nCSEe8GqVRH1worwBXxyvwL3T2ZSP+s4x1uAYhyistkAIwBTtvM33+VZuy4W1y/RWa5sNK7fler5Q\nL+j1CkJKeU0Pf3YQwEGPVuQCIcR8APPT09N9/dKkImuz8qHXafDLaalKl0Iq5lgLUVxjxeQhHVNe\nk6L0CAvp/ZeO0lrnnWC7O+5vXL4pK4T4ry7fa4UQb3m+pN6xFxP1pqyuGZsPl+C2TDNi2ZSP+uHn\n1dQdVxBFNb23+XYwRTvv9trdcX/jzqidWQixAgCEEGEANgE44ZWqiPrp1d0FsNklHpzJthrUP3qd\nFgkDws6113BlHwiH5XMzYNBdeKVh0GmxfG6Gx+v0BncC4lcAxnaGxBYAO6SUT3ulKqJ+ONeUb5yJ\nTfnII1JiOtZCNLfZUF7f4vIiuQUTk/HMwrFI7rxiCA/V4pmFY1UxQA24EBBCiElCiEkAJgJ4EcBt\n6LhyyOo8TuRX3t5X2NGUbxavHsgzOvaFsJxbUW3uZaOg8y2YmIzdT12J8SlGTB4So5pwAFwYpAbw\nXJfvawCM6jwuAVzp6aKI+qql3YZXdxdg5vB4jEnmGBV5hjnWgH8cPYNTVR39lFy9gjhfktGAExXq\nWk3tyiymK3xRCJEnOJryvbB4gtKlUAAxx4TDZpfYf6q64/u+BES0HlknKiGlVE3LF1duMd0tevjb\nCCGGCSFmeLYsIvfZ7RJrsvIx2hSFy9PZlI88xxEIe/KqoNdpkBAZ5vY5kqMNsLTaUG9t93R5XuPK\nLaY4AIeFEI41D46V1OkAZgOoAvCU1yokctEXx8uRX9mEv9wxUTW/oZE6OBbL/VBaj/SEyD79fCUZ\nO8YtSmqtMIb3vo+EP+j1CkJK+SKASQDeAZAA4KrO70sA3COlvEVKyemupLg1WflIiTHgejblIw9L\nitZDIwAp+3Z7yXEOADhTp45FcoBrVxCQUtoAbO/8IvI72aeqcfB0Df7zxtFsykcep9NqkGQ0oKTW\n2qcBagAwdV5BlNY1e7I0r3KlWd9/9PDHUkr5ew/WQ+SWzTklWLktFyW1VmgEYNAxHMjzHJMfAOCj\nnGJMMEe7PV01YUAYQjQCZ1TSZgNw7QrC2T55EQAeQMf4BAOCFOHolOlohmaXwO8+OYbQEK2q5pqT\nf3P8nLXaOjqy1lnb+9SRVasRGBilV00fJsC1MYjnHF8A1gIwALgfwLsAuBKJFKP2TpmkDp78OTNF\n61V1i8ml63EhRKwQ4g8AvkPHVcckKeVvpZQVXq2u+3rmCyHW1tXVKfHy5CfU3imT1MGTP2dJRoOq\nBqldWQexEsABAA0Axkopn5ZS1ni9sh6wmysB6u+USergyZ8zU7QBZXXNsNtlf8vyCVeuIJ4EYALw\n7wBKhRD1nV8NQoh675ZH1L2lc4ZddExNnTJJHTzZkdUUrUebTaKqscVT5XmVK602OC2E/FJFfce9\n3MQBYahsaIEp2oDlczM4QE0e5fh5WrktF6W11n79nCWdN9U1Mar3HemU5tI6CCJ/09Dchg17TmHu\n6IFYc0+m0uVQgFswMdkjv3gkGTsXy9VaMcEc3e/zeRuvDkiV3t5XiPrmdiybw21nST0c4xZqmcnE\ngCDVaW6zYf03BZiRHo/xKvgtjMghJlwHvU6jmsVyDAhSnY0Hi1HZ0IJlTgapifyZEAImowGlKpnq\nyoAgVWm32bEmKw8TzNGYNowtvUl9kqL1KK3lLSYij/v0uzMoqrZi2ZxhbOlNqqSmxXIMCFINu13i\npa/zMGJgJK6+ZKDS5RD1iSnagIqGFrR19nbyZwwIUo0vf6xAbnkDls4ZBo2GVw+kTiajHlICZSqY\nycSAIFWQUmLVjpNIiTFg/jiT0uUQ9VlS51TXMwwIIs/Ym38Wh4tq8fDsYdwQiFTNZFTPznJ8p5Eq\nvPR1HuIjw3Dr5BSlSyHqF8cVhBpmMjEgyO99V1yLXSeq8ODModB3aZpGpDaRYSGI0oeooi09A4L8\n3uodeYjSh+CuSwcrXQqRR5ii1THVlQFBfu1kRQO2HSvDvdNTMUCvU7ocIo9IMqpjsRwDgvzaS1/n\nIyxEg/umpypdCpHH8ArCi7jlaHAorrHg48MluGPqYMRFhildDpHHmKINqLG0wdpq6/3BClJlQHDL\n0eCwLisfQgBLZqYpXQqRRzn2hfD3pn2qDAgKfFWNLXj3QBEWTEjmHtMUcBw7y53x83EIBgT5pVe/\nKUCrzY5H2NKbAlDyuY2DeAVB5Jb65jb8fe9pXDdmEIYlRCpdDpHHDTR2jKn5+1oIBgT5nb/vPY2G\nFm4nSoErLESL+Mgw3mIickdzmw2v7S7ArBEJGJPMSQgUuEzRet5iInLH+9lFqGpsxaMce6AAZzIa\n/L6jKwOC/EabzY41O/MxeUgMpg6NVbocIq/q2HrUCiml0qV0iwFBfuOTw6UoqbXi0Su4nSgFPpPR\nAEurDfXWdqVL6RYDgvyC3S7x0s48jBw0AFdkJCpdDpHXJUX7/2I5BgT5hc+PleNkRSOWzuHVAwUH\n07md5RgQRN2SUmL11ycxJC4cN4xNUrocIp8wda6mLvHjqa4MCFLc7pNn8V1xHR6exe1EKXgkDAhD\niEbgjB8vluO7kRS3+uuTSBwQhlsmJytdCpHPaDUCA6P0fj3VlQFBisoprMGevLNYMjMNYSHcTpSC\ni6lzqqu/YkCQolZ/nQejQYc7uZ0oBaEkP18sx4AgxeSWNWD7sXLcNz0VEWEhSpdD5HNJ0XqcqbPC\nbvfPxXIMCFLMyzvzEB6q5XaiFLRMRgPabBJVTS1Kl+IUA4IUUVRtwSdHSnHn1MGIiQhVuhwiRZxb\nC+GnU10ZEKSINVl50AjgQW4nSkHMsfWovy6WY0CQz1U0NOP97GIsmpyCQZ1vEKJg5LiC8NfFcgwI\n8rlXvilAu82Oh2expTcFt5hwHcJCNH67WC7opo5szinBym25KK21whRtwPK5GVgwkQu0fKXO0oa3\nvi3EDeNMSI2PULocIkV9fLgU7XaJ9d8U4LPvy/zu8yioAmJzTglWbDoKa5sNAFBSa8WKTUcBwK/+\nUwLZG3tPobGlHUtn8+qBgpvj88jWOcXVHz+PguoW08ptuefCwcHaZsPKbbkKVRRcLK3teG3PKVw5\nMhGjTFFKl0OkKDV8HqkyIIQQ84UQa+vq6tx6XndL2v15qXsgeXd/EaqbWrGM24kSqeLzSJUBIaXc\nIqV8yGh0b1N7x4yBrnRaDX4qb/BEadSN1nY71u3Kx9TUWGSmcjtRou4+j7o7rgRVBkRfLZ+bAYPu\nwoZwOq2AVgNc/+IuPPvZj7C0+u/2f2q2OacEZ+qasewKXj0QAc4/jww6LZbPzVCooosFVUAsmJiM\nZxaORXK0AQJAcrQBKxeNx+6nrsLCScl4eWcernk+C18eL1e61IBis0u8vDMPo01RmD0iQelyiPyC\ns8+jZxaO9ZsBagAQUvpnkyhXZGZmyuzsbI+db39BNf5981H8VN6IuaMH4nfzR/vV5Z5abT16Bsve\nOoRVd07CDeO4YxyR0oQQB6WUmb09LqiuIHozdWgsPn1sJn47byR2/lSJq5/fiXVZ+Wiz2ZUuTbWk\nlFi14yTS4iMwb8wgpcshIjcwILoIDdFg6Zxh2P74bExLi8N/bz2O+X/9BgdP1yhdmiplnajCD6X1\neGT2MGg1QulyiMgNDIhumGPDsf7eTKy5ZzLqrG245aU9WLHpKGotrUqXpiqrdpxEklHvV/dVicg1\nDIgeCCGbXxquAAALtUlEQVQwd/QgfPHEbCyZORTvZxfhqud2YtOhYqh57MZXsk9VY39BNZbMTENo\nCH/UiNSG71oXRISF4N9uGIVPH5uBIXHheOL9I7hj3bc4WdGodGl+bfXXeYgJ1+H2qWalSyGiPmBA\nuOGSpChsfGQ6nlk4FsfPNOC6F7Pwp225aO6yXJ6A42fq8dWPFbj/8qEIDw2qll9EAYMB4SaNRuCO\nqYPx5ZOzMX+8CX/bcRLXvpCFr3MrlC7Nr7z0dR4iQrW4d1qq0qUQUR8xIPooPjIMzy+egLeXXIoQ\nrcB9rx3Ao28dQlmdf2784Uunqprw6XeluPuyITCG65Quh4j6iAHRT9OHxeOz38zEv1w7Al8cL8fV\nz+/Ea7s7NsQJVmuy8hGi1eCBGUOVLoWI+oEB4QFhIVr8+srh2P74bEweEoP/3HIMN63ajcNFtUqX\n5nPl9c348GAxbp2cgsQobidKpGYMCA8aHBeODfdPweq7JqGqsQU3r96N/7f5e9RZ25QuzWfW78qH\nTUpuJ0oUABgQHiaEwPVjk/DFE7Nx3/RUvLXvNK56bic+PlwS8Gsnappa8da+Qswfl4TBceFKl0NE\n/cSA8JIBeh1+N380Pvn1DCRH6/Gbdw/jnlf2o6CqSenSvOb1vadgabVh6Zx0pUshIg9gQHjZmGQj\nNi27HL+/aTSOFNdi7gtZeGH7TwG3dqKppR0b9pzC1ZcMRMagAUqXQ0QewIDwAa1G4J5pqfjyydm4\nbuwgvPjlCcz7cxZ2nahUujSPeWd/IWotbdwQiCiAMCB8KHGAHi/ePhFvPnAphBC455X9+Kd3clDR\noO61Ey3tNqzblY9paXGYNDhG6XKIyEMYEAqYMbxj7cQ/Xz0c//d9Ga760068sfcUbHZ1DmJvOlSC\n8voWXj0QBRgGhEL0Oi3++eoR2Pb4LIw3R+M/Pv4BN6/eje9L6pQuzS3tNjte3pmHcSlGzEiPV7oc\nIvIgBoTChsZH4O8PTMVf7piI0tpm3Pi3b/D0Jz+goVkdaye2fl+G02ctWDZnGITghkBEgYQB4QeE\nELhxvAlfPjkbd182BK/vPYWrntuJT78r9eu1E1JKrN5xEsMSInDtKG4nShRoGBB+xGjQ4b9uGoPN\nyy5HYlQYfv12Du597QBOn/XPtRM7civwY1kDls5Jh4bbiRIFHAaEHxpvjsbHj87A0/NH4dDpGlz7\nQhb++uUJtLT719qJ1TvykBxtwE0TTEqXQkRewIDwU1qNwH2XD8WXT87G1aMG4rntP+G6F3dhT16V\n0qUBAPYXVCP7dA0empUGnZY/RkSBiO9sPzcwSo9Vd07ChvunoN0mcee6fXjivcOoamxRtK5VO04i\nLiIUizO5nShRoGJAqMScjER8/vgsPHZlOrZ8V4or//Q13tp3GnYF1k58X1KHnT9V4lczhsIQqvX5\n6xORbzAgVESv0+LJazPw2W9mYbTJiH/76Hvc8vIeHCut92kdL32dhwFhIbhn2hCfvi4R+RYDQoXS\nEyPx9pJL8cJt41F41oL5f/sGf/j0GBpb2r3+2vmVjdj6/RncM20IovTcTpQokDEgVEoIgZsnpuCr\nJ+fgtilmrP+mANc8vxP/932ZV9dOvLwzD6FaDX7F7USJAh4DQuWM4Tr88eax+HDpdBgNOjzy5kE8\n+Ho2iqotHn+tM3VWfJRTgtunmBEfGebx8xORf2FABIjJQ2Lw6WMz8O83XIK9+WdxzQs7sfrrk2ht\nt3vsNdZlFUBKYMmsNI+dk4j8FwMigIRoNXhwZhq+eGI25oxIxP/+Xy5u+Msu7Ms/2+9zVze14p39\nhbhxggkpMdxOlCgYMCACkCnagJfvmYxX7s2EpdWG29Z+i+UfHEF1U2ufz7lhdwGa221YNoctvYmC\nBQMigF11yUB88cRsLJ0zDB/llODK577GewcK3V470dDchg17TuHaUQORnsjtRImCBQMiwBlCtfjt\nvJHY+puZGJE4AL/98CgWr9mL3LIGl8/x9r5C1De3Y9mcdC9WSkT+hgERJEYMHID3Hr4MKxeNQ15l\nI274yy4889lxWFp7XjvR3GbD+m8KMCM9HuPN0T6qloj8AQMiiAghcGumGV89OQe3TErBmp35uOb5\nLGw/Vt7tczYeLEZlA7cTJQpGfhMQQogIIcTrQoh1Qoi7lK4nkMVEhOJ/Fo3DxkemITIsBEveyMaS\nN7JRUmu94HHtNjvWZOVhgjka09LiFKqWiJQS4s2TCyFeBfALABVSyjHnHZ8H4EUAWgDrpZTPAlgI\nYKOUcosQ4j0Ab3mzNgIyU2Px6T/NwKvfFODPX5zA1c/txOPXDEdcRCie337iXGBce8lAbidK5Cc2\n55Rg5bZclNZaYYo2YPncDCyYmOyV1/JqQADYAOBvAN5wHBBCaAGsAnANgGIAB4QQnwBIAXC082H+\ntTNOANNpNXh49jDcMC4JT3/yA/649UcIAOfPc3p7fyHGpkR77YeQiFyzOacEKzYdhbWt4yOypNaK\nFZs6Pja98f706i0mKWUWgOouh6cCOCmlzJdStgJ4F8BN6AiLFF/URRdLiQnH+nunIDYiFF0nwVrb\n7Fi5LVeRuojoZyu35Z4LBwdrm81r708lPoiTARSd931x57FNAG4RQrwEYEt3TxZCPCSEyBZCZFdW\nVnq30iBU081iutIu4xNE5HvdvQ+99f709i0ml0kpmwDc78Lj1gJYCwCZmZm+3y0nwJmiDRcNVjuO\nE5GyfP3+VOIKogTA+ftUpnQeIz+wfG4GDLoLd4kz6LRYPjdDoYqIyMHX708lriAOABguhBiKjmC4\nHcCdCtRBTjgGunw1S4KIXOfr96fw5uYyQoh3AMwBEA+gHMDvpJSvCCGuB/BndExzfVVK+d99OX9m\nZqbMzs72VLlEREFBCHFQSpnZ2+O8egUhpbyjm+NbAWz15msTEVH/cDopERE5pcqAEELMF0Ksraur\nU7oUIqKApcqAkFJukVI+ZDQalS6FiChgqTIgiIjI+7w6i8nbhBCVAE6fd8gIwNX7TvEAqjxeVGBx\n599TaUrV6u3X9eT5+3uuvj6/L8/je9mzuv57DpFSJvT6LCllwHwBWOvGY7OVrtffv9z591T6S6la\nvf26njx/f8/V1+f35Xl8L3v2q6//d4F2i6nbHk7UJ2r691SqVm+/rifP399z9fX5fXmemn721KBP\n/56qvsXUH0KIbOnCQhEi8m98L3tPoF1BuGOt0gUQkUfwvewlQXsFQUREPQvmKwgiIuoBA4KIiJxi\nQBARkVMMiE5CiAghxOtCiHVCiLuUroeI3CeESBNCvCKE2Kh0LYEgoANCCPGqEKJCCPF9l+PzhBC5\nQoiTQoinOg8vBLBRSrkEwI0+L5aInHLnfSylzJdSPqBMpYEnoAMCwAYA884/IITQAlgF4DoAowDc\nIYQYhY6tT4s6H2bzYY1E1LMNcP19TB4U0AEhpcwCUN3l8FQAJzt/02gF8C6AmwAUoyMkgAD/dyFS\nEzffx+RBwfhBmIyfrxSAjmBIBrAJwC1CiJfAZf5E/s7p+1gIESeEeBnARCHECmVKCxxe3XJUTaSU\nTQDuV7oOIuo7KeVZAI8oXUegCMYriBIA5vO+T+k8RkTqwfexDwRjQBwAMFwIMVQIEQrgdgCfKFwT\nEbmH72MfCOiAEEK8A2AvgAwhRLEQ4gEpZTuAXwPYBuA4gPellD8oWScRdY/vY+WwWR8RETkV0FcQ\nRETUdwwIIiJyigFBREROMSCIiMgpBgQRETnFgCAiIqcYEEQeJIRI7dqWmkitGBBEROQUA4LISzp3\nN8sRQkxRuhaivmA3VyIvEEJkoGOPgvuklEeUroeoLxgQRJ6XAOBjAAullMeULoaor3iLicjz6gAU\nApihdCFE/cErCCLPawVwM4BtQohGKeXbShdE1BcMCCIvkFI2CSF+AWB7Z0hwrwJSHbb7JiIipzgG\nQURETjEgiIjIKQYEERE5xYAgIiKnGBBEROQUA4KIiJxiQBARkVMMCCIicur/A109LqPfv2FjAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4a0bea82e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "degree_dist = final_network.degree_distribution()\n",
    "xs, ys = zip(*[(left, count) for left, _, count in degree_dist.bins()])\n",
    "matplotlib.pyplot.loglog(xs, ys, marker=\"o\")\n",
    "pylab.xlabel(\"k\")\n",
    "pylab.ylabel(\"N(k)\")\n",
    "pylab.show()"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}