tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configuration Models of Random Hypergraphs\n",
    "*Phil Chodrow, MIT Operations Research Center and Laboratory for Information and Decision Systems*\n",
    "\n",
    "This notebook is a basic tutorial for conducting null hypothesis tests by randomizing polyadic data under hypergraph configuration models. The `hypergraph` class is implemented in `hypergraph.py`, and includes methods for some basic descriptives as well as the relevant Monte Carlo algorithms. \n",
    "\n",
    "We also illustrate the triadic closure analysis shown in the main text. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from hypergraph import *\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_email_data():\n",
    "    '''\n",
    "    Read the email data from the supplied .txt, returning it as a list of tuples. \n",
    "    '''\n",
    "    with open('email-enron.txt', 'r') as f:\n",
    "        F = f.read()\n",
    "    L = F.split('\\n')\n",
    "    L = [tuple([int(v) for v in g.split(' ') if v is not '']) for g in L]\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = read_email_data() # read the data\n",
    "G = hypergraph(L)     # construct a hypergraph from the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean node degree is 188.0, and the standard deviation is 238.0.\n",
      "The mean edge dimension is 2.5, and the standard deviation is 1.4.\n"
     ]
    }
   ],
   "source": [
    "def print_moments(G):\n",
    "    D = G.node_degrees()\n",
    "    K = G.edge_dimensions()\n",
    "    print('The mean node degree is ' + str(round(D.mean())) + ', and the standard deviation is '+ str(round(np.sqrt(np.var(D)))) + '.')\n",
    "    print('The mean edge dimension is ' + str(round(K.mean(), 1)) + ', and the standard deviation is '+ str(round(np.sqrt(np.var(K)), 1)) + '.')\n",
    "    \n",
    "print_moments(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Sampling and Triadic Closure\n",
    "\n",
    "We'll now illustrate the use of Monte Carlo sampling to conduct a simple hypothesis-test using the `hypergraph` class. We'll replicate the triadic closure study shown in the main text. Recall that our metric of interest is the proportion of 3-cycles in the projected graph that are \"filled in\" by a higher-dimensional edge. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_triangles(C):\n",
    "    '''\n",
    "    Count the number of triangles in the (unweighted) projected graph of hypergraph C\n",
    "    '''\n",
    "    G = projected_graph(C, as_hyper=False)\n",
    "    G_ = nx.Graph(G)\n",
    "    triangles = nx.triangles(G_)\n",
    "    n_triangles = sum(triangles.values())/3\n",
    "    return(n_triangles)\n",
    "\n",
    "def count_filled_triangles(C):\n",
    "    '''\n",
    "    Count the number of distinct sub-edges of size 3 in C. \n",
    "    '''\n",
    "    container = set()\n",
    "    for f in C.C:\n",
    "        for t in itertools.combinations(f, 3):\n",
    "            container.add(tuple(sorted(t)))\n",
    "    return(len(container))\n",
    "\n",
    "def analyze_closure(G):\n",
    "    tri = count_triangles(G)\n",
    "    filled = count_filled_triangles(G)\n",
    "    print('The rate of simplicial closure in triangles is ' + str(round(filled / tri, 2)) + '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rate of simplicial closure in triangles is 0.81.\n"
     ]
    }
   ],
   "source": [
    "analyze_closure(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compare this value to a sample from a null model. Let's start with a stub-labeled null. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 steps completed.\n",
      "Task completed in 3.24 seconds\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "tick = time.perf_counter()\n",
    "G.MH(n_steps = 100000, label = 'stub')\n",
    "tock = time.perf_counter()\n",
    "\n",
    "print('Task completed in ' + str(round(tock - tick, 2)) + ' seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that the randomized graph has the same degree sequence and edge-dimension sequence as the old: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean node degree is 188.0, and the standard deviation is 238.0.\n",
      "The mean edge dimension is 2.5, and the standard deviation is 1.4.\n"
     ]
    }
   ],
   "source": [
    "print_moments(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the null value for simplicial closure in triangles compare to the observed value? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rate of simplicial closure in triangles is 0.25.\n"
     ]
    }
   ],
   "source": [
    "analyze_closure(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simplicial closure in triangles occurs in the observed data three times as often as would be expected by the stub-labeled configuration model. Formally, we should repeatedly randomize and pull many samples, but the spread is not large. Take my word for it, or try it yourself. \n",
    "\n",
    "What about the vertex-labeled model? \n",
    "\n",
    "The vertex-labeled MH algorithm accepts an additional parameter, `n_clash`. Tracking the counts of each edge is very expensive, so instead we implement an approximation that resamples the edge counts only when a \"clash\" requires an update. We recommend setting `n_clash = 1` when using vertex-labeled randomization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2373 epochs completed, 100010 steps taken, 477951 steps rejected.\n",
      "Task completed in 25.98 seconds\n"
     ]
    }
   ],
   "source": [
    "tick = time.perf_counter()\n",
    "G.MH(n_steps = 100000, label = 'vertex', n_clash = 1)\n",
    "tock = time.perf_counter()\n",
    "print('Task completed in ' + str(round(tock - tick, 2)) + ' seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The degree and dimension sequences are again preserved: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean node degree is 188.0, and the standard deviation is 238.0.\n",
      "The mean edge dimension is 2.5, and the standard deviation is 1.4.\n"
     ]
    }
   ],
   "source": [
    "print_moments(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about simplicial closure? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rate of simplicial closure in triangles is 0.21.\n"
     ]
    }
   ],
   "source": [
    "analyze_closure(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preserving Node-Edge Dimension Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases, it may be desirable to consider a more restrictive class of nulls that preserve more detailed information about the data. A fairly natural definition is to require that each node preserve not only its total degree, but also the number of edges of each dimension incident to it. Class `hypergraph` implements this as \"detailed\" MCMC, as we demonstrate here. \n",
    "\n",
    "First, we'll read in the data and take a look at the \"node-dimension matrix,\" whose $ij$th entry is the number of edges of dimension $j$ incident on node $i$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = read_email_data() # read the data\n",
    "G = hypergraph(L)     # construct a hypergraph from the data\n",
    "A0 = G.node_dimension_matrix().T[:10,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll do MCMC, finding that this matrix is preserved when we specify the `detailed = True` in the MCMC call. You can do this with either stub- or vertex-labeling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3132 epochs completed, 100011 steps taken, 1330856 steps rejected.\n"
     ]
    }
   ],
   "source": [
    "G.MH(n_steps = 100000, label = 'vertex', detailed = True, n_clash = 1)\n",
    "A1 = G.node_dimension_matrix().T[:10,]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(A0 == A1) # check that all entries are equal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Closing Thoughts\n",
    "\n",
    "While our particular illustration is quite simple, the recipe is extremely general. To conduct null hypothesis tests, \n",
    "\n",
    "1. Measure the empirical value of your quantity of interest on the observed data. \n",
    "2. Repeatedly randomize under the appropriate configuration model, measuring the value at appropriate intervals. \n",
    "3. Compare the empirical value to the null distribution.  \n",
    "\n",
    "Some care is required in choosing $n_\\mathrm{steps}$ to ensure that the chain is appropriately mixed (i.e. that the hypergraph is close to fully-randomized). To my knowledge, there are no general results on this. A practical approach is to track your quantity of interest over many iterations, and only use samples after they have \"leveled off,\" which suggests stationarity. "
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}