Implementing an ERGM from scratch

ERGM_from_scratch.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Implementing a toy ERGM from scratch\n",
      "\n",
      "### David Masad \n",
      "\n",
      "Code to go along with blog post:\n",
      "http://davidmasad.com/blog/ergms-from-scratch/"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import networkx as nx"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "G = nx.florentine_families_graph()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx.draw(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAE+CAYAAADMGBp/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4TOfbwPHvTPZJJgRJSCJB0NROVSxVW4tEFaX21lJL\nUW1V0dYvqiiqomipUkIstbxpLSGoEKIItW+1L0EtIc0iiSzzvH9MpImlgmQm4v5c11wmc55zzn0k\nM/c859k0SimFEEIIIUxGa+4AhBBCiOeNJF8hhBDCxCT5CiGEECYmyVcIIYQwMUm+QgghhIlJ8hVC\nCCFMTJKvEEIIYWKSfIUQQggTk+QrhBBCmJgkXyGEEMLEJPkKIYQQJibJVwghhDAxSb5CCCGEiUny\nFUIIIUxMkq8QQghhYpJ8hRBCCBOT5CuEEEKYmCRfIYQQwsQk+QohhBAmJslXCCGEMDFJvkIIIYSJ\nSfIVQgghTEySrxBCCGFiknyFEEIIE5PkK4QQQpiYJF8hhBDCxCT5CiGEECYmyVcIIYQwMUm+Qggh\nhIlJ8hVCCCFMTJKvEEIIYWKSfIUQQggTk+QrhBBCmJgkXyGEEMLEJPkKIYQQJibJVwghhDAxSb5C\nCCGEiUnyFUIIIUxMkq8QQghhYpJ8hRBCCBOT5CuEEEKYmCRfIYQQwsQk+QohhBAmJslXCCGEMDFJ\nvkIIIYSJSfIVQgghTEySrxBCCGFiknyFEEIIE5PkK4QQQpiYJF8hhBDCxCT5CiGEECYmyVcIIYQw\nMUm+QgghhIlJ8hVCCCFMTJKvEEIIYWKSfIUQQggTk+QrhBBCmJgkXyGEEMLEJPkKIYQQJmZp7gCE\nEIVXRkYGu3bt4urVq9y5c4eiRYtSu3ZtXFxczB1avouJiSFo7ly2rV1LbGwsVpaWuLi58XavXrRp\n0wYrKytzhyjMSKOUUuYOQghRuNy4cYN5c+Ywa+pUHFNS8NZosFaKW1otu1NSaOXnx8BPP6V+/fpo\nNBpzh5unTp48ydcjR7JqzRraaLW0TU6mOJAOXACC9HrOWFjQd9Aghn3+Ofb29maOWJiDJF8hRJ5a\ntHAhg/v1o51Gw4DkZF6+Z/stYIFGw486HRVefpmlq1ej1+vNEWqei4iIoGPr1nyUlER/g4ESDyl3\nGBhja8sFb29CN29+Lu4EiJwk+Yp8oZQiIyMDS0tp2XiezJg+nW8//5y1SUlUfkTZdGCgjQ37y5Vj\nc1TUM5+Ad+/eTasmTViWlETTXJRXwEgrK9aXLcu2vXtxcHDI7xBFASIdrkSeOXPmDMM++giPYsWw\nsrTExtoava0tzevVY+XKlaSnp5s7RJGP1q5dy/jPPiMiF4kXjB1Ofrpzh1rnztG5dWue5XrAnTt3\neMvPj7m5TLwAGuDrtDRqXbjAkPffz9U+AwYMYNy4cU8cpyhAlBBP6cKFC8qvYUNVwtZWfWplpY6D\nSgVlABULaiGoenq98ihWTM2bM8fc4Yon9P7776uxY8c+cJvBYFA1vL1VKCj1kEcQqFce8HoaqIr2\n9mr79u2PFc+iRYtU8+bN8+LS7uPl5aXs7OyUvb29cnV1VT179lSJiYkPLb948WLVTK9/6LX/1+MG\nqKK2tiomJiZfrkUUTJJ8xVM5fPiwci9WTI23sFBJj/iQ2QOqvE6nRg4bpgwGg7lDFw9wN+k4ODio\nsmXLqsGDB6vo6OhH7rdz505Vzt5eZTxB8lWgvtNoVNc2bUxwhblTpkwZFR4erpRS6sCBA6p06dLq\ns88+e2j5BtWqqZBs12PIfOQ2Ab9jZ6e+/eYbU12eKADktrN4YpcuXcKvcWO+uXWLzzMysHtE+drA\nzqQkVs6YwdTAQFOEKB6TRqMhNDSUhIQEdu7cSXR0NIMHD37kfjMnT2ZAUtITt2P1UIp1GzZw/fr1\nXJXPyMh4wjM9vurVq+Pn58eRI0fYtWsX9evXp2jRotSoUYOtW7dy6tQpzpw6xTRgPNAcKAqcBcoA\n4dmONRp4J/P5B4A+8xGSnMywESP46quvWL58OXq9PuthY2NDkyZNAOjZsycBAQEmuW6RvyT5iifW\n+JVXqHDrFt0eY58SwLqkJMYFBHDx4sWHllu8eDEtWrR46hjFk3N1daVTp06cOXMGyPnBn5SURJ8+\nfShTpgzFixdnxerVvJHZZhsNvAW4YPx935u6xwDuQDlgfeZrTkAta2u6d+9OqVKlcHd3JyAgAIPB\nAMD8+fN55ZVXGDVqFF5eXowePZr58+fTsGHDfLt+g8FARkYGe/bsISwsDBcXF1q0aEGPHj04duwY\nAwYM4M0332Tnzp28YGWFBpgBfATcBDwxtutmH0iV/fkPQELmY0PmtrZt29KxY0cSEhJISEjgypUr\neHt707VrV+P+Gk2hG5r1vJKuqOI/bdiwgQkTJnD8+HEMBgPe3t706NGDdu3acfPqVfY+QScZT6C7\nwcCcmTMZO3HiA8t069aNbt0eJ62LvKIyf6eXL19m8eLFvPbaa0DOD/6goCCSk5M5dOgQ9vb2ONnb\nUwLIAN4AGmMc06oF9mY79m7gdeAIMA14D7icue1YfDyJO3dSoUIF0tLSCAwMJDg4mOLFixMTE8Ol\nS5c4ceKEMdmvWEFcXBz//PMPFSpUwGAwZD2UUo98nptyr7/+elbcWq2WBQsWoJTik08+YejQoWi1\nWm7fvk3fvn2pk5qKBdASaPWY/983gO6Zz6tVq5b1usFgoEuXLjRp0oS+ffve9/sRzzZJvuKh5syZ\nQ0BAADNnzqRFixZYWVmxf/9+JkyYwLUrV+io0VDkCY89IDWVJrNmETBmDNbW1jm2ZWRkYGFh8fQX\nIB6bUoq2bdtiaWlJQkICPj4+LFiw4L5yBoOBmJgYLl++zIsvvoi9jQ2pd+6wGzgKbIasZoj62faz\nB/6HMSn3x1gLvo5x2M014JP+/Wnfvj0ajYbly5cTGRnJTz/9xJo1a5g4cSIbNmzA2toarVbLqlWr\nCAkJ4ZdffkGr1aLVatFoNI98npty5cuXZ+7cuTRt+m/f5YEDBxIUFJTj71Wn09G9e3eO/PILmrg4\nfB/1/3vPz2lAB6A9EKTT5ajVjhw5ktu3bzN9+vRHHFU8iyT5igdKSEhgxIgRzJ07l3bt2mW97uvr\ny8qVK6lRrhzuKSkEAGOBCIzf3kcBE4EUjO1fPTP3i8N4+3E9oAP6Al4GA5GRkURHR/Pzzz/TtGlT\nFixYwLvvvou3tzdz584lMjLSRFcswFi7XbVqFU2bNiUhIYEFCxZQqVIlLly4APxb63rvvfdITEyk\ndeYQIa2lJWeAv4GKQPGHHL8y/7Z1lcr8NxGIwZiYfvrpJ+bOnQsYE7ynpye1atXi0KFD1KhRg1q1\namUda/fu3djZ2VGhQoW8/C94KE9PTxo3bkxYWFiO169fv84L8+dThfs/UO2Bq9l+3g85vrAOxtg+\nXB/484UXsl5funQpy5YtY8+ePfJFtJCSNl/xQEeOHCEhIYE2bdo8cPv1W7dwIGcb1nXgABCFsUYz\nCGPSBeOHTCxwDtgKBAEZaWncuHEDMH6Qpqenc+jQIUaOHJkPVyQel16vZ8CAASQkJLB58+Yc23Q6\nHZ9//jmnT59m9erV3Lp9m0AbG0oDJzG2eT6O2Mx/MzIyqFatGp9//jl//PEHhw4dyipj7glbunfv\nzs6dOwkODiY2NpaUlBQiIiJIS0uj+WuvcY37a7Y1gF8wfrlYifFv/66fgG3AImCmXs+A4cMB2L9/\nP4MHD+a3336jePGcX2PklnPhIclXPFB0dDQlSpRAq/33T6R+/fo4OTmh0+lISk1FQ84PGwPGpOuM\nscarBU5gbAdchjEB2wNeQGfgYlpa1sQblpaWjB49miJFimBra5vv1yce7u4HfGJiIrNnzyY1NRUv\nL68cH/yhoaGcPn0ag8GARqPBxtaWcKUoD1QBPgFOY7wDsiMX5/w/Ozte9PGhbdu2dO3alQsXLtCy\nZUtcXV15//332b9/v0l7OD+Ih4cHGzduZN68eXh7e+Pp6UlgYCBKKQZ++imXH1BDHYuxTbcisBTo\nl23bUoxfRksCWxIS6Nu3LxMmTGD16tX8888/vPLKK1k9nlu1MrYkS4erwkNuO4sHKl26NDExMRgM\nhqwEvGPHjqxtJCWRkpycY59SkDWXrWXm87u3FNOAl7KVfQn4TimcnJy4ceMG1atXv6/tV5hH69at\nsbCwQClFuXLlCA4OplKlSjk++E+fPs3gwYOJiYmhWrVqjB49mv1//ME3q1axJi2NDwFfjF/AumG8\nrXpvz18yfz4NhAB7QkMJDAxk3LhxxMfH4+3tzTvvvIPBYGDu3LmcOHGC5s2b4+/vj7+/f74monPn\nzj3w9Tp16hAREXHf6+7u7tSsWpXzR49CWlrW62WBXQ85xxbgNtDIzo7Oo0bx6WefZW378ssvH7hP\nUFBQ7i5AFHymHlgsng3x8fHKyclJhYSE3LfNw8NDtW7RQtXQaNT/MicJ2ALK456JA8qACgeVDsoa\n1Pps2z4GZanVqpiYGBUUFKReeeWVHOd40GuiYLt+/bryLlVK/aDV5npyiQugyuh0an5Q0COPHxcX\np3799VfVp08f5ebmpsqXL68+/PBDtX79epWcnJz/F/gI169fV+Xd3NQYS8tcTbARC8oXlN7SUu3b\nt8/c4QsTk9vO4oH0ej0TJ05k4MCBhISEkJCQgMFg4MCBA9y+fZs2b7/NSa0WQy6OZQF0BGZirAlf\nABYCvi+9dF+blnh2OTs7syEyku9cXfnUyopb/1HWgLHzXX07O4aMG0ePnj0feXxHR0fatWvHnDlz\nuHTpEitWrMDV1ZWxY8fi4uJC69at+fHHH7M6h5mas7MzkXv3sqp8eV63t2cNxiaXe/0DTNdoqGVv\nT8WOHbFwcKBx48b88ccfJo5YmJW5s78o2MLCwlSjRo2Us7OzcnZ2Vr6+vmrOnDkqNTVVuTg5qbey\n1XxLP6Tme/dbfndQzpk1ZCcrKxUVFaWUUmr+/PmqYcOGOc77oNfEs+FuDdDB0lL1srNT20FdAXUT\n1ClQUzQaVcHeXlUvV06tXr06T8558+ZN9csvv6h33nlHlShRQlWqVEkNGzZMbdmyRaWmpubJOXIr\nOTlZLViwQPlWqqQ8dTo12NpajQY1UqNR3XU6VdTWVnVu3Vpt27ZNGQwGdfDgQVW0aFHl6OioNmzY\nYNJYhfnIkoLiiUVGRtK+RQs2JydTJZf7GICudnbw+ussXbUqP8MTZhIZGUn37t2JjIxk6eLFLJkz\nh6s3b3InLY2i9vY0aNiQgZ9+Sr169fKlzTYjI4O9e/eybt061q5dy6lTp3jttdfw9/fHz8+PUqVK\nPfogeWTv3r1s3bqV2Fu3sLaxwcXFhTZt2lCyZMkc5Xbv3p01o9vcuXN56623TBajMA9JvuKp/LJk\nCZ/06cOvycnUe0TZZOBdOztuVq3Kuq1bpVdzIWQwGPD19eXjjz8uMDOUXbt2jfXr17Nu3To2btxI\n2bJlszpt+fr6FphxtNu2baNt27ZoNBqmTJlCjx49zB2SyEeSfMVTCwsL492OHakLDExMpAU5x7Bd\nBH7UaAiytaW5vz9zFi/GxsbGPMGKfLV48WKmTp1KVFRUjmFqBUV6ejo7d+5k3bp1rFu3jsuXL9Oi\nRQv8/f1p0aIFJUqUePRB8tGGDRvo2rUrVlZWfPHFF3z44YdmjUfkH0m+Ik8kJSWxdOlSZkycyPXL\nlyljZYUNcEspzqelkWEwsH3PHqpWrWruUEU+SU5OxsfHh0WLFuXrggd56dKlS1mJeMuWLVSqVCmr\nVlyzZk2zfIH47bff6N+/P7a2tvTt25f//e9/Mra3EJLkK/KUUopTp05x9epV7ty5Q9GiRalUqRJv\nvPEGH3zwAe3btzd3iCKfjB8/nr179xISEmLuUJ7InTt3iIyMzErGcXFx+Pn54e/vz+uvv06RIk86\nk/njW7RoEcOHD0ev19O6dWu+/fZbScCFjCRfYRLBwcEsW7aMtWvXmjsUkQ+uXbtG5cqV2bVrF+XL\nlzd3OHnizJkzhIWFsW7dOrZv306tWrWyasWVK1fO92Q4e/Zsxo0bR4kSJXjppZeYNWtWgWmfFk9P\nkq8widu3b1O6dGmOHDmCm5ubucMReez9999Hp9MxZcoUc4eSL5KSktiyZUtWD2qDwZCViJs1a4a9\nvX2+nPe7775j5syZuLq64u7uzsKFC2UmuEJCkq8wmb59++Lt7c1n2abRE8++o0eP0qRJE/766y+K\nFStm7nDynVKKv/76K+v29O7du6lfv35WMs7rVZbGjBnD8uXL8fT0BOD//u//0Ol0eXoOYXqSfIXJ\n7Ny5k549e/LXX39J+1Uh4ufnR4sWLfj444/NHYpZxMfHs2nTpqxkbG9vn5WIGzVq9NRD6pRSjBgx\ngs2bN+Pt7c3Vq1dZs2YNjo6OeXQFwhwk+QqTUUrx4osvMnfuXBo0aGDucEQe2LhxI4MGDeLo0aNy\nOxTj3/jBgwezEvGhQ4do3Lhx1gQfXl5eT3zcQYMGceTIEXx8fNi3bx/r16+/b2iUUoqoqCi2bdvG\nrRs3sLKxwbVkSdq1a4e7u3teXKLII5J8hUlNmjSJkydP8vPPP5s7FPGUMjIyqFGjBmPGjKFdu3bm\nDqdAunXrFhs3bmTdunWEhYXh6uqaVStu0KABVlZWuT6WwWCgV69eXL16lerVqxMaGsrvv/+Ou7s7\nSUlJLFm8mJmTJhH399+0vXOH4unppAPnbW1ZCTRt3JiBw4bRpEkTufNUAEjyFSb1999/U6lSJS5d\nupRvnVSEacyZM4eFCxeydetW+TDPhYyMDP7888+sWvGTTHuZnp5O586dSU9Pp06dOvz8888sWrSI\nAe++i9vVq3x0+zbNuX+h9gRgETDF3p7XO3Zk+uzZWFrKirLmJMlXmFzr1q3p0KGDTJ/3DEtISOCF\nF15g9erV1K5d29zhPJPunfayXLlyWbXiOnXqPHRYUWpqKm3btsXJyYmqVavy9ciRjNRoGJGRcd96\nyfeKAzrodJT08yN4xQr50mRGknyFyf36669Mnz79gYuSi2dDQEAA58+fZ+HCheYOpVBIS0tj165d\nWUOZrly58p/TXiYnJ+Pn58eZQ4d4Lz6e0RkPWrzwwZKAZvb2vPXFFwz74os8vhKRW5J8hcmlpqbi\n4eHBzp078fb2Nnc44jFFR0dTo0YNDhw4QOnSpc0dTqEUHR2dNcFH9mkvW7VqRY0aNdBqtaxdu5bB\nb77JGYPhkTXee50AGjk6cvHGjf/sKHfx4kUqV65MfHy81JLzmCRfYRZDhgzBwcGBsWPHmjsU8Zje\nffddSpcuzddff23uUAql7du3M3z4cI4ePYqlpSU+Pj5069aNVatWERkZSZEiRfDz8+Pk/v10OXiQ\nQU/4Ef6aXs97P/1Ely5d8vgKRG5I8hVmcfjwYfz9/Tl//rxMmfcM2bt3L2+88QYnT55Er9ebO5xC\nJz4+Hk9PT0aPHk3v3r2xtrYmMjKSkiVLsnfvXubOnUtQUBDLli1jwqhRXDEYyD7a18D9na0e5ldg\nWs2abN23L+8vRDySJF9hNrVr12b8+PE0b97c3KGIXFBK0aRJE7p27Uq/fv3MHU6hdPjwYerUqUN8\nfHyOYUh//fUXNWvWJC0tDTs7OzQaDTW0WsrFxeEEXAbCgRVAKWAAcBBwByYArYErwAvZzpWBcY1t\npRTVq1fn7NmzWdtu375NREQEnp6elCtXjvT09AK5ROSzTP43hdn06tWLefPmmTsMkUurVq0iJiaG\n3r17mzuUQqtq1aq4u7vTpUsXQkNDiY2NBcDHx4dZs2ZRr149EhISWLRoEY6Z9aY5QDvgGvAyxkT7\nOnAD+B7oBJwE3DAOObr7aAtZbcUHDx4kISGBhIQEAgMD8fHxoVatWqa67OeSJF9hNl26dGH9+vXc\nunXL3KGIR0hNTWX48OEEBgbK+NB8Fhoail6vZ8CAAbi4uNCmTRtiYmLIfpPS3t6epMwOUL5AF8AS\nOABcBD7O/LkJUAv45Z5zfAP8BTjY2OR4ffv27QQEBLB69WocHBzy5fqEkSRfYTbFihWjZcuW/PLL\nvR8NoqCZNWsW5cqVo0WLFuYOpdDz8fEhKCiI6Ohotm/fzunTp5k+fXqO3sZeXl4cT01FYUy+d10B\nKgLZp6+pjfG29F1hwHTgf0CZbBN7REdH06lTJ4KDgwvNspAFmSRfYVa9e/cmKCjI3GGI/xAbG8u4\nceOYPHmyuUN57vj6+tKlSxfWrl2LRqPJqv2WL1+e8hUrEg1k767ohvEW8+1sr+0BPDKfnwB6Ymwb\nDtHp6D14MGAcN9y2bVuGDBkiX7BMRJKvMKtmzZpx7do1Dh06ZO5QxEOMGzeOdu3aUaVKFXOHUuid\nOHGCKVOmcPnyZTIyMti3bx+LFi2idu3auLq6curUKRITEwEYOGIEx+9pAqgLlAamAWlABLAf6AzE\nA22Ar4HywDqDgR69egHGL8Evvvgin376qUmuU0hvZ1EABAQEkJiYyHfffWfuUMQ9zpw5g6+vL0eP\nHsXV1dXc4RR6V65cYciQIfzxxx/Exsbi5ORE69at+fbbb7GxsaFdu3Zs374da2trLl26hJOjI63u\n3GF5tmMcAwZibP/1wJhs22BMxE0x3pJOATSWltja2REfH49Wq0Wn0+W4tb1+/Xrc3d3x9vYmLS1N\nejvnMUm+wuzOnDlDvXr1uHTpkixLV8B06NCBmjVrMnLkSHOHIh7gzz//xK9RI5YkJfF6LsobgE+0\nWv544QUi9uyRxU3MSL7KCLPz9vamUqVKhIaGmjsUkc327dvZvXs3n3zyiblDEQ9Ru3Ztfl2/nu56\nPWMsLLj+H2X3Ae1sbFhoZcXbPXtK4jUzSb6iQOjVq5d0vCpADAYDQ4cOZfz48djZ2Zk7HPEfGjZs\nyB/79hHduTMv2NrS0dqaFcBmYAMwG6in1+Nnb88fdnYsX7uW7777jnXr1pk38Oec3HYWBcLt27fx\n8PDg2LFjuVrXVOSvX375hSlTphAVFSVtfc+Q2NhY/Fu2JP2ff7DX6bC2tsbFzY1OvXvj7+/PtGnT\n+OmnnwgMDKRXr15s2rSJ6tWrmzvs55IkX1Fg9OnTh4oVKzJ8+HBzh/JcS05OxsfHh4ULF/Lqq6+a\nOxzxmKpUqUJwcPBDZ6j66quvCAkJ4aOPPuKrr75i165duLm5mThKIclXFBg7duygd+/eHD9+XJYv\nM6OJEyeye/dufv31V3OHIh5TfHw8bm5uxMbG5pgbOjulFMOHD2fr1q20bNmS0NBQtm3bJjNamZjc\nTxIFRr169QDYtWuXmSN5fl2/fp3JkyfzzTffmDsU8QT27NlDzZo1H5p4ATQaDZMmTeKll14iIiKC\nKlWq0LVrVzIyMkwYqZDkKwoMjUZDz549ZbEFM/ryyy/p3r07FSpUMHco4gns2rWLunXrPrKcRqNh\nxowZlClThmvXrhEfHy8TbJiYJF9RoLz77ruEhIRw+/btRxcWeerYsWP83//9H6NGjTJ3KOIJ5Tb5\nAmi1WubNm4der8fBwYGwsDBmzJiRzxGKuyT5igLFzc2NevXqERISYu5QnjvDhg3jiy++oFixYuYO\nRTwBpRS7du3C19f30YUzWVpasmTJEgwGAy+88ALjxo2TIUgmIslXFDiy2ILp/f7775w4cYJBgwaZ\nOxTxhM6dO4eNjQ0eHh6PLpyNtbU1ISEhxMfHU6dOHd59910OHDiQT1GKuyT5igKndevWHDlyhLNn\nz5o7lOdCRkYGQ4cOZdKkSTK95zPscW4538vOzo7Vq1dz7do16tatS+vWrbl8+fKjdxRPTJKvKHCs\nra3p1q0b8+fPN3coz4X58+dTtGhR2rVrZ+5QxFN4muQLoNfrCQsL49KlS1SsWJHWrVtnraAk8p4k\nX1Eg9erVi/nz58vwh3yWmJhIQEAAgYGBMrb6Gfe47b0P4uTkxMaNG7l8+TKWlpZ06dJF3oP5RJKv\nKJCqV6+Os7MzmzdvNncohdqkSZNo2rQpL7/8srlDEU8hJSWFo0eP8tJLLz31sVxcXAgPDycmJoZT\np07Jwhr5RJKvKLBksYX8denSJWbMmMH48ePNHYp4Svv378fHxwedTpcnx3N3dyc8PJzExERWrFjB\n999/nyfHFf+S5CsKrK5du7Ju3TpiY2PNHUqh9L///Y/+/fvj6elp7lDEU3ra9t4HKVu2LJs2bSIj\nI4NRo0bJkp95TOZ2FgVap06dePXVV6lSpQrnz5/n9u3bODo6Uq1aNapVq2bu8J5Z+/bto1WrVpw4\ncQJHR0dzhyOeUqdOnWjVqhXvvvtunh/74MGDNGnShIyMDCIiIqhZs+YDy6WmpmJpaSmrYOWSJF9R\nYN28eZPPhw9nRXAwXnZ2VFMKXUYG8ZaWbM/IwKNcOQaOGEGHDh2wtbU1d7jPDKUUTZs2pXPnzvTv\n39/c4Yg8UKZMGTZu3EjFihXz5fi7du2iRYsWWFtbs3//fjw8PEhNTeW3335j5jffsPvwYdINBhTg\n6uhI95496T94MOXKlcuXeAoDSb6iQNqwYQPdO3TA32BgYFISdYDsfXHTgXXADAcHzjo6snbLlnz7\n4ClsVq9ezeeff87BgwextLQ0dzjiKf39999UqVKFmJiYfO2xHhERwRtvvEGpUqUY0Lcv344fz4sG\nAwMTEmgJ2AMZwGlgjrU1C7Ra6tWty6xFi3B3d8+3uJ5VknxFgbNy5Ure79qVX5OTqZ+L8j9rNAQ4\nOrJl1y58fHzyPb5nWVpaGlWqVGHatGm0bNnS3OGIPLBy5Upmz55tkmkh165dS8c2bXBTit8MBqr8\nR9kkYLKlJXOLFmVDZKS8N+8hN+dFgXLgwAH6duvG2lwmXoA+SjEhPh7/xo2Ji4vL1/iedbNmzaJM\nmTKSeAuRvBjfm1t7duygoqUlfz4i8QLogFHp6Xx58yZ+jRpx9epVU4T4zJDkK8yuZ8+eBAQEAPDJ\nwIGkJyXxsNGKA4BxDzqGUtRJSGD+I4Ym+fv7s3DhwqeK91kVGxvLuHHjmDx5srlDEXkoKioqz3s6\nP8hff/35pQVfAAAgAElEQVTFj999x/o7dyiSy31GA1uUouOtWwzu1w+tVovBYHjqWAYMGMC4cQ/6\nJHh2yG1nYVKNGzfm0KFDXL16NWse4V69elG6dGkGDhxIBU9PiqSlcekJjr0d6OPuzvHoaJmt6QE+\n/fRT4uLimDNnjrlDEXkkPT0dJycnLl68iJOT033bw8LCmDhxIsePH0er1VKrVi1GjhxJgwYNHvtc\nHw0YwJaffqKNUozN5T5fYWwDngKUt7EhITWV9PR06RGN1HyFCZ0/f57du3fj4uLC6tWrc2xTSjF3\n9mwaaTQ8adpsAFjHxbFly5anjrWwOXv2LPPnz2fs2Nx+bIpnwdGjR3F3d39g4p0zZw49e/ZkyJAh\nXLhwgVOnTtG7d2+WLVv22Oe5ffs2i4KDeeEx62p3SzsDTeULcQ6SfIXJBAcH89prr/HOO++wYMGC\n+7bv3ryZOqmpAPwAlM18LMlWpicQkPm8NaDP9rAESt++zbfffoter896WFlZ0bt3b8BY8547d26+\nXF9B9tlnn/Hxxx9TsmRJc4ci8tDDJtdISEhgxIgRzJw5k7Zt22JnZ4der6dDhw5Mnz6d3bt3U69e\nPYoWLUqpUqUYPHgwaWlpWftPnDiRGjVqUKRIEapVq0ZQUBAuBgOrgEkY329tMsteAdpjTLBlgYfN\nhdUhJYXsN1qDgoKoVKkSer0eb29vZs+enbUtIiICDw8PZs+eTbly5XBzc8ux0Er2pqpnlSRfYTLB\nwcF06tSJjh07smHDBm7cuJFj+z///IMjcA04A5wEFgDvAacyy2j4d8jRGiAh87EcKAVUVYpqlSuT\nkJBAQkICx48fx8XFhU6dOhn312ieu1vSO3bsYOfOnTJHbyH0sPbeI0eOkJCQ8NCVqiwtLZk2bRo3\nb94kJCSENWvWMHPmTMBYm54/fz7r1q0jLi6OFStWkJqaSj0LC7oBIzC+51YBBoxfgj2AI0AwMBHY\n+IBzlrrnZ1dXV9auXcutW7f4+uuvGThwIPv378/afv36dQ4cOEBUVBRjxoxh0KBBWR0qC8P7WJKv\nMInt27dz+fJl3nzzTSpUqEClSpVYvHhxjjK2trakYhzDOxiwAl4FagLZb5Tde+PrJMYa8XLABrC1\ntwcgOTmZNm3a8PHHH9OiRYt8uKqCTynFJ598wvjx4/Ns3l9RcDys5hsdHU2JEiUe2rZaq1Yt6tSp\ng4WFBfXr16d79+5s3boVMK7vnJKSwqlTpzAYDLzwwgvY29tjkdlRKvv7bw9wFmPCdQUaAm8DSx9w\nzntHlPv7+1O2bFmsrKzo3Lkz9erVIzIyMmu7wWBgzJgxODs707NnT7RaLSdOnMja/qx3V5LkK0xi\nwYIFNG/eHL1eD8Dbb7+ddev57pvIzcuLaMAByD4vTi3g74ccNw7j7a+vgfrAKZ0ONzc3AN577z1e\nfPFFhg0blufX86xYtmwZ6enpdOvWzdyhiDz2zz//cPHiRapUuX/QT+nSpYmJiXloz+KTJ09mTZhR\npEgRvvvuO27evAlAtWrV+Prrr/nss89wdXWlf//+6HQ6Yh4wIcsFIBFwA5wyH/OA6w845617fg4L\nC6Nu3boUL14cJycndu/enRUDQKlSpShRogRgrKmXKFGiUK0vLMlX5Lvk5GSWL1/O5s2bKVWqFKVK\nlSIwMJBDhw5x6NChrNtH3fv1Y62dHYkYbzvftRfjm/teBqAr0AzoA8QAYQYDHTp0YOLEiZw+ffq5\nbN+9KyUlhc8++4zAwEDpXVoI7d69m5deeumBs5RVqVIFR0dHfvvttwfuO2DAAEqWLMnp06eJi4tj\nyJAhORJ1t27d2LlzJ8eOHePQoUNER0cTkZpKGjlrvqUxtv9eA2IzH/HA3SUYst8YDssW5507d2jf\nvj19+vTh+vXrxMbGUqdOnWe+Nvs4ZG45ke9WrlyJpaUlBw8ezBpepJSiY8eOBAcHZ5Vr2rQp6fb2\naJOTmQF8A+wEDmBs+4Wcb/yRGGfRmZr5889AGS8vNm3axPfff09UVBQ2Njb3xfO8vMGnTZtGzZo1\nadSokblDEfngv8b36vV6vvnmGz744AMsLS15/fXXSU9PJzw8nIiICBwcHHByckKr1bJlyxaCg4Mp\nU6YMAH/++ScZGRnUqlULpRQWFhakpaXxerNmXF+3jniMTUOWgC/GTlZfAB9gTMbHgRSgNv++X5OB\nFRYWaDIysmJ0cHCgePHipKSksHz5cqKiomjSpEmurr0wvIfl67DId8HBwfTu3RsPDw9cXFxwcXHB\n1dWVDz74gMWLF5ORkZHVgaJt585YajR4ARWBdzEm1buzNmfvcLUUiMJ4q8sB4weAlaMjPXr04Nq1\na7z44otZPZ4HDhyYFc+z3lEjN27cuMG3337LN998Y+5QRD551DKCffr0ISgoiClTpuDl5UXFihWZ\nN28eXbp0YfTo0Rw4cAAPDw++/fZbPvjgg6z94uPj6devH8WKFaNJkya8/PLLDBs2jIHDhnFGp+Mi\nUBx4C2MCCcXY47kuxh7P/TDWfuHf9+tioFrVqlnvPRsbG6ZNm8aYMWMoX748f/75J507d84R/3+9\nTwtDhyuZZEMUKOnp6bRr0YKiO3YwLyUFq1zscxNoZGFBqqcne/bv58aNG4waNYotW7bwxRdf0K9f\nvwfWgAuzQYMGZfVoFYWPUooSJUpw+PDhrD4Opjhn2+bNcdq+nXkpKbmuue0B/HU61m/bxksvPWzu\nuueP1HxFgWJpacnS1auJqVOHN3Q6jvxHWQVsBurrdPh/8AHN/f1p0KAB1tbWLFmyhLCwMNavX4+P\njw/BwcFkZLvlVZgdP36c5cuXM2rUKHOHIvLJ6dOnsbe3N1niBWNtc8nKlZz28eEdW1vu5GKfcKCV\nnR3zli6VxHsvJUQBlJaWpsZ++aUqVbSoelWvV0tAHQZ1FtR+UNM0GvWCg4Oq7OWlFi1cqJRSymAw\nqClTpih3d3e1d+/erGNt27ZNNWjQQFWuXFn99ttvymAwmOuyTKJVq1YqMDDQ3GGIfLRw4UL19ttv\nm+XcSUlJqvObb6qSdnYqwNJSRYNS2R5poH4F9ZqDgypZtKiKiIgwS5wFnSRfUaClpqaq5cuXq5YN\nGqhKHh7Kq3hxVdXLS3V58021bdu2BybSkJAQ5ezsrEJDQ7NeMxgMKjQ0VFWrVk35+vqqzZs3m/Iy\nTGbTpk2qXLlyKiUlxdyhiHw0aNAgs3/BOnr0qPqgb19VxMZGlbewUI2LFFF1ixRRrnZ2qkG1amrx\n4sXyd/gfpM1XFEq7du2iXbt2jBo1igEDBmS9bjAYWLp0KaNGjcLb25vx48cXmtthGRkZvPTSSwQE\nBNC+fXtzhyPyUe3atZk2bdoTLZCQ17799luioqIYMGAAtra2uLm5UbZsWXOHVeBJm68olOrWrcv2\n7duZOnUqw4cPzxrDqNVq6dq1K8eOHaNt27a8+eabvP322zlmznlWLViwAL1ez1tvvWXuUEQ+Sk5O\n5tixY9SqVcvcoQDGCTuaNm1Ks2bNaNCggSTeXJLkKwotb2/vrHmNO3fuTEpKStY2a2trBgwYwKlT\np6hduzavvPIKffr0ITo62owRP7nExEQCAgIIDAx85odgiP+2b98+KleujJ2dnblDAeDw4cNUrVrV\n3GE8cyT5ikKtePHi/P7771hYWNCsWTNiYmJybNfpdIwYMYKTJ0/i4uJCjRo1GDp06H3lCrrJkyfT\nqFEj6tSpY+5QRD571PheUzIYDBw9evSBU1yK/ybJVxR6tra2LF68mFdffZV69epx6tSp+8o4OTkx\nfvx4jhw5QkpKCj4+PowZM4aEhAQzRPx4Ll++zPfff8+ECRPMHYowgV27duHr62vuMADjGt1FixZ9\n4HrC4r9J8hXPBa1Wy4QJExg2bBgNGzZkx44dDyxXqlQpZsyYQVRUFCdPnqRChQpMmzaNO3dyM6rR\nPP73v//Rr18/vLy8zB2KMIGCVPOVW85PTpKvCSml2L59Ox/2708nf3/eeu01enfuTHBwcI72SJF/\n+vXrx/z582nbti0rVqx4aDlvb28WLVrExo0b2bRpExUrViQoKIj09HQTRvto+/fvJywsjM8//9zc\noQgTuHz5MsnJyXh7e5s7FMCYfKtVq2buMJ5NZh7q9FxITU1Vs378UVUtU0ZVtLdXEzQa9QuoFaBm\ngfJzcFDODg5q2EcfqcuXL5s73OfC/v37lYeHh5o0aVKuJt3Yvn27atiwoXrxxRdVSEhIgZiow2Aw\nqCZNmqgff/zR3KEIEwkJCVH+/v7mDiNLx44d1aJFi8wdxjNJkm8+i4+PVy0aNlSNdTq1CZThntlg\n7j5Og/rYykq5OznlmJ1J5J/o6GhVrVo1NWDAAJWWlvbI8gaDQa1bt07VqFFD1alTR23atMkEUT7c\n6tWrVaVKlXIVuygchg0bpsaMGWPuMLL4+PiogwcPmjuMZ5JMspGP7ty5Q8uGDSl36BA/3bmTq/Ub\nQ4CBej1bdu2iUqVK+R3icy8+Pp4OHTpgbW3N0qVLcXBweOQ+BoOB5cuXExAQgJeXFxMmTODll182\nQbT/SktLo2rVqkyZMgV/f3+TnluYz6uvvkpAQACvv/66uUMhJSUFJycn4uLispYKFbknbb756ItP\nPqHokSPMzmXiBWgPTE5M5M3XXsvRvhgZGYmPj0++xPk8c3R0ZO3atZQsWZJXX32VK1euPHIfrVZL\n586dOXbsGG+//Tbt2rWjffv2HD9+3AQRG82ePRtPT0/8/PxMdk5hXunp6ezbt6/ADCc7fvw43t7e\nknifkCTfPFCmTBl0Oh0ODg6ULFmSXr16cfXqVYKCgvg+ORmLxzzeO0pRKjGRNWvWZL3WsGFD/vrr\nr7wNXABgZWXFnDlzaN++PfXq1ePIkf9aSynnfv379+fkyZP4+vrSqFEjevfuzcWLF/M13n/++Ycx\nY8YwefJkmVDjOXL48GE8PT0pUqSIuUMBpKfz05Lkmwc0Gg2hoaEkJiayYcMGwsPD6dmjB020Wjye\n8JgDExKYKQuhm4xGo2HkyJGMHz+epk2bEh4enut9dTodw4cP5+TJk5QqVYqaNWsyZMgQbty4kS+x\njh8/ntatW0sv0+dMQRpiBJJ8n5Yk3zxWvXp1/P392bV9O5dv38YFcAJaA5czy+wE9NketsDd2VCL\nZr7WDwiPikKr1XLx4kUiIiIoXbq0Sa/ledStWzdWrFhB165dmT9//mPtW7RoUb7++muOHj1Keno6\nPj4+jB49mvj4+Mc6zpUrV5g4fjx9u3Wja+vW9H/3XaZNncqtW7c4d+4c8+bNY+zYsY91TPHsk+Rb\nuEjyzSN3+63t37+fdevWkXznDoOAi8CfQBrwQWbZekBC5iMWqAt0zdz2T7ZtpW1sqFatGu7u7ia7\nDgGNGjVi69atjBkzhi+//JLH7ZNYsmRJvv/+e/7880/Onj1LhQoV+O677x45lnv79u287e9PFW9v\nzo8ZQ+0lS/APDaXGwoXsGTkSb3d3WjZqRKdOnShVqtTTXKJ4BkVFRUnyLUzM29m6cPDy8lIODg6q\naNGiysvLSw0aNEhpNRp1J9tQot9BOT1giNH7oFo/4PWloOw0GhUcHKyUUmrLli3Kw8PDzFf6fLl6\n9aqqU6eOeuedd9SdO3ee+DiHDh1Sb775pipdurSaO3fufUODDAaD+nbCBOWu06kZGo2Ke8hwtOug\nvgblotOpxTK28rly8+ZN5eDgoNLT080dilLKGI9ery8Q492fVVLzzQMajYZVq1YRGxvL+fPn+eGH\nH3C0taUvUAYogrEXcxyQvQ71E7ANWHLP8fYDg4FK9vZS6zUjV1dXtmzZQnx8PC1btuSff/55ouNU\nrVqVVatWsWzZMoKDg6latSohISFZNerJEycSNHYsO5OSGKgUjtn2jQDuNjY4A18A4UlJjOjXj6W/\n/PLkF3ePixcvotfrH7uWL0xj9+7d1K5dGwuLx+2+mT8OHz5MlSpVpMPfU5Dkm0/0ej3bgCiMSTcE\nY+K9+9EWCYwCVgHZR5ZeB9oB3wHn09OpUKGCyWIW99PpdISEhFCtWjUaNGjA+fPnAdiwYQONGzfG\n1dUVZ2dn6taty48//vifx6pXrx5btmxh6tSpfP3119SpU4fJkyczfdw4NiQlkdsW/SrAuqQkPuzT\nJ8+GN3l6epKQkCAfpgWUtPcWPpJ880mDpk2J02qxAY4B2fstRwMdgYVA+WyvpwMdgO7AHaBe3brS\nyaoAsLCwYOrUqfTr148GDRowcuRIevTowYcffsjZs2e5fPky06ZNY8OGDaSmpv7nsTQaDS1atODP\nP/9k2LBhTBo1ijFJSY/dK74qMCA1lRmBgU98XeLZURDbe6W3/VMy933vwqBMmTIqPDw8x2t///23\nsrawUMVA1QK1BJQWVAaooMznDtkeVUCdB6UBZZ+53c7OTun1ehUdHa22bNmiSpcubaYrFHctWbJE\naTQa9fnnnz+0TKNGjdTPP/+c9XNQUJB65ZVXsn7+8MMPVenSpZWDg4Oy0GjU79nadZNA9QBVHJQv\nqJmgPLJtnwDKO/NvpiIoe2trFR8fr5RSKiMjQ33yySfK2dlZ1ahRQ82fP19pNBqVkZGRFdeECRNU\n8+bNVdGiRVXz5s1VTEyMUkqpc+fO5SgrCo6MjAzl5OSk/v77b3OHkqVevXpq69at5g7jmSY13zxw\n7tw5mjZtmuO1kiVL8vnw4dTQ6dgJdAEyMN5q6Jn5PCHb4zDgBRiA8RYWVC5XjsTEROLj4/Hw8KBx\n48b5PnmDeLQyZcqg1WoJCgpixowZDyyj0Wj+8/ZtnTp1OHjwIO/36UNNjYZ3gLv15a+A05mP/wO+\nB7IfqTywHWPzxPtASmoqc+bMAYyzXoWGhnLgwAG2bNnC/Pnz74tj5syZDB8+nCNHjhAbG8vkyZMf\n/z9BmNSpU6dwdHSkZMmS5g4FME6veuTIEbnt/JQk+eajgLFj0TdoQGc7O5Jzuc8crZZvHB1ZHR6O\nViu/noImOjoaZ2dnduzYwffff8/QoUOpX78+Tk5O6HQ6IiMjH3mMbt264eTkxMUTJxhgMBAHnMjc\nthzojXG8twfGZojsXaA6ACUBO+BjjP0F/ti2zbjv8uV06tQJNzc3ihYtyjvvvJOjA5VGo6Ft27Y0\na9YMd3d33nrrLQ4cOPC0/yUinxW09t4LFy7g6OiIk5OTuUN5psmnez6ysLBg2Zo12DRvTn17e37h\n3xpOdgrYDXSzs+MbV1cioqIoU6aMSWMVuVO6dGliYmLw8vJix44d7N27Fzc3N65cuULx4sUxGAyP\nPMbkyZOpVKkSv27cyIcY2/djMrf9DdTIVrbmPfsGZ253AoqROVY8xrj333//Tc2a/+5Rq1at+85d\no8a/Ry9ZsiSJiYmPjFeYV1RUFL6+vuYOI4t0tsobknzzmY2NDYt//ZX/LVjAzy+/jKedHUMtLZkO\n/AiMA152cKCZRoP7wIHsOXpUejgXYFWqVEGv17Ny5UqKFSvGhg0bsLGxoWnTplmJ197enqtXr2bt\ns3///qznkZGRfPnllyxatIi3mjdnOsZa7N36aSmMQ83u2pft+QWMteIJwK3MhyNgY29v3LdUqRzn\n2rcv+97iWVXQar6SfPNGbhfbEU9Bq9XSvn37rJVv/m/5ck5dukRqSgpOJUsypkkTIiMjuX37ttzK\nKeD0ej0TJ05k4MCBKKVo3rw5wcHBvP/+++zevZtLly5Ro0YNVq5cSY8ePThx4gSrVq3C09MTMA5d\nsra2pkiRIriVL890rZY72WrLHYH5GMeFJwK/ZTu3TeajOMbhazMz/3XNnO2qY8eOBAYG4ufnh16v\nZ968efe1+SoZx1tgKaXYtWsXP333HccOHiQuIQF7OzvOnD/P9evXMRgMBaIp6vDhw7Rq1crcYTz7\nzNvfS9x18eJF5eTkpOLi4swdisiFsLAw1ahRI+Xs7KycnZ2Vr6+v6tGjh3JxcVFr1qxRzZs3V0WK\nFFEtWrRQkyZNUg0bNsza94svvlCenp6qbNmyyt7SUpUBFZ6tt/O7oIpl9nb+EVTpbL2dp2X2cvYE\nNQSUpVarZsyYoZRSKj09XQ0ZMkQVL15cVatWTc2cOVNZWVllnbdx48Zq7ty5WT/Pnz8/K65z584p\nrVYrvZ3N5JclS1QNb29V3t5eTdZo1C5Qx0DtATUHVC0HB1XO1VVNnTLF7L+jSpUqqf3795s1hsJA\nkm8B0qFDBzV9+nRzhyGewoYNG5Szs7NaunRprsq/0bixmveQ6SQf9RhlaakGvffeQ4/9ww8/qDJl\nyuTVpYl8YDAY1IghQ1QFnU6tzxyK+KDftQHULlD1dTrVsXVrlZKSYpZ4U1JSlK2trdnOX5iY/x6G\nyDJ48GB++OGHXHXaEQVT8+bN2bRpE8OGDeObb7555G3eT0aNYpROx6XHPM8hYEpGBq82b571WkpK\nCuvWrSM9PZ2IiAhmz55Nu3btHv8ihMmMCQjg959+YkdSEi14eCccDeCLcWrRlE2b6N2li1k+J/76\n6y/Kli2LjY2Nyc9d6Jg7+4t/GQwGVb16dbV+/XpzhyKe0qVLl1T16tVVv3797ltI4V5jR49WZbVa\ndTGXNd5DoDx0OvXxRx8pZ2dnFRgYqAwGg0pKSlIvv/yycnBwUC4uLqp3794qISHBRFcscuP9999X\nY8eOVUoptWPHDlVap1NXs/1uNaDOPOL3nwSqrLW1qlChQtZxo6KiVMuWLZWjo6NauXKl8vPzUwsW\nLMjz+BcuXKg6duyY58d9HknyLWB+/vln1apVK3OHIfJAfHy8atmypWrZsmXWLFT3Sk5OVk2bNlX1\nXn5Zuel06nuNRv3zkA/da6DGa7XKOduqRmfPnlW+vr6qVatW6saNG6a8PPEQXl5eys7OTtnb2ysX\nFxfVvXv3B/bl6NKmjZqq0eT4Hecm+arM9n57W9usVYWaNm2qRo8ene/XNnz4cDVu3Lh8P8/zQG47\nFzBdu3YlKiqK06dPmzsU8ZT0ej1r1qzB09OThg0bcvny5Rzb09PT6dKlC8WLFydy506Wb9zINj8/\nytjY0N/WllnAIoy9mrvpdFS0teV0p05s2rmTrt26AVC2bFkiIyOpVKkSNWvWZOvWrSa/TpGTRqMh\nNDSUxMREwsPD2b59O+PGjctR5tq1a4Rt2ECPJ+x9XhVIT01lz549gHFVqvr16z9t6I8kw4zyjiTf\nAsbOzo733nvvoVMXimeLpaUls2bNokuXLtSrV49Dhw4Bxin6+vbtS1JSEosWLcLCwoIGDRqwfO1a\njp49S7lRo9jfvTthrVtzuEcPfMeP59yVK8xdsuS+Ce2trKyYNGkSc+bMoXPnznz11VdkZGSY43LF\nPapUqUKLFi2yvkz37NmTgIAAli1bRhuNhjmAG1AJWHHPvnHAu4ALxqVJv+bf8eAaoKTBQPDs2Xh7\ne3P27FnatGmDo6MjqampNG7cmLlz5+b59UjyzTsyzrcAGjBgALVq1WLs2LE4ODg8egdRoGk0GkaM\nGEGZMmV47bXXWLhwIRs2bODEiRP8/vvvWFtb5yjv5ubGiM8/f+zztGzZkr1799K9e3ciIiJYvHgx\nbm5ueXUZ4jGozBrtgQMHWL9+PUOHDgX+nff78sWLaJKTmQRsxZhce99zjMFALHAO4wxozTBOwnK3\nnB1w+dw5zpw5Q9myZZk7d27WHPOPml/8ScTGxhIXF4eXl1eeHvd5JTXfAsjLy4tGjRoRHBxs7lBE\nHurUqRO//vor7du3Z8WKFaxduxb7zNmp8oqbmxu///47TZs2pVatWqxbty5Pjy8eTSlF27Zt0ev1\n1KpVixYtWjB48OAcZVISEzkG+GGs9eqAftm2ZwDLMCZge4yLrnTGuAzpXRZAclJS/l3IPQ4fPkzl\nypULxEQfhYH8LxZQH374IT/88IPMSFTIHDx4kOLFi2NhYUFgYGC+/H4tLCwICAhgxYoVvP/++3z6\n6aePXGdY5B2NRsOqVauIi4tj5cqVLFy48L4FLIo4O3OLh8/jHQOkAS9le+0lIHuvgXSgaLFieRv8\nf5BbznlLkm8B1ahRIywsLAgPDzd3KCKPLF68mIkTJxIREcGePXsIDw+ne/fu3LlzJ1/O17BhQ/bv\n38/Jkyd55ZVXOHv2bL6cRzyYVqvlzTff5MMPP8y67ayU4sqVK5w8dYqbPHwe7xKAFfBnttf+xLjS\n1V2xWi21GjbMp+jvJ8k3b0nyLaA0Gg2DBw9m+vTp5g5F5IHQ0FCGDh3K+vXrKVu2LM7OzmzevJmU\nlBSaN2/OrVu38uW8xYsXZ9WqVXTr1o26deuyfPnyfDmPeLgBAwbwxx9/0KZNG5YvX87KlSvx8PAA\nnY61wHEgCfg52z4WGOf5nolxju8LGJeb7J65PRG4qRS9+/R56Hnz+q6KJN+8Jcm3AOvWrRs7duyQ\nGsszbuvWrfTu3ZvVq1dTuXLlrNft7OxYsWIFL7/8Mg0aNODcuXP5cn6NRsNHH31EWFgYX3zxBf37\n9yfJhG2Fz6Pr168za9YsWrVqRdWqVSlRogTnzp3Dz8+PgQMHMnnyZIaOGEF5S0uaYryl3B5jL+a7\nvse4alU5oCHQA+iVuS1So6F48eKUKFHioTHkZYcrpRRHjhyR5JuXzDfEWOTGp59+qoYOHWruMMQT\n2rt3r3J2dlbh4eH/We6HH35QpUqVUlFRUfkaT1xcnOrSpYuqXLmyOnr0aL6e63mSkZGhoqKiVEBA\ngKpRo4YqXry46t69u1q2bJmKjY194D7Xr19Xbk5OKuQx5/TeA6qETqcOHjxosus7f/68KlWqlMnO\n9zyQ5FvAnT17VhUrVkwlJiaaOxTxmI4fP65Kliypfvvtt1yVX716tXJ2ds51+SdlMBjU3LlzVYkS\nJdTPP/+cNUuSeDyJiYnqt99+U++9954qWbKkevHFF9Xw4cNVZGTkI6cUvWvv3r3K2cFBLctl4v0D\nlMqigoAAACAASURBVKudXb7/jdzr7kpdIu9I8n0GtGnTRs2aNcvcYYjHcOHCBeXp6anmz5//WPvt\n2bNHubm5qalTp+ZTZP86evSoqlKliurSpYssZZlLFy5cUDNmzFB+fn5Kr9erZs2aqalTp6rTp08/\n8TG3bt2q7LVa1cbWVm3KXMHo3qT7J6jetraqhL29CgsLy8Mryp3x48fLHbg8plFKxrIUdOHh4Xz0\n0UccPnw4zwfOi7x3/fp1GjZsyMCBA/noo48ee//z58/j7+9P8+bNCQwMxMLCIh+iNEpOTubjjz8m\nPDycpUuXUrt27Xw717PIYDCwZ88e1qxZQ2hoKJcvX8bPz4/WrVv/f3v3HR1ltTVw+DfpmXRIAUNI\ngjSpgaCItFAEKREUMPSOgIAFERVFUREuiAJXL8K9ggkYkIiCNKkRBGmiCCoo1QAJoCAJ6XV/fyTM\nl1ADJDMp+1lr1pp5yzl7JpPZc86c9xw6dOiAm5vbPdfRr18/XFxcCGrYkHkzZ5Jx8SLNc3Jwzcoi\nycaGg9bW/GVvz6jnn2foiBF4e3sXwTO7M3379uWxxx5j4MCBZq+7rNLkWwqICPXq1ePDDz80zWCj\nSqaEhATatGlDaGgob7311l2XEx8fz5NPPombmxuRkZEYjcYijPJ6UVFRjB07lkmTJvHcc8+V6y95\nSUlJbN68mTVr1rBu3Tq8vLzo2rUroaGhPPzww0X6ZWjFihW89tprHDhwAKPRiIiwe/duDh8+zJUr\nV3ByciIwMJB27doV65ew26lfvz6LFy+mUaNGtz9YFYom31Ji/vz5bNy4kZUrV1o6FHUTKSkpPPbY\nYwQFBTF37tx7TmAZGRkMHz6c33//nTVr1uDj41NEkd7YyZMn6d27Nz4+PoSHh1OxYsVira8kiYmJ\nYe3ataxZs4Zdu3bx8MMP07VrV7p27Uq1atWKpc7z588TFBTE119/TdOmTYuljqKQkZGBm5sbly9f\nxsHBwdLhlBmafEuJpKQk/P39+fHHHwkICLB0OOoamZmZdO/enQoVKhAREVFkU/CJCG+99RaLFy9m\n/fr11K5du0jKvZmMjAxee+01li9fTmRkJC3NOImDOWVnZ5u6k9esWcO5c+fo3LmzqTvZ1dW1WOsX\nER5//HEaNmx43YpHJc2hQ4cICwvjyJEjlg6lTNHkW4q8+OKLWFtbM3PmTEuHovLJycmhf//+JCUl\n8eWXX2Jra1vkdYSHh/Pyyy8TFRVF69ati7z8a61fv56hQ4cyZswYJk2aZNEuz6KSmJho6k5ev349\n3t7epu7kpk2bmvU5Llq0iA8//JC9e/det7BGSRMZGcnXX3+tE7QUMU2+pcjJkydp2rQpMTExxf4b\noCocEWHMmDEcPnyYb775BkdHx2Kra+vWrfTp04c5c+bQt2/fYqvnqtjYWPr374/BYOCzzz4rlSsk\n/fnnn6bu5N27d9OsWTNTd3JgYKDFYnrwwQeJjo4uFZNWvPLKKzg5OTF58mRLh1Km6AxXpUi1atV4\n+OGHWbp0qaVDUXkmT57Mvn37WL16dbEmXoB27doRHR3Nq6++yrvvvlvsi274+vqyZcsWQkJCCA4O\nZsOGDcVaX1HIzs5m165dTJo0ifr16/PQQw/x448/MnLkSGJjY9m4cSPjxo2zWOLNyclhyJAhvPTS\nS6Ui8YJOK1lszH1tk7o3mzZtkvr16+vECCXArFmzpHbt2vLXX3+Ztd7Y2Fhp3LixDBs2TDIyMsxS\n57Zt26RKlSry0ksvma3OwkpISJAvvvhCBg0aJF5eXlK/fn2ZNGmS7N69W7KysiwdXgGzZ8+W5s2b\nl7i4bsXPz++ermNWN6bdzqWMiFCnTh3mz59vlt/+1I0tWrSIt99+mx07duDn52f2+pOSkujduzcZ\nGRmsWLGi2AcIAVy8eJHBgwdz8eJFli1bZrHWI8CpU6dM197u2bOHRx55hNDQULp27VpiF3v//fff\nadGiBXv27KF69eqWDqdQ4uPjqVKlCleuXNF1fIuahZO/ugsfffSR9OjRw9JhlFsrVqyQypUryx9/\n/GHRODIzM2X06NFSv359OXPmjFnqzMnJkQ8++EC8vLzkiy++MEudIiJZWVmyc+dOefnll6Vu3bri\n7e0tQ4YMka+++kquXLlitjjuVmZmpjz44IMyb948S4dyR3bs2CEPPfSQpcMok7TlWwolJiYSEBDA\ngQMHqFq1qqXDKVc2b95Mv3792LRpE0FBQbc/oZiJCO+//z5z585lzZo1Zotp//79hIWF0aFDBz74\n4INi+b37ypUrbNy4kTVr1vDNN99w3333ERoaSmhoKA8++GCpaom988477Ny5kw0bNpT4CUyys7OJ\nj48nOTmZL774gsOHD7Nw4UJLh1X2WDj5q7v03HPPySuvvGLpMMqV3bt3i5eXl+zYscPSoVwnKipK\nPD09zTrvb3x8vISFhUn9+vXl8OHDRVLmiRMnZM6cOdKuXTtxcXGRTp06yX/+8x+JiYkpkvItYf/+\n/eLl5WW23om79ccff8gLzzwjFZycxN3eXnyNRnGyshJPo1HemjxZ4uLiLB1imaLJt5Q6evSoeHl5\nSUpKiqVDKRcOHTok3t7esm7dOkuHclPff/+9VKpUSRYsWGC2OnNycuR///ufeHp6yqJFi+54IGBm\nZqbs2LFDJk6cKA888ID4+PjIsGHDZOXKlZKYmFhMUZtPamqq1KlTRz777DNLh3JTf/31l3QJCRFv\nR0d51dZWTl2zqMNBkFEODuJuby/D+vSR1NRUS4dcJmjyLcU6d+4sCxcutHQYZd6JEyfE19dXli1b\nZulQbuvo0aNSvXp1eeWVVyQ7O/uWx2ZnZ0t8fHyRfIH79ddfpW7dutK3b9/b/gYbHx8vy5cvl/79\n+0vFihUlKChIXn/9ddm7d+9tYy5tJkyYID169CixVyfExMRI9fvuk1dtbSXtNssZJoD0dnSU5kFB\npeJ39pJOk28p9s0330hQUFCJ/ccuC2JjY6VatWry8ccfWzqUQvv777/lkUcekd69e1/XSklLS5PI\nyEhp3qCBWFtZibOtrdhZW4urg4OMHDTonhZoT05OlhEjRkj16tXlxx9/LLDv2LFjMnv2bGnbtq24\nuLhI586dZd68eXL69Om7rq+k++6776Ry5cpmvxStsOLj46VuQIC8Z21dqLWEBSQbZLi9vTzWqlWh\n1yxWN6bJtxTLzs6WatWqydgxY6R+QIC4OjiItZWVeBiN0qpRI1m6dKmkp6dbOsxS69KlS1KvXj2Z\nNm2apUO5Y6mpqdKrVy9p0aKFXLx4UURE5n34ofi4uko7Fxf5EiQ97wM1B+QsyNvW1uJrNEqLoKB7\nGsn9+eefi6enp4wbN05eeuklqV27tlSqVEmGDx8uq1atkqSkpKJ6miXWlStXJDAwUL7++mtLh3Kd\nU6dOicFgkDcmTZL+9vY3XD/4VrdMkIednEw9Qa1bt5ZPPvnkrmKJiYkRZ2fnctmA0ORbSiUlJcmo\nQYPExcZGelhby3aQy3kfqH+DfAnSzsVFfFxdZcbUqeXyzX0vEhMTpWnTpjJhwoRS+9plZ2fLSy+9\nJDVq1JARAwfKA0aj/FqID9b/GAzi4+oqe/fuvaP6Ll++LMuWLZN+/fqJh4eHGI1GqVmzpmzevLnM\ndSffztNPPy1DhgwpsvL8/f3F0dFRnJ2dJTAwUMaNG3fXA7iuJt9Kbm5y6A4T79XbFyCtGjUSEZGQ\nkBD9+esuaPIthS5evChNHnhABjo4yIXb/JMcAWlqNMqAnj1L1aw6lpSWlibt27eXYcOGldrEm98T\noaFSA+SfO/hw/RqkkpubHDt27JZlHz16VN5//31p06aNuLi4SJcuXWT+/Ply5swZSU9Pl/Hjx4uf\nn1+JHCFeXNavXy/+/v6SkJBQZGUGBATI1q1bRUTk/Pnz0r17d+nevfsdl5OZmWlKvi2cnO4q8QpI\nBsh9jo7yyy+/3HXyLe/d1pp8S5mUlBR5pGFDedHWttDdRckg7YxGeaaMJJPilJmZKU8++aT0LCNf\nVmJiYqSCg4OcLcT7xAByIu/+KJCOBoOEtm1boLzMzEzZtm2bvPjii1KrVi2xt7eXkJAQWb16tSQn\nJ98whrVr14qPj49MnTq1TLymt3Lp0iXx9fWV6OjoIi03f/IVEVm2bJnUr19fRHJf36CgIHFxcRE/\nPz+ZMmWK6biriTYqKkrq1q0rISEh8ueffwogw0Fqg1QGmZXvfbAX5GEQN5BKIGPzku3V/ZtAaoE4\ngbRu3VpCQkJM3c45OTnyzjvviL+/v3h7e8vAgQNNX0JuFovBYCh3PSMimnxLnXffflu6OTqKP8iW\nQnygfgbSASQepLqT0x1/KHTq1EkWL15cTM+mZMnJyZGhQ4fKo48+KmlpaZYO55b8/f3Fzs7O9Hvu\nVUFBQWIwGEzXxb42caI8a2dXqC9p+ZOvgCSBVHBwkIMHD8rSpUulT58+4uHhIY0bN5Y333xT9u/f\nX+gPzTNnzkirVq2kXbt2Zfp60d69e8uzzz5b5OUGBATIli1bRETk7Nmz0rVrV3nhhRdEJHfe7V9+\n+UVSU1Nl1apV4uLiIhEREXLmzBnZvn27GAwGadu2raxatUq2bdsmn3/+uQDSCiQF5BeQCvk+T37M\nS8BZIN+D+IPMydv3N4gzSBTIv0BaNW8udnZ2ppbvwoULJTAwUE6dOiVJSUnSrVs3GTBggIj8f/Lt\n3r27nDhxQtLS0kzbymPy1RmuLGznzp1MnDiR3377DRsbGx544AHmzJlDkyZNrjs2OzubQB8fVl26\nRA9gIdD2DuqaB3zbqRNfrF9/w/1TpkzhxIkTLFmy5G6eSqkmIkyYMIFdu3axefNmnJ2dLR3SLQUG\nBuLg4MCYMWMYO3YskLv6TK9evTh27BinTp2iUqVKVPXyYvuVK9QqRJlWwHGgWr5tzxgMRNjY0LZj\nR0JDQ+nSpQu+vr53FXN2djbvvPMOCxYsIDw8nI4dO95VOeYgImRlZZGRkVHglpmZedNtW7duZcmS\nJUyfPh0rK6vrjrvd+bfadvDgQbKzs4HclZGsra1xdnYmKyvLdKyNjQ12dnZkZmZiY2NDhQoVMBgM\nnD17lsDAQDw8PLCzsyMnJ4d9+/YxFXgt7/n2BYzAJzd4LV4HDgNfAYuByUAMMAs4PWoU//30U+bN\nm8fQoUNp164dTZs2Zdq0aQBs2bKFzp07k5aWxunTp6lWrRrbt2+nZcuWQO7yitWqVSMrK6tUzVhW\nFGwsHUB5duXKFbp27cqUKVMYOnQodnZ27NixA3t7+xsev27dOu7LyKDxXdbXH3j922+Ji4srlWuz\nFqfp06ezefNmtm/fXuIT71X9+/dn8eLFpuQbERHBwIEDef311wHYtm0bhrQ0OgBpwJPAbMAh7/z3\n8h67A29dU/ZgwA8YLkK0lxfDhw/nzTff5KWXXsLLy4t58+bRoUMHQkJCGDBgAMOGDbtlrCJCdnY2\nEyZMIDg4mMGDB9O9e3eee+45RKTQSeleEtidlJmZmYmtrS22trbY2dkVuN1oW05ODnv37qVZs2as\nW7fupsfZ2tpiNBpxd3e/ZXnXbuvevTtTp06lZcuWZGRksGrVKubOncvx48c5cuQIb7zxBr/99hsZ\nGRnk5OTQq1cvIiIiTMnt999/x87ODshNeIGBgbjk+/s0Ar7Nu38UGA/8CKQAWcDVpkAc0CDv/gVr\na7x9faldu7apnHPnzhEcHGx6HBwcTFZWFhcuXDBta9q06S3fK+WFJl8LiomJIT09nTFjxmBrawvA\no48+CsCJEycYMWIEhw4dwmAw0LFjRzLi4xmSmFigjDigOhALeORtOwB0AM4DS8htIe8AXIHKmZnU\nq1ePrKwsatasyZw5c2jRogUbNmxg+vTpiAirVq2ievXqHDhwoNAfrqXZxx9/zKJFi9ixYwceHh43\nPEZyf6IhJyfnprfb7S/KW2pqKlZWVpw7d4758+fj4+NDeHg47777LgBr1qzhk08+ISsriy3kJtzh\nwNvANGADMBPYDgQAQ695voa8my9w5sIFevbsSfPmzWnbti2JiYm89dZbTJs2jQMHDvDnn38yd+7c\nWya7q8nsajKxsbEhPDycTz/9lCpVquDk5FToRHTtNmdn5zs6tzDH2traFnoOZhEhNDSUiRMn8s47\n79zhu69w7Ozs8Pf3N62GVKtWLd59911++OEHxo4dS//+/dm4cSN2dnb079+fazs0bWwKftQbDAYi\n7O15Ni0NgJ/I/VsDjAYCye0FcSK35bs9b19l4CCQDUTZ27OsbVumTp1qKve+++5j//799OjRA8id\nA9zGxgYfHx9Onz59w1jKK30VLKh+/fr4+vrSp08fBg8eTPPmzQt8+L/22mu0atWKI0eOMGTIEC7G\nxjL8mjLuA5oBX4Jp31KgF2B9gzprZWfjGxxMh44d2bBhA126dGHy5MlYWVkREhLCxYsX6dWrFzk5\nOUydOpWYmBhWr15NTEzMXScKcyalO60zNTWV5ORknJ2dqVWr1i3LMxgMWFlZ3fB2q33Fcbty5Qpf\nffUVjo6OzJo1C09PT6ysrNiwYQMiwvr16zl06BAPGwzUzfvbP0dui2YaEAV0Aurk7Xs6b9u1rIFM\nERo0bMjo0aOvS1jjx4+nS5cu9OzZ85aJzcbG5rpklpOTw+zZs5kxYwYzZswwfWCXNosWLSIuLo6v\nvvqqWOu5mlCTkpJYsmQJGRkZ+Pv74+zsTIUKFcjOzmb16tWsXr2a7t2737a837KyOEhuEtgMfJ63\n3ZncL/JW5LaGFwP+efu6AGPJ7a6uHBDA7t27CyT6Pn36MG3aNEaOHImnpyfz5s2jd+/e5a5LuTA0\n+VrY2rVrmTFjBqNHj+b8+fN07tyZhQsXcv/993P//fcD0KBBA0aPHs0Lzz57wz9YX3IT7nBAgOV5\nj2+kCbDir784efKk6feXn3/+GW9vbxITE0lLSyM+Pt70IX/1H8vGxuaOE4S5E9Kd1rl161bGjx/P\nhg0bqFu37m3LLUmr0QQGBjJjxgyqV69Oy5YtCQwMZNy4cYSFhWFnZ8eUKVP45ptv2JeTY+oRESAn\n7/454NF85TW6ST1/A7Y2NvTu3ZuwsLDr9ru6uhIQEECDBg2uP/k2rKysePHFF2nZsiW9e/dm69at\nfPDBBzg4ONz+5BLi1KlTvPzyy2zbts3UrVtcQkNDsba2RkSoVq0aixcvpk6dOrz//vu8+eabTJ06\nlfbt2zNy5MgC3bw3et8aDAYee+wxWq9diyPwCtA+b98UYCJQBWhKbrJdm7fPk9wvaT0MBuzj4vjn\nn39o3ry5qdyhQ4cSFxdHq1atSEtL47HHHuPDDz+8bSzlknnHd6lb2bNnj9SpU0cmT54s58+fl7Cw\nMPH19RVXV1dxdnYWR3t7icobdRgAsjXv/j8gjiDnQLbljU68OmL1U5AW+R63srISb29vcXNzE3d3\nd7GysjKNgH7zzTelf//+BWIqqxfQb9++Xby8vGTPnj2WDuWu5L/0JCQkRNzc3CQlJUUyMzPFYDDI\nyZMnxWg0iou9vcTeYGTzEJD++R5vuWa082CQySCTDQap5u8vjz/++A3jKKr3R3x8vPTq1UsaNGgg\nR44cuefyzCE7O1tat24tM2fOtHQod+Xvv/+W+ytXlgUGQ6Gv780BecHOTloFB+vsefdI+wJKkKZN\nm9KnTx/WrVvHa6+9Rnx8PL/++isJCQnMnj0bewcHVhmN153nQe5vvFdbvH1uUv52YEdODjNnziQ+\nPp7Tp0/j6Ohoat0aDIbrfisqi3766Sd69uzJ0qVLy8Tgj4ULFxIdHV1gTV0rKytGjBiBu6cnb1tb\nk0PuuIBNefufAjYCR8gdVHPtKFch93e9+dbWJKaksHbtWsLCwjh16hSxsbH88ccf/39sEbxn3Nzc\nWL58OWPGjKFly5ZERETcc5nFbe7cuWRnZzN+/HhLh3JXPD09+Wb7dqZWrMgMa2uyb3N8GjDKwYFv\n/f1ZuWlTsbf0yzpNvhb0xx9/8MEHHxAbG0t2djY//fQTkZGRNGnSBCcnJzw8PLC1tWXfvn189NFH\nGJ2cWJ+Tw183KKsvEEHub799b1LfL+R+KDdv3pwLFy4wadIk0tPTTfsrVarE4cOHC2yDovlwLSn+\n+OMPunTpwvz582nfvv3tTygFqlWrRuPG/z8G/moX+YwZM+jYqROf5OTgTm4389G8Yx4DXiT3UrVg\noAe5A6xMZZCbmGvVrcv58+d5++232bx5M9WqVaNmzZpERkaa3idF1W1oMBh4+umniY6OZsaMGQwc\nOJDEawYYlhSHDx9m2rRphIeHY219o9EVpUONGjXYdeAAa+rX536jkX9ZWV33+XIKeNnWlqoODlwK\nCeG7H3+kQoUKlgi3bLFsw7t8i42Nlaeeekp8fX3FaDSKr6+vjBo1ShITE+XUqVPSrVs3qVixorRo\n0UIWLFggfn5+MqxfP5libV2g21lAUkFcQOpd000UDtIyr7uok5OTdO7cWapWrSo1a9aU+fPnS2Bg\noKn78tKlS9KiRQtxdXWV4OBgESlb3c6nT5+WqlWryqJFiywdilmNGDBAOjs6Fpil6Ha3AyBeRqN8\n//33BcpKTEyUiIgIadOmjXh6esqYMWPkhx9+KPKZ05KSkmTYsGFSo0YN+emnn4q07HuVkZEhwcHB\npWqlq8LYt2+fDOndW1zt7cXfyUnqurpKFaNRKjo5yfixY+Xo0aOWDrFM0eRbyhw7dky8XVxk0x18\nkArIdBsbaVi9+k2nACzrLly4ILVq1ZIPPvjA0qGYXUZGhoS2aycdjUaJL8R75VsQb0dHWfHFF7cs\n99SpU/L2229LtWrVpG7dujJz5swin71q6dKl4unpKf/+979LzNSoU6ZMkY4dO5aYeIpacnKynDx5\nUg4dOiQxMTElfra30kqTbym0fft28TQaZWUhB0hMsbWVQB8fOXv2rKVDt4j4+Hhp1KiRvP7665YO\nxWIyMzNlzLBh4mFvL+Ps7OTINe+TDJAVIG1dXKSSu7tpKsPCyMnJke+++06GDRsm7u7u0rlzZ4mK\nirpuLeG7dezYMQkODpbu3bvLpUuXiqTMu7V//37x8vIqt/9Lquho8i2lfvjhB6lSsaK0dXGRFXBd\nl2ICuUvD1XV2lqb16sm5c+csHbJFpKSkSKtWrWTMmDFltqVyJ2JiYuS1iRPFx9VVqjk6Sl1rawly\ndZWKDg7SMihIli1bdk+jWJOSkmTJkiXSrl07qVixoowePVr27t17z699WlqaPP/881K1alXZuXPn\nPZV1t1JSUuSBBx6QyMhIi9Svyhad27kUS09P56uvvmLejBmcOHaM2ra2OImQYDDwa3o6bUNCeOal\nl2jTpk25vJYuMzOTJ554Ajc3N5YsWaIX+ueTkZHBggULiIqKYs6cOfj4+FClSpUireP06dMsWbKE\n8PBwbG1tGTRoEP3797/ruaEhd+au4cOH89xzz/HKK6+Y9W/64osvcvr0aaKiosrl/5MqWpp8y4hj\nx45x+vRpkpOTcXV1pVatWlSuXNnSYVlMTk4OAwYMICEhgZUrV5qm71T/Lzw8nG3bthEeHl6s9YgI\nu3btIjw8nC+//JKmTZsyePBgunXrdlcTapw9e5a+fftib2/PkiVLqFSpUjFEXdD27dvp06cPhw4d\nwtPTs9jrU2WfNgXKiBo1atCuXTsef/xxQkJCynXiFRGeffZZzp49yxdffKGJ9yYuX76Mu7t7sddj\nMBho3rw5//vf/zh79iwDBgxg4cKF+Pr6MmrUKPbs2XNHl7NVqVKF6OhomjVrRuPGjdm8eXMxRg+J\niYkMGTKEBQsWaOJVRUaTrypz3nzzTXbv3s3q1asLTDyhCoqPjzdL8s3PaDTSt29fNm3axM8//4y/\nvz+DBg2idu3aTJ8+nbNnzxaqHBsbG95++20+++wzBg8ezKuvvkpmZmaxxPziiy/Spk0bQkNDi6V8\nVT5p8lVlyuzZs4mKimLDhg24ublZOpwSLT4+/qarOJmDn58fr776Kr///jvh4eH8+eefNGzYkI4d\nO7J06VJSUlJuW0bbtm05cOAABw4coHXr1sTExBRpjOvWrWPTpk3Mnj27SMtVSpOvKjM+/fRT5syZ\nw6ZNm/Dy8rJ0OCWeubqdb8dgMNCsWTMWLFjA2bNnGTJkCIsXL6ZKlSo8/fTTfP/997fslvb29mb9\n+vU88cQTPPjgg6xcubJI4rp06RJPP/00n376Ka6urkVSplJX6YArVSasXLmSMWPG8O2331KrVi1L\nh1MqPP744wwbNoxu3bpZOpQbio2N5bPPPiM8PJzs7GwGDRrEgAEDqFq16k3P2bt3L3369KFz587M\nmjXrpgO6EhMT+WzJEravW8flS5ewsbHB29eXpwYPpmPHjlhZWdG7d28qVarEnDlziuspqnJMk68q\n9bZs2ULfvn3ZsGFDgTmO1a21bNmSd999l1atWlk6lFsSEfbt20dERATLly+nUaNGDB48mCeffBLj\nDRYaiY+PZ8SIERw/fpzPP/+8wJexmJgYZr7zDsuWLqWNlRVPJCdTEcgCYoBPnZ2Jd3KiWdu2/PDj\nj/z88886bkAVD4tcXaxUEdmzZ494enrK9u3bLR1KqVOvXj05ePCgpcO4I6mpqbJ8+XLp1KmTuLu7\ny7Bhw2THjh3XTeKRk5MjH3/8sXh6ekpERISI5L5XKrm5yas2NnLmFjPC7QVpB9K4Vi25ePGiJZ6m\nKge05atKrV9//ZX27duzcOFCunTpYulwSh0/Pz927dqFn5+fpUO5K3FxcURGRhIeHk56ejqDBg1i\n4MCB+Pv7m445dOgQYWFh1KhRgz3R0SxKTqZrIcrOAV60s2NXjRp8u2/fDVvYSt0LTb6qVDp58iSt\nWrXivffeo0+fm61grG7F2dmZc+fO4eLiYulQ7omIsH//fsLDw1m+fDkNGjRg8ODB9OjRAycnyW1I\nmAAAFCdJREFUJ+Lj47m/cmU+Sku76VrX4cBCYEf+coFBDg44hYXxcTFNRLJjxw5GjBjB77//Xizl\nq5JLk68qdc6dO0eLFi2YMGECo0ePtnQ4pVJmZiaOjo5kZmZafKpEZ2dnUwzJyck4ODiY1sj973//\ne0dfrtLT01mzZg3h4eF8//33dO/encDAQLbMmsV3t1gbOJzrky/AX0AtBwdOxsVZ9LIsVfbopUaq\nVPnnn3/o0KEDQ4cO1cR7D65OsGHpxAuQlJREYmIiiYmJ+Pv7s3btWtPjwiberKwsAOzt7enZsydr\n167lyJEj1KtXj/9Mn84zt0i8t+INdLGyIuLTT+/qfKVuRpOvKjWSk5Pp0qULHTp0YNKkSZYOp1Sz\nxOxWdyo9PZ3nn38eX19ffH19eeGFF8jIyABg27ZtVKlShfnz51OzZk2GDh3KW2+9RZ8+fRg7dize\n3t60b9+emjVrIsCTwL+A6oALUBdYdYM6Xwd8gQeA6Lxto1JSmDhxIlu3bjUdN2XKFAYMGGB6vHr1\naurWrYuHhwdt2rQp0I0cEBDA/PnzadasGR4eHvTu3Zv09HTT8yitv7mre6PJV5UK6enpPPHEE9Sp\nU4dZs2aViBZbaXb58uUS34367rvv8v3333Pw4EEOHjzI7t27mTp1qmn/hQsX+OGHH/juu+/473//\ni4jw1Vdf0bBhQ44ePUqzZs145ZVXeMDeHjtyE+9OcruShwP9gAv56tsLeACngTfITdjxQEMgKzu7\nwHsu//2jR48SFhbGv//9by5evEjHjh0JDQ01tcYNBgPz58/n3//+N/v27WPv3r3FvpiFKvk0+aoS\nLzs7m379+uHi4sKCBQs08RaB0tDyjYyMZMSIEXh6euLp6cnw4cNZsmSJaX92djZTpkyhUqVKpsk0\natWqxYgRI3B3d2fo0KEcP34cp7zjewKVAEfgeaAisC9ffVbAGMAa6AM4AesAI7mDr242PObqtcft\n2rXD2tqaZ599ljNnzrBr1y7TMQMHDuTBBx+kRo0adOzYkZ9//rkIXiFVmtlYOgBVfv3666+cOnWK\npKQk3NzceOCBBwgMDCxwjIgwcuRI4uPjWbduHTY2+pYtCqUh+Z47d47g4GDT4+DgYOLi4kyPfXx8\nruuybdiwoel+5cqVyczMJCEvMS8GPiB3Mg2AJOBivnNrAPnnw2oExAGJ5LZe83/py5+I4+LiCkzu\nYjQaqV27NrGxsaZtQUFBBeI6fvz4rZ+8KvP0k0yZVWpqKsuXL2fejBmcP32aBra2OOXkcMXKih/T\n02ncuDHPvPwyXbp0wcrKiokTJ/Lrr7+yZcsW7O3tLR1+mVEaup3vu+8+9u/fb0rA+/fvx9fX17T/\n2i9iN+sR+S0tjV+BocAa4LG87YHktmivOgakktsyBviJ3BbwTsDBzo7z58+bjj1w4IBp4Q5fX1/W\nrVtn2pecnMzvv/9eINb89AITBdrtrMxo3759VPf1JWrcON74/XdOpaSwNiGB5YmJfJOQwOm0NPrt\n2sW0fv1oWKMGL7/8Mhs2bGD9+vU4OztbOvwypTS0fPv06cOiRYu4ePEiFy9eZNGiRfTv3/+mx98s\nqbVv25a1gD25Xc0JwHRyW7X55QAfA5nA50AK0BmY5+xMUOPGLFu2jIsXL7Jq1Sq2b99uOq9Xr178\n/PPPREdHk5mZyUcffUSVKlV45JFH7vq5q7JPW77KLLZv307Pzp1ZlJLCzVZFdQAGAAOSkpidksLk\n999n044dVKhQwYyRlg+loeX7+uuvc+XKFRo0aADAU089xeuvv27af21L99qu4avbRk+YwLM7dzIt\nKYkBQBowBGiR/zjgYeAS4A+4AiuAf4AfgB2LFjF48GBq1qxJhw4dePrpp00t4Vq1arFs2TLGjRtH\nbGwsjRo1Ys2aNTf9ieTaOHUMQ/mkk2yoYnf8+HGaN27M0sRE2t3BefOtrPigcmX2Hz6sS7oVsVGj\nRtGwYcNyca20iPBw/fqE/vEHr+eNQC6MdKCD0Uj78eOZ/M47xRegKpe021kVCSsrK06ePHnDfe+9\n8w6jU1LuKPECjMrJIejyZcKLaYKDW8UcGRlJx44di6XekqA0dDsXFYPBwKrNm1lYoQKzrK0pTGsj\nGehhNOLTpg2vvfVWcYeoyiFNvuVAQEBAgQkCzCkhIYGoqChGZWff0XlWwEng2ZQU5r33ntkHqfTr\n14+NGzeatU5zKg3dzkWpcuXKfLd/PxH+/nQzGtkKN0zCKcAi4EEnJ7y7dSNy5UqsrPRjUhU9fVeV\nAzf6LcxcFkdE0NHKikqFPD5/ihagOWCXkEB0dPRNzlB3ozy1fK/y8/Njz6FDdJ41ixcCAnjA2ZlX\nrKyYRe4ArNF2dlR1cGBVmzZ8sGIFCyMjsbW1tXTYqozS5FtOxcfH07VrV7y9vfHw8CA0NLTAdYkh\nISFMnz6dRx99FC8vLwYNGkRGRgajRo3Cx8eHhx56iDNnzhQo87vvviMoKAhvb28mTpyIiPD9xo2E\npqSwiNwp+9zJvdTjdL7zrMi9BrMRUBNonbe9IbkDX6onJfH9zp2sXbuWoKAg3NzcaN68Ob/88oup\njIULF9KsWTPc3NyoXbu2KVlnZ2czbdo0qlevjqurK02aNCnwPHfv3k1QUBAVKlRg7Nixpu3h4eG0\nbNnyXl/mEqu8tXyvcnJyYtTo0Rw8eZL/rV+P05QpxI0dS/yLL1J92jT2HznC6uhoHnvsMR0IpYqX\nuRcQVuYXEBAgW7duLbDt0qVL8tVXX0lqaqocP35cOnbsKN27dzftb926tfj7+8u2bdvkxIkTEhgY\nKHXr1pWIiAi5ePGiPP744/L222+bjjcYDBISEiKXL1+W06dPy/333y+ffPKJdGjaVN4AqQKyGeQf\nkHEgj+RbwNwA0hLkZ5C0fNtO5N2fBdK3Vy9xdXWVVatWSUJCgkybNk0CAgIkIyND/v77b6lSpYoc\nPXpURERiYmLkxIkTIiIyc+ZMqV27tmnfwYMH5dKlS6aY27ZtK0eOHJEff/xRXFxcZMOGDSIi8umn\nn0qLFi2K749iYV5eXnL+/HlLh6FUuaUt33KqQoUKPPHEEzg4OHD//fczYcKEAtcuGgwGunXrRuvW\nralWrRrt27fH3t6egQMHUrFiRcLCwq77HXnAgAG4u7vj5+fHU089xbJly7Czs2MtMBpoT+7cuW+S\nO61f/nbzCHJbujeaRiMDOJw3f263bt1wdXVl4sSJJCUlsWfPHgwGA6mpqRw9epTMzEyqVq1KtWrV\nAPjkk08YOXIkNWrUAKBBgwYFLl0aM2YMtWvXpnHjxjzyyCPlYto/ESmX3c5KlSSafMuplJQURo4c\nSUBAAG5ubvTo0YOEhIQCA5vyT9Xn4+Njut4SwNvbu0D3LRScQq9Ro0bExcVRqWpVYsn9Tc0j71Yd\nsAPyn930FrGesrcnPSODyMhIPDw88PDwwNPTk9TUVOLi4qhYsSJLlixh9uzZ+Pj40KdPHy5cyJ0y\n/+zZszRv3vymZV877V9ycvItIikbUlJSsLGx0RnDlLIgTb7l1Pvvv8+ePXvYu3cvCQkJfPnll4jI\nPY0qPnDggOn+Tz/9hK+vL72HDCHZ2prXgcv5bsnkTmpw1Y2mIxBy59/90mAgODiYgQMHcvnyZdMt\nKSmJsLAwADp16sSWLVuIiYkhLS2NWbNmAbmDbHbu3HnXz6ks0lavUpanybecyMjIIC0tzXSzt7fH\n3d0de3t7Dh8+zIwZM647J38iLkxSXrp0KfHx8Zw5c4YVK1YQFhZG27ZtcXVzYy6widwu5ATgi9uU\n5QPsByKB1i1b8vzzzxMVFcWqVatITk4mOTmZdevWkZSUxNGjR4mOjiY9PZ3s7GxsbGy4eDF3yvzh\nw4ezYMECvvnmG7Kysjh06BD//PPPDeu81y8fpUV5HWylVEmiybec6Ny5M0aj0XQ7e/Ysvr6+1KxZ\nkwEDBjB06NAbTs2X//6t9kPusmmtW7emcePGdO/e3VTmy5Mn42FnxwTAE6gP5L+C9kZjSqcAzwLP\nAPUefpjg4GCWLl3KtGnT8PX1pUaNGixevBjIXev31VdfxcvLiyZNmuDh4cH06dMBGD9+PAMGDGDM\nmDF4eHgwYsQI0tLSbhh//udoycuzipu2fJWyPJ1eUhW7zMxMurRpg//+/cxPT8e6EOckA48bjTzQ\nty8f/ve/ZTYRWsLatWv5+OOPC6zEo5QyL235qmJna2vLivXrOd6gAb0cHU3rqd7MIaCN0UjV0FDm\nzp+vibeIabezUpanyVeZhaurKxt27KDG008T7OTE405OrCd3xHMCuZcdLQNaurjQ2cODvm+9xaJl\ny7C2Lkw7Wd0J7XZWyvI0+Sqzsbe3Z8acOcRcuEC3uXOZWrcuTdzc8LO3p5mHB+HNmjE+IoI///qL\n5ydM0BZvMdGWr1KWp+v5KrNzcnJi2LBhDBs2zNKhlEvx8fH4+vpaOgylyjVt+SpVzly+fFm7nZWy\nMG35KlUOJCQksGbNGs6dO8eePXuwsbEhKCiI4OBgS4emVLmklxopVYYdPHiQeR98QFRUFK1tbKie\nloZtVhb/2Nmx0cYG76pVeebllwkLC8PR0dHS4SpVbmjyVaoMEhGmv/MOH/7rXzyTkcHw7GwqX3NM\nNvAN8KGTE3E+Pqzftg0/Pz8LRKtU+aPJV6kyaNKECXwzfz7rkpO57zbHCjDb2pq5FSrw/U8/UaVK\nFXOEqFS5pslXqTIm8rPPeGvkSHalpOB5B+fNsLFhWWAg+w8fxsZGh4MoVZx0tLNSZYiIMPXVV8lK\nSeFwIc9xAf4EJmZl4XD+POvXr7+rukNCQli4cCEAkZGRdOzY8a7KUao80OSrVAkSEhJChQoVyMjI\nuKvzv/32W2zi4zkBtCrkOYlAALkLXIxJTOQ///rXXdWdfzGKfv36sXHjxtucoVT5pclXqRLizz//\nZN++fXh7e7N69eq7KmP+++/zTFLSDVeKKoxe5K7LfOLEibssQSlVGJp8lSohFi9eTPv27RkwYAAR\nERGm7WfOnOHJJ5/E29sbT09Pxo0bZ9r3v//9jzp16uDq6krdunXZs2cPbcltyUbnHbMPaAa4A5WB\ncUBmvnqtgJN59x8DfIFDhw4BEB4eTsuWLU3H7tq1iwcffBB3d3ceeughdu/efcPncu15SqmCNPkq\nVUIsXryYsLAwnnrqKTZu3Mjff/9NdnY2Xbt2xc/Pj5iYGGJjY+nduzcAX3zxBZMnT2bJkiVcuXKF\n1atXk5yWhhsF10i2AeYCl4AvgTXAvJvEYACMOTnEx8dft++ff/6hS5cujBs3jn/++YfnnnuOLl26\ncPny5SJ8FZQqHzT5KlUC7Ny5k9jYWB5//HFq1KhBnTp1iIyMZN++ffz222+88cYbODo6Ym9vT/Pm\nzQH45JNP6N27t2mWqvvvvx8Xo5HUa8puDDwEWAOPAP2B7beIJd3KCqPReN32devWYTQaGThwIFZW\nVvTr1w+j0ciaNWvu/QVQqpzR5KtUCRAREUGHDh1wcXEBoFevXkRERHD27Flq1qxJxYoVrzvn7Nmz\npkR8ld999103yvko0JXcLmc3YDa5reCbiRO54WQbcXFxNGrUqMC2Jk2aEBsbe7unp5S6hiZfpSws\nNTWVqKgooqOjqVy5MpUrV+b999/n0KFD+Pj4cPToUS5duj5d+vn5sXPnzgLbBj/7LAucnQtsGw1U\nAo6Tu3byC0DOTWJJBzLs7Hj44YeB3MFXV0cw33ffffz0008Fjt+/f7+ukKTUXdDkq5SFrVq1Chsb\nG44cOcLBgwc5ePAgR44coUWLFqxcuZJ69eoxfvx4jh8/TlpaGrt27QJg+PDhREVFERUVRUZGBseP\nH6dFixZ8n5NDVr7ynQEPcv/ZvwUW3yKWKzY2OLm7c+7cOaKjo/n666+5Og9P586dSU1NZcmSJWRl\nZbF06VJSUlLo2rVr8bwwSpVhmnyVsrDFixczdOhQqlSpgre3N97e3vj4+DB27FiWL1/O2rVruXLl\nCk2bNsXPz4+oqCgAevbsyZQpU3jjjTeoWLEiTz75JGlpaQwdNozLBoOpdTsF+BmoArwHjKXggKyr\n93cBcXZ2VK9Rg7p16zJz5kzGjBljavlWrFiRtWvXMnfuXDw9PZk9ezZr166lQoUK1z2n/Nf8KqWu\np9NLKlXGpKam0r5ZM4KPHGFORkahvmH/DHQyGlm0YgWdOnUq7hCVKve05atUGePo6Mia6GgO1K3L\nk46O/HKLY5OA+UAHo5GPIiI08SplJtryVaqMSktL473p01nw4YcEZmUxPDGR+wF74DKw1t6eSIOB\nVi1a8Nr06TRp0sTCEStVfmjyVaqMy8rKYs2aNUTOn8+52FjS09Nxd3enWfv2PP3MM7qGr1IWoMlX\nKaWUMjP9zVcppZQyM02+SimllJlp8lVKKaXMTJOvUkopZWaafJVSSikz0+SrlFJKmZkmX6WUUsrM\nNPkqpZRSZqbJVymllDIzTb5KKaWUmWnyVUoppcxMk69SSillZpp8lVJKKTPT5KuUUkqZmSZfpZRS\nysw0+SqllFJmpslXKaWUMjNNvkoppZSZafJVSimlzEyTr1JKKWVmmnyVUkopM9Pkq5RSSpmZJl+l\nlFLKzDT5KqWUUmamyVcppZQyM02+SimllJlp8lVKKaXMTJOvUkopZWaafJVSSikz0+SrlFJKmZkm\nX6WUUsrMNPkqpZRSZqbJVymllDIzTb5KKaWUmWnyVUoppcxMk69SSillZpp8lVJKKTPT5KuUUkqZ\nmSZfpZRSysw0+SqllFJmpslXKaWUMjNNvkoppZSZafJVSimlzEyTr1JKKWVmmnyVUkopM9Pkq5RS\nSpmZJl+llFLKzDT5KqWUUmamyVcppZQyM02+SimllJn9HwDoU9EPwcPoAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1075fc350>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_weight(G, edge_coeff, tri_coeff):\n",
      "    '''\n",
      "    Compute the probability weight on graph G\n",
      "    '''\n",
      "    edge_count = len(G.edges())\n",
      "    triangles = sum(nx.triangles(G).values())\n",
      "    return np.exp(edge_count * edge_coeff + triangles * tri_coeff)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def permute_graph(G):\n",
      "    '''\n",
      "    Return a new graph with an edge randomly added or subtracted from G\n",
      "    '''\n",
      "    G1 = nx.copy.deepcopy(G)\n",
      "    d = nx.density(G1)\n",
      "    r = random.random()\n",
      "    if (r < 0.5 or d == 0) and d != 1:\n",
      "        # Add an edge\n",
      "        nodes = G.nodes()\n",
      "        n1 = random.choice(nodes)\n",
      "        n2 = random.choice(nodes)\n",
      "        G1.add_edge(n1, n2)\n",
      "    else:\n",
      "        # Remove an edge\n",
      "        n1, n2 = random.choice(G1.edges())\n",
      "        G1.remove_edge(n1, n2)\n",
      "    return G1\n",
      " "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mcmc(G, edge_coeff, triangle_coeff, n):\n",
      "    '''\n",
      "    Use MCMC to generate a sample of networks from an ERG distribution.\n",
      " \n",
      "    Args:\n",
      "        G: The observed network, to seed the graph with\n",
      "        edge_coeff: The coefficient on the number of edges\n",
      "        triangle_coeff: The coefficient on number of triangles\n",
      "        n: The number of samples to generate\n",
      "    Returns:\n",
      "        A list of graph objects\n",
      "    '''\n",
      " \n",
      "    #current_graph = G\n",
      "    v = len(G.nodes())\n",
      "    p = nx.density(G)\n",
      "    current_graph = nx.erdos_renyi_graph(v, p)\n",
      "    current_w = compute_weight(G, edge_coeff, triangle_coeff)\n",
      "    graphs = []\n",
      "    #for i in range(n):\n",
      "    while len(graphs) < n:\n",
      "        new_graph = permute_graph(current_graph)\n",
      "        new_w = compute_weight(new_graph, edge_coeff, triangle_coeff)\n",
      "        if new_w > current_w or random.random() < (new_w/current_w):\n",
      "            graphs.append(new_graph)\n",
      "            current_w = new_w\n",
      "    return graphs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "graphs = mcmc(G, -1.25, 0.15, 10000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 5.92 s per loop\n"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sum_weights(graphs, edge_coeff, tri_coeff):\n",
      "    '''\n",
      "    Sum the probability weights on every graph in graphs\n",
      "    '''\n",
      "    total = 0.0\n",
      "    for g in graphs:\n",
      "        total += compute_weight(g, edge_coeff, tri_coeff)\n",
      "    return total\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fit_ergm(G, coeff_samples=100, graph_samples=1000, return_all=False):\n",
      "    '''\n",
      "    Use MCMC to sample possible coefficients, and return the best fits.\n",
      "\n",
      "    Args:\n",
      "        G: The observed graph to fit\n",
      "        coeff_samples: The number of coefficient combinations to sample\n",
      "        graph_samples: The number of graphs to sample for each set of coeffs\n",
      "        return_all: If True, return all sampled values. Otherwise, only best.\n",
      "    Returns:\n",
      "        If return_all=False, returns a tuple of values,\n",
      "            (best_edge_coeff, best_triangle_coeff, best_p)\n",
      "        where p is the estimated probability of observing the graph G with\n",
      "        the fitted parameters.\n",
      "\n",
      "        Otherwise, return a tuple of lists:\n",
      "            (edge_coeffs, triangle_coeffs, probs)\n",
      "    '''\n",
      "    edge_coeffs = [0]\n",
      "    triangle_coeffs = [0]\n",
      "    probs = [None]\n",
      "\n",
      "    while len(probs) < coeff_samples:\n",
      "        # Make the jump size larger early on, and smaller toward the end\n",
      "        w = coeff_samples/50.0 \n",
      "        s = np.sqrt(w/len(probs)) \n",
      "        # Pick new coefficients to try:\n",
      "        edge_coeff = edge_coeffs[-1] +  random.normalvariate(0, s)\n",
      "        triangle_coeff = triangle_coeffs[-1] + random.normalvariate(0, s)\n",
      "        # Check how likely the observed graph is under this distribution:\n",
      "        graphs = mcmc(G, edge_coeff, triangle_coeff, graph_samples)\n",
      "        sum_weight = sum_weights(graphs, edge_coeff, triangle_coeff)\n",
      "        p = compute_weight(G, edge_coeff, triangle_coeff) / sum_weight\n",
      "        # Decide whether to accept the jump:\n",
      "        if p > probs[-1] or random.random() < (p / probs[-1]):\n",
      "            edge_coeffs.append(edge_coeff)\n",
      "            triangle_coeffs.append(triangle_coeff)\n",
      "            probs.append(p)\n",
      "        else:\n",
      "            edge_coeffs.append(edge_coeffs[-1])\n",
      "            triangle_coeffs.append(triangle_coeffs[-1])\n",
      "            probs.append(probs[1])\n",
      "    # Return either the best values, or all of them:\n",
      "    if not return_all:\n",
      "        i = np.argmax(probs)\n",
      "        best_p = probs[i]\n",
      "        best_edge_coeff = edge_coeffs[i]\n",
      "        best_triangle_coeff = triangle_coeffs[i]\n",
      "        return (best_edge_coeff, best_triangle_coeff, best_p)\n",
      "    else:\n",
      "        return (edge_coeffs, triangle_coeffs, probs)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "edge_coeffs, triangle_coeffs, probs = fit_ergm(G, 100, 10000, True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print edge_coeffs[-1], triangle_coeffs[-1], probs[-1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-1.66727441786 -0.25833462772 1.91109415929e-05\n"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "i = np.argmax(probs)\n",
      "print max(probs)\n",
      "print edge_coeffs[i]\n",
      "print triangle_coeffs[i]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "50.5102270216\n",
        "-2.02335420406\n",
        "-0.0789307290118\n"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(8,4))\n",
      "ax = fig.add_subplot(111)\n",
      "p1, = ax.plot(edge_coeffs)\n",
      "p2, = ax.plot(triangle_coeffs)\n",
      "ax.set_ylabel(\"Coefficient value\")\n",
      "ax.set_xlabel(\"Iteration\")\n",
      "l = ax.legend([p1, p2], [\"Edge coefficient\", \"Triangle coefficient\"])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAEKCAYAAAD3mecXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYU2f7B/BvAFEElClbQVBABATELaItqDhBtCgO1KJ9\nbWut1Vrft61aW0dbV4d7Uve2dSsVcaIIiopbZA8BUVBACOf3x/kRiRCSkHGScH+uK5fm5Ml5biBw\nn+c5z+AxDMOAEEIIIRpNi+sACCGEEKJ4lPAJIYSQRoASPiGEENIIUMInhBBCGgFK+IQQQkgjQAmf\nEEIIaQSUnvAnTZoECwsLuLu7iywzd+5ceHh4oFu3brh//74SoyOEEEI0k9IT/sSJE3Hy5EmRrx8/\nfhy3bt1CUlISVq1ahYiICOUFRwghhGgopSf83r17w9jYWOTrf//9NyZMmAAA6Nq1K4qKipCbm6us\n8AghhBCNpMN1AO/LzMyEnZ2d4LmtrS0yMjJgYWEhVI7H4yk7NEIIIYRTsiyOq5KD9t7/gkQld4Zh\n6KHAx7x58ziPQdMf9D2m77OmPOh7rPiHrFQu4dvY2CA9PV3wPCMjAzY2NhxGRAghhKg/lUv4Q4cO\nRVRUFADg6tWrMDIyqtWdTwghhBDpKP0e/ujRo3H+/Hnk5+fDzs4OCxYsQEVFBQBg6tSpCAoKQmxs\nLNzd3aGvr48tW7YoO0Ty//z9/bkOQePR91g56PusePQ9Vn08Rh43BjjA4/Gkvqdx8N5BjD80Xmy5\n7/t8j697ft3Q0AghhBC5a0jeq0nlRukr0p28O5jmOw3f9/leZJnD9w/j4L2DSoyKEKLKTExM8OLF\nC67DII2IsbExCgsL5X7eRpXws0uy4WbuBgNdA5Fl3Fu5Y8nFJUqMihCiyl68eCGXEdKESEpR085V\nbtCeImUXZ8PKwKreMk4mTnjy4gmqmColRUUIIYQoXuNK+CXZsDKsP+Hr6+rDVM8U6S/T6y1HCCGE\nqJPGlfAlaOEDQHvT9nhY8FAJERFCCCHK0WgSPsMwyCnJEdvCByjhE0I007Nnz6ClpYWqKvW4ZZmX\nl4eIiAiYm5tj9uzZAIBvv/0WDg4O6NatGy5evAgXFxex59mxYwf69++v6HBVXqMZtFdQWgB9XX00\n02kmtmw7k3Z4VPhICVERQohs7O3tkZeXB21tbcGxiRMn4rfffuMwKvlYt24dcnJykJmZCV1dXVy4\ncAFr167FgwcPYGpqCgASbaEeHh6O8PBwucSkpaWFx48fo23btnI5nzI1moQvaXc+wLbwo1OiFRwR\nIYTIjsfj4ejRo+jXrx/XochdamoqfHx8oKurK3ju7OwsSPZcUddZG42mS1+SAXvVqEufEKIJqqqq\nMGvWLLRq1QpeXl6IjY0Vej0lJQV+fn4wMzNDcHAwZs2ahXHjxglev3r1Knr06AEjIyN06tQJ58+f\nF1lXeno6QkJC0KpVK5iZmeHzzz8XxPDjjz/C3t4eFhYWmDBhAl69eiW2joiICERFReHXX3+FoaEh\n1q9fj8jISMTFxcHQ0BALFixATEyM0O6qomLYunUrevfuLSh3//59BAQEwMTEBC4uLti3b5/gtYiI\nCMycORNhYWEwMTFBt27d8PTpUwCAn58fAMDT0xOGhoZC71MLjJqSNvStiVuZ8APhEpUtryxnmi5s\nypRXljckNEKIBlH1P5P29vbM2bNn63xtzZo1TLt27ZiMjAymsLCQ8fPzY7S0tBg+n88wDMN069aN\n+fLLL5mKigrm4sWLjIGBATNu3DiGYRgmIyODadGiBbNlyxbm1atXzLZt25iWLVsyz58/r1VPZWUl\n4+HhwUyfPp158+YNU1ZWxly6dIlhGIbZtGkT4+DgwKSkpDAlJSXMsGHDxNaRn5/PMAzDREREMN99\n952gnq1btzK9evUSPD937hxja2srNoYtW7YI3ldSUsLY2NgwixYtYl68eMEcPXqUMTQ0ZJKTkxmG\nYZgJEyYwLVq0YA4ePMg8f/6cGTx4MBMWFiaok8fjMU+ePJH2xyQVUZ85WT+L1MKvg662Lmxb2CLl\nRYrE53/64mlDQyOEqDkeT/ZHQzEMg+HDh8PY2Fjw2LRpEwBg7969CAsLg42NDYyNjTFhwgRBd3Ra\nWhri4uLw2WefQUdHBz179oSXl5fgvNu3b0fXrl0REREBQ0NDjB8/Hg4ODjh+/HitGK5du4a7d+/i\n+++/h56eHpo2bYoePXoAYAfMhYWFwd7eHvr6+vjss8+we/du8Pl8kXUcO3ZM6Our6//SxFDT0aNH\nYWBggLlz58LIyAiDBg1Cv379hFrrffr0QXBwMMzMzBAeHo6bN29K+uNQaY3nHn5JNtq0bCNx+epu\nfWczZ7Fl7+ffh+ufrjj80WEMcxkmS5iEEDXE5S1dHo+HI0eO1HkPPzs7G506dRI8r5nQs7KyYGBg\nIDT4zNvbW7Cka2pqKi5cuABjY2PB65WVlcjJyalVT3p6Otq3b1/nvfXs7Gz4+PgInvv4+KCyshK5\nublS1SFOfTHUlJqaipSUFKE6+Xy+YBt2Ho8n9D2ztLRESUmJ1PGoosaT8Iuz0c2mm8Tl25u2l3ik\nfmxqLHytfRH5TyScTJzg1sqtoWESQojcWFlZITExESEhIQCAhIQEoddKSkrw9OlTQdJPSEiAvb09\nAKB169bw9/fHiRMnxNZjZ2eHhw8foqCgoFbCtba2Rnx8PEaMGAEAiI+Ph46ODiwtLaWqQ5YY3i/n\n6OiI5ORkkWXq60mQ5HVVRV36IrQzaSfxwL0LaRcQ6R2JXwN/xfA9w1FYKv9NDwghRBRRCWjUqFHY\nu3cvMjMz8eLFC2zfvl3wWps2bdC1a1csXLgQOTk52Llzp1DX9dixY3HlyhVERUXhxYsXKCsrQ0xM\nDDIzM2vV07VrV3Ts2BEzZ87E48ePUVZWhsuXLwNgt0Tfu3cvnj17hpKSEqxevRphYWHQ0tISW4c0\nibVLly4iY6hp8ODBKCkpwa+//oqcnBxUVFTg+vXrgul94uq0sLBAfHy8xHGpksaT8KWYlgdIN1L/\nQuoF9G7TG+M9x2Oo81CE7Q9DZVVlQ0MlhBCpDBkyBIaGhoJHdWs6MjISgwYNQqdOndC3b19MmjRJ\naGOWHTt24MmTJ+jQoQP279+P8PBwwRQ4W1tbnD59Gps3b4ajoyNat26NZcuW1bloj5aWFv755x+8\nevUKXbt2hZ2dHfbu3QsAmDRpEiIiIuDn54e2bdvC0NAQv//+e711VCddHo8nFO/7z6uPAYC2trbI\nGGq+z9DQEGfPnkVMTAzc3d1hZWWFuXPn4u3bt2LrAID58+fjiy++gJGREfbv3y/Vz4lrPEZN+yak\n3RfYYJEBsr7KQoumLSQqn1qUip6beyJjZka95dJepqHz+s7InZULHo+HyqpKBO0IgoeFB34N/FXi\n+AiRRlZxFqwMrBS2qxZ5R9Y9yNWJr68vBg8ejHnz5nEdSqMm6jMn62exUdzDLy4vBgMGhrqGEr/H\nrqUdCkoL8Prta+jr6ossV926r/7Dq6Olg92hu9FlQxfEPIuBtta71a88LDywYciGhn8hGuJJ4ROM\nOTgGZZVlgmM6WjrYNHQTOll2quedpFr3Td2xYcgGBDoGch0KUWPx8fEwNjaGpaUloqKicPPmTWzc\nuJHrsIiCNIqEn12SLXVrSIunBUdjRzwufAxPS0+R5S6kXUDv1r2FjpnomSDu4zg8LnwsdDx4TzAe\n5D+QaOS/Jvvj+h/obN0ZU7ynCI7tvLMTq6+vxvoh6zmMTD3kvc5D2ss0bE/aTgmfyCQnJwchISEo\nKCiAnZ0dNm7cCE9P0X/viHprHAm/WLoBe9Wq7+OLS/iR3pG1jps2N4Vpc+GRoqPcRmHP3T34vs/3\nUseiKUorSvHXrb8QPyUe9kb2guPGesbwWueF3wf+jqY6TbkLUA0kZieiY6uO+PvB32J7oAipz+DB\ngzF48GCuwyBK0igG7VW38KXVzrT+TXTy3+Qj41VGvRcENYV1DMOuO7sazf3Auuy5uwddbbsKJXsA\naN2yNTq26ogTj2WfnqPpEnMSEegYiJ6te+Lw/cNch0MIURONI+E3tIVvUv9I/YtpF9Hdtjt0tCTr\nKOlq0xWlFaW4nXdb6lg0xdr4tfjE55M6XxvrPhY7bu9QckTqJyE7Ad6W3hjrPhbbb28X/wZCCEFj\nSfgNbOGLm5pX1/37+vB4PIR1DMPuO7uljkUTJGYnIqs4C0Htgup8PbRDKE4/OY2XZS+VHJl6SchO\ngJeVF4Y6D8WV9CvILcnlOiRCiBqghF8PcV36F1IvwK+Nn1Tn/MjtI+y5u6dRduuvu7EOkd6RQjMX\najLWM0Y/h344cO+AkiNTHy/LXiKnJAfOps7Q19XHUOehjfYCkhAiHU4SfmxsLLy9veHh4SFYgKGm\nmJgYtGzZEl5eXvDy8sKPP/4oU30N7dK30LdAeWV5nSvnlbwtQfLzZPja+Ep1zk6WnaCjpYP4LPVc\nqamhisuLsefuHnzs/XG95cLdw6lbvx43c27Cw8JDcNE01oO69QkhklF6wufz+Zg0aRIOHjyIGzdu\nYNOmTbh3716tcn369EFiYiISExPx7bffylRnQ1v4PB6PXVO/oHYr/0r6FXhZeaGZTjOpzxnWMQy7\n7zauVtn2pO34wOEDsRdeg9sPRmJ2IjJf1V6+k7AD9rys3m2A0s+hHzJeZeB+/n0OoyLq4MKFC3Bx\ncVF4Pc+ePYOWlladK/JxZeXKlXBxcYG1tTUA4NSpU/Dz80OLFi1w8+ZNdOzYEbGxsWLPY2hoiGfP\nnik4WsVR+rS8a9euwcnJSbBBQ1hYGI4cOQJXV1ehcpJ0ec+fP1/wf39/f/j7+9dZrqEtfOBdt35X\n265Cxy+kSd+dX+0jt48Q+Fcgfgn4BVo8zb+rwjAM1sSvwfL+y8WWbabTDMGuwdh9Zze+6vGVEqJT\nLwnZCejTpo/guY6WDkZ3HI0dt3dgYd+FHEZGuGBgYCBYX+T169do1qwZtLXZ3p/169dj9OjRgrK9\ne/cWrBffmKSlpWHWrFl49OgRHBwcAADTpk3Dd999h4iICADAnTt3JDpXcXGxXGKKiIiAnZ0dFi6s\n/3c2JiYGMTExcqkT4CDhZ2Zmws7OTvDc1tYWcXFxQmV4PB4uX74MNzc3tGnTBr/++is6dOhQ61w1\nE74oZZVleF3xGqZ69W+ZKIqogXuxqbH4ptc3DTpnB/MOMNEzwaW0S+jdRvJBf+rqasZVlFaWop9D\n7e076xLuHo5Zp2epfcIvryzHhoQNyCzORMGbAhSUFqCorAhLP1yKztadG3TOhOwEfNntS6FjYz3G\nYsTeEfjB/wdaareRqbltq4ODAzZt2lTnNrmVlZXQ0WkUy67UkpaWBgsLC0GyZxgGaWlp6NmzJ8eR\nifd+Q3bBggUynU/pzUtJ/iB5e3sjPT0dN2/eREhICIYOHdrg+nJKcmChb9HgP4R1Tc0rryxHfFY8\netj1aHBcYR3DsOfunga/Xxk2J25G+9/by/wYsmsIpnWeJnFvRp82fZD7OhfJz0VvX6kONiZsxLZb\n22CoawhvK2+EuYWhY6uOWBO/pkHne1PxBk9fPK21/bKXpRf0dPRwOb32zmCkcYqJiYGtrS3Wrl2L\n9u3bY/LkyYiJiRFqbC1ZsgROTk4wNDSEm5sbDh9+t6bD1q1b0atXL/zwww+wsbFB27ZtcfLkScHr\nKSkp8PPzg5mZGYKDgzFr1iyMGzeuzlhevnyJyZMnw8rKCjY2Nvjuu+9EdvdXVVVh0aJFcHJyQosW\nLdC5c2dkZLD7mVy+fBm+vr4wMjJCly5dcOXKFbF1nD17FoGBgcjJyYGhoSHGjBmDFi1agM/nw8vL\nC+3atQMA2NvbIzo6GgB72/n9GKp379PS0sLTp08BAOXl5Zg1axbatGkDCwsL/Oc//0FZWZnQ93/9\n+vVo27YtrK2tsXXrVgBsz8vOnTvx888/w9DQEMOGDZP8BysjpV/y2djYID09XfA8PT0dtra2QmUM\nDd+teT958mTMmTMHhYWFMDExkbo+WbrzgbpH6t/IvgFnM2eJN+Kpy0duH6HH5h5YOWClxPP4lamK\nqcJPF37Cyv4rZV4KmAcenEycJC6vraWN0R1HY92NdVgWuEwlvz/iVPAr8PPln7Fv5D50sekiON7D\nrgfc17jjz6A/pR7/cTv3NlzMXKCrrSt0nMfjYbzneHy0/yNYGlgKjutq62KC5wRM9JpY6z1E8+Xm\n5uL69euIjY2FkZERrl69KvS6k5MTLl68iJYtW2Lt2rUIDw/H06dPYWFhAQC4fv06AgMDcefOHaxa\ntQqTJ08WJL4xY8age/fu+PfffxEXF4cBAwYgODi4zjgiIiLQrFkzXLlyBS9fvsTYsWNhZ2eHKVOm\n1Cq7bNky/PXXXzhx4gTatWuH27dvo3nz5igsLMSgQYOwatUqjB07Frt27cKgQYPw5MkTGBsb11vH\niRMnMHbsWKG8o6WlhaSkJLRt2xaA8A55y5cvF4ohKSkJenp6tWL95ptvcPv2bRw/fhzNmjXDxx9/\njB9++AGLFi0CAOTl5eHmzZuIi4vDkSNH8OmnnyI4OBhTpkzBlStXYGdnhx9++EHaH6tsGCWrqKhg\n2rZty6SkpDDl5eWMp6cnk5ycLFQmJyeHqaqqYhiGYY4cOcLY2NjUOo+koR9IPsAM2zWswfG+KH3B\nNPuxGRN+IFzw6LqhK/PFiS8afM5qXTZ0YbzXeTM9N/UUPGadniXzeeUhJiWGcfvTTfBzULZHBY8Y\nzzWejPESY2bk3pHM5oTNTOarTE5iaYjNCZuZD6M+rPO1ftv6Mfvv7pf6nKuvrWYmH5lc52uV/Eom\nPjOeuZ55XfA49fgUExAVwDisdGA2J2xmKvgVUtepDir4FUxpRanQg1/Fl9v5Jflbg/mQ+SEre3t7\nJjo6mmEYhjl37hzD4/GYtLQ0wevnzp1jbG1t63xvVVUVY2dnx/z9998MwzDMli1bGGNjY4bPZ7+P\nWVlZDI/HY3Jzc5nU1FSGx+MxT548Eby/d+/ezLhx4xiGYZiUlBSGx+MxfD6fycnJYXR0dITiWLFi\nBdO3b98642jfvj2zYsWKWsejoqIYa2troWM2NjbMtm3bxNZR19f9fvw1v3eiYqj5vqqqKkZfX5+5\ndOmS4LVDhw4xDg4Ogjq1tbWZ58+fMwzD5j0DAwMmLi6OYRiGiYiIYL799ts662AY0Z85WVO20ptO\nOjo62Lx5M4KDg1FZWYnIyEi4urpi3bp1AICpU6di//79WLNmDXR0dODh4YEjR440uD5ZW/hGzYyw\nb+Q+FJUVCR3va9+3weesdnT0UTwoeCB4XlxejPCD4fj5w585vxe75eYWTOw0kbM4nEyccPOTm8gq\nzsKpx6dw8slJzDg1A7ERsRIvZcwVfhUfiy8uFrkRUPXUwxEdRkh13sScRHhZetX5mraWNnysfWod\nD3QMRGxqLL4/9z0WXVyElf1XYlD7QVLVq8rKK8vh+qcrsoqzBMcYMPCy9MKxMcdq7WehKMw81VtX\nw8LCQqgL/31RUVFYvnw5UlNTAbDjAfLz8wWvu7m5QUuLvQ1nZWUlKJOXlwcDAwNB6xgAfHx8hN5b\nLTU1FVVVVfDw8BAcq6qqQuvWreuMKSMjo85761lZWfDyEv7sV3e1p6WlSVWHOKJiqOn58+d48+YN\nBg1697vEMIzQrQorKyuYmZkBYPOemZmZ0JgLLnDSV1o95a6mqVOnCv7/6aef4tNPP5VLXQ2dklfT\n4PaK2VzCXN8c5vrmQsd0tXWRWZwJ2xa2It6leMXlxTh8/zCWfriUsxiqWRtaY6LXREz0mogp/0zB\nhbQLKp/w9yfvh7m+udBo+ppGuI7Al6e+xIvSFzDWM5b4vAnZCZjYaaLU8fi18cO5Cedw9ulZRByJ\nwKzuszCj2wzOLyrlYe/dvXAyccLTL54KjjEMg7nRc9F7S2+cHnea098lLtU3SC81NRWTJk3CP//8\ngwEDBgBgB/0xEsyOsrKyQklJCZ4+fSpI+jdu3BDMvKrJzs4OWlpaSE5OFlw01MfOzg4XL16Er6/w\n+iY2NjZISEgQOhYfH4+QkBCp62hoDDWZmZlBT08Pp06dQpcuXUSWE0XWfe0bSuPnhMkj4SuTh4UH\nknKTOI1hX/I++Nv7w8LAgtM43tfZurPKL1hUPfbhf73/JzKhtmzWEgFtA6RaUbCCX4Hk58nwsPAQ\nX7gOPB4PAY4BuDzpMjYlbsJnJz5DZVVlg86lKhiGwaq4VZjedbrQcR6PhyUfLkFEpwj03tK7znU0\nGrumTZuiadOmMDU1xcuXL7F48WJkZWWJfyOANm3aoGvXrli4cCFycnKwc+dO3Lx5s86yVlZWGDRo\nEObMmYN79+6hqqoKT548ETnn/eOPP8a6detw4sQJVFZWIikpCYWFhQgKCkJpaSn++usvVFZWYufO\nnXjz5g0GDx4MS0tLqeoQR1QMNWlpaSEyMhLff/89EhISUFVVhczMTJw+fVqiOiwsLJCUlITKSuX+\nDmp+wpexS1/ZVCHhV3fnqxp1SPhHHx6FjpYOBjoNrLectCsKJj9Phr2Rvcxb4bYxaoNLky7hUcEj\nDNs9DMXl8plXzIWrGVdRVFYkcm+Gr3t+jf/1/h/6bO2DxOzEOstosrouOKuPWVpaYvHixRg3bhw8\nPT3x9u1b9OrVS6jc+++v+XzHjh148uQJOnTogP379yM8PBy6urp1lo2KioKBgQECAwNhbGyMkSNH\nIicnp86YZ86ciXHjxuHTTz+FsbExIiMjUVZWBhMTExw9ehSrVq2CmZkZVqxYgaNHjwoGcouro76v\nRdIY3n/f0qVL4ePjg9DQUBgZGSEgIAAPH76b0VVfHZMnT0Z6ejpMTU0REhIispy88Rgu+hXkQNIu\nEa91Xtg4ZGOd9zdVUdStKJx8fBI7R+zkpP6HBQ/ht8UP6V+mo4l2E05iEOUt/y2Mlhjh+eznKrkH\nPMMw6LapG2b3mI3QDqH1li2vLIf1cmvc+uSWRF3OWxK3IDolGttD5LOMbgW/Ap8e/xTnU8/Dzdyt\n3rJ2Le2wsv9KlbsFMPrAaHS16YoZ3WbUW+5A8gH859h/cGPKDdi1FH1PWxSuul/Via+vLwYPHox5\n8+ZxHYpGEPWZk/WzqH7znaSkji38ny/9zFn9W29uRbhHuMole4Ad3+DWyg03c26iZ2v5LZqxP3m/\nRN9zEz0TOJk4oZ1JO7QzbQd7I3uhKYMJ2Ql4WfYSwS51T0+qqalOU4S4hmDX7V2Y3XO22PL1Ddhr\niCbaTbBu8Dr8m/IvXpbXvzvhV6e/QoRnhNCSvlzLfJWJU49PYe2gtWLLjugwArdyb2HO2TmcXUhr\nmvj4eBgbG8PS0hJRUVG4efMmNm7cyHVYRAyNTviVVZUoKC1AK/1WXIciMVczVzx58QTlleVoqtNU\nqXXzq/hsD8PYk+ILc6S6W1+eCX/rza0Y2WEk+tjXPcgOYFvv+W/y8ajwER4WPsSxR8eQ9jINVcy7\nUbk8Hg+/BPwicjfA94W7h+OLk19IlPATshMkupCQBo/HwwdtPxBb7krGFRy6f0ilEv7aG2sR7hGO\nls1aSlR+Ts85cPnTBZfSLsn1s9NY5eTkICQkBAUFBbCzs8PGjRvh6anag2mJhif8vNd5MNUzVauF\nW5rqNIWjsSPu5d9DJ8tOSq37zNMzsDK0QsdWHZVarzQ6W3VGTGqM3M5Xwa/AhbQL2Dp8K8yam8nt\nvJLwa+OHwtJC3Mm7U+/3nF/Fx63cW5wl3GCXYEw9OhU/9FXyIiEilFWWYf2N9YiNkHxQlr6uPpZ8\nsARfnPwC1yKvNYo9LBRp8ODBGDxYMbOXiOJo9Kde3brzq3E1cE9VB+vV5GPtgxtZN+R2vmuZ1+Bo\n7Kj0ZA8AWjwtjO44Gr9f+x3XMq+JfPz94G+YNzeHUTMjpccIAN1suyH/TT4eFz7mpP737bmzB95W\n3lKvADnGfQx0tXWx7eY2BUVGiGpTn6ZvA6jblLxqXCT8NxVvcPzRcYnuiXLJzdwNqS9TUVxeDMOm\nhuLfIMbZp2cl6tZWlCk+UxBxOAKfHf+s3nLjPOtep1wZtHhaGOY8DIfuHZLo9oMiVU/F+6nfT1K/\nl8fjYdWAVRi2exhGdBgh09LYhKgjzU74atzCX3l1pVLrvJR2CZ4WnlItBMOFJtpN4N7KHYk5iQ3e\nnrimsyln8W3vb+UQWcM4mTjh4qSLnNUvqWCXYCw4v0DpCX/JxSW4lXtL8Lx698v+Tv0bdD5fG1/0\nd+qPny78pBILSxGiTJqd8KmFL7HolGhOW7rSqB64J2vCL3lbgsTsRPRq3Ut84Uaur0NfhB0IQ1Zx\nFqwNrZVS55X0K/jz+p/4+UPhGRQ+1j4y3YNf1G8R3Ne4Y1C7QTBv/m6lSwNdgzqn7RkbG6vclESi\n2YyNFdPw0viE39FcdQegiWJjaIO3/LfILclV2mp3/6b8i18CflFKXbLysfLB2ZSzMp8nNjUWvja+\nKjmnX9XoausiqF0Qjtw/gv/4/kfh9TEMg2+iv8EC/wUY7T5arue2MrTCog8WYerRqULH016mIePL\njFq9XO+vskaIuqJBeyqIx+PBw8IDt/NuK6W+orIi3Mu/h2623ZRSn6zkteLe2adn8aHDh3KIqHEI\ndgnGofuHlFLXyccnkfc6D+M9xyvk/FN8puDep/eEHn3a9EF0SrRC6iNEFWh2wlfTLn1Aud36Mc9i\n0N22u9Ln/TeUq7krMl9l4mVZ/QvGiMP1gD11M8BpAK5mXMWL0hcKraeKqcLc6LlY1G+RUqfUDnAa\ngBOPTyitPkKUTWMTfsGbAjwqeARHE0euQ2kQZSb8f1P+xQcO6pP4dLR04GnpiYTsBPGFRcgtyUXa\nyzR0tu4sx8g0m4GuAfzt/XHs0TGF1rP7zm4002mG4S7DFVrP+wY6DcTJxydpGV2isTQ24a+4ugKh\nHULVapW9mpSZ8NVpwF41Hysf3Mhu+Hz8f1P+hb+9v1otyqQKFN2t/5b/Ft+d+w5LPlyi9IFyTiZO\naKbTDHfXHW3yAAAgAElEQVTy7ii1XkKURSMTfsGbAqyJX4P/9v4v16E0mJu5G+7n31f4FqY5JTnI\nKs6S6zrtyiDrffyzKWfxYVu6fy+tIc5DcPbpWZRWlCrk/BsTNqKdSTv42/sr5Pz14fF41K1PNJpG\nNm9WXF2BENcQ2BvZcx1Kg+nr6sO2hS0e5D+AW6v6dzOTRXVLV9L131VFZ+vO+OF8w5Z6ZRgGZ56c\nwewe3C4io47MmpvB28obk/+eDCcTJxjoGsBA1wDNmzQHD/W3yHk8HnS0dNBEqwmaaDeBrrau0PQ6\nhmHwY+yPODZGsbcM6jPQaSCWX1mOr3t+zVkMotx7fg95r/OEjlVUVaDkbQlev33N/lvxutYtif5O\n/VV6uWyiPBqX8AtLC7Emfg3iI1V733RJeFp6Iik3SeEJv599P4WdX1GcTZ2R+zoXL0pfSL1Y0OPC\nx+AzfDibSrc0K2GtDlqNE49PoORtCfLf5ONZ0TO8rngt9n1VTBUqqypRwa9ARVUF3vLf1kpO03yn\ncbpJj7+9P0YfGC23lRzlgV/Fx8LYhVgTvwYuZi5Cr+lo6QguuqovvGpeRL0ofYHNNzcj6ZMktbuo\nJ/KnMQn/+nXg5k0gzYlt3TsYO3Adksw8WnkgKS8JoyHfecg1RadEY2b3mQo7v6Joa2mjk2UnJGQn\nSD3+IDolGh+2/ZAWU2kgV3NXuJq7ch2GQhjoGqCrTVf8m/IvhrkME1s+MTsRl9Mviy3na+MLX2tf\nqT9zuSW5CD8YDj7Dx61PbsHSwFKq9zMMgx6be2Bf8j6EdQyT6r1E82hEwq+sBD7+GKhqWois0NUa\n0boH2IF76xPWK+z8T188RVllGVzN1POPd/V9fGkT/tmnZzHMWfwfc9I4DXQaiJNPTopN+IWlhQja\nGYQh7YegiXYTkeUqqyqxKm4V+AwfYR3DMLrjaIm62GOexSD8YDgme03GvD7zGtRC5/F4+MH/B3x+\n4nOM7DCSWvmNnEYk/HXrAH194LrBCox3DtaI1j2g+JH6/6b8i34O/dS2pdvDtgfGHhqLpZekWxP9\nTcUb/DbwNwVFRdTdAKcBGLxrMBiGqfd3Y/aZ2RjZYaREnyWGYZCQnYDdd3dj4I6B4FfxoautW+97\nyvnl2DZ8GwIdA6X+Gmr6sO2HMNc3x647uzDWY6xM5yLqjceo6aRTHo8HhmGQnw906AAcPFGIPvvb\n48yoa+jn1Zbr8OSiiqmC0RIjPJvxDCZ6JnI//+gDoxHQNgCTvCbJ/dzKwDAMisqKwEC6j3ATrSYq\nc3+WqB6GYdB6ZWucGXem1j3zaudSzmHC4Qm4O+2u1J+lKqYK6S/TxX5uTfVM5fY5PZdyDlOPTkXy\np8k0FVWNVee9huLkJx8bG4sZM2agsrISkZGR+Pzzz2uVmTt3Lo4dO4bmzZtj69atcHGp+xfvu++A\n4aNfYM6tIWhTNB6vUtsC6jXDTCQtnha62HSB9TLrejcL0dbSxnCX4ZjVfRY8LT0lOjfDMPg35V8s\n/mCxvMJVOh6Pp/K7+xH1w+PxBIvw1JXwSytKMeXoFPwR9EeDErIWTwttjNrII1SJ9XXoC2tDa+xI\n2oEJnSYotW6iOpSe8Pl8PiZNmoSzZ8/CxsYGvr6++PDDD+Hq+u4+8vHjx3Hr1i0kJSUhLi4OERER\nuHr1aq1zJSYC+09lw3xmIAbYBEDnwa+4excYrtwFuhTq5NiTeMt/W2+Z129fY3PiZgzcMRDuFu6Y\n3WM2PnD4oN7uyLvP78JA10Ctpy4SoigDnAZg/Y31mNFtRq3XfrrwEzwtPDHUeSgHkTXcAv8FmPT3\nJIxxH1PvmAOiuZSe8K9duwYnJyfY29sDAMLCwnDkyBGhhP/3339jwgT2KrRr164oKipCbm4uLCyE\nd477+Osn4E0KxFjPyZjbay7+yuXh5EmlfSlKoaOlI7YLrnmT5pjTaw5mdJuBHbd3YPqJ6cguya73\nfeWV5Qj3CJd3uIRohA8cPsCEwxPwpuINmjdpLjh+J+8O1t1Yh1uf3OIwuobpY98HDkYOiLoVhcne\nk7kOp05FZUWYHzMfs3vMhk0LG67D0ThKT/iZmZmws3u357StrS3i4uLElsnIyKiV8G+WemFgWh+8\njX6L8/zzcHPzx7Jlio1flTXVaYpJXpMQ0SkCBW8KxJan7nBC6tayWUt4W3lj0YVFaG/aXnD8j2t/\nYGHfhbA2tOYwuoZb4L8AYw+NxXjP8SrXyo/PiseofaPQRJtdmEldtutWpJiYGMTExMjtfEpP+JKO\nCH9/YEJd75s/cyO+CxkleP7mDfDoETtNT6cRj0vR4mnBXN+c6zAIUWtf9/gae+7uQdrLNMGxD9p+\ngCk+UziMSjY9W/dEW+O22H1nN8Z5juM6HADs3/o/r/+JH87/gNWDVqOzdWd0Xt8Z8/rMg4GuAdfh\nccrf3x/+/v6C5wsWLJDpfEpPizY2NkhPTxc8T09Ph62tbb1lMjIyYGNTu3unZrIHgObNASsr4MkT\nwJkWUSOEyGBQ+0EY1H4Q12HI3ZyeczDz1EyEe4TXOxhYVqUVpTj5+CT4DL/ecnvv7sWjwke4PPky\nnEycALCDDLckbsHnXWsP6CYNp/SE37lzZzx69AjPnj2DtbU19uzZg127dgmVGTp0KP744w+EhYXh\n6tWrMDIyqtWdL4qbG3D3LiV8QgipS0DbAOhq6+L4o+MY3H6wQupIf5mO4D3B0GuiBwv9+v92tzNt\nh6jgKDTTaSY49mW3LzH+0HhM851GiwXJkdITvo6ODjZv3ozg4GDBtDxXV1esW7cOADB16lQEBQUh\nNjYW7u7u0NfXx5YtWyQ+f3XCDwlR1FdACCHqi8fjYU7POVhycYlCEv7FtIsYtW8UZnSbgdk9Zjdo\nYa/utt1homeCow+PSrTEMZGM2i+8877t24GjR4HduzkIihBC1EBlVSWc/3BG1PAo9GzdU2z5t/y3\nKCwtFDqmzdOGiZ6JUAt8Xfw6fB/zPbYN34YBTgNkinHX7V1Yn7Ae5yack+k8mkTWhXc0LuEnJgLj\nxwO3b3MQFCGEqIm18Wtx/NFx/D3673rLpb9Mx4AdA/D89XOhe/4VVRV4WfYSxnrGaKXfCno6enhT\n8QZHwo6gnWk7meOr4Feg7W9t8c/of9DJspPM59MElPDfU1oKmJgAr14BTVRr1gkhhKiM0opStP2t\nLc6MOyNyM5+7eXcxcMdAfNH1C3zV46tar/Or+CgoLUDe6zzkv8lHZ+vOch1Zv/TiUiTnJ2Pb8G1y\nO6c6kzXhix2imZ+fj6VLl2LoUHZVqeTkZGzatKnBFSqanh5gaws8fsx1JIQQorr0muhhepfp+OVy\n3fPdL6ZdRL+oflj8weI6kz3ALuvdSr8VOrbqCH97f7lPo5viMwX/PPgH2cXZcj1vYyW2hR8YGIjh\nw4dj7dq1SEpKQkVFBby8vHDnzh1lxVin+q50hg8Hxo4FQkOVHBQhhKiRorIiOP7miNiIWFgZWgmO\n/5vyL6Ydm4YdITsQ4BjAYYTAp8c/xf38+3Azd5Pqfbraugh0DEQ/h34as2GQwrv0XVxccP/+fXh5\neSExMREMw8DZ2RkPHz5scKXyUN8X/r//Abq6wLx5Sg6KEELUzOILi2u18k30TLAndA98rH04iuqd\ngjcF2HVnF6qYKqneV1xejL8f/o1nRc8Q2iEUozuORmfrzuDh3awBbS1ttboYUPhueUZGRsjIyBA8\nP3jwIMzNVXsVNzc34PBhrqMghBDVN7f3XMztPZfrMEQybW6Kz7p81qD3/s/vf3hS+AR77u7Bf479\nBw8LhBuqLZq2QPK05EazMqnYFv6NGzcwceJEpKamwsSE3ZP98OHD8PSUbBtWRanvSufWLWD0aCA5\nWclBEUIIURtT/pkCK0MrLPCXbclaZVHKKP2Kigo8ePAAAODs7IwmKjD8vb4vvKwMMDJiR+rr6io5\nMEIIIWrhYcFD9NrcCylfpEBfV5/rcMRSeMLftm2boJKaKyaNHz++wZXKg7gv3MUF2L8f6Fj3bBNC\nCCEEoXtD0bt1b3zR7QuuQxFL4ffwr1+/Lkj0BQUFOH36NAIDAzlP+OJUL7FLCZ8QQogoc3rOQei+\nUEzznaZyWwbLm9iE/8cffwg9z8zMxKRJkxQWkLxUJ3xCCCFEFF8bXzgaO6rUlsGKIvV8hJYtWyIz\nM1MRsciVmxuwdy/XURBCCFF1c3rOwVenv1LYlsGX0y9j3Y11QlMLGYYBA0bQRc+g7q76ZYHLYG1o\nLZc4xCb8IUOGCP5fXl6O5ORkfP3113KpXJGohU8IIUQSgY6BaKLdRO5bBldWVWLRhUVYfX01/tv7\nvzBuZiz0Oo/HAw88wW3zmmsEVJPn6oViB+3FxMQI/t+sWTN06tQJzZo1E/0GJRE3eOHtW6BlS6Co\nCGjaVImBEUIIUTu7bu/C6vjVuDDxglzOl1qUirGHxqKpdlNEBUfJpZVOm+fUw8sL+O03oHdvJQVF\nCCFELVVWVaL97+0xznMcTPVMZTpXcXkxVsWtwuwes/FVj6/kdptAYQnfwMBAaBre+5W+evWqwZXK\ngyRf+OLFQFoasGaNkoIihBCitmKexeDgvYMyn4fH42Gcxzh0tu4sh6iEz0stfBFSUwEfHyArixbg\nIYQQot4Uvj1utby8PKSlpQke6qBNG6BDB+DkSa4jIYQQQrglNuFHR0fDzs4ODg4O8Pb2hr29PYKC\ngpQRm1yMHQts3851FIQQQgi3xCb8ZcuW4dKlS3ByckJeXh62b98OPz8/ZcQmF6GhwKlTwMuXXEdC\nCCGEcEdsws/JyUHr1q2hr6+P169fY8yYMTh37pwyYpMLExOgXz/goOzjMAghhBC1JTbhm5iYoLi4\nGEFBQQgNDcWkSZPg6uqqjNjkJjycuvUJIYQ0bmJH6b9+/RrNmjWDtrY2YmJikJmZieHDh0Nfn9ut\nBKUZrVhWBlhbA7dvAzY2Cg6MEEIIUQCFT8tbtmwZwsLCYCOHTFlcXIxx48bh6dOncHR0xF9//QUD\ng9rLBtrb26NFixbQ1tZGkyZNcO3atdqBS/mFT54MuLoCs2bJ9CUQQgghnFD4tLzi4mIEBgaiV69e\n+OOPP5Cbm9vgyhYuXIgePXogKSkJ3bp1w48//lhnOR6Ph5iYGCQmJtaZ7Bti7Fhgxw65nIoQQghR\nOxIvvHPr1i3s3bsX+/fvh62tLaKjo6WuzMXFBefPn4eFhQVycnLg7++P+/fv1yrn4OCA+Ph4mJqK\nXt5Q2iudqiqgdWt2xL6bm9ShE0IIIZyStYUv8fa4rVq1gqWlJUxNTfH8+fMGVZabmwsLCwsAgIWF\nhcjeAh6Ph379+kFLSwvTpk1DZGRkneXmz58v+L+/vz/8/f1F1q2lBYwZw7b0HRzqj3PSJGCw/DZM\nIoQQQqQWExMjtIGdrMS28FevXo29e/ciLy8PI0eOxEcffYQOHTqILB8QEICcnJxax3/66SdMmDAB\nL168EBwzMTFBYWFhrbLZ2dmwsrLCvXv3EBQUhKioKPR+bwechlzpvH4NnD4N1Pe206fZf9eulerU\nhBBCiEIpvIWfnp6OlStXolOnThKd8MyZMyJfq+7Kt7S0RHZ2Nlq1alVnOSsrKwCAq6srgoODce3a\ntVoJvyH09YHg4PrL8PnA7t0yV0UIIYSoFLGD9hYvXixxshdn6NCh2LZtGwBg27ZtGD58eK0yb968\nQXFxMQDg+fPnOH78ONzd3eVSvyTMzYG8PKVVRwghhCiFUnfLEzUtLysrC5GRkTh27BiePn2KkJAQ\nAICpqSlGjRqFqVOn1g5cxq4NUZKTgZAQoI6xhIQQQghnaHtcOXv+HHBxAQoK5H5qQgghpMEUPg9/\nzpw5Eh3TFCYm7EY7lZVcR0IIIYTIj9iEf7p62HoN9Q3MU3fa2oCxMZCfz3UkhBBCiPyITPhr1qyB\nu7s7Hjx4AHd3d8HD0tISPXr0UGaMSmduznbtE0IIIZpC5LS8MWPGYODAgfjmm2+wdOlSwX0DQ0PD\nelfA0wStWlHCJ4QQollEJvyWLVuiZcuW2P3/k9Lz8vJQVlaG169f4/Xr12jdurXSglQ2auETQgjR\nNGIX3tm9eze+/fZbaGtrQ1dXV3D89u3bCg2MS5TwCSGEaBqxCX/RokU4d+4c7OzslBGPSqDFdwgh\nhGgasaP0TU1NYWhoqIxYVAbdwyeEEKJpxLbwXVxc4Ofnh2HDhsHIyAgAO/l/5syZCg+OK9SlTwgh\nRNOITfgWFhYICQkBj8dDSUmJMmLiHCV8QgghmkbipXVLS0uhp6en6HgkpqildQHgzh1g1Ch2XX1C\nCCFEFSh8ad1r165h0KBB6NChAwDg1q1bmDZtWoMrVAc0aI8QQoimEZvwIyMjsXTpUsH9e09PT5w/\nf17hgXHJ1BQoKgL4fK4jIYQQQuRDbMKvqKhAx44dBc/Ly8tRWlqq0KC4pqMDGBnRjnmEEEI0h9hB\ne+3bt8eRI0cAAGlpafj999/h6+ur8MC4Vj1wr1UrriMhhBBCZCe2hb9hwwbs2bMHGRkZ8Pb2RlZW\nFlavXq2M2DhF9/EJIYRoErEtfHNzc+zcuVMZsagUWnyHEEKIJhGZ8JcuXYo5c+bg888/r/Uaj8fD\nb7/9ptDAuEZz8QkhhGgSkQm/ehqej48PeDye4DjDMELPNRUlfEIIIZpE4oV3VI0iF94BgN9/B+7f\nB/78U2FVEEIIIRJT+MI7AQEBKCoqEjwvLCxE//79G1yhuqBBe4QQQjSJ2ISflZUlWHQHAExMTJCR\nkaHQoFQBDdojhBCiScQmfB6Phxs3bgiex8fHK7QrXVXQPXxCCCGaRGzC//XXX9G/f3/BY8CAAVi+\nfHmDKtu3bx/c3Nygra2NhIQEkeViY2Ph7e0NDw8P/P777w2qS1aU8AkhhGgSiQbtPX/+HFevXgWP\nx0O3bt1gZmbWoMru378PLS0tTJ06FcuWLYO3t3etMnw+H87Ozjh79ixsbGzg6+uLXbt2wdXVVThw\nBQ/aq6gA9PSA8nJAW1th1RBCCCESUdigvXv37gEAbty4gfT0dFhbW8PKygppaWn1ts7r4+Ligvbt\n29db5tq1a3BycoK9vT2aNGmCsLAwwdK+ytSkCdCyJVBYqPSqCSGEELkTOQ9/+fLl2LBhA7766qs6\n592fO3dOIQFlZmbCzs5O8NzW1hZxcXF1lp0/f77g//7+/vD395drLNXd+ubmcj0tIYQQIlZMTAxi\nYmLkdj6RCT8gIAAAsHnzZrRt21biEwYEBCAnJ6fW8UWLFmHIkCFi3y/Noj41E74i0H18QgghXHm/\nIbtgwQKZzicy4f/8888YNWoUQkNDperCP3PmjEwB2djYID09XfA8PT0dtra2Mp2zoSjhE0II0RQi\nE76zszP8/f3x9OnTWi1zHo+Hv//+W6aKRQ086Ny5Mx49eoRnz57B2toae/bswa5du2Sqq6Fo8R1C\nCCGaQmTC37FjB5KSkjBixAjMmjVLKEE3dC39Q4cOYfr06cjPz8egQYPg5eWFEydOICsrC5GRkTh2\n7Bh0dHSwefNmBAcHo7KyEpGRkbVG6CsLLb5DCCFEU4icljdu3Dj89ddf+Pnnn/H1118rOy6xFD0t\nDwB++w14+BD44w+FVkMIIYSIpbBpeQ8ePEBqaiq2b9+OwsLCWo/GgO7hE0II0RQiu/SnT5+OkJAQ\nPHjwAD4+PrVeT0lJUWhgqoDu4RNCCNEUYlfa++STT7B27VplxSMxZXTp37oFhIcDd+4otBpCCCFE\nLIVvj7t27VpkZ2djx44dANhldhtD6x6gQXuEEEI0h9gW/sKFCxEdHY2srCw8fPgQmZmZGDlyJC5f\nvqysGOukjBZ+RQXQvDm7nr6W2EsjQgghRHEU3sKPiorC6dOnoa+vD4BdGOfVq1cNrlCdNGkCGBgA\nL15wHQkhhBAiG7EJv0mTJtCusV1cWloaSktLFRqUKqGBe4QQQjSB2IQ/ceJEjB49GkVFRViwYAGC\ngoLwxRdfKCM2lUD38QkhhGgCkdPyqs2ePRunT59G69atUVxcjBUrVgg21mkMaC4+IYQQTSA24QOA\nh4cHysvLAQCenp4KDUjVUMInhBCiCcR26e/duxceHh5Yu3Yt1q5di44dO2Lfvn3KiE0lUMInhBCi\nCcS28P/73//iwoULcHZ2BgA8fPgQQUFBGDlypMKDUwXm5sCTJ1xHQQghhMhGbAu/qqoKlpaWgucW\nFhaoqqpSaFCqhAbtEUII0QRiW/gjRozAwIEDERoaCoZhcOjQIYSGhiojNpVAXfqEEEI0gdiV9gAg\nKSkJR48eBY/Hw+DBg+Hu7q6M2OqljJX2AODmTWD8eCApSeFVEUIIISLJmvdEJvxHjx4hKysLffr0\nEToeGxsLGxsbODo6NrhSeVBWws/MBHx8gJwchVdFCCGEiKSwpXVnzJiBpk2b1jqup6eHGTNmNLhC\ndWNuzi6t6+4u/DhwgOvICCGEEMmJbOG7ubnh7t27db6pY8eOuMPxnrHKauEDQEoKUFLy7vnu3cDL\nl8AffyilekIIIUTmvCdy0F6TJk1Evqnm2vqNgYOD8PO0NOD337mJhRAASE8HWrQAWrbkOhJCiLoQ\nmfC9vLywfv16TJkyRej4hg0b4OXlpfDAVJmjI83NJ9xhGKBbN6CoCLC1BXx9gS5dACcngMd7V05P\nD+jTR/gYIaTxEtml/+LFCwQFBUFXVxc+Pj4AgBs3bqC8vBzHjx+HiYmJUgN9nzK79N9XXs62rl6/\nBnQkWpyYEPl5/Bjw9weePQOSk4Hr19lHaqpwuUePgI8+An76iYsoCSHyprBR+gDAMAzOnTuHO3fu\ngMfjwc3NDf369WtwZfLEZcIHgDZtgHPngLZtOQuBNFKbNwNnzwI7d9Zf7vlzoGdPYMYMYNo05cRG\nCFEchd3Drz55v379VCbJq5Lqbn1K+ETZLlwAevcWX87cHDh5ki1rZQUEBys+NkKI6hK7tK487du3\nD25ubtDW1kZCQoLIcvb29vDw8ICXlxe6dOmixAglR/fxCVckTfgAe0H6zz/A1KnAxYuKjYsQotqU\negfa3d0dhw4dwtSpU+stx+PxEBMTw/k4gfo4ObH3UglRpqwsdl2IDh0kf4+3N7B9OzBiBLB8OcDn\nA69esY/iYnYQYH1atQJmzpQtbkII95Sa8F1cXCQuy+X9eUk4OgJXr3IdBeFKbi7www9Aaem7Yzwe\nMGsW4OqquHovXAB69QK0pOybCwwE/vwT+Ouvd9P5WrQADA3Fn+uHH4DRo9nbAoQQ9aWSY8yrxw5o\naWlh2rRpiIyMrLPc/PnzBf/39/eHv7+/cgIEdek3ZunpwIcfAgMGsMm32oULwLJlwMaNiqtbmu78\n94WGsg9pXb7M1jtqVMPqJYQ0TExMDGJiYuR2Pok2z5FGQEAAcupYeH7RokUYMmQIAKBv375YtmwZ\nvL296zxHdnY2rKyscO/ePQQFBSEqKgq93/srx/Uo/Vev2BZPSQnNc25Mnjxhk/1nnwFffSX8Wmoq\nu+9CdjZQz7pVMvHwYC8olDm0ZdkydrVJWlmSEG4pdJR+Q5w5c0bmc1j9f9+hq6srgoODce3atVoJ\nn2stWgD6+uymOtTV2Tjcu8d2jf/vf8Ann9R+vU0bdpDc+fPsRYG8vXjBJl5lr3vl5wds26bcOgkh\n8sdZl76oq5Q3b96Az+fD0NAQz58/x/Hjx/Hbb78pOTrJVHfrU8JXb7dvA7du1V+mrAz47jvg55+B\nceNElwsNBfbvV0zCv3QJ6NpVcb0Honh5sYv8FBYCKjyOlhAihlIT/qFDhzB9+nTk5+dj0KBB8PLy\nwokTJ5CVlYXIyEgcO3YMOTk5CAkJAQCYmpriyy+/RGBgoDLDlFh1wq95H5eoh6oq4PhxYMUK4P59\nduU6cbdm1q8H/v+ulEgjRgA9erAD5OS95URsLNvaVjYdHaB7d/Y+/rBhyq+fECIfcr+Hryxc38MH\ngPnz2SlOCxdyGkajVVgIHD4sflrZ+168ADZsYG/JzJzJDkbT1ZVfXF5ewMqV7Dr28tS9O7BoEdC3\nr3zPK4mffmK/38uWKb9uQghL5e7hNyaOjsCJE1xH0Xht28Y+/n+rB4np6rKtdT8/xQy4rO7Wl2fC\nf/0aSEpiu/S54OdHc/EJUXeU8GVAU/O4df488M03QFgY15EICw1lW+GrVkk/X16UuDjA0xNo3lw+\n55OWry87aLG4mJ27TwhRP0pdWlfTODrSantcqapi7ynLu9tcHpydAVNT4MoV+Z1Tlvn38tCsGduT\ncvkydzEQQmRDCV8GrVoBb9+y+5IT5bpzh02qqjpDorpbX15iY7lN+ADbrR8by20MhJCGoy59GfB4\n77r137+PHBcHnDolfKxtW2DsWOXFp8nOn1fN1n210FBg4EB27XpZxwlUVADXrrFb3XLJz49dZpcQ\nop4o4cuoroTPMMDkyWxCMjV9d+yTT9hE0KwZN7FqkvPngeHDuY5CtA4d2Pvt16+zq+K9fs2uY79m\nDdsjZGrKPszM2HXt67soePmSvVg0NlZe/HXp0QNITGT3D9DT4zYWQjRJVRXw++9AVFTtWUcHDwL2\n9vKphxK+jOraNS8xEXjzhv0B1hy0deoU2/JX5ZapOmAYtmt5xQquIxGNx2Mv7tasYX9hN21iE+aK\nFWzyzs8HCgrYx8uX4s/32WeKj1kcfX2gY0f2M6zEbSsI0WgZGUBEBJszVq5kf89qsrSUX12U8GXk\n6Mj+Aaxp2zZg/PjaI7T79gXOnaOEL6t79wADA8DOjutI6hcWxibGcePYnRUdHd+9Jq8rdmWrvo9P\nCZ8Q2e3dC3z+OXtBP3cuu8iVIlHCl5GjI7Bz57vnb98Cu3bVPUK7b19gyRLlxaapVP3+fbWOHdmW\nvCbx82NbIYQQySUksH/7+fx3xwoLgcxM4J9/lLcZFiV8Gb3fpX/yJDstq2ZrrlrPnkB8PN0DlVVs\nLHo7wjsAABXwSURBVNC/P9dRNE49ewKjR7MXtvJcnZAQTZWYyA7gnTtXuFdSWxsICKjdha9IlPBl\nZGvLtuKqk/i2bcCECXWXNTQE3N3Z7l0ulkfVBAzDtvAXLeI6ksbJ2Ji9mE1IALp14zoaQlTbrVts\nsl+9mt1ng2uU8GWkrc1ui5qSAlhYANHRwObNosv7+7P38SnhN8zjx+x9LnW9B64JAgKAAwfUO+G/\nesX2ttVkZgZ4eHATD9E8t28DAwawg7dVIdkDlPDlonrFvXPn2Ku5li1Fl+3blzbbkUX1/XtFrIFP\nJPP550CnTuyyxtXTTtVNeDh7/9TI6N2xhATg4kV27AVpnEpL2a2gZZWby976WrkSGDlS9vPJCyV8\nOXByYufi794NLFhQf9mePd9N2+NqXXR1pi4D9jRZ69ZASAi7V4A6LsTz779AcjL7aNr03fElS4DF\ni4EdO7iLjXCHzwcCA9kLwZqfi4bQ0mIX3froI/nEJi+0Pa4crFrFjrRMTgbS0sRPrejZk70w+PBD\n5cSnKRiGvX0SHQ20a8d1NI3bkyfszn1PntTfo6VqqqrYRbL++9/aLa9Xr9g1Eq5eZS/iSeOyYgVw\n6BAQEyO/Ta/kTda8p6JflnpxdGSTUHi4ZPMo+/ZlP1REOs+escvM0h9j7jk6AkFBwB9/cB2JdP76\nix1cGxpa+7UWLYBp02jqbGP08CHw00/s+CtVTfbyQC18Obh/H3B1ZQdpSHL/7+xZYN484NIlxcem\nSbZuZac97t7NdSQEYD/3fn7A06fsQkjSeP0a+OIL9hw1desG/Pqr/GKs6c0bdsrs3r1A9+51lyko\nYHuPbt1S/YWdiHzw+ezn+KOPgOnTuY6mftTCVwFOTmxLR9LBPj16sH9QSkoUG5emofv3qsXFhe2t\nWrNGuvdlZbF/YCsr2dZ09WPxYvbW2JEjiol3+XI20YtK9gA7CHHyZOCXXxQTA1E9q1axPbOqsHy1\nolELnyN+fsC337KDRBqbigp2edyqqnfHtLQANzd2mmNdnj5l7xlfuwY4OCgnTiJeUhK7CNLTp5It\nJnXzJjB0KPCf/7Cj/N+fbXHuHLuOxd277LoV8pKTw36+rl9n79OLK9uhA/sZtbCQXwxE9Tx4wI6p\nioure7E0VSNr3qOEz5Hvv2cT3+LFXEeiHBUV7OjovXuBw4eBVq2ER8IWFrJzVtetq50EKivZveA/\n+giYMUO5cRPxhg8H+vUT3x169CgwaRLw55/1T1WaOJEdCCjPJXw/+YRd0WzZMsnKf/YZW37p0vrL\n5ecDz5/XX0Zbm71NQFNJG6aiovYGZZLg89mNqQoL2UdBAVBWJlzmwAF245rPP5dLqApHCV9NnTvH\nLrV49SrXkYh29Ciwfr3s5+Hz2Svodu3YP/ShoezUrpqKi9mkHh4OzJ4t/Nq8eez36cQJzR5Qo65u\n3GB7rMT9OhoZsaOgu3atv1xBAdsaP3oU6NxZ9vgOHGDHC9y+LfkWw2lpgJcX8OgRYGLy7virV+zS\nztHR7AVsaipgZVX/uQoK2N/1L79s+NfQGL19y47bWbSI7XJv0kS69/N47GfO1JT9GZqY1O6FMjNj\nGxHq8neFEr6aKi0FzM2B7Gz5dl3K07Bh7LgEcX+gJeHpyU6pq09GBjto67ff2HneADuwccQIdu0C\ncX9YCXfKy4U3BqmLrq7ku4FFRbEt/GvXZNtB7OpVYMgQdmtqb2/p3jt9OrtKWs2Wua4uOwbngw/Y\nabU+PuLju3OHLfv4sfSDG9VZYiKwZYv4z4WeHtud3q4dOx6qVSt2ifIlS9jB0PPm1T/uojGhhK/G\n+vVjV2SquQBPt27sHxmuVVWxV7937gDW1sqr98YNtmv/+HGgfXt2RbdVq9j7vqTxYBh2Cd+gIGDm\nzIad48kToFcvYONGYNCghsXw/p8YHq9hXfPh4WyvxX//K/171c2rV+wty1272IummqsZ1qWkhP1Z\nPXrEXhRlZbHjQubNk09jQ5OoVcKfPXs2jh49Cj09Pfj5+WHBggVoWceqHbGxsZgxYwYqKysRGRmJ\nz+u4waIJCT83l+06rMYwbGv2n3/YRMel27fZVvajR8qv+8gRdj60tze7OZG0o8CJZnj0iG3ZBQQI\nH7eyYsdzdOkiOvkWFrLv/eIL9rPEtYcP2cFhjx6JT4DqimGA/fvZWxf9+7PjH8zMpD9PZaXi94VX\nV2qV8M+cOYMPPvgAADB16lSYmppiyXurXPD5fDg7O+Ps2bOwsbGBr68vdu3aBVdXV+HANSDh1+Xn\nn9lk+9df3Mbx55/s2uKbNnFT/6pVbHfg5cu0BHFjdv167QFbDx8CO3eyiWH0aGDUKMDS8t3rfD57\nQdC1q2pNr5s8GbCxUc/liMUpLmZnVzx8CKxdy/asEPlTq4Rf0/79+3H48GFs375d6PiVK1ewYMEC\nnDx5EgAEFwTffPONUDlNTfhFRey0oaQktnXLlY8+YrtTRW31qwwMQyObSd0Yhr1HvHMnOxCwuFj4\n9YED2QtGVRqM9ewZOwjx/v2GtXxVVWoqO06iWzd2PRJdXa4j0lyy5j3OOk42bNiAjz/+uNbxzMxM\n2NVY4srW1hZxcXF1nmP+/PmC//v7+8Pf31/eYSqdkREwfjx7H1/clCBFYRh2JDLXS4xSsiei8Hjs\nLR9vb8WtzCdv9vbshfTPP7MPSVVWAhs2sD0D+fnCr/Xpw94ClGQNBEW4fJmddTNnDnu/nn5n5Ssm\nJgYxclyHXe4t/ICAAOTk5NQ6vmjRIgwZMgQA8NNPPyExMRH79++vVe7AgQM4efIkNmzYAADYvn07\n4uLi8Pt7I9k0tYUPACkpgK8v+y8XI/gfPWIHFKal0S8wIfKUlQW4u7ODYcXNOmEYdmri11+ztwJ+\n+UV4Nc+qKvY2QVER28sh7bQ1aVRU1L7YOHWKjW3bNrZHhSieyrXwz5w5U+/rW7duxfHjxxEdHV3n\n6zY2NkhPTxc8T09Phy2XfdsccHBgE+7mzeygI2kwDDtvVZKBbp9+ys4Pfl9sLDuvmpI9IfJlbc0u\n9BIRIX6qWWwsO7B3+XJ25kpdv49btrADfSdMYMf9iFqpUhZ8Pjtw8u5d4cF0lpbseiJubvKvkyiG\nUu/hnzx5El999RViY2NhampaZ5nKyko4OzsjOjoa1tbW6NKlS6MatFctLg4IC2Nb25KOWOXz2RXC\nrl5l723W1zuQl8f+Eqem1p4bPGECO9d46tSGx08IqVtRETsotqKi/nJt2wJjxoj//S8tZVvYrq7A\n6tXyv1D/9Vf2tsG5c6o1JqIxUqtBe+3atcPbt29h8v9LV3Xv3h2rV69GVlYWIiMjcezYMQDA+fPn\nhablTa9jzU5NT/gAO9J1+nR2FLI4ZWXsXN/q7r0WLcS/JziYbTm8n9gdHNh58O9dYxFCVNSrV+xi\nQH36sDMXanJ1bfhMlzt32A2SaA8L1aBWCV+eGkPCP3SIHTh39Wr9V+1FReyqeFZW7P20mmvU1yc6\nml1WMinp/9q7+5gqyzcO4N9zDMsN1oLxYqCoUbx00POQIS8SJSVbA2SSDTdlWLIKW5nFYtP+yqDN\nPwxXOnvBmGZUTothoi48xkDkaKYNIWgeBAGRxjxIJIRcvz+ecQwRiB/nxXPO9/Pf83Kec3ONs+u5\nn+e+r/v29Vtb1eph167xkT6RM/nzT3XAb1fX7X29vep4oP37p369wUF1auPrr6tjBcjxmPBd2K1b\n6vuxgYGJH6Vdvw6sWQNs3z61R24i6qpgu3er7+wB4KuvgIMH1frjROTcenrUVwOXL6sLEk3Fe++p\nqxuWlfHm/15xzw3aI+uZMUMtPHLt2sTneXgAc+ZM/Uep0agD9z7++HbCHxmwR0TOz9tbHQB84MDU\neum1tepUwF9/ZbJ3Jezhu7neXnV+8G+/qVN/wsPVGtiOLu1LRNbx/ffqQkTjTec+eVL9zf/b0aPq\nYL2MDJs3j6aAj/Rp2jZsUCt/bdgAhIaq7wJtMb2HiOxvcFC9mTca1Zv7O4899hiwfv3o6n/+/uqg\nXrq38JE+TVturrp8Z1iYusAHkz2R65g5U53ps28fsGXL6GOff67+7u/cT66JPXwCoL7na2kBXn1V\nrZ5FRK7j9Gl1BH9j4+138v396hr0P/yg1vine9908x7LKBAAdeqNycQBe0SuKDpanZVTV3d7386d\n6rQ7Jnv3wR4+AVAX6Hj3XaCwkKtdEbmirVuBq1fVWTm9vWrvvrKSpXGdCQftERHRpFpa1CI87e1q\nQa+mJvW9PjkPDtojIqJJzZunFtratw/YsUOda0/uhe/wiYjcRFaWOitn5UogJMTRrSF74yN9IiI3\nYTarA/WOHQPmznV0a2iq+A6fiIj+MxGWy3VWnJZHRET/GZO9+2LCJyIicgNM+ERERG6ACZ+IiMgN\nMOETERG5ASZ8IiIiN8CET0RE5AaY8ImIiNwAEz4REZEbYMKncRkMBkc3weUxxvbBONseY3zvs2vC\nz8vLQ3h4OKKiorBx40aYzea7njdv3jwsXLgQiqIgOjrank2kf+EP2PYYY/tgnG2PMb732TXhL1++\nHPX19Thz5gz++usvFBYW3vU8jUYDg8GAc+fOoa6uzp5NJCIickl2TfjPPfcctFottFotkpOTceXK\nlXHP5cI4RERE1uOw1fKSk5Oxfv16rFq1asyxBQsWwMvLC1qtFrm5ucjJyRlzjoYrQBARkZuZTsq+\nz4rtAKD24q9evTpmf0FBAVJTUwEAH3zwAby8vO6a7AGguroas2fPRkNDA55//nmEhYUhISFh1Dl8\nAkBERPTf2b2H/+WXX+Kzzz7DTz/9hAceeGDS8zdt2oTAwEC8/fbbdmgdERGRa7LrO/yKigps27YN\nZWVl4yb7/v5+3LhxAwDQ3d2NH3/8EZGRkfZsJhERkcuxaw//0UcfxeDgILy9vQEAsbGx2LlzJzo6\nOpCTk4PDhw/j0qVLWLlyJQDAx8cHL774Il555RV7NZGIiMg1iRM6efKkKIoikZGRsmPHDkc3xyW0\ntrbK008/LREREZKYmCh79uwREZHe3l5ZsWKFREZGSnp6uty4ccOxDXUBQ0NDotfrJSUlRUQYY2vr\n6+uTrKws0ev1Eh4eLrW1tYyxDXz66acSGxsrUVFR8uabb4oI/5ena926deLn5yc6nc6yb6KYFhUV\nSWRkpCiKIlVVVZNe3+kq7d26dQsvvfQSDh48iLNnz+KLL75AQ0ODo5vl9Dw8PLB9+3bU19fjwIED\nyM/PR0NDA95//33ExcXhwoULiImJwdatWx3dVKdXVFSEiIgIy0wTxti6cnNzkZiYiHPnzuHChQsI\nCwtjjK2sp6cHBQUFOH78OIxGI5qamnD06FHGeZrWrVuHioqKUfvGi+nFixdRXFyMs2fP4uDBg8jO\nzsbw8PDEX2CT2xQbqqmpkeTkZMt2YWGhFBYWOrBFriklJUWOHz8uoaGhcvXqVRER6ezslNDQUAe3\nzLm1tbVJUlKSVFZWWnr4jLH1XL9+XebPnz9mP2NsXf39/RIcHCzt7e3S19cniYmJUltbyzhbgclk\nGtXDHy+mBQUF8uGHH1rOS05OllOnTk14bafr4be3t2POnDmW7aCgILS3tzuwRa7njz/+QH19PWJi\nYtDV1QV/f38AgL+/P7q6uhzcOuf21ltvYdu2bdBqb//0GGPrMZlM8PX1RXZ2NnQ6HXJyctDf388Y\nW9msWbOwa9cuzJs3DwEBAYiPj8eSJUsYZxsYL6YdHR0ICgqynPdfcqHTJXwW3LGtvr4+ZGZmYvv2\n7fD09Bx1TKPRMP7TUF5eDj8/PyiKMm4dCcZ4eoaGhmA0GpGRkQGj0YiBgQF89913o85hjKevu7sb\nr732Gi5evIiWlhacOnUK5eXlo85hnK1vsphOFm+nS/iBgYFoa2uzbLe1tY26y6H/3z///IOMjAys\nWbMGK1asAKDeUY4UUurs7ISfn58jm+jUampqUFZWhvnz52P16tWorKzE2rVrGWMrCgoKgo+PD1JT\nUzFr1iysXr0aFRUVCAgIYIytqK6uDjExMQgJCYGPjw9WrVqFqqoq/i/bwHgxvTMXXrlyBYGBgRNe\ny+kS/uLFi9Hc3IyWlhYMDg7im2++QVpamqOb5fREBC+//DIef/xxbNy40bI/LS0NJSUlAICSkhKk\np6c7qolOr6CgAG1tbTCZTCgtLcWyZcuwd+9extiKAgICEBISgtOnT2N4eBiHDx9GUlISUlNTGWMr\nSkhIwJkzZ9DT04OBgQEcOXIEy5cv5/+yDYwX07S0NJSWlmJwcBAmkwnNzc2Try5r9REHdmAwGESv\n14tOp5OioiJHN8clVFVViUajkUWLFolerxe9Xi9HjhzhNBsbMRgMkpqaKiKcymRtv//+uyxZskQe\neeQRSU9Pl76+PsbYBvbs2SNPPfWULF68WLZs2SK3bt1inKcpMzNTZs+eLTNnzpSgoCApLi6eMKYf\nffSR6HQ60ev18vPPP096fYctnkNERET243SP9ImIiGjqmPCJiIjcABM+ERGRG2DCJyIicgNM+EQu\nbqSA0uXLl/H1119b9doFBQWjtuPj4616fSKyHiZ8Ihc3Un3LZDJh//79U/rs0NDQhMcLCwtHbVdX\nV0+tcURkN0z4RG4iPz8fVVVVUBQFRUVFGB4eRl5eHiIiIhASEoLi4mIAgMFgwDPPPIOMjAwsXLgQ\nAJCeno4nnngCy5Yts5Sqzc/Px99//w1FUbB27VoAGFWOOS8vDwsWLEBISAi+/fZby7WTkpKQmZmJ\niIgIbN682Z4hIHJvtishQET3Ak9PTxFRi/2MrNAnIrJ7927ZtGmTDA8PS19fnyiKIgMDA3LixAnR\narXyyy+/WM7t6ekRERGz2Sxz5861bI9c+87vOn/+vMTGxsr169elsbFRfH19pbOzU06cOCEeHh7S\n2NgoN2/eFJ1OJ62trTb9+4lIxR4+kZuQO2psHTt2DOXl5YiKisLSpUthNptRW1sLANDr9VAUxXJu\naWkpkpKSEB8fD7PZjObm5gm/q7y8HC+88AIefPBBhIaG4sknn4TRaIRGo0F0dDRCQ0Nx//33Iy4u\nDjU1Ndb/Y4lojPsc3QAicpzNmzcjKytr1D6DwYCHH37Ysn3p0iXs2rULBoMB3t7eUBQFN2/enPC6\nGo1mzA3GyFiChx56yLJv5syZk16LiKyDPXwiN+Hl5YXu7m7LdnJyMkpKSiz7mpqa0N/fP+ZzHR0d\n8PX1hbe3N6qrq3H+/HnLsRkzZtz1MykpKTh06BDMZjOamppgNBoRHR097rLARGR77OETubiRnvWi\nRYsQHBwMRVGQnZ2NN954Ay0tLYiOjoanpyf8/Pxw6NChMWtuL126FMHBwQgPD4dOp8Ozzz5rOfbO\nO+8gISEBERER2Lt3r+VzkZGRiIuLg6Io0Gq1+OSTT+Dn54eGhoYxa3ZzzXQi++DiOURERG6Aj/SJ\niIjcABM+ERGRG2DCJyIicgNM+ERERG6ACZ+IiMgNMOETERG5gf8BCTVIC7QQW9IAAAAASUVORK5C\nYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109ec2810>"
       ]
      }
     ],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "weighted_edge_coeffs = np.array(edge_coeffs[1:]) * np.array(probs[1:])\n",
      "print np.sum(weighted_edge_coeffs)/np.sum(probs[1:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-2.00472506195\n"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "weighted_tri_coeffs = np.array(triangle_coeffs[1:]) * np.array(probs[1:])\n",
      "print np.sum(weighted_tri_coeffs)/np.sum(probs[1:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0.0171693500275\n"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(12,4))\n",
      "ax1 = fig.add_subplot(121)\n",
      "ax1.hist(edge_coeffs[1:], weights=-np.log(probs[1:]))\n",
      "ax1.set_title(\"Edge coefficient\")\n",
      "ax1.set_xlabel(\"Coefficient value\")\n",
      "\n",
      "ax2 = fig.add_subplot(122)\n",
      "ax2.hist(triangle_coeffs[1:], weights=-np.log(probs[1:]))\n",
      "ax2.set_title(\"Triangle coefficient\")\n",
      "ax2.set_xlabel(\"Coefficient value\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 79,
       "text": [
        "<matplotlib.text.Text at 0x10a74d050>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEVCAYAAAD94xMaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VOW9//HPhAQtTUKAmAETNEBUSEjIJALBIzhQwJZr\nED2HoDFcpOtQq0dp0bpWqYlYxbaUCq1aLTYRpBylCpSCWnsY5KINIgglWC4GSEIS7iYBuSR5fn/w\nc7pDbpOQyeTyfq2VtZLZl+f77D3znW9mnv1smzHGCAAAAIAkyc/XAQAAAAAtCQUyAAAAYEGBDAAA\nAFhQIAMAAAAWFMgAAACABQUyAAAAYEGBjGt2+PBh+fn5qbKy0teheOT48eOaNm2abrjhBs2dO1eS\n9NOf/lS9evVSUlKStmzZor59+9a7nzfffFN33323t8MFAI9s3rzZo9x1rVpizv/Nb36jvn376sYb\nb5Qkvf/++xo2bJiCg4O1a9cu9e/fXx999FG9+wkKCtLhw4e9HC1aAxvzIKMmkZGROn78uDp06OB+\nbPr06Vq8eHG1dQ8fPqzevXurvLxcfn4t/3+u+fPna+vWrVq7dq06duyozZs3a9KkSfrXv/6lbt26\n+SQmPz8/HTx4UL179/ZJ+wBapsDAQNlsNknSuXPndP3117vz8quvvqqUlJRmj6ml5fyjR4+qd+/e\nOnDggHr16iVJ6tOnj+bNm6dp06b5JKZp06apZ8+emj9/vk/ax7Xz93UAaJlsNpvWrVunESNG+DqU\nJnfkyBElJiaqY8eO7r9vu+02nxXH3+B/VQBXKysrc//eq1cvLV26tMa8XF5eLn//9vmWfvToUdnt\ndndxbIzR0aNH9R//8R8+jgytme//9UOrU1lZqR//+McKCwuTw+Go9rVVbm6uhg0bptDQUE2aNEk/\n/vGPlZqa6l7+ySef6I477lBISIji4+O1adOmWtvKy8vTPffco7CwMIWGhuqRRx5xx/Dss88qMjJS\ndrtdaWlpKikpqbeNadOm6Y033tCvfvUrBQUF6dVXX9WsWbP0j3/8Q0FBQcrIyJDL5VLPnj3rjSEz\nM1NDhw51r/fFF19o1KhR6tq1q/r27au3337bvWzatGmaM2eOpkyZoq5duyopKUlffvmlJGnYsGGS\npAEDBigoKKjKdgBQE5fLpYiICL3yyiu69dZbNXPmzGq5a8GCBYqKilJQUJBiYmK0evVq97LMzEzd\neeedeuaZZxQeHq7evXvrvffecy+vL49bffXVV5o5c6Z69Oih8PBwzZs3r9bhF5WVlXruuecUFRWl\n4OBg3X777crPz5ckbdu2TQMHDlRISIgGDRqkjz/+uN42PvzwQ40ePVpFRUUKCgrS1KlTFRwcrIqK\nCjkcDt1yyy2Srnwr+ve//12SVFFRUS2GgoICSVe+zfsmN1+8eFE//vGPdfPNN8tut2v27Nm6cOFC\nleP/6quvqnfv3rrxxhuVmZkp6con+ytWrNAvfvELBQUFaeLEiZ6fWLQcBqhBZGSk+fDDD2tc9vLL\nL5tbbrnF5Ofnm9OnT5thw4YZPz8/U1FRYYwxJikpyTz++OPm8uXLZsuWLSYwMNCkpqYaY4zJz883\nwcHB5o9//KMpKSkxWVlZpnPnzubEiRPV2ikvLzdxcXHm0UcfNefPnzcXLlwwW7duNcYYs3TpUtOr\nVy+Tm5trysrKzMSJE+tt4+TJk8YYY6ZNm2bmzZvnbiczM9Pceeed7r83btxoIiIi6o3hj3/8o3u7\nsrIyEx4ebp577jlz5swZs27dOhMUFGRycnKMMcakpaWZ4OBg884775gTJ06YcePGmSlTprjbtNls\n5tChQw09TQDakcjISPP3v//dGHMlT/n7+5sZM2aYwsJC8/XXX1fJXcYY8/bbb5vCwkJz/vx58+tf\n/9p06tTJFBUVGWOu5K+OHTuajIwMc/r0afP000+bG2+80b1tXXk8NzfX2Gw2d85PTk42U6ZMMbm5\nuWbXrl2mf//+5ve//32NffjFL35h+vbta/bv32+MMWb37t3m1KlT5tSpUyYkJMRkZWWZiooKs3z5\nctOlSxdz+vTpettwuVxV+m1M9ZxqPXZXx/D555+bU6dOVdvuscceM9/5znfMP//5T3Pw4EHjdDrN\nU0895T7+AQEBZvbs2eb48ePmtddeM506dTJnz541xlR/n0HrQ4GMGt18880mMDDQhISEuH/+8Ic/\nGGOMGT58eJUX/tKlS93J8siRI9US09ChQ92JdcGCBWbUqFFV2oqPjzdZWVnVYti2bZvp0KGDu7C1\nGjFihDtRGWPM3/72NxMQEGDKy8vrbWPatGnmpz/9qXuZtdA1pmqBXFcM1u1WrlxpbrvttirLJ06c\naDIyMowxVwrk8ePHu5f96U9/Mn379nX/TYEMoD5XF8g2m80cPXrUvfzqAtmqsrLS9OzZ06xdu9YY\ncyV/denSxV3kHjt2zNhsNlNcXFxvHrcWyEVFRcbf379KHIsWLTLDhw+vMY5bb73VLFq0qNrjb7zx\nRpUC3RhjwsPDTVZWVr1t1NTvugrk2mKwbldZWWm+/e1vuz8QMcaYd9991/Tq1cvdZocOHdwf7ly+\nfNkEBgaaf/zjH8aY6u8zaH3a54Al1Mtms2nNmjU1jnUrLCxUfHy8+2+Hw+H+/dixYwoMDKxysVlC\nQoJOnz4t6cp4382bN6tLly7u5eXl5SoqKqrWTl5enm699dYaxwYXFhYqMTHR/XdiYqLKy8tVXFzc\noDbqU1cMVkeOHFFubm6VNisqKhQeHi7pyvG0HrPu3btXGVsIAA1lt9urDKm42htvvKFf//rXOnLk\niKQr45lPnjzpXh4TE+O+yK5Hjx7udY4fP14tjycmJlbZ9htHjhxRZWWl4uLi3I9VVlbqpptuqjGm\n/Pz8GscGHzt2rMp7iST30IejR482qI361BaD1YkTJ3T+/HmNHTvW/ZgxpsrQkR49eig0NFSS5O/v\nr9DQUPJ6G0KBjAbr0aOHdu7cqXvuuUeS9Nlnn1VZVlZWpi+//NKdXD/77DNFRkZKkm666SY5nU5t\n2LCh3nZ69uyp/fv369SpU9UK1BtvvFGffvqpJk+eLEn69NNP5e/vr+7duzeojWuJ4er1+vTpo5yc\nnFrXMfVchFffcgCwquuivCNHjmjGjBn6y1/+ou9+97uSrlzk50meqSmP79ixw53HrXr27Ck/Pz/l\n5OS4i+y69OzZU1u2bNHAgQOrPB4eHl7lvUS6ktfvueeeBrfR2BisQkND9a1vfUvvv/++Bg0a1OA2\nbDYbOb2V4yI91Kq2F/d//ud/6q233lJBQYHOnDmj5cuXu5fdfPPNGjx4sObPn6+ioiKtWLFCu3bt\nci9/4IEH9PHHH+uNN97QmTNndOHCBblcLvcFElaDBw9W//79NWfOHB08eFAXLlzQtm3bJEkpKSl6\n6623dPjwYZWVlemll17SlClT5OfnV28bDUlagwYNqjUGq3HjxqmsrEy/+tWvVFRUpMuXL2v79u36\n4osvPGrTbrfr008/9TguAKjLddddp+uuu07dunXTV199peeff17Hjh3zaNv68rhVjx49NHbsWD35\n5JPat2+fKisrdejQoVrnHH7ooYf0+9//Xhs2bFB5ebl2796t06dPa8yYMfr666+1bNkylZeXa8WK\nFTp//rzGjRun7t27N6iN+tQWg5Wfn59mzZqln/3sZ/rss89UWVmpgoICffDBBx61YbfbtXv3bpWX\nlzcqRvgeBTJqNX78eAUFBbl/vvm0dtasWRo7dqzi4+M1fPhwzZgxwz1Pp3TlBhqHDh1SdHS0Vq1a\npfvvv989pVpERIQ++OADvf766+rTp49uuukmLVy4sMYrnv38/PSXv/xFJSUlGjx4sHr27Km33npL\nkjRjxgxNmzZNw4YNU+/evRUUFKQlS5bU2cY3RarNZqsS79V/f/OYJHXo0KHWGKzbBQUF6cMPP5TL\n5VJsbKx69Oihp556SpcuXaq3DUlKT0/X//zP/ygkJESrVq1q0HkC0D5dnVOsj3Xv3l3PP/+8UlNT\nNWDAAF26dEl33nlnlfXqykl15fGr133jjTcUGBio0aNHq0uXLrrvvvtqHdI2Z84cpaam6uGHH1aX\nLl00a9YsXbhwQV27dtW6dev04osvKjQ0VIsWLdK6devUtWtXj9qoqy+exnD1di+88IISExN17733\nKiQkRKNGjdL+/fs9amPmzJnKy8tTt27d3N+2onXhRiHwuoEDB2rcuHF6+umnfR0KAKARyONobzz6\nBPmb+QTHjx8vSSotLVVycrLi4uI0adKkKoPSFy9erLi4OCUkJGjLli3eiRot2qeffqpDhw7p3Llz\nevnll7Vr1y4lJyf7OiygzcvLy9Pw4cMVExMjp9PpnpeVnI2GIo+jvfOoQH7xxRcVHR3t/jph/vz5\nuuOOO7R7924lJSXp2WeflSTl5OTo9ddf144dO/TOO+9o2rRpLepe7WgeRUVFGj58uMLCwvTiiy/q\nD3/4gwYMGODrsIA2LyAgQIsWLdLevXu1atUq/eQnP9G+ffvI2Wgw8jjau3oL5Pz8fK1fv14PPfSQ\newzn2rVrlZaWJklKS0tz351nzZo1SklJUUBAgCIjIxUVFaXs7Gwvho+WaNy4cTp69KjOnTunL774\nwv1cAeBd3bt3d08nGBoaqoEDB6qgoICcjQYjj6O9q3eat8cff1y//OUvq9zGt7i4WHa7XdKVKzWL\ni4slXZnHMCkpyb1eREREtdkJ6hrUDgCtQWu4dOPgwYPau3evkpKSyNkA2rXG5Ow6P0Fet26dwsLC\n5HA4at15TVfCXr38aubKHfxa/c/TTz/t8xjoS9vsB31puT+tQVlZmaZMmaJFixYpMDCwyrK2nrNb\n03OtNcXa2uJtTbG2tnhbU6zGND5n1/kJ8rZt27R27VqtX79eFy5cUElJiVJTU2W321VUVKTu3bur\nsLBQYWFhkq5M9J2Xl+fePj8/330nMQCA912+fFmTJ0/WAw88oIkTJ0oSORsAGqjOT5Cfe+455eXl\nKTc3VytXrtSIESO0bNkyTZgwQVlZWZKkrKws95WtEyZM0MqVK3Xp0iXl5ubqwIEDjboDTUsXHNxV\nNptNGRkZ7k9jmvInOLirr7sIoBUyxmjmzJmKiYnRY4895n68vedsAGioBt1q+puv3ubNm6fU1FTF\nxcWpT58+WrZsmSQpOjpa06dPV2Jiovz9/ZWZmdkmx6+Vlp6RZCS5JDm9sP/mP2ZOp7PZ2/SGttIP\nib6g4bZu3arly5crLi5ODodDkvT888+3q5zdmp5rrSlWqXXF25pilVpXvK0p1mvR7DcKaQv3J7/y\nBuLNPrT+YwS0VW0hhzVEe+svgLalsTmMW00DAAAAFhTIAAAAgAUFMgAAAGBBgQwAAABYUCADAAAA\nFhTIAAAAgAUFMgAAAGBBgQwAAABYUCADAAAAFhTIAAC0EsHBXWWz2ZrlJzi4q6+7C/gMt5puBG41\nDbRfbSGHNUR7629L5/33nyqtce7R6nGraQAAAKAJUCADAAAAFhTIAAAAgAUFMgAAAGBBgQwAAABY\nUCADAAAAFhTIAAAAgAUFMgAAAGBRZ4F84cIFDR48WPHx8UpKStKiRYskSenp6YqIiJDD4ZDD4dCG\nDRvc2yxevFhxcXFKSEjQli1bvBs9AAAA0MTqvZPe+fPn1alTJ128eFGJiYl69913tWLFCgUFBWnO\nnDlV1s3JydHUqVO1fft2FRQUaOTIkdq/f7/8/P5dh7eFuzJxJz2g/WoLOawh2lt/WzrupAc0jNfu\npNepUydJUllZmcrLy3XddddJUo2NrVmzRikpKQoICFBkZKSioqKUnZ3d4KAAAAAAX/Gvb4XKyko5\nHA7t3btXv/nNb3TTTTdJkpYsWaKlS5dqyJAhWrhwoUJCQnTs2DElJSW5t42IiFBBQUG1faanp7t/\ndzqdcjqd194TAPACl8sll8vl6zAAAM2o3iEW3zh8+LDGjBmjN998U+Hh4brhhhtUUlKiuXPnqqKi\nQkuXLtUjjzyipKQk3X///ZKkhx56SGPGjNE999zz7wbbwNd1DLEA2q+2kMMaor31t6VjiAXQMF4b\nYvGNyMhIjRkzRps2bVJYWJhsNps6d+6shx9+2D2MIjw8XHl5ee5t8vPzFR4e3uCgAAAAAF+ps0A+\nefKkzp49K0k6deqUNmzYoNjYWBUVFUmSysvLtWLFCsXGxkqSJkyYoJUrV+rSpUvKzc3VgQMHNGjQ\nIC93AQAAAGg6dY5BLiwsVFpamioqKtS9e3fNmTNH3/nOd/Tggw9q165d6tixo4YNG+ae/i06OlrT\np09XYmKi/P39lZmZ+f+/DgIAAABaB4/HIDdZg21gPBtjkIH2qy3ksIZob/1t6RiDDDSM18cgAwAA\nAO0BBTIAAABgQYEMAAAAWFAgAwAAABYUyAAAAIAFBTIAAABgQYEMAAAAWFAgAwAAABYUyAAAAIAF\nBTIAAABgQYEMAAAAWFAgAwAAABYUyAAAAIAFBTIAAGg3goO7ymazNctPcHBXX3cXjWQzxphmbdBm\nUzM32eRsNpskb/ah9R8joK1qCzmsIdpbf1s677//VGmtTZ57jmH70tgcxifIAAAAgAUFMgAAAGBB\ngQwAAABYUCADAAAAFnUWyBcuXNDgwYMVHx+vpKQkLVq0SJJUWlqq5ORkxcXFadKkSSorK3Nvs3jx\nYsXFxSkhIUFbtmzxbvQAgCpmzJghu92u2NhY92Pp6emKiIiQw+GQw+HQhg0b3MvI2QBQXb2zWJw/\nf16dOnXSxYsXlZiYqHfffVevvfaaQkND9cQTT+iFF17QmTNntGDBAuXk5Gjq1Knavn27CgoKNHLk\nSO3fv19+fv+uw9vCFdHMYgG0Xy09h23evFmBgYF68MEHtWfPHklSRkaGgoKCNGfOnCrrtpec3ZYw\nA8O14xi2L43NYf71rdCpUydJUllZmSoqKnTddddp7dq12rRpkyQpLS1NTqdTCxYs0Jo1a5SSkqKA\ngABFRkYqKipK2dnZSkpKqrLP9PR09+9Op1NOp7PBgQNAc3C5XHK5XL4Ow2NDhw7V4cOHqz1e0xsE\nORtAW9NUObveArmyslIOh0N79+7Vb37zG910000qLi6W3W6XJNntdhUXF0uSjh07ViWxRkREqKCg\noNo+rckWAFqyqwvCjIwM3wVzDZYsWaKlS5dqyJAhWrhwoUJCQsjZANqcpsrZ9V6k5+fnp88//1wH\nDx7USy+9pJ07d1ZZ/s3dYmpT1zIAgPfNnj1bubm5+vjjj9WhQwf96Ec/qnVdcjYANGAWi8jISI0Z\nM0abNm2S3W5XUVGRJKmwsFBhYWGSpPDwcOXl5bm3yc/PV3h4eBOHDABoiLCwMNlsNnXu3FkPP/yw\nsrOzJZGzAaA2dRbIJ0+e1NmzZyVJp06d0oYNGxQbG6sJEyYoKytLkpSVlaXk5GRJ0oQJE7Ry5Upd\nunRJubm5OnDggAYNGuTlLgAA6lJYWChJKi8v14oVK9wzXJCzAaBmdY5BLiwsVFpamioqKtS9e3fN\nmTNH3/nOdzRo0CClpqYqLi5Offr00bJlyyRJ0dHRmj59uhITE+Xv76/MzEy+rgOAZpSSkqJNmzbp\n5MmT6tmzpzIyMuRyubRr1y517NhRw4YNc0/ZSc4GgJrVO81bkzfYBqYMYpo3oP1qCzmsIdpbf1s6\npii7dhzD9qWxOYw76QEAAAAWFMgAAACABQUyAAAAYEGBDAAAAFhQIAMAAAAWFMgAAACABQUyAAAA\nYEGBDAAAAFhQIAMAAAAWFMgAAACABQUyAAAAYEGBDAAAAFhQIAMAAAAWFMgAAACABQUyAAAAYEGB\nDAAAAFhQIAMAAAAWFMgAAACARZ0Fcl5enoYPH66YmBg5nU5lZmZKktLT0xURESGHwyGHw6ENGza4\nt1m8eLHi4uKUkJCgLVu2eDV4AAAAoKnZjDGmtoVFRUUqKipSfHy8Tp48qf79+2vjxo166623FBQU\npDlz5lRZPycnR1OnTtX27dtVUFCgkSNHav/+/fLz+3cdbrPZVEeTrYLNZpPkzT60/mMEtFVtIYc1\nRHvrb0vn/fefKq21yXPPMWxfGpvD6vwEuXv37oqPj5ckhYaGauDAgSooKJCkGhtbs2aNUlJSFBAQ\noMjISEVFRSk7O7vBQQEAAAC+4u/pigcPHtTevXs1ZMgQbd26VUuWLNHSpUs1ZMgQLVy4UCEhITp2\n7JiSkpLc20RERLgLaqv09HT3706nU06n85o6AQDe4nK55HK5fB0GAKAZ1TnE4htlZWVyOp2aN2+e\nJk6cqOPHj+uGG25QSUmJ5s6dq4qKCi1dulSPPPKIkpKSdP/990uSHnroIY0ZM0b33HPPvxtsA1/X\nMcQCaL/aQg5riPbW35aO4QHXjmPYvnhliIUkXb58WZMnT9YDDzygiRMnSpLCwsJks9nUuXNnPfzw\nw+5hFOHh4crLy3Nvm5+fr/Dw8AYHBQAAAPhKnQWyMUYzZ85UTEyMHnvsMffjhYWFkqTy8nKtWLFC\nsbGxkqQJEyZo5cqVunTpknJzc3XgwAENGjTIi+EDAAAATavOMchbt27V8uXLFRcXJ4fDIUl67rnn\n9Kc//Um7du1Sx44dNWzYMC1atEiSFB0drenTpysxMVH+/v7KzMz8/19lAAAAAK2DR2OQm7TBNjCe\njTHIQPvVFnJYQ7S3/rZ0jJ+9dhzD9sVrY5ABAACA9oQCGQAAALCgQAYAAAAsKJABAAAACwpkAAAA\nwIICGQAAALCgQAYAAAAsKJABAAAACwpkAAAAwIICGQAAALCgQAYAAAAsKJABAAAAC39fBwAAQGsX\nHNxVpaVnfB0GgCZCgQwAwDW6UhybZmjJ1gxtAGCIBQAAAGBBgQwAAABYUCADAAAAFhTIAAAAgAUF\nMgAAAGBRZ4Gcl5en4cOHKyYmRk6nU5mZmZKk0tJSJScnKy4uTpMmTVJZWZl7m8WLFysuLk4JCQna\nsmWLV4MHAFQ1Y8YM2e12xcbGuh8jZwNAw9iMMbXOS1NUVKSioiLFx8fr5MmT6t+/vzZu3Kg//vGP\nCg0N1RNPPKEXXnhBZ86c0YIFC5STk6OpU6dq+/btKigo0MiRI7V//375+f27DrfZbKqjyVbBZrPJ\nu9P5tP5jBLRVLT2Hbd68WYGBgXrwwQe1Z88eSdITTzzRrnN2c/D++4K7pWZq50pbbfHcN9+5ktrq\nMWxNGpvD6vwEuXv37oqPj5ckhYaGauDAgSooKNDatWuVlpYmSUpLS9Pq1aslSWvWrFFKSooCAgIU\nGRmpqKgoZWdnNzgoAEDjDB06VF26dKnyGDkbLV1wcFfZbLZm+QE84fGNQg4ePKi9e/cqKSlJxcXF\nstvtkiS73a7i4mJJ0rFjx5SUlOTeJiIiQgUFBdX2lZ6e7v7d6XTK6XQ2MnwA8C6XyyWXy+XrMK4J\nORstXfPdaEXiZittW1PlbI8K5LKyMk2ZMkWLFi1SYGBglWX1/UdW0zJrsgWAluzqgjAjI8N3wTQB\ncjaAtqypcna9s1hcvnxZkydP1gMPPKCJEydKuvIJRFFRkSSpsLBQYWFhkqTw8HDl5eW5t83Pz1d4\neHijAgMANA1yNgA0TJ0FsjFGM2fOVExMjB577DH34xMmTFBWVpYkKSsrS8nJye7HV65cqUuXLik3\nN1cHDhzQoEGDvBg+AKA+5GwAaJg6h1hs3bpVy5cvV1xcnBwOhyTp+eef17x585Samqq4uDj16dNH\ny5YtkyRFR0dr+vTpSkxMlL+/vzIzMxkQDwDNKCUlRZs2bdKpU6fUs2dPPfPMM+RsAGigOqd580qD\nbWDKIKZ5A9qvtpDDGqK99bexmObtGltq5qnX2uIxRM0am8M8nsUCzcnfq5/iBAV1UUnJaa/tHwAA\noDWjQG6RyuXN/25LS/kKFQAAoDb1zmIBAAAAtCcUyAAAAIAFBTIAAABgQYEMAAAAWFAgAwAAABYU\nyAAAAIAFBTIAAABgQYEMAAAAWFAgAwAAABYUyAAAAIAFt5oGAAA18JfNZvN1EIBPUCADAIAalEsy\nzdQWhThaFoZYAAAAABYUyAAAAIAFBTIAAABgQYEMAAAAWFAgAwAAABZ1FsgzZsyQ3W5XbGys+7H0\n9HRFRETI4XDI4XBow4YN7mWLFy9WXFycEhIStGXLFu9FDQAAAHiJzRhT6xwumzdvVmBgoB588EHt\n2bNHkpSRkaGgoCDNmTOnyro5OTmaOnWqtm/froKCAo0cOVL79++Xn1/VGtxms6mOJluFK/NCerMP\n3t9/az8HgK+0hRzWEO2tv43l/fcFd0vN1A5tNVVbvH58q7E5rM55kIcOHarDhw9Xe7ymhtasWaOU\nlBQFBAQoMjJSUVFRys7OVlJSUoODAgC0TWVlZXrvvfeapWh48MGZunCh1OvtAGh7GnWjkCVLlmjp\n0qUaMmSIFi5cqJCQEB07dqxKMRwREaGCgoIat09PT3f/7nQ65XQ6GxMGAHidy+WSy+XydRhtxocf\nfqjU1EcVEPAfXm/rSnHMjS4ANFyDC+TZs2frZz/7mUpKSjR37lz96Ec/0tKlS2tct7ZbVFoLZABo\nya7+Jz4jI8N3wbQBxhh17DhYJSVvN0NrFK0AGqfBs1iEhYXJZrOpc+fOevjhh5WdnS1JCg8PV15e\nnnu9/Px8hYeHN12kAAAAQDNocIFcWFgoSSovL9eKFSvcM1xMmDBBK1eu1KVLl5Sbm6sDBw5o0KBB\nTRstAAAA4GV1DrFISUnRpk2bdPLkSfXs2VMZGRlyuVzatWuXOnbsqGHDhmnRokWSpOjoaE2fPl2J\niYny9/dXZmZmrUMsAAAAgJaqzmnevNJgG5gyiGnegParLeSwhmjq/r777ruaNu0NlZS822T7rF1b\nnDqsLfapbbfVnvJFS9TYHMad9AAAAAALCmQAAADAggIZAAAAsKBABgAAACwokAEAAACLRt1quqWr\nrKzkqlEAAAA0SpsrkC9fvqzOnbvqwoXzvg4FAAAArVCbK5DLy8t1+XK5jKnwYivcAAUAAKCtYgwy\nAAAAYEE5au0gAAAWiklEQVSBDAAAAFhQIAMAAAAWFMgAAACABQUyAAAAYEGBDAAAAFhQIAMAAAAW\nbW4eZHjCXzab9+ZyDgrqopKS017bPwAAgDdRILdL5ZK8dyvu0lJupAIAAFovhlgAAAAAFhTIAAAA\ngEWdBfKMGTNkt9sVGxvrfqy0tFTJycmKi4vTpEmTVFZW5l62ePFixcXFKSEhQVu2bPFe1ACABouM\njFRcXJwcDocGDRokqe6cDgDtVZ0F8vTp0/Xee+9VeWz+/Pm64447tHv3biUlJenZZ5+VJOXk5Oj1\n11/Xjh079M4772jatGmqrKz0XuQAgAax2WxyuVzauXOnsrOzJdWe0wGgPauzQB46dKi6dOlS5bG1\na9cqLS1NkpSWlqbVq1dLktasWaOUlBQFBAQoMjJSUVFR7gQMAGgZjKl6gW5tOR0A2rMGz2JRXFws\nu90uSbLb7SouLpYkHTt2TElJSe71IiIiVFBQUOM+0tPT3b87nU45nc6GhgEAzcLlcsnlcvk6jCZh\ns9k0YsQI+fn56Qc/+IFmzZpVa063ImcDaC2aKmdf0zRvNputzvl0a1tmTbYA0JJdXRBmZGT4Lphr\ntHXrVvXo0UP79u3TmDFj1Ldv3yrLa8vp5GwArUVT5ewGz2Jht9tVVFQkSSosLFRYWJgkKTw8XHl5\nee718vPzFR4e3qigAABNr0ePHpKkfv36adKkScrOzq41pwNoCv7ufzy9+RMc3NXXHW1zGlwgT5gw\nQVlZWZKkrKwsJScnux9fuXKlLl26pNzcXB04cMB9lTQAwLfOnz+v0tJSSdKJEye0fv16xcbG1prT\nATSFb27M5d2f0tIzzdaj9qLOIRYpKSnatGmTTp06pZ49e+qZZ57RvHnzlJqaqri4OPXp00fLli2T\nJEVHR2v69OlKTEyUv7+/MjMzvXo7YwCA54qLizVp0iRJUrdu3fT4449r9OjRGjJkSI05HQDaM5u5\n+pJmbzdos1W7iropff311woO7qry8q+91oZkkzdv1dwW9t/MTyug2Xg7h7U0Td3fd999V9OmvaGS\nknebbJ+183au80VbbbFPtNUU7bSnvNQQjc1h3EkPAAAAsKBABgAAACwokAEAAAALCmQAAADAggIZ\nAAAAsKBABgAAACwokAEAAAALCmQAAADAggIZAAAAsKBABgAAACwokAEAAAALCmQAAADAggIZAAAA\nsKBABgAAACwokAEAAAALCmQAAADAggIZAAAAsKBABgAAACwokAEAAAAL/8ZuGBkZqeDgYHXo0EEB\nAQHKzs5WaWmpUlNT9eWXX6pPnz5atmyZAgMDmzJeAAAAwKsa/QmyzWaTy+XSzp07lZ2dLUmaP3++\n7rjjDu3evVtJSUl69tlnmyxQAAAAoDlc0xALY0yVv9euXau0tDRJUlpamlavXn0tuwcAAACaXaOH\nWNhsNo0YMUJ+fn76wQ9+oFmzZqm4uFh2u12SZLfbVVxcXOO26enp7t+dTqecTmdjwwAAr3K5XHK5\nXL4OAwBahODgriotPdMsbQUFdVFJyelmaetqjS6Qt27dqh49emjfvn0aM2aM+vbtW2W5zWaTzWar\ncVtrgQwALdnV/8RnZGT4LhgA8LErxbGpd72maavmOrI5NHqIRY8ePSRJ/fr106RJk5SdnS273a6i\noiJJUmFhocLCwpomSgAAAKCZNKpAPn/+vEpLSyVJJ06c0Pr16xUbG6sJEyYoKytLkpSVlaXk5OSm\nixQAAABoBo0aYlFcXKxJkyZJkrp166bHH39co0eP1pAhQ5Samqq4uDj3NG8AAABAa9KoArlXr17a\ntWtXtceDgoKYuQIAAACtGnfSAwAAACwokAEAAAALCmQAAADAggIZAAAAsKBABgAAACwokAEAAAAL\nCmQAAADAggIZAAAAsKBABppZcHBX2Ww2r/4EB3f1dTcBAGi1GnUnPQCNV1p6RpLxchs2r+4fAIC2\njE+QAQAAAAsKZAAAAMCCAhkAAACwoEAGAAAALLhIDwAAoFXzl83GxdlNiQIZAACgVSuXt2dH+rf2\nUYhTIANtkrc/TQiQdNlrew8K6qKSktNe2z8AAHWhQAbaJG9/mmDz6v6ZxxkA4EtcpHdNXL4OoAm5\nfB1Ak3C5XL4OoQm5fB1AE3L5OgC0Gy5fB9AALl8H0EAuXwfQAC5fB9BALl8HgKt4pUD+6KOPlJCQ\noLi4OC1ZssQbTbQQLl8H0IRcvg6gSVAgt1QuXweAOrStnO3ydQAN4PJ1AA3k8nUADeDydQAN5PJ1\nALhKkw+xqKio0IwZM/Thhx8qPDxcAwcO1MiRI9WvX7+mbgoAcI3I2QBQXZN/gpydna2oqChFRkYq\nICBAU6ZM0Zo1a5q6GQBAEyBnA0B1Tf4JckFBgXr27On+OyIiQv/4xz+qrNM8c/V5u41v9p/h5f17\nS037b7q++HI+xowMb52TpuTp8bmWvvjiOVSXhvWFOT2bh+9ytjfPr/W51pzPo8a01ZjXuC/75M38\n2tT9qivWlvi8uNZj2xL71AQt+ei9oMkL5Po6YkxzzdMHAKgPORsAqmvyIRbh4eHKy8tz/52Xl6eI\niIimbgYA0ATI2QBQXZMXyLfffrsOHDigw4cP69KlS/rf//1fTZgwoambAQA0AXI2AFTX5EMs/P39\n9frrr2vSpEkqLy/XrFmzuBoaAFoocjYAVOeVeZDvuusu7dy5U3v27FFeXp769eunhIQEPfbYY/rq\nq69q3CYyMlJxcXFyOBwaNGiQN8K6ZnPnzvWoL61hTtG3335bMTEx6tChgz777LNa12sN58XTvrT0\n81JaWqrk5GTFxcVp0qRJKisrq3G9lnxOPDnGTz31lOLi4pSUlKQvvviimSP0XH19cblc6ty5sxwO\nhxwOh5599lkfRNk0rDn70Ucf9fi5KF2ZJs7hcGj8+PHNGPG/eRJrXl6ehg8frpiYGDmdTmVmZjZ7\nnK3ptVFfrG+++aYGDBigAQMGaOrUqfrnP//pgyiv8DSvb9++Xf7+/nrnnXeaMbrqPIl3+/btuvPO\nOzVgwAA5nc7mDdCivli//vprpaWlyeFw6K677vLp7DczZsyQ3W5XbGxsres0+PVlvOyDDz4wFRUV\npqKiwjz00EPmySefrHG9yMhIc+rUKW+Hc0086Ut5ebnp06ePyc3NNZcuXTIDBgwwOTk5Poi2bvv2\n7TP/+te/jNPpNDt27Kh1vdZwXjzpS2s4L3PnzjUvvPCCMcaYBQsWtLrXiifH+K9//av53ve+Z4wx\n5pNPPjGDBw/2Raj18qQvGzduNOPHj/dRhN7l6XPRGGMWLlxopk6d6rNj4UmshYWFZufOncYYY06c\nOGHsdnuzvv5b02vDk1i3bdtmzp49a4wxJjMzs0XH+s16w4cPN2PHjjWrVq3yQaT/jqO+eM+cOWOi\no6NNXl6eMebK89UXPIn15ZdfNrNnzzbGGHP48GHTu3dvU1lZ6YtwzUcffWQ+++wz079//xqXN+b1\n5fVbTY8aNUp+fn7y8/PT3Xffrfz8/LqKdW+Hc0086UtrmVO0b9++uvXWWz1at6WfF0/60hrOy9q1\na5WWliZJSktL0+rVq2tdtyWeE0+OsbWPgwcP1tmzZ1VcXOyLcOvk6fOlJZ6HpuDpczE/P1/r16/X\nQw895LNj4Ums3bt3V3x8vCQpNDRUAwcO1LFjx5otxtb02vAk1iFDhqhz586SpLFjx9b5vu5Nnr5O\nlyxZonvvvVc33HCDD6L8N0/iXbFihSZPnuy+UDY0NNQXoXoUa+fOnVVaWqrLly/r9OnT6tSpk8+m\nZBs6dKi6dOlS6/LGvL68XiBbvfbaa5o4cWKNy2w2m0aMGCGHw6HXXnutOcNqlNr6UtOcogUFBc0Z\nWpNqbeelNq3hvBQXF8tut0uS7HZ7rS/elnpOPDnGNa3jqzfXunjSF5vNpm3btikmJkZjxoxRTk5O\nc4fpNZ4+Fx9//HH98pe/lJ9fs76VVOFprN84ePCg9u7dq6SkpOYIT1Lrem00NFe++uqrtb6ve5un\nx3XNmjWaPXu2JN/Or+5JvAcOHNDp06c1dOhQORwOvfnmm80dpiTPYk1JSVFFRYVCQ0N15513+ixW\nTzTm9dUkF+mNGjVKRUVF1R5/7rnn3OPSfv7znysoKEj33XdfjfvYunWrevTooX379mnMmDHq27ev\nhg4d2hThNci19qUl3dzAk77UpzWdl7q0lPNSWz9+/vOfV/nbZrPVGnNLOSdX8/QYX/1JY0s5N1ae\nxJSQkKC8vDwFBAQoKytLEyZM0MGDB5shuqZxrc/FdevWKSwsTA6HQy6Xy1thSmqa140klZWVacqU\nKVq0aJG+/e1vN3mctWlNr42GtLlx40YtX75c27Zt82JEtfMk1scee0wLFiyQzWaTMcan3/p4Eu/l\ny5flcrn04Ycf6vz58xo1apTuuecefetb32qGCP/Nk1h/+9vfyt/fX4WFhdqzZ4/Gjh2rI0eO+PQf\n5ro09PXVJAXy3/72tzqXZ2Zmav369fr73/9e6zo9evSQJPXr10+TJk1Sdna2T970r7UvLWlO0fr6\n4onWcl7q01LOS139sNvtKioqUvfu3VVYWKiwsLAa12sp5+Rqnhzjq9fJz89XeHh4s8XoKU/6EhQU\n5P595syZevLJJ3X69Gl17dq12eK8Ftf6XNy2bZvWrl2r9evX68KFCyopKdGDDz6oN954o8XFKl0p\nPCZPnqwHHnig2T/xbE2vDU9z5e7du/X9739fGzZsUEhISHOG6OZJrDt27NCUKVMkSSdPntSGDRsU\nEBDgk6kMPYm3Z8+e+t73vqfu3btLujIN40cffaS77767xcX60UcfaebMmerUqZMGDx6sG2+8Ufv3\n71ffvn2bNVZPNOr11QRjo+u0YcMGEx0dbU6ePFnrOufOnTMlJSXGGGOOHz9ubrvtNvP+++97O7QG\n86Qvly9fNr179za5ubnm4sWLLfJiMCun02k+/fTTGpe1lvPyjbr60hrOy9y5c82CBQuMMcY8//zz\nNV5s1JLPiSfH2HqhxMcff9xiL9LzpC9FRUXuC1LWrFljwsPDfRGqV3jyXLRyuVxm3LhxzRFaNZ7E\nWllZaVJTU83jjz/e3OEZY1rXa8OTWI8cOWKioqLMJ5984pMYv9HQvD5t2jTz5z//uRkjrMqTePft\n22cGDhxozp07Z06dOmVuueUWU1pa2iJjfeWVV8zDDz9sKioqzKFDh0xUVFSzx2mVm5vr0UV6nr6+\nvF4gR0VFmZtuusnEx8eb+Ph49xWPBQUFZsyYMcYYYw4dOmQGDBhgBgwYYEaMGGFeeeUVb4fVKJ70\nxZgrbxbx8fGmf//+5sUXX/RVuHV65513TEREhLn++uuN3W433/3ud40xrfO8eNIXY1r+eSkpKTET\nJ040sbGxJjk52Z0UW9M5qekYv/LKK1XifPLJJ03//v3N4MGDW9w/KVb19eW3v/2tiYmJMQMGDDCp\nqam1/nPWGnnyXLRyuVw+m8XCk1g3b95sbDabGTBggDt/b9iwoVnjbE2vjfpinTlzpunatav7WA4c\nOLDFxmrl6wLZGM/ifemll0y/fv1Mr169zO9+9ztfhVpvrGfPnjWPPvqocTgcZvTo0eavf/2rz2Kd\nMmWK6dGjhwkICDARERFm6dKl1/z6shnTRi/DBgAAABqhZY6kBgAAAHyEAhkAAACwoEAGAAAALCiQ\nAQAAAAsKZDTImTNnNH36dMXHxysmJkZjx47VgQMHGrWv5cuXa9CgQUpNTZUk3X///br99tv11ltv\nadasWdq3b1+t2/7lL3/RCy+80Kh2v/rqK7388suN2rYmgYGBTbYvAGhK5OzqyNnwBLNYwGPGGN1x\nxx2aPn26vv/970uSPvnkE128eFF33XVXg/eXmJiod955RzfffLN2796tRx991Ot35JKkw4cPa/z4\n8dqzZ0+T7C8oKEilpaVNsi8AaCrk7JqRs+EJPkGGxzZu3KiOHTu6E60kJSUluRPt3Llz1bt3b0VF\nRemtt95yr/PLX/5S0dHRuuWWWzR//nxJ0n//93/rn//8p8aPH69f/OIXmjx5snbt2qWEhAR9+eWX\ncjqd2rFjhyTpvffeU//+/RUTE6NRo0ZJunJHw0ceeUSSdOLECSUnJys6OloDBgxQdna2JCk9PV2z\nZ8/W8OHDFRcXp5UrV0qSfvKTn+jQoUNyOBx68sknq/Txqaee0ksvveT+Oz09XQsXLtS5c+c0cuRI\nJSQkaMyYMTXeSdHlclW57fUPf/hDZWVlSbpyN6ekpCT17dtXo0eP1smTJxtzCgDAY+RscjaugXem\nbEZb9OKLL9Z6J6rPP//cDBkyxJw9e9Z88cUX5oYbbjCFhYXm/fffN/fdd58pLy83Fy9eNHfddZc5\nduyYMcaYyMhIc+rUKWNM9TtxOZ1Os2PHDnP8+HHTrVs3900Yzpw5Y4wxJjMz0/zwhz80xhiTkpJi\nVq9ebYwxZs+ePe6bAzz99NMmNjbWnDlzxhw9etT06dPHGGPM4cOHa73bzs6dO81dd93l/js6Otrk\n5+eb8vJy9x3sjhw5YiIiItzrBAYGGmOM2bhxY5U+/PCHPzRZWVnm0qVLJjY21uTl5RljjPnd737n\nvvMXAHgLOZucjcbz93WBjtbDZrPVumzdunW699571blzZ3Xu3FkDBw7U9u3btWnTJm3fvl233367\nJOncuXP6v//7P91///1Vtjc1jPQxxuiTTz5RYmKiEhMTJUkhISHV1vvb3/6mffv2KT09XZJ09uxZ\nff3117LZbJo4caJCQkIUEhKiDh066Pjx4zW29Y34+HgdP35chYWFOn78uLp06eK+X/uLL76o9evX\n69y5czp9+rSOHz+usLCwOo+ZMUb/+te/dOTIEfcnFRUVFYqMjKxzOwC4VuRscjYajwIZHouJidGq\nVatqXGaz2WpNYtOmTdPTTz/dqDbrSvDfMMZozZo1uummm6otsybnjh076sKFC/Xu77777tOqVatU\nVFSkKVOmSLryVdzmzZv1/vvv69vf/rbCwsKq7ev666/XxYsX3X+fOnXKHX/Xrl21c+fOetsGgKZC\nziZno/EYgwyPjRgxQhcvXtRrr73mfmz79u366KOPNG7cOL377rv66quvtH//fm3fvl2DBw/W3Xff\nrbfffltHjx6VJBUUFOjEiRMetWez2ZSUlKQdO3bo008/lSSdPn1aUtVPL+6++24tWbLEneg+//zz\nOvcbFBRUZwz/9V//pT/96U9atWqV7rvvPnfc4eHhCgoK0sqVK91xWMXHxysnJ0dlZWUqKCjQBx98\nIEm67bbbJEl//vOfZYzR5cuXlZOT49ExAIDGImeTs9F4FMhokPXr12vr1q0aMGCA+vfvr2eeeUbh\n4eGKjY3VHXfcIYfDoTFjxuh3v/udwsLCNGrUKH3/+9/X8OHD1b9/f913330qKyurtl+bzVbjJw+h\noaFavny5HnzwQcXExCglJaXa+osXL9bRo0fVr18/xcTE6Pe//32V/V6tW7duuvfee5WQkFDtgg9J\nio6OVllZmSIiImS32yVJycnJOnv2rPr166ctW7YoOjq6WhvXX3+9nnjiCSUlJWnGjBkaPXq0JCkg\nIECrV6/Wr3/9a912221yOBz6+OOPPT7mANBY5GxyNhqHad4AAAAACz5BBgAAACwokAEAAAALCmQA\nAADAggIZAAAAsKBABgAAACwokAEAAACL/wcG0trFWFvHowAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a725410>"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}