class26_bayesnet_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 26:  Bayesian Networks\n",
    "## Infer a Bayesian network from a matrix of discretized phospho-flow cytometry data.\n",
    "### Based on supplementary data from the 2005 article by Karen Sachs et al. (Science v308, 2005). \n",
    "\n",
    "In this class exercise, we will use the fundamental theorem for the likelihood of a Bayesian network structure for categorical variables, in order to score the posterior probability of the network shown in the Sachs et al. article (Figure 3A) vs. the phospho-flow cytometry data that the same authors provided in their supplementary data. The phospho-flow cytometry data have been already discretized for you (see \"class26_bayesnet_dataprep_R.ipynb\").  We will need to implement a single-vertex log-likelihood function using Theorem 1 from the article by Cooper & Herskovits in *Machine Learning* (volume 9, pages 309-347, 1992)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the tab-delimited data file of discretized phosphoprotein expression data (12 columns; first 11 columns are the expression levels -- \"low\", \"medium\", \"high\"; last column is the experiment identifier for the row; there are nine experiments). Print out the first six lines of the data frame, so you can see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>praf</th>\n",
       "      <th>pmek</th>\n",
       "      <th>plcg</th>\n",
       "      <th>PIP2</th>\n",
       "      <th>PIP3</th>\n",
       "      <th>p44.42</th>\n",
       "      <th>pakts473</th>\n",
       "      <th>PKA</th>\n",
       "      <th>PKC</th>\n",
       "      <th>P38</th>\n",
       "      <th>pjnk</th>\n",
       "      <th>expt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>high</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     praf    pmek    plcg PIP2    PIP3  p44.42 pakts473     PKA     PKC  \\\n",
       "0     low     low     low  low    high     low      low  medium  medium   \n",
       "1     low     low  medium  low     low  medium   medium  medium     low   \n",
       "2  medium    high  medium  low  medium  medium   medium  medium  medium   \n",
       "3  medium    high  medium  low     low     low      low  medium  medium   \n",
       "4     low  medium     low  low  medium  medium   medium     low     low   \n",
       "5     low     low  medium  low     low  medium   medium  medium  medium   \n",
       "\n",
       "      P38    pjnk       expt  \n",
       "0    high    high  1_cd3cd28  \n",
       "1     low    high  1_cd3cd28  \n",
       "2  medium  medium  1_cd3cd28  \n",
       "3  medium  medium  1_cd3cd28  \n",
       "4  medium    high  1_cd3cd28  \n",
       "5    high    high  1_cd3cd28  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas\n",
    "g_discret_data = pandas.read_csv(\"shared/sachs_data_discretized.txt\",\n",
    "                                 sep=\"\\t\")\n",
    "g_discret_data.head(n=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log likelihood function for the vertex degree distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "\n",
    "def log_prob_network_prior(network):\n",
    "    degarray = numpy.sum(network, axis=0) + numpy.sum(network, axis=1)\n",
    "    degarray_float = numpy.zeros(len(degarray))\n",
    "    degarray_float[:] = degarray\n",
    "    return numpy.sum(numpy.log(numpy.power(1.0 + degarray_float, -2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a vectorized log-factorial function, since it doesn't seem to be a builtin in python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.special\n",
    "\n",
    "def lfactorial(n):\n",
    "    return scipy.special.gammaln(n+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log likelihood function for a single vertex, based on Theorem 1 in the 1992 article in *Machine Learning* by Cooper & Herskovits. Note:  we are using igraph's adjacency matrix format which is the transpose of Newman's adjacency matrix definition!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 373,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def log_likelihood_network_vertex(network, vertex, discret_data):\n",
    "    # network is a NxN numpy matrix (N = 11 vertices)\n",
    "    # vertex is an integer\n",
    "    # discret_data is \"g_discret_data\" (N = 11 vertices, M = 7466 samples)\n",
    "    parents_vertex = numpy.where(network[:,vertex]==1)[0].tolist()\n",
    "    all_vertices = parents_vertex.copy()\n",
    "    all_vertices.append(vertex)\n",
    "    df1 = discret_data[all_vertices]\n",
    "\n",
    "    # change the name of the vertex column to \"vertex\"\n",
    "    df1_column_names = df1.columns.tolist()\n",
    "    df1_column_names[len(parents_vertex)] = \"vertex\"\n",
    "    df1.columns = df1_column_names\n",
    "\n",
    "    # count the data, grouped by all columns (parents & vertex)\n",
    "    df1 = df1.groupby(df1.columns.values.tolist()).size().reset_index(name='count')\n",
    "    \n",
    "    # make a new column called \"count factorial\" that is the log of the factorial of the count\n",
    "    df1[\"countfactorial\"] = lfactorial(df1[\"count\"].values)\n",
    "    \n",
    "    # drop the \"count\" column, as we no longer need it\n",
    "    df1 = df1.drop(\"count\", 1)\n",
    "    \n",
    "    if len(parents_vertex) > 0:\n",
    "\n",
    "        # sum up log-factorial-counts values for all possible states of \"vertex\" and its parent vertices, \n",
    "        # for each possible combination of parent vertices\n",
    "        nijkdf = df1.groupby(by=df1.columns[list(range(0,len(parents_vertex)))].tolist(), \n",
    "                    as_index=False).sum()\n",
    "        \n",
    "        # count number of cells with each possible combination of states for its parent vertices\n",
    "        df3 = discret_data[parents_vertex]\n",
    "        nijdf = df3.groupby(df3.columns.values.tolist(), as_index=False).size().reset_index()\n",
    "        nijdf_col_names = nijdf.columns.values\n",
    "        nijdf_col_names[len(nijdf_col_names)-1] = \"count\"\n",
    "        nijdf.columns = nijdf_col_names\n",
    "        \n",
    "        # compute the log factorial of the counts\n",
    "        nijdf[\"countfactorial\"] = math.log(2) - lfactorial(2 + nijdf[\"count\"])\n",
    "        \n",
    "        # drop the \"count\" column as we no longer need it\n",
    "        nijdf = nijdf.drop(\"count\", 1)\n",
    "        \n",
    "        # merge the two log-factorial-count values from nijdf and nijkdf, into two columns in a single dataframe\n",
    "        nmerge = nijdf.merge(nijkdf, how=\"outer\", \n",
    "                     on=nijkdf.columns[0:(len(nijkdf.columns)-1)].values.tolist(),\n",
    "                     copy=False)\n",
    "        \n",
    "        # sum the log-factorial-count values from nijdf and nijkdf\n",
    "        llh_res = numpy.sum(nmerge[\"countfactorial_x\"]+nmerge[\"countfactorial_y\"])\n",
    "    else:\n",
    "        # we handle the case of no parent vertices specially, to simplify the code\n",
    "        M = discret_data.shape[0]\n",
    "        llh_res = math.log(2) - lfactorial(M + 2) + numpy.sum(df1[\"countfactorial\"].values)\n",
    "    return llh_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log-posterior-probability function for the whole graph, using the per-vertex likelihood and the network structure prior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_posterior_prob_network(network, discret_data):\n",
    "    Nvert = network.shape[1]\n",
    "    lpvert_values = []\n",
    "    for i in range(0, Nvert):\n",
    "        lpvert_value = log_likelihood_network_vertex(network, i, discret_data)\n",
    "        lpvert_values.append(lpvert_value)\n",
    "    return log_prob_network_prior(network) + numpy.sum(numpy.array(lpvert_values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define an adjacency matrix for the \"real\" network shown in Fig. 3A of the Sachs et al. article (not including the \"missed\" edges which are the dotted arcs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          praf  pmek  plcg  PIP2  PIP3  p44.42  pakts473  PKA  PKC  P38  pjnk\n",
      "praf       0.0   1.0   0.0   0.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "pmek       0.0   0.0   0.0   0.0   0.0     1.0       0.0  0.0  0.0  0.0   0.0\n",
      "plcg       0.0   0.0   0.0   1.0   1.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "PIP2       0.0   0.0   0.0   0.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "PIP3       0.0   0.0   0.0   1.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "p44.42     0.0   0.0   0.0   0.0   0.0     0.0       1.0  0.0  0.0  0.0   0.0\n",
      "pakts473   0.0   0.0   0.0   0.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "PKA        1.0   1.0   0.0   0.0   0.0     1.0       1.0  0.0  0.0  1.0   1.0\n",
      "PKC        1.0   1.0   0.0   0.0   0.0     0.0       0.0  1.0  0.0  1.0   1.0\n",
      "P38        0.0   0.0   0.0   0.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n",
      "pjnk       0.0   0.0   0.0   0.0   0.0     0.0       0.0  0.0  0.0  0.0   0.0\n"
     ]
    }
   ],
   "source": [
    "real_network_adj = numpy.zeros(shape=[11,11])\n",
    "molec_names = g_discret_data[list(range(0,11))].columns.values    \n",
    "real_network_adj = pandas.DataFrame(real_network_adj, index=molec_names, columns=molec_names)\n",
    "real_network_adj[\"PKA\"][\"PKC\"] = 1\n",
    "real_network_adj[\"praf\"][\"PKC\"] = 1\n",
    "real_network_adj[\"pjnk\"][\"PKC\"] = 1\n",
    "real_network_adj[\"P38\"][\"PKC\"] = 1\n",
    "real_network_adj[\"pjnk\"][\"PKA\"] = 1\n",
    "real_network_adj[\"P38\"][\"PKA\"] = 1\n",
    "real_network_adj[\"praf\"][\"PKA\"] = 1\n",
    "real_network_adj[\"pmek\"][\"PKA\"] = 1\n",
    "real_network_adj[\"p44.42\"][\"PKA\"] = 1  # p44.42 = ERK\n",
    "real_network_adj[\"pakts473\"][\"PKA\"] = 1\n",
    "real_network_adj[\"pakts473\"][\"p44.42\"] = 1\n",
    "real_network_adj[\"pmek\"][\"PKC\"] = 1\n",
    "real_network_adj[\"pmek\"][\"praf\"] = 1\n",
    "real_network_adj[\"p44.42\"][\"pmek\"] = 1\n",
    "real_network_adj[\"PIP2\"][\"plcg\"] = 1\n",
    "real_network_adj[\"PIP3\"][\"plcg\"] = 1\n",
    "real_network_adj[\"PIP2\"][\"PIP3\"] = 1\n",
    "print(real_network_adj)\n",
    "real_network_adj = real_network_adj.as_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make an igraph network out of the adjacency matrix that you just created, and print the network summary and plot the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH D--- 11 17 -- '"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import igraph\n",
    "real_network_igraph = igraph.Graph.Adjacency(real_network_adj.tolist())\n",
    "real_network_igraph.summary()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the log posterior probability of the real network from Sachs *et al.* Figure 3A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-84131.0311113\n"
     ]
    }
   ],
   "source": [
    "lp_real = log_posterior_prob_network(real_network_adj, g_discret_data)\n",
    "print(lp_real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate 10000 random rewirings of the network -- eliminating any rewired digraphs that contain cycles -- and for each randomly rewired DAG, histogram the log ratio of the \"real\" network's posterior probability to the posterior probabilities of each of the random networks. Does it appear that the published network is pretty close to the maximum *a posteriori* (MAP) estimate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 369,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "lprobs_rand = []\n",
    "for _ in itertools.repeat(None, 10000):\n",
    "    graphcopy = real_network_igraph.copy()\n",
    "    graphcopy.rewire_edges(prob=1, loops=False, multiple=False)\n",
    "    if graphcopy.is_dag():\n",
    "        lprobs_rand.append(log_posterior_prob_network(numpy.array(graphcopy.get_adjacency().data), \n",
    "                                                      g_discret_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 370,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAF5CAYAAABX68ozAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcHVWZ//HPl81okABGSVwCMkRoFpFu8AcRZEAQdGRx\nGbQ1I4giiI6YwYVB0RlmmBEHaBlEBXXYMrQ6LsMyiCKLMBBxSCMYbMIqDYEITSBKQrPl+f1xqkl1\n5fZ2033r1s33/XrdV1KnTtV9Tt2+9z731DlVigjMzMzMqmK9sgMwMzMzGw8nL2ZmZlYpTl7MzMys\nUpy8mJmZWaU4eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlVipMXMzMzqxQnL2ZmZlYpTZe8SDpB\n0ipJZ+TKrsvKBh8vSLq4sN2mki6S9KSkJyRdKGla41tgZmZmk6mpkhdJuwEfB24rrArgXGALYAYw\nEzi6UKcbmA3sBewNbAdcOJnxmpmZWeNtUHYAgyRtDMwHPgacVKPKyoh4bJhttwMOAN4YEYuysqOA\nWyXNjoi7JylsMzMza7Bm6nk5G7gsIq4ZZv1cSY9Lul/S2YVTQnsA/YOJC0BE3AYsA+ZMXshmZmbW\naE3R8yLpA8CbgF2HqXI+cD+wFGgDTgN2BvbM1s/I1hUtzdbVes5XkHpr/gAM1Be5mZnZOmkKsBXw\n84h4vNFPXnryIum1wNeB/SLiuVp1IuKC3OJdku4GFkmaExE3jbD7GGHdAcB/jjtgMzMzG/Qh4OJR\na02w0pMXoAN4JdAjSVnZ+sBbJX0KeElEDElCIuL3kp4GZgE3MXwPy0xq98hA6nFh/vz5tLW1rXUj\nmtm8efPo6uoqO4xJ53a2nnWlrW5na1kX2tnb28vcuXMh+y5ttGZIXn4J7FQoOx/oBb5aTFwAJM0G\nXgo8lBUtAF4hacfcgN2dgc1IyU0tAwBtbW20t7evbRua2rRp01q+jeB2tqJ1pa1uZ2tZV9qZKWXY\nRenJS0SsAH6fL5O0Ang8InolbQ28D/gZ8DhpzEsXsBC4MdvHnZKuBM6VdAwg4BzSAGDPNDIzM2sh\nzTTbKC/f2/IscBDwK1L31IWknpYDC70yHwTuAa7P6i4GPtyIYM3MzKxxSu95qSUi9s39/yHShedG\n22Y5TlbMzMxaXrP2vNgE6uzsLDuEhnA7W8+60la3s7WsK+0sk2qMh10nSGoHFi5cuHBdGlhlZma2\n1np6eujo6ADoiIieRj+/e17MzMysUpy8mJmZWaU4eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlV\nipMXMzMzqxQnL2ZmZlYpTl7MzMysUpy8mJmZWaU4eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlV\nipMXMzMzq5QNyg7AzKzV9PX10d/fX3YYa5g+fTqzZs0qOwyztebkxcxsAvX19bHttm0MDKwsO5Q1\nTJnyMhYv7nUCY5Xn5MXMbAL19/dnict8oK3scHJ6GRiYS39/v5MXqzwnL2Zmk6INaC87CLOW5AG7\nZmZmVilOXszMzKxSmi55kXSCpFWSzsiVbSTpLEmPSXpK0iWSXlPY7nWSLsvWPyrpTEk+LWZmZtZi\nmip5kbQb8HHgtsKqM4EDgXcBHcD6wOWSlG23HnBFVrcdOAT4K+D0BoRtZmZmDdQ0yYukjUnD8z8G\nPJkr3wQ4EpgXETdHxGLgCGAHYL+s2gHAbODwiLgrIhYAxwNHZfs1MzOzFtE0yQtwNnBZRFxTKN+V\nNCvq2sGCiOgHbgfmZEW7A7dHxLLcdtcAU0g9NWZmZtYimmJMiKQPAG8iJSpFWwArImJFoXwpMCP7\n/4xs+UUR8WdJK3N1zMzMrAWUnrxIei3wdWC/iHhuHJvGBNUxMzOzCik9eSGd1nkl0DM4AJc0IPet\nkj5FGqi7saSphd6XmcAt2f+XArvkd5qNdZlKoUemaN68eUybNm1IWWdnJ52dnXU2x8zMrHV0d3fT\n3d09pGz58uUlRZM0Q/LyS2CnQtn5QC/wVWAJ8BywD3A5gKTp2TYnZPUXAF+QtHlu3MvbgAFg4UhP\n3tXVRXu7r4JpZmZWS60f9D09PXR0lDektPTkJetN+X2+TNIK4PGI6M2WvwecLqmfNBPpNGARcHW2\nyS+AxcB5kj4PbJ7VOTcinmpIQ8zMzKwhmmm2UV5xrMpxwJWknpdbgOeBgyMiACJiFem6LiL1tFya\n1f1cowI2MzOzxii956WWiNi3sPwcKYE5boRtHgIOnuTQzMzMrGTN2vNiZmZmVpOTFzMzM6sUJy9m\nZmZWKU5ezMzMrFKcvJiZmVmlOHkxMzOzSnHyYmZmZpXi5MXMzMwqxcmLmZmZVYqTFzMzM6sUJy9m\nZmZWKU5ezMzMrFKa8saMZmaj6evro7+/v+ww1tDb21t2CGYtz8mLmVVOX18f227bxsDAyrJDMbMS\nOHkxs8rp7+/PEpf5QFvZ4RRcAZxUdhBmLc3Ji5lVWBvQXnYQBT5tZDbZPGDXzMzMKsXJi5mZmVWK\nkxczMzOrFCcvZmZmVilOXszMzKxSnLyYmZlZpTh5MTMzs0px8mJmZmaVUnryIukTkhZJWiFpQNIt\nkg7Nrb9O0qrc4wVJFxf2samkiyQ9KekJSRdKmtb41piZmdlkKz15AR4A5pEulbkj8D/AjyXtkq0P\n4FxgC2AGMBM4urCPbmA2sBewN7AdcOGkR25mZmYNV/rtASLiikLRVyQdDXQAt2ZlKyPisVrbS9oO\nOAB4Y0QsysqOAm6VNDsi7p6k0M3MzKwEzdDz8iJJ60n6a2Bj4NrcqrmSHpd0v6SzC6eE9gD6BxMX\ngIi4DVgGzGlI4GZmZtYwpfe8AEjaEVgATAFWAodFxL3Z6vOB+4GlpFNLpwE7A3tm62dk64qWZuvM\nzMyshTRF8gLcSUpIpgLvBrolvS0ibomIC3L17pJ0N7BI0pyIuGmEfcZYnnjevHlMmzZ0bG9nZyed\nnZ3ja4GZmVkL6u7upru7e0jZ8uXLS4omaYrkJSKeB+7LFn8naXfgWODIGnV/L+lpYBZwE8P3sMyk\ndo/MEF1dXbS3t9cbupmZWUur9YO+p6eHjo6OkiJqsjEvBTV7TiTNBl4KPJQVLQBekZ16GqyzM7AZ\nKbkxMzOzFlJ6z4ukrwA/Bx4kDdQ9DNgPOEDS1sD7gJ8Bj5PGvHQBC4EbASLiTklXAudKOgYQcA5w\nmWcamZmZtZ5m6HnZCvgBaVDub0iJyyERcS3wLHAQ8CvgD6RrtywADoyIfM/MB4F7gOuzuouBDzcm\nfDMzM2uk0nteIuIjI6x7iHThudH2sRwnK2ZmZuuEZuh5MTMzMxszJy9mZmZWKU5ezMzMrFKcvJiZ\nmVmlOHkxMzOzSnHyYmZmZpXi5MXMzMwqxcmLmZmZVYqTFzMzM6sUJy9mZmZWKU5ezMzMrFKcvJiZ\nmVmlOHkxMzOzSnHyYmZmZpXi5MXMzMwqxcmLmZmZVYqTFzMzM6sUJy9mZmZWKU5ezMzMrFKcvJiZ\nmVmlOHkxMzOzSnHyYmZmZpXi5MXMzMwqpfTkRdInJC2StELSgKRbJB2aW7+RpLMkPSbpKUmXSHpN\nYR+vk3RZtv5RSWdK2qDxrTEzM7PJVnryAjwAzAPagB2B/wF+LGmXbP2ZwIHAu4AOYH3gckkCkLQe\ncEVWtx04BPgr4PRGNcDMzMwap/TeiYi4olD0FUlHAx2S7gWOBN4bETcDSDoCeBjYD7gKOACYDewd\nEcuyOscD3ZK+GBFPNaYlZmbNr7e3t+wQapo+fTqzZs0qOwyriNKTl7ysF+W9wMbAdaSelg2Aawfr\nRES/pNuBOaTkZXfg9sHEJXMNMCXb/lcNCd7MrKk9AqzH3Llzyw6kpilTXsbixb1OYGxMmiJ5kbQj\nsICUcKwEDouIeyTtBqyIiBWFTZYCM7L/z8iWXxQRf5a0MlfHzGwd9ySwCphPOkvfTHoZGJhLf3+/\nkxcbk6ZIXoA7gZ2BqcC7Sad83jZC/RjDPsdSx8xsHdNGGh5oVl1NkbxExPPAfdni7yTtDhwLXARM\nlTS10PsyE7gl+/9SYJfcOiRtTEqEhvTI1DJv3jymTZs2pKyzs5POzs56mmJmZtZSuru76e7uHlK2\nfPnykqJJmiJ5GUYAC4HngX2AywEkTQd2Ak7I6i0AviBp89y4l7cBA9n2I+rq6qK93b9CzMzMaqn1\ng76np4eOjo6SImqC5EXSV4CfAw+SBuoeRppJdEBE/EnS94DTJfWTTtqeBiwCrs528QtgMXCepM8D\nm2d1zvVMIzMzs9bTDNd52Qr4AXA/8BtS4nJIRAzOMDoOuJLU83ILqSfm4IgIgIhYRbqui0g9LZdm\ndT/XuCaYmZlZo5Te8xIRHxll/XOkBOa4Eeo8BBw8waGZmZlZE2qGnhczMzOzMXPyYmZmZpXi5MXM\nzMwqxcmLmZmZVYqTFzMzM6sUJy9mZmZWKU5ezMzMrFKcvJiZmVmlOHkxMzOzSnHyYmZmZpXi5MXM\nzMwqxcmLmZmZVYqTFzMzM6uUupIXSX8jacpEB2NmZmY2mnp7Xk4Hlko6R9KbJzIgMzMzs5HUm7y8\nGjgSmAn8r6RFko6X9MqJC83MzMxsTXUlLxHxfET8JCIOBmYBFwIfBR6S9BNJ75Lk8TRmZmY24dY6\nwYiIpcDVwLVAALsC/wncI2mftd2/mZmZWV7dyYuk6ZI+I+k24EbgVcChwJak00k/Av5jQqI0MzMz\ny2xQz0aSfgq8E7gf+C5wQUQ8lquyUtJpwGfXPkQzMzOz1epKXoA/AftFxA0j1HkMmF3n/s3MzMxq\nqit5iYjDx1AngHvr2b+ZmZnZcOq9SN1Zkj5eo/xoSWeufVhmZmZmtdU7YPcw4KYa5QuA99cfjpmZ\nmdnI6k1eNgWW1yh/AthsPDuS9GVJCyU9LWm5pCsktRXqXCdpVe7xgqSLC3U2lXSRpCclPSHpQknT\nxt0yMzMza2r1Ji/3APvVKN8fuG+c+2oH/g3YHtgdeBq4RtLUXJ0AzgW2AGaQpmIfXdhPN2mA8F7A\n3sB2pIvnmZmZWQupd7bR14F/lzSddHE6gH2Bk4DPjWdHEXFoflnSR4FlwB7AL3OrVhamY+e32Q44\nAHhjRCzKyo4CbpU0OyLuHk9MZmZm1rzqnW30neyu0icCp2bFjwJfiohvr2VM00k9LcsK5XMlHU6a\npn0FcGJEDJ662gPoH0xcshhvk7QMmAM4eTEzM2sR9fa8EBFnAWdJmpktPzJBMXUBN0RET67sfNIF\n8ZYCbcBpwM7Antn6Gdm6oqXZOjMzM2sRdScvgyYwaUHS2cAOwFsKz3FBbvEuSXcDiyTNiYhas55e\n3HS055w3bx7Tpg0d19vZ2UlnZ+fYAzczM2tR3d3ddHd3DylbvrzWnJ3Gqff2ADNIPST7kmYeKb8+\nIjaqY59nAe8C9hotIYqI30t6mnRH65sYvodlJrV7ZF7U1dVFe3v7eMM1MzNbJ9T6Qd/T00NHR0dJ\nEdXf83Ix6UaM/0Aa67JWJH0DOATYOyL6xlB/NvBS4KGsaAHwCkk75gbs7kyatj1Sz4yZmZlVTL3J\nyxxg1/wA2XpJ+ibQCRwMrJC0RbZqeUQMSNoaeB/wM+Bx0piXLmAh6W7WRMSdkq4EzpV0DKkn6Bzg\nMs80MjMzay31XuflfsYwnmSMjgE2Aa4DHs49DsvWPwscBPwK+APp2i0LgAOz+ycN+iDp+jPXZ3UX\nAx+eoBjNzMysSdTb8zIP+FdJR0XEH9cmgIgYMYGKiIdIF54bbT/LcbJiZmbW8upNXi4CXg4skfQk\nsCq/MiJetbaBmZmZmdVSb/JywoRGYWZmZjZG9V5h93sTHYiZmZnZWNQ7YBdJbZJOlXSxpFdlZftL\n2n7iwjMzMzMbqq7kRdIBpKnKuwDvBTbOVu0EnDIxoZmZmZmtqd6el38B/j4i3k6ayjzoatI1YMzM\nzMwmRb3Jyw7AT2qU95NuF2BmZmY2KepNXvpJtwco2onVl+w3MzMzm3D1Ji/fB06VtFm2HJLeSLps\n//wJiczMzMyshnqTlxNJ9xl6FJgK/A64Ffgt8E8TE5qZmZnZmuq9zsuzwPslvQFoJyVBPRFx50QG\nZ2ZmZlZU7xV2AYiIu4C7JigWMzMzs1HVlbxIOnek9RHx8frCMbNm09fXR39/f9lhDNHb21t2CGZW\nonp7XmYWltcjTZ/eFLh+rSIys6bR19fHttu2MTCwsuxQzMxeVO+Yl4OKZZIEfAO4d22DMrPm0N/f\nnyUu84G2ssPJuQI4qewgzKwkazXmJS8iQtJpwE3AGRO1XzNrBm2ksfnNwqeNzNZldd+YcRhbsfo+\nR2ZmZmYTrt4Bu1+rUfxK4FDgx2sVkZmZmdkI6j1ttEdhOYDlpAvUfWOtIjIzMzMbQb0Ddvea6EDM\nzMzMxmKix7yYmZmZTap6x7z8H+lU0agi4s31PIeZmZlZLfWOebkWOBpYBCzIynYHdgLOAZ5Z+9DM\nzMzM1lRv8rIpcEZE/GO+UNJJwJYR8bGx7kjSl4FDgO2BZ4EbgeMjojdXZyPgdOADwEuBq4FjI2JJ\nrs7rgG8C+wArge5sP8/X1UIzMzNrSvWOeTkMOK9G+QXAX49zX+3Av5GSl92Bp4FrJE3N1TkTOBB4\nF9ABrA9cnl3VF0nrkS65Obi/Q4C/IiU8ZmZm1kLq7Xl5hpRo9BXKd2ecp4wi4tD8sqSPAstI07F/\nKWkT4EjgvRFxc1bnCOBhYD/gKuAAYDawd0Qsy+ocD3RL+mJEPDWu1pmZmVnTqjd5+XfgPEm7AjeT\nBu/uDhwL/OtaxjQ929+ybLkji/PawQoR0S/pdmAOKXnZHbh9MHHJXANMybb/1VrGZGZmZk2i3uu8\nnCLpj8Angc9kxb8H5kXEd9Yypi7ghojoyZZnACsiYkWh3tJs3WCdpYUY/yxpZa6OmZmZtYC6b8wY\nEd8FvjuBsSDpbGAHYM+xhDBBdczMzKxC6k5eJE0nDc59LXBaRDwhaWfgsYh4uI79nUUakLtXYful\nwFRJUwu9LzOBW3J1dinsb2NgKoUemaJ58+Yxbdq0IWWdnZ10dnaOtwlmZmYtp7u7m+7u7iFly5cv\nLymapN6L1O0K/Bx4HHg98D3gCdIsn22BD41zf9/Itt07IoqDgBcCz5OmQF+e1Z9OuqbMCVmdBcAX\nJG2eG/fyNmAg235YXV1dtLe3jydcMzOzdUatH/Q9PT10dHSUFFH9U6XPAM6NiDeQEoRBVwB7j2dH\nkr5JSnY+CKyQtEX2mAIQEX8iJUenS9pd0nbA+aQL5F2d7eYXwGLSIOJtJe0BnJbF6JlGZmZmLaTe\n5GU30pV0ix4BXjXOfR0DbAJcR5r+PPg4LFfnOOBKUs/LLaSemIMjIgAiYhXpui4i9bRcmtX93Dhj\nMTMzsyZX75iXPwMb1yj/C+DR8ewoIkZNoCLiOVICc9wIdR4CDh7Pc5uZmVn11NvzcilwoqT1s+XI\nxqF8FfjJhERmZmZmVkO9ycvxwFbAEtK9hn4GPAisAk6ckMjMzMzMaqj3InXLJb0F2J90L6H1gB7g\n54PjUMzMzMwmw7iTF0kbkHpaPhERvyDN9DEzMzNriHGfNoqI50kXhHt24sMxMzMzG1m9Y17mA0dM\nYBxmZmZmY1L37QGA4yXtT7quypBemIj4/FpFZWZmZjaMepOXDuD23P/zPGDXzMzMJs24khdJWwP3\nR8RekxSPmZmZ2YjGO+blbuCVgwuSfiBpi4kNyczMzGx4401eVFh+JzB1gmIxMzMzG1W9s43MzMzM\nSjHe5CVYc0CuB+iamZlZw4x3tpGA8yU9ky1PAb4taUW+UkS8ZyKCMzMzMysab/JyQWF5/kQFYmZm\nZjYW40peIuIjkxWImZmZ2Vh4wK6ZmZlVipMXMzMzqxQnL2ZmZlYpa3NjRjMzswnT29tbdghrmD59\nOrNmzSo7DCtw8mJmZiV7BFiPuXPnlh3IGqZMeRmLF/c6gWkyTl7MzKxkTwKrSFffaCs5lrxeBgbm\n0t/f7+SlyTh5MTOzJtEGtJcdhFWAB+yamZlZpTRF8iJpL0mXSloiaZWkgwvrz8vK84+bCnU2knSW\npMckPSXpEkmvaWxLzMzMbLI1RfICTAV+CxzL8Dd6/BmwBTAje7yzsP5M4EDgXUAHsD5wuSRNRsBm\nZmZWjqYY8xIRVwJXAoyQbDwTEY/VWiFpE+BI4L0RcXNWdgTwMLAfcNVEx2xmZmblaJael7H4S0mP\nSnpI0kWFU0IdpETs2sGCiOgHbgfmNDhOMzMzm0RVSV4uAz4A7Al8BNgW+JWkKdn6GcCKiFhR2G5p\nts7MzMxaRFOcNhpNRPwkt3iXpB7gQeA9wMUjbTravufNm8e0adOGlHV2dtLZ2VlPqGZmZi2lu7ub\n7u7uIWXLly8vKZqkEslLUUQ8Lul+YPCqQUuBqZKmFnpfZgK3jLSvrq4u2tt9XQEzM7Naav2g7+np\noaOjo6SIqnPaaAhJ04AtgYeyooXA88A+uTrTgZ2AGxseoJmZmU2apuh5kTQV2AYYnGm0taSdgWXZ\n4wTgv0k3wNgS+CrwOPBTgIj4k6TvAadL6idda/o0YBFwdQObYmZmZpOsWXpedgVuJfWgBHA60AP8\nI/ACsDvpOi8PAJcAS4C3Fk4RHUeabn056VTR88DBETHquBczMzOrjqboeYmIXzFyIrX/GPbxHCmB\nOW6i4jIzM7Pm0yw9L2ZmZmZj4uTFzMzMKsXJi5mZmVWKkxczMzOrFCcvZmZmVilOXszMzKxSnLyY\nmZlZpTh5MTMzs0px8mJmZmaV4uTFzMzMKsXJi5mZmVVKU9zbyGxd19fXR39/f9lhrKG3t7fsEMzM\n1uDkxaxkfX19bLttGwMDK8sOxcysEpy8mJWsv78/S1zmA21lh1NwBXBS2UGYmQ3h5MWsabQB7WUH\nUeDTRmbWfDxg18zMzCrFyYuZmZlVipMXMzMzqxQnL2ZmZlYpTl7MzMysUpy8mJmZWaU4eTEzM7NK\ncfJiZmZmldIUyYukvSRdKmmJpFWSDq5R5x+y9SslXSNp+8L6TSVdJOlJSU9IulDStMa1wszMzBqh\nKZIXYCrwW+BYIIorJX0hW3cEsCPwB+AqSVNz1bqB2cBewN7AdsCFkxm0mZmZNV5T3B4gIq4ErgSQ\npBpVjgNOjoirsjpHAw8BHwS+I6kNOAB4Y0QsyuocBdwqaXZE3N2AZpiZmVkDNEvPy7AkvR6YAVwz\nWBYRzwE3AHOyot2B/sHEJatzG7AsV8fMzMxaQNMnL6TEJYClhfKl2brBOsX1xTpmZmbWApritFGd\n1hgbU0+defPmMW3a0HG9nZ2ddHZ21huXmZlZy+ju7qa7u3tI2fLly0uKJqlC8rIUEKkHZVmufCar\ne1uG62HJ16mpq6uL9vb2CQjTzMys9dT6Qd/T00NHR0dJEVXgtFFE3E9KQPYdLJO0IbAncGNWtAB4\nhaQdc3V2BjYDbmpctGZmZjbZmqLnJZvyvA2phwVg6yz5WBYRDwJfB06StBi4DzgReJ40PZqIuFPS\nlcC5ko7J9nMOcJlnGpmZmbWWpkhegF2Ba0ljVAI4PSu/ADgyIr4maUq2vClwM/D2iFiR28cHgbOA\n67PlS4C/bUDsZmZm1kBNkbxExK8Y5RRWRJwMnDzC+uXAhyc4NDMzM2syTT/mxczMzCzPyYuZmZlV\nipMXMzMzqxQnL2ZmZlYpTl7MzMysUpy8mJmZWaU4eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlV\nipMXMzMzqxQnL2ZmZlYpTl7MzMysUpy8mJmZWaU4eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlV\nipMXMzMzqxQnL2ZmZlYpTl7MzMysUpy8mJmZWaU4eTEzM7NKcfJiZmZmlVKJ5EXSVyStKjweLtT5\nB0lLJK2UdI2k7cuK18zMzCZPJZKXzCJgC2BG9thpcIWkLwDHAkcAOwJ/AK6SNLXhUZqZmdmkqlLy\n8nxEPBYRj2aPx3PrjgNOjoirIuI+4GhgA+CDpURqZmZmk6ZKyctsSY9IWirpp5LaACS9ntQTc81g\nxYh4DrgBmFNOqGZmZjZZqpK83ETqRflL4H3AS4HrJQ2eQgpgaWGbpdk6MzMzayEblB3AWETEVbnF\nxZIOBu4FjgSuHW6zsex73rx5TJs2bUhZZ2cnnZ2d9YRqZmbWUrq7u+nu7h5Stnz58pKiSSqRvBRF\nxLOSbgdmkXpYROplWZarNpM1e2PW0NXVRXt7+6TEac2lr6+P/v7+ssNYQ29vb9khmJkNq9YP+p6e\nHjo6OkqKqKLJi6QNgO2BmyLifklLgX2B32frNwT2BL5UXpTWTPr6+th22zYGBlaWHYqZma2lSiQv\nkv4ZuAJ4AHgVKSnZDLgwq/J14CRJi4H7gBOB54HuNfdm66L+/v4scZkPtJUdTsEVwEllB2FmVhmV\nSF6AN5DGt0wHngJuBPaOiAcBIuJrkqYAFwCbAjcDb4+IFSXFa02rDWi204Q+bWRmNh6VSF4i4rAx\n1DkZOLkB4ZiZmVmJqjJV2szMzAxw8mJmZmYV4+TFzMzMKsXJi5mZmVWKkxczMzOrFCcvZmZmVilO\nXszMzKxSnLyYmZlZpTh5MTMzs0px8mJmZmaV4uTFzMzMKsXJi5mZmVWKkxczMzOrlErcVdrMzKws\nvb29ZYdQ0/Tp05k1a1bZYZTCyYuZmVlNjwDrMXfu3LIDqWnKlJexeHHvOpnAOHkxMzOr6UlgFTAf\naCs5lqJeBgbm0t/f7+TFzMzMitqA9rKDsBwP2DUzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlVigfs\n2oTr6+ujv7+/7DCGaNbrNJiZ2fg5eVkHdHd309nZ2ZDn6uvrY9tt2xgYWNmQ51s3dQONeT3Lt660\n1e1sLetKO8vTcsmLpGOBzwIzgUXAvIj433KjKlcjk5f+/v4scSnjugjzgK5h1l0BnNTAWCbTuvTB\nuK601e1sLetKO8vTUsmLpPcDpwMfBm4GjgZ+JqktIh4qNbh1ThnXRZg2wnP6tJGZWatoqeSF9NP7\n7Ij4r2wFZaozAAAT0klEQVT5i5LeCXwC+GJ5YU28yy+/nCVLloyp7gMPPMA555wzyRElfX19DXke\nMzNbd7VM8iJpQ6ADOLmw6mpgTuMjmjy33XYbBx10EKDsMZpVHHPMsZMc1ernMjMzm0wtk7wA04H1\ngaWF8qXAjBr1p0A1Z6Hcfvvt2f+OA7YYwxYXks6kNcK3gD7SGJNGH9uHgP8cZt2N2b9lxDWa8cY2\nUjsnWtnHbbi2lh3XSOqJrRGvaTMcs6q9nvXG1YjX836gvO+w3PNOKeP5FRFlPO+EkzQTWALsGhE9\nufLjgY9FRFuh/gdp3DeAmZlZK/pQRFzc6CdtpZ6XfuAF1uxlmcmavTEAPwc+BPwBGJjUyMzMzFrL\nFGAr0ndpw7VMzwuApF8D/xsRn82VLQSujIiWGrBrZma2rmqlnheAM4DzJN0M/AY4CpgNHFpqVGZm\nZjZhWip5iYgfStoc+Brp9NEi4B0R8WC5kZmZmdlEaanTRmZmZtb6fFdpMzMzqxQnL2ZmZlYpLZW8\nSNpK0oWS+iQ9m/17qqSNCvV2knSdpJWSHpS0xh37JL1X0h2SBiQtkrTGoF9J/yBpSbafayRtP5nt\nq4ekYyXdJ+lpSf8nac+yYxqOpC9LWpjFulzSFZKK1+fZSNJZkh6T9JSkSyS9plDndZIuy9Y/KulM\nSRsU6uwt6Zbsue6RdHQj2liLpBMkrZJ0Rq6sJdop6dWSLpLUn72XbpPUXqgz4vtI0qbZPp6U9ET2\nHp9WqDPqe3qySHqJpDOyNjybPf8pktYr1KtUOyXtJenSLOZVkg6uUachbRrL5/FktVXSBtnre4ek\nZyQ9LukHNd6PTd/Wsbymubrfzup8uinbGREt8wDeDpwN7AW8JlteAnwjV+flwCPAd0hz1N8OPEa6\n+/RgnT2AZ0n3RJoFfCpb3i1X5wvAo8D+wNbAf2TPNbXs45CL8f3A08BfZ+04Bfgz8NqyYxsm3v8G\nPgC8nnRnxx9nr9XUXJ1vAXcD/w/YFrgcuJXV47fWA34HXAa8IXst7wHOzO1jK+Ap0q0kZuWO07tL\naPNuwH1ZG85opXYCm5IuA3oWsH32ntwL2CpXZ9T3EfAz4NfATsAbSTMJL82tH/U9PcntPDV7/gOz\n4/xu4Eng81VuZ9aek4FDSNfQOriwviFtYgyfx5PZVmAT4JJs3SxgF+B64PbCPpq+raO9prl67yZ9\n3jwIfLoZ2znpb+yyH8BngIdyy5/IDuwGubLjgAdzy98HflzYz0+B/8wtPwx8Kre8IfBH4Kiy25yL\n6dfAaYWyW4FTyo5tjPFvSrpZ0n7Z8ibAM8C7cnWmZ3/0+2fL7yBddHDzXJ1DgJXAxtnyqcDCwnN1\nATc2uH0bA4uBfYFryZKXVmkn8FXSNZZGqjPi+4iUxK4CdszV2Tkrm50tj/qenuR2/g/wzUJZNzC/\nVdqZxVFMXhrSJsbweTzZba1R501ZvW2q2tbh2kn6kdGXtel+cskLsF2ztLOlThsNYzqwLLe8O3BD\nRDyfK7saeLWkLbPlPbIyCnXmAEjamjQV+5rBlRHxHHADTXITSK2+UeU1hVVVulHldCBY/fp1kKb3\nXztYISL6gdtZ3abdSb+I8q/5NaSrQXbk6tR6fXeVtP5ENmAUZwOXRUTxNdqV1mjnQcAtkn6odPrr\nTkmfGVwp6fWM/j7aHeiPiEW5OreR/ibydUZ7T0+mS4H9JM0GkLQzqYfp0my5Vdr5oga3acTP45JM\nJ31hP5ktt0RbJYl0M7yvRUStmybtQZO0s6WTF0l/QeqO6soVz6D2zRvF6lsLDFdncP0WpC/Vsd4E\nsgzjvVFlM+oivQkG71U1A1gRESsK9fJtWuO1i4g/k3okRnt9NyAdt0kn6QOkX29/X2P1FrRGO7cG\n/o50euutpNOWp0g6JhffaO+jWm0YS53ie3rSRMQ5wI+AxZKeBRYCZ0XED3PxVb6dBY1s02ifxw0l\n6SXAv5J6Cfqz4lZp6wnAsxHxjWHWN007K3GROklfAb4yQpUgnSvL35Dx1aRzcz+IiPNGeYqxXOxm\nouqUrQoxIulsYAdgLAOMJ+K1adhxkfRa4Ouk02HPjWPTSrWT9OPohoj4p2y5V9IbgY8D3x5hu0q1\nU9LfAX9DGifQS0pKz5L0x4g4f4RNK9XOMWpUm0ppt9KA+B+Qvog/OYZNKtNWSR3Ap0ljesar4e2s\nSs/LWaRzbcM92khX0wVeTFyuIZ3bL86uqJXdzSQduEdGqbM0t77Wr53hbgJZhvHeqLJpSDoLeBfw\nlxHxcG7VUmCqpKmFTYqvzZA2S9oYmMror+/zpOM22TqAVwI9kp6T9BywN3Bc9sv9j8DGLdDOR0hj\nevLuJJ1Th7G9j4b7NTZancH3dCP+1r8I/GNEXBIRd2U9Lqdn5YPxtUI78xrRprF+HjdElrj8kDQI\ndb+IeCq3uhXauifpc+nB3OfSlsAZku4bIb5ijA1pZyWSl4hYln0ojPR4FiCbvnYtcAtwZI3dLQD2\nLJzzfxvwcET05ersW9huX+CmLJ77SQf5xTrZGJM9gRvXusETIPtFv5A127EPWTuakaRvkO5FtU/u\n9Ri0kPTFu0+u/nTSqPfB474AeKPSbSIGvY00uLUnV6d4XN4G3BIRL0xEO0bxS1LMbyINdtuZ9Pc6\nP/f/56h+O28k3Vss7w3AQzDm99EC4BWSdszV2RnYjNV/xyO9px+YsNYMbyPSD4W8VVl5K7XzRQ1q\n05g+jxshS1z+C9gG2DcinixUaYW2XkiaPbRz7vEw6XY7B+Tia452TvQI5jIfpMztbuAq0q+7LQYf\nuTqbZC/IuaRz8m8nTff7TK7OHqTZHoPTuD5J+lLYNVfn86RfyPsDfwF8j/Sh3ExTpQ8DVpCmSm8J\n/DPwJ+B1Zcc2TLzfBJ4gDXbcIveYUqizmDQobDvSFOKFDJ1CfBtpauO22Wt5N/D13D62Ik0Z/8fs\n9R08ToeW2PYXZxu1SjtJA48HSONeXge8N/v7+2iuzqjvI9JsnptY/cH6a+C/c+tHfU9Pcjvnk6a7\n70f63HkX6Rfm16vcTlIv3s6snlnzmWz5dY1sE2P4PJ7MtpLGDl4CPED6AZH/bNqwSm0d7TWtUX/I\nbKNmauekv7Eb+QAOJ/0Cyj9WAS8U6u0AXEca3LgE+FKNfb2HdP56ALgDOKRGnS9nL9JK0pfP9mUf\ngxoxHpP9AT4N/B/wlrJjGiHWVTVevxeAD+fqbAicSTrt8RTp2jCvKezntaSZHk+Rri/Qlf+Qyers\nRUoGngbupeQp7qTTnPnkpSXaCbyTNEtqoNYHYVZnxPcRMI30q/DJ7HEBsEmhzqjv6Uls48akmWMP\nk6azPwj8G7BRldtJOpVZ6z35H41uE2P4PJ6stpJ++BXXDS6/tUptHctrWqh/H2smL03RTt+Y0czM\nzCqlEmNezMzMzAY5eTEzM7NKcfJiZmZmleLkxczMzCrFyYuZmZlVipMXMzMzqxQnL2ZmZlYpTl7M\nzMysUpy8mJmZWaU4eTHLkXStpDMmYb/XS/pQo7Yri6TDJT1RKFslaVYd+3qVpP7sZqsTEduGku6R\n9JaJ2F9ZJO2dHdNNcmX3S3prHfvaSNJDkjoK5adL+veJiNdsMjh5MZtkkg4k3bDv4lzZH7IvoFWS\nnpW0WNJnRtsuK99F0vclPZxt+6ikn0n6QCPaMwbFe44MWZZ0Xa7tA9mX56WS3j1ko4hHgW7gi7We\nRNJ5kv4l+/+qwj5vk/TBwiYfBZZExI2F/ewj6bLsOD6bHdefZMe/WY14X5fs2KyS9ELWpqWSfiHp\nI5L04k4iniXdE+ufCrs4DfjIRCWOZhPNyYvZ5PtbYH4MvZFYAF8CZgDbABcBZ0j68EjbSTqEdDv5\n50h3C3896eaLFwOfyv8aH4/C7esnmgrLQbrj7Baku86+h3Rjtu9L+nah7oXA4ZJePmSH0nrAX5Fu\nWDnocNLx3B74X2C+pL1z6/8WOK+wn2NJd6G/l3QTya1Id0n+BXDieBqZ26fyCUKJfkY6HlsCB5Ju\n/nkmcFl2/AbNB/aXtM1gQUQ8AlxPurGrWfOZrDuS+uFHFR+kO+Oekf1/U9KX5zJgBXAFsE2h/lFA\nH+nOzj8k3WL+idz6aaQ7DbcXtqt1q/lFwE9yzz1kO+BlpLtH17wD7DjaOHhn2beT7jQ+AOydrTsI\nuIV0F+qHgVOBDXLbfpaUaDxDuuP1ecDLc+sPB5YVnm8VMKvWMS7UOyKru2+h/GGgs1C2F/BI4TkO\nzi2vR7rj7eBruSPp7rmb5+q8LmvHl9fyeB4OPEFKpu7IXrctgd1ICcOT2fFcAMypcWw+Cvwk+xvr\nAw4r1HknsJh0h96rs+d7gdydfLO/p/wdjs8b/Fsq7Guf7DmPLJTfBPx9oezjwOKy35N++FHr4Z4X\ns+H9mPSltz/QTurt+KWkKQCS3g58C/hXoI2U3HyRoV36byF90dw+hud7BthwhO0OADYHvl5fc9Zw\nCinZegNwm6QDSF96p5F6g94PHMrQUwpPAR8j9Zi8D9gVmKixEReQkoD3FMp/Q0q48g4CLhtuRxGx\nCnie1cfzrcAfImJZrtr7gPWZmPhfBvwdMJf0t/BoVvZNYGfgTcBvgcslTStsexKp7dsBPwDOl/QK\nAEnbkhKbH2f7PQv4l3qDjIhrgdsY2zH+DbCNpBn1Pp/ZZHHyYlaDpB1Iv1KPjIiFEbEY+BtgOjA4\nNuNo4KcR8a2IeDAizgcuLexqNqmH4PkRnkuSOklfcFdnxdvU2G529m9vbttdJf0593jnOJr5xYi4\nMSL6IuJJ0mmSUyLi+xGxJCJuIH2xfnxwg4j4dkQsyNZfl61/7ziec1gREcBdpFM3eQ+QEqy8Q4BL\nau1H0vqSPgdsxurjOTvbT95s4I9Z2we3fU/heO4wxvA3AD4REbdGxL0R8XRE/CoifhQRD2R/P58k\nJbb7FbY9NyIuiYgHScdzQ2BOtu4ooDciTsz289/Ad8cY03DuZGzH+AHSKb9iuVnpNig7ALMmtQMw\nEBG/HSyIiD9JuoPUGwNpbMV5he1+w+rkBuDlpN6KWk6VdArpy+oF4GxW9wKMtF3ebaRf9gD3MPb3\ndAA9hbIOYFdJJ+fK1gc2kvSyiFgp6R3AF0htn5I930skTY2IFWN87pGINQej/ol0PFIFqQ14NfDL\nQr1uSauAl5BOwXwl+7KHsR/PK0nH87Wk01tjHQv0dETclS/IeixOJfVobEb6sfhS0jiUvN8P/ici\nBiQ9RTptCOk4/7pQf8EYYxrOqMc4V0aNcrPSOXkxm1xPARsPs+7fgPNJp6OKvSy1trs7+3d7UtJC\nRDwH3AdQxxjRZwvL65ESk/8pVswSl78g9XacBhxPGssxJ2vDWg/4zQaRziYlgHmbMDTxOAi4KiKe\nKdT7DKmn5Rng4awnZ9BTpFNdeXcDW0jaLCKegNRO4D5JL7DmQOORFI8lwPeBjUhjWh4knca6kTWP\n1aoa2yr374gzi+rQRhojk1c8xoNl1Cg3K51PG5nVdgcwRdIugwXZWIUdgN9lRb8nDcrMe3Nh+W5g\n5jCzefoj4r7slFPxtFKt7X5BGhMyb3xNGbMeYHYW05BHtv7NwPPZKYyFEXEva/YirI0jSD0OPyqU\nb8nqxA3SKaPi6TlIp4Duy05pFb/w7waK15r5ESmhOK7uiEe2B3BmRFyd9cqsBF45zn3cAexeY791\nkbQvsBOjH+PBsqhRblY697yY1RARd0i6BviupI+TTkOcSprtM3gq4tvAFZI+AVwO7AsczNBfyoO/\ntN8I3DqOENbYLiJWSPoYaUrxRqTBoPeTemjekT3vC2Pcf61ehZOBSyUtIX25vQDsAuwUESeTphO/\nRNKnScdgd9L043q8TNIWpM+g15IGkH4G+GZEXF+ouxvweUgXriMNEj5knM93PfDvhV6WByV9FuiS\ntDlpuvoS0qDouYzveNZyL+laKbeRThd9jTSzazy+A3wyO714Lmng+EfHuO1LsmO8Pmla+juAE0iJ\n30WFum9mzR63NwP3RcTD44zZbNK558VsqHzi8T5SL8tVpOnDGwD7RcQAQERcRboOxomkQbQHk04F\nvfgFlQ0GvSJbN9zzrBnEMNtlYzjmkL6Q/ot0yugG0rTn90fEGqd8xtDOwX3/gnRK5iDSaak7gL8H\nlmbrf0OaTfUlUq/TEdn6ehxFmgJ9D2kmzXbAX0fEkGRI0m6k3pjBmUUHATdHRP9o7RmyMuJ3pNey\neDy/QZrFtTXpeN9P6uGaDRwQEXeMu2WrfZjUM/Vb0sX2vkWahTRa3C+WZQN930v6W+wl9RLVvGhf\nDQeSjvH9pGu+7A18KiIOzfdMZQnhbqS/p7yDSdP/zZqO1uxdNbN6SfoGqadi71zZ/qRemm1qnM4Y\naV91bddsskG0W0VEXx3bngVsGBHHZMuXADdExGl17OsoYG7+tWkVku4HDq/RazWWbY8H9o+IA3Nl\nM0kzv7bPZkGZNRX3vJitBUmflrS9pFnZ6aWPkgawvijroXkAKF6ufkT1btcqsh6BDwD/nCu+gdSL\nUY/zgBmS5oxacx2RnX6s1ZvzWeA8Jy7WrNzzYrYWJP2I1B0/jTSj5JsRcXrJMX2LNGajKEi3Gzi2\nwfG8ALy+np6XZiDpCtIVfYsC+JeI+GqDQ1qDpPuAI+rpeTGrIicvZi1G0nRWT3Mt+lON8SI2guwU\nykuHWb0sf5E7M2sMJy9mZmZWKR7zYmZmZpXi5MXMzMwqxcmLmZmZVYqTFzMzM6sUJy9mZmZWKU5e\nzMzMrFKcvJiZmVml/H8SNwXVpgzWpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff28bfff0b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot\n",
    "matplotlib.pyplot.hist(lp_real - lprobs_rand)\n",
    "matplotlib.pyplot.xlabel(\"log(P(G_real|D)/P(G_rand|D))\")\n",
    "matplotlib.pyplot.ylabel(\"Frequency\")\n",
    "matplotlib.pyplot.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
}