ExponentialExpanderHopfield_v009.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Implementing the network structure from the paper https://arxiv.org/abs/1704.02019\n",
    "# Associative Content Addressable Memory with Exponentially Many Robust Stable States\n",
    "# By Rishidev Chaudhuri and Ila Fiete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Seems to work nicely now!\n",
    "\n",
    "# Added a draft cell to count the number of satisfied states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/adammarblestone/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import pprint\n",
    "import time\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.cm as cm # color maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "global_num_timesteps = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def matching_assignment(N, N_c, z, z_c):\n",
    "    \n",
    "    start = time.time()\n",
    "    \n",
    "    if N_c * z_c != N * z:\n",
    "        print \"MISMATCHED SIZES\"\n",
    "\n",
    "    outgoing_edges_left = []\n",
    "    for i in range(N):\n",
    "        for j in range(z):\n",
    "            outgoing_edges_left.append((i,\"?\"))\n",
    "\n",
    "    incoming_edges_right = []\n",
    "    for k in range(N_c):\n",
    "        for l in range(z_c):\n",
    "            incoming_edges_right.append((\"?\", k))\n",
    "\n",
    "    matches = []\n",
    "    while len(matches) < N*z:\n",
    "        random_index_left = np.random.randint(0, len(outgoing_edges_left), size=1)[0]\n",
    "        random_index_right = np.random.randint(0, len(incoming_edges_right), size=1)[0]\n",
    "        s = (outgoing_edges_left[random_index_left][0], incoming_edges_right[random_index_right][1])\n",
    "        if not s in matches:\n",
    "            matches.append(s)\n",
    "            outgoing_edges_left.pop(random_index_left)\n",
    "            incoming_edges_right.pop(random_index_right)\n",
    "        if time.time() - start > 10: # This is to re-start if it doesn't find a good assignment on this round\n",
    "            print \"Trying another assignment...\"\n",
    "            return matching_assignment(N, N_c, z, z_c)\n",
    "        \n",
    "    input_to_constraint_node_connectivity = []\n",
    "    for n in range(N):\n",
    "        outgoing = []\n",
    "        for o in matches:\n",
    "            if o[0] == n:\n",
    "                outgoing.append(o[1])\n",
    "        input_to_constraint_node_connectivity.append(outgoing)\n",
    "        \n",
    "    return input_to_constraint_node_connectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 60\n",
      "N_c = 25\n",
      "z = 5\n",
      "z_c_avg = 12\n",
      "r = 8\n",
      "M_typical = 16\n",
      "Trying another assignment...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEZCAYAAACXRVJOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGcxJREFUeJzt3XmUbWV9p/Hne7mIgoIXlYuIoJKWQaLIUpuIkTK2A9gC\nEUOcQVuNbVRsEwW1s7iahFZXXAlRo40mBFBjQIKiRkGFcmoHREAFRFEhSOAig4zqQvn1H/stOBQ1\nnBrOUHWfz1pn1Z737z2nzv6ePZx9UlVIkrRm1AVIksaDgSBJAgwESVJjIEiSAANBktQYCJIkwEDQ\nLJKcneRls4x7aJKbkmTYdY2jJPsluWJAy57rddg5yR1JFvw+TnJYkq/09N+R5BELrUGri4GgBauq\nK6pq65rnSyzTNzqr3Ki+0LOU9dYs3QOTZDLJL9sHipuTXDyM9ao/BoIGKQxwQ+MeyrIa1nNZwKvb\nB4r7VdXuQ1qv+mAgbEKSHNo+ld3UHr9KctYcszwsyVfbtJ9Lsm1bzt0OVSQ5PMmP23Q/TvL8JLsB\n7wd+r63z+jbt1klOTHJNkp8meWtPfWuSvDvJz9ty/nTaes5O8letpluBh7d1X9TWfWmSV/Ysb78k\nVyR5Y1vflUkOTrJ/kh8muTbJUQt4/g5IcmFb1xVJ3nD30XlDko1tPYf3jJirzUcnOamnf9bDQO35\n+Zv2/FwKPGueeo9sz8lNSb6f5OB+2zqD30nyzSQ3Jjktyf171rNPkq8luSHJeUn2m2dZBvm4qiof\nm+ADuB9wIfDyWcafDfwI2AXYovUf08btDPyW7gPFlsCNwO+0ceuB3Vv3YcCXpy33ROC0Nt/OwCXA\nS9u4VwHfBx4MbAN8fmo9PTVdBuzW1r0W2B94WBv/+8CtwF6tfz/gduCtwGbAy4GfAx9p698DuA3Y\nuU2/L3D9HM/ZfwJPbN3bzLCeo9t69m91bNNHm48GTuxZx53PbU+bX9bz/FwE7ADcHzird9oZ6j0E\nWN+6/wi4paf/bq8NcAfwiDn+F64AdgfuA3wcOKmNewhwLfCM1v/U1v+AOZa1EbgG+Aqw36jfCz56\nXp9RF+BjBC969wnt08B755jmbOAtPf3/E/j31j09EK4H/hC497RlTN/orAF+DezaM+yVwFmt+4vA\nK3rGPXWGjeOGedp2GvDa1r1f2zCn9d+3bfge1zP9t4ED+3zeLgNeAdxv2vCp9azpGbYReEIfbV5I\nIHwReGXPtE+bKxBmqP884NmzvDbzBcIxPf27A79q/0dvAk6YNv3ngBfPsqzHA1sBmwMvAW4CHj7q\n94SP7uEho03TMXRvyiPmme7qnu7b6Daod1NVtwF/TBcYVyX5VJJdZ1neA+k+1f9Hz7DL6T5lQvfJ\nt/dqnZmu3LnbsHb45+tJrktyA92n8wf2THJdtS0R8Mv295qe8b+cqV2zOITuMM3l7fDVPtPWc0dP\n/9TzNV+bF2L683P5XBMneUk7hHNDe24exd2fm4WYvt7N27J2Bg5Ncn173EC3p/XgmRZSVedU1a1V\ndXtVnQh8DThgkTVpmRkIm5gkz6PbgB9SVb9djmVW1eer6unA9nSHQ46bGjVt0mvpDq3s3DNsZ+DK\n1n0VsGPPuJ1mWt1UR5J70R2+eBfwoKpaB3yWAR2jrqpzq+pg4EHAJ4GT+5htvjbfSreXNWXGDWlz\nFfDQacuZUZKd6F6HV1fVuvbcXMjin5vp672drm1X0O3hbNse66o7WfyuPpdbS6hJy8xA2IQkeSzw\n98DBVXX9UhfXlrldkgOTbEm3kbiF7vADdIdNdkyyOUD7BH0y8NdJ7ptkZ+B/AVMnVU8GjkiyQztp\n+aZ5arhXe1xbVXck2R94+hLbNaMkmyd5QZKtW5DeTHe4Zk6tzacwe5vPB56c7rsd2wBzneQ+GXhd\nkockWQccOce0W9G9Dte2k9EvBfacr945vCjJbu11fhtwStvz+jDw7CRPb+u5dzuZv8P0BSTZpk23\nRZLNkryQ7rzP55ZQl5aRgbBpOZDuZORXc9fVRp+ZZdr5LhedGr8GeAPdJ95rgSfTHT6C7qTnhcDV\nSaYO07yO7nDKT4AvAx+uquPbuA8CZwLfBc4FPgP8pudQzN1qqqpb2vJOaVcxPY/uk3s/dd+jP8mT\nktw0x7wvBn6a5Bd05wFe0Od6Xsssba6qLwD/Stfmc4BPzbGcDwJnABfQnfs4ddaVV10MvBv4Bt2h\nv0cBX+2z3pnGnQScQHdi/V60w41V9TPgIOAtdCfsLwf+nJm3LZsDf0V3yO7nwJ8CB1XVpXOsW0M0\ndbJtcCvoPvV8iO7TyR3Ay4Af0r0JdqY7UXdoVd040EK04iR5JvD+qnr4qGuRNgXD2EM4lu7qlN2B\nxwA/oNst/kJV7Ur3KfLNQ6hDY64dbti/HU54CN0VOP826rqkTcVA9xCSbA2cV1W7TBv+A7rrjzcm\n2R6YrKrdBlaIVoQk9wG+BOxKd/XPp4HXt0NDkgZs0IHwGLorHS6i2zv4NvB64Mp21cPUdNdX1bYD\nK0SSNK9BHzJaC+wNvK+q9qa7xO4o5jixJ0kajbUDXv7PgCuq6tut/1S6QNiYZH3PIaNrZpo5iUEh\nSYtQVQv+fsdA9xCqaiNwRZJHtkFPpbsM8XTg8DbsMOa4VHDUX+Ue5OPoo48eeQ22b3GP9t85xMdw\n3wur+bXbFNq3WIPeQ4DuOvGPtC8n/QR4Kd0NwE5O96MblwOHDqEOSdIcBh4IVXUB3Q2tpvtvg163\nJKl/flN5hCYmJkZdwkCt9vatZqv9tVvt7VusgX9TeSmS1DjXp01XMtAfg5tpjUs6NqxNSxJq3E4q\nS5JWDgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkx\nECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgTA\n2kGvIMllwI3AHcDtVfWEJOuAfwV2Bi4DDq2qGwddiyRpdsPYQ7gDmKiqx1bVE9qwo4AvVNWuwFnA\nm4dQhyRpDsMIhMywnoOAE1r3CcDBQ6hDkjSHYQRCAZ9Pck6Sl7dh66tqI0BVXQ1sN4Q6JElzGPg5\nBGDfqroqyYOAM5NcQhcSvab3S5KGbOCBUFVXtb8/T/IJ4AnAxiTrq2pjku2Ba2abf8OGDXd2T0xM\nMDExMdiCJWmFmZycZHJycsnLSdXgPpwn2RJYU1W3JNkKOBN4G/BU4PqqemeSI4F1VXXUDPPXIOuT\nFisJw92xDb4X1K8kVFUWPN+AA+HhwGl075y1wEeq6h1JtgVOBh4KXE532ekvZpjfQNBYMhA0zsYy\nEJbKQNC4MhA0zhYbCH5TWZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKk\nxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkS\nYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEnAkAIhyZok30lyeutfl+TMJJckOSPJNsOoQ5I0u2Ht\nIRwBXNTTfxTwharaFTgLePOQ6pAkzWLggZBkR+AA4EM9gw8CTmjdJwAHD7oOSdLchrGH8LfAG4Hq\nGba+qjYCVNXVwHZDqEOSNIe1g1x4kmcBG6vq/CQTc0xas43YsGHDnd0TExNMTMy1GEna9ExOTjI5\nObnk5aRq1m3x0heeHAO8CPgNcB/gfsBpwOOAiaramGR74Oyq2n2G+WuQ9UmLlYQ5PscMYo34XlC/\nklBVWeh8Az1kVFVvqaqdquoRwPOAs6rqxcCngMPbZIcBnxxkHZKk+Y3qewjvAJ6W5BLgqa1fkjRC\nAz1ktFQeMtK48pCRxtlYHjKSJK0cBoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUG\ngiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElS01cgJPli\nP8MkSSvX2rlGJrk3sCXwwCTrgLRRWwMPGXBtkqQhmjMQgD8BXg/sAJzLXYFwE/DeAdYlSRqyVNX8\nEyWvrar3DKGe6eutfuqThi0JMMz/zeB7Qf1KQlVl/imnzdfvP1mSJwIPo2evoqpOXOgKF8JA0Lgy\nEDTOFhsI8x0ymlr4ScAuwPnAb9vgAgYaCJKk4ekrEIDHAXv4cV2SVq9+v4fwfWD7hS48yRZJvpnk\nvCQXJjmmDV+X5MwklyQ5I8k2C122JGl59XtS+WxgL+BbwK+nhlfVgX3Mu2VV3ZZkM+BrwJ8BBwLX\nVdW7khwJrKuqo2aY150SjSXPIWicDfQcArBhoQueUlW3tc4t6PZIbgAOAvZrw08AJoF7BIIkaXj6\nCoSq+tJiV5BkDd13GHYBPlBVFyVZX1Ub27KvTrLdYpcvSVoe/V5ldDN37R/fC9gcuLWqtp5v3qq6\nA3hskq2BM5JMcM997Vn3hTds2HBn98TEBBMTE/2ULEmbjMnJSSYnJ5e8nL6/h3DnDN3B04OAfWY6\n7j/PvH8B/BL4H8BEVW1Msj1wdlXtPsP0nkPQWPIcgsbZYs8hLPhup9X5BPCMPop64NQVREnuAzwN\nOA84HTi8TXYY8MmF1iFJWl79HjJ6Tk/vGrrvJfyqj1kfDJzQ9irWACdV1ReTnAecnORlwOXAoQsr\nW5K03Pq97PT4nt7fAJcBH6yqawZU19R6PWSkseQhI42zgd/LaBQMBI0rA0HjbKDnEJLsmOS0JNe0\nx6lJdlx4mZKkcdXvSeXj6U4E79Aen2rDJEmrRL/nEM6vqr3mG7bcPGSkceUhI42zQV92el2SFyXZ\nrD1eBFy30JVJksZXv4HwMrpLQ68GrgKey13fI5AkrQL93tzu7cBhVXUDQJJtgb+hCwpJ0irQ7x7C\no6fCAKCqrgceO5iSJEmj0G8grEmybqqn7SH0u3chSVoB+t2ovxv4epJTWv8fAX89mJIkSaPQ9zeV\nk+wB/EHrPauqLhpYVXet08tONZa87FTjzFtXSENkIGicDe3215Kk1clAkCQBBoIkqTEQJEmAgSBJ\nagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBAw4EJLsmOSs\nJBcm+V6S17Xh65KcmeSSJGck2WaQdUiS5jfQX0xLsj2wfVWdn+S+wLnAQcBLgeuq6l1JjgTWVdVR\nM8zvL6ZpLPmLaRpnY/mLaVV1dVWd37pvAS4GdqQLhRPaZCcABw+yDknS/IZ2DiHJw4C9gG8A66tq\nI3ShAWw3rDokSTNbO4yVtMNFHweOqKpbkkzf9511X3jDhg13dk9MTDAxMTGIEiVpxZqcnGRycnLJ\nyxnoOQSAJGuBTwOfrapj27CLgYmq2tjOM5xdVbvPMK/nEDSWPIegcTaW5xCafwIumgqD5nTg8NZ9\nGPDJIdQhSZrDoK8y2hf4MvA9uo9TBbwF+BZwMvBQ4HLg0Kr6xQzzu4egseQegsbZYvcQBn7IaCkM\nBI0rA0HjbJwPGUmSVgADQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKk\nxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkS\nYCBIkhoDQZIEGAiSpMZAkCQBAw6EJP+YZGOS7/YMW5fkzCSXJDkjyTaDrEGS1J9B7yEcDzxj2rCj\ngC9U1a7AWcCbB1yDJKkPAw2EqvoqcMO0wQcBJ7TuE4CDB1mDJKk/oziHsF1VbQSoqquB7UZQgyRp\nmnE4qVyjLkCSBGtHsM6NSdZX1cYk2wPXzDXxhg0b7uyemJhgYmJisNVJ0gozOTnJ5OTkkpeTqsF+\nQE/yMOBTVfW7rf+dwPVV9c4kRwLrquqoWeatQdcnLUYShrtzG3wvqF9JqKoseL5B/pMl+SgwATwA\n2AgcDXwCOAV4KHA5cGhV/WKW+Q0EjSUDQeNsLANhqQwEjSsDQeNssYEwDieVJUljwECQJAEGgiSp\nMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIE\nGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkScAIAyHJ\nM5P8IMkPkxw5qjokSZ2RBEKSNcB7gWcAjwKen2S3UdQySpOTk6MuYaBWe/tWs9X+2q329i3WqPYQ\nngD8qKour6rbgY8BB42olpFZ7f+Uq719q9lqf+1We/sWa1SB8BDgip7+n7VhkqQRWTvqAsbFlVde\nyate9aqhrW/PPfdkiy22GNr6JGk+qarhrzTZB9hQVc9s/UcBVVXvnDbd8IuTpFWgqrLQeUYVCJsB\nlwBPBa4CvgU8v6ouHnoxkiRgRIeMquq3SV4DnEl3HuMfDQNJGq2R7CFIksbPWHxTOcmaJN9Jcvoc\n0zw+ye1JnjPM2pbDfO1LMpHkvCTfT3L2sOtbqrnal+QBST6b5Pwk30ty+AhKXLQklyW5oL0+35pl\nmr9P8qPWxr2GXeNSzNe+JC9o4y9I8tUkvzuKOhern9evTbcity99/n/2vX0Zl6uMjgAuAraeaWT7\nIts7gDOGWdQymrV9SbYB3gc8vaquTPLAYRe3DOZ6/V4DnF9V+7e2XZLkw1X1m6FWuHh3ABNVdcNM\nI5PsD+xSVf8lyX8FPgDsM8wCl2jO9gE/AZ5cVTcmeSbwQVZX+1b69mW+/88FbV9GvoeQZEfgAOBD\nc0z2WuDjwDVDKWoZ9dG+FwCnVtWVAFV17bBqWw59tO9q4H6t+37AdSsoDADC3O+Tg4ATAarqm8A2\nSdYPo7BlMmf7quobVXVj6/0GK+/7QvO9frCCty/M374FbV9GHgjA3wJvBGY8mZFkB+Dgqno/XeNX\nmjnbBzwS2DbJ2UnOSfLi4ZW2LOZr3weBRyX5T+ACur2JlaSAz7fX5hUzjJ/+JcsrWVkbzfna1+vl\nwGeHUNNymrN9q2D7Mt/rt6Dty0gPGSV5FrCxqs5PMsHML8jfAb03v1sxL1qf7VsL7A38AbAV8PUk\nX6+qS4dX6eL02b43AxdU1VOS7EL3z/voqrplmLUuwb5VdVWSB9HVfnFVfXXURS2jvtqX5CnAS4En\nDb3CpZmvfSt2+9LM174FbV9GvYewL3Bgkp8A/wI8JcmJ06Z5HPCxJD8Fngu8L8mBQ65zsfpp38+A\nM6rqV1V1HfBl4DFDrnOx+mnfvsApAFX1Y+CnwIq5kWFVXdX+/hw4je4+XL2uBB7a079jG7Yi9NE+\nkjwaOA44cK5j8eOoj/at5O1LP+1b2PalqsbiAewHnD7PNMcDzxl1rcvZPrqN4+eBzYAtge8Be4y6\n3mVs37uBo1v3errDK9uOut4+27QlcN/WvRXwNbqTc73THAB8pnXvA3xj1HUvc/t2An4E7DPqegfR\nvmnTr6jtS5+v34K2L+NyldHdJPkTultZHDdt1Kr40kRv+6rqB0nOAL4L/BY4rqouGm2FSzPt9fs/\nwPFJLqDbHX9TVV0/0gL7tx44rd1CZS3wkao6c9rr9+9JDkhyKXAr3WGVlWLe9gF/AWwL/EOSALdX\n1T32IsZUP+3rtdK2L/38fy5o++IX0yRJwOjPIUiSxoSBIEkCDARJUmMgSJIAA0GS1BgIkiTAQNAq\nlOSwJNsvYf5nJ3nTctbUxzpvHub6pJn4PQStOu2e739eVeeOupZ+Jbmpqma8/bs0LO4haOwleUP7\ncZ3vJjmiDds5yfd6pvmzJEcnOYTu/jQfTvejPVu0bxJf3O72eGyST7V51iU5rf3AyP9LsmcbfliS\n97Tu49s8X0ty6dQPqKTzD0kuSnJGks9khh9XSfKIdD8QdE6SLyV5ZBv+sLbOC5L8Zc/0sy43yd5J\nJtuyPpt2m+0kr0tyYbof6PnoYF4FbQoMBI21JHsDhwGPB34PeEWSqZtzTd+9rao6Ffg28IKq2rsN\n/wDwjKp6PPCgnvneBnynqh4DvBU4qXdZPd3bV9W+wLOBd7ZhhwA7VdUewEtabTM5DnhNW/cbgfe3\n4ccC72vrvqpn+hmXm2Qt8B7gkLas44Fj2jxHAntV1V7Aq2apQ5qXgaBx9yTgtOru1ngr8G/A7/cx\n39RtjHcDflxV/9H6/2Xask8CqKqz6e4bf98ZlvWJNs3FwHZtWO9dXDcC9/hpwiRbAU8ETklyHvB/\n6e4/MzX/x1p3bxDNttxdgT3pbnF8Hl2A7dDGXQB8NMkL6e5XIy3KWN7cTurDb+ju4Djl3nNM2+89\n7meb7teLWBZ0H7hu6NlT6VXctRfSzzIDfL/tqUz3LODJwIHAW5PsWVV3LKBOCXAPQePvK8DBSe7d\nPnH/Id093TcCD2rnAbYA/nvPPDdz1+87XwI8PMlOrf+Ppy37RdD9EDnw85r/h3umNt5fAw5px/zX\nAxPTJ6yqm4GfJnnunTN3vy0wNf/zW/cLe2abbbmXtPbu05azNskebdxOVfUl4KjW7pn2cqR5uYeg\nsVZV5yX5Z+Acuk/Ux1XVdwGSvL0N/xlwcc9s/wx8IMltdMfgXw2ckeSWnuUAbAD+qd2a+1a6Y/b3\nKGGW/lPpfoXqQrrfeDgXuJF7emGr5X/Tvd8+Rncr4tfTHeZ5E/DJnulnXG5V3d6C5T3pfjh9M+Dv\nkvyQ7gT61nRhdWxV3TRDHdK8vOxUq16Srdr5B5K8D/hhVR27XMtNsi3wTbqfM1zyD7UParnSfNxD\n0KbgFUkOA+4FfIfu5O5y+HSS+wObA29fxo32oJYrzck9BEkS4EllSVJjIEiSAANBktQYCJIkwECQ\nJDUGgiQJgP8PKYjlZds1vJMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103f89ad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "total outgoing edges\n",
      "300\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEZCAYAAACNebLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHAFJREFUeJzt3Xu0HGWd7vHvEwJBCAnhGiAQiAOooAKOHhCVLR6FdRgB\n0UEFhotnOOANEC9cPJ4EXToDIsqI4FIggyAygGYkKJoRslWQayAkCCKKIJckCAQIIgyE3/njfTtp\nmrf37t27b7t5Pmvtle7quvze6k49VW91VykiMDMzqzWu2wWYmVlvckCYmVmRA8LMzIocEGZmVuSA\nMDOzIgeEmZkVOSCsLkl7SHpgiNfPkfT5TtbUyyTNl/SRNsx3uPdhtqQvNjnvP0naMz+eKenCZmqw\n/uSAsOHU/aFMRHw0Ir483AyqN0LWtE79YGmo5bSkBklTJf1Y0kOSXpS0Vc3rX5X0e0lPSrpT0j+1\nYrk2cg4IG9MkrdHtGmzEXgSuAg6gHDpPA/tExGTgcOBMSbt2rjyrcECMcZIOlLRC0lP571lJ1wwz\nzdqSvibpPknLJf1K0oT6o+t4ScvyHt/hVS+s6tqQtKGkuXl+j0n6ZR7+PWArYG6u7zN5+L6S7pD0\nuKRrJL2mar67SLo170FeKumSquXsIekBSZ+TtAQ4X9L6edmP5GXPlbRF1fzmS/qSpOvyuvpxrvei\nvIwba/dih1h3EyRdKOnR3NYbJW1cNcrWkq7Nbf2ZpA2qph2qzS9KmlFat4Uadpa0INd+CbD2EPXO\nkHR1rveR3OZJjbS1PDudJOkvku6VdFDVC2tJOl3S/ZKWSDq73mcqIh6JiG8DtwAqvH5KRNyTH98E\n/BrYrcmabRQcEGNcRFwaEetFxCRgC+CPwMXDTPY1YGdgV2AD4HOkvbqSqcB6wObAPwPfkjS5MN6n\ngQeADYFNgJNzfYcCfwb+ISImRcTpkrbLNR4DbEzam5wrabykNYEfAefn2n4AvK9Q0/qk4Pk/pM/x\n+cCWedgzwFk103wQODi34++A64HzgCnA74CZlRFzwHyuzvo4DKis6w2Ao4G/Vb3+4TzOxsAEoBKI\ndducp2uo+yavnznABXn5lwHvH2oS4CukdfZaYBowq5FlFUzNy9yctGf/HUnb5tdOJa3XN+R/twD+\nX5PLWUXSq4A3A78d7bxs5BwQfUKSSBvT+RFx7jDjHQEcExFLI7khIp6vM8l/A1+KiJURcRXp8H/7\nwnjPA5sB2+Rxr6tddNXjA4ErI+KaiFgJnE7aC34rKbTWiIiz8nzmADfVzGslMDMino+I5yLi8YiY\nkx//FfgX4B0108yOiPsiYgVp43xPRMyPiBdJG9mdKyNGxHsj4rQ66+N5Ughul9fdbRHxdM1y/hgR\nzwGXAjsN0eZX5TbXrp+h7AaMj4h/y+vnh8DN9UbOtVwdES9ExGPA14E9GlzWy2YHfCGv918BP8nt\nAjgS+FREPJnfg38lheVofRu4LSLmtWBeNkLjhx/FxoivAOsCxw4z3kakPdt7G5zvY3kjWvEMMLEw\n3ldJe6bzJAXw3Yg4tc48NwfurzyJiJD0IGmv80XgoZrxa78985fqQMt7md8A9iIdWQiYKEmx+mqU\ny6qm/1vhealNJd8j7YVfko+kvg+cnDf6AEurxq1eV6U2P0Bq80hsxsvXz/2lEQEkbQKcCbw917IG\n8PgIl1mxPCKerVnu5rmLbR1gQdr/ANLOZ6OhVyTpq8DrgHeOZj7WPB9B9AFJHyJ1oby/akNVz6PA\ns8CrW1lDRDwdEZ+JiFcD+wLHS6r8x67tPnkYmF4zbEvShm8JaQNc+9pLFlfz/NPAtsCbI2J9Vh89\njGoDVZL32r8UETuQ9v7/ATi0gUnrtfnB/PgZ0ka2Ymqd+Szh5aEy1PmTr5BCd4e8bg6h+fUyJYdx\n9XIfJn2mnsnL2CD/rZ9PMjdF0imkwH93zRGadZADYoyTtDPwb8D+ETHsnmHeoz4fOEPSZpLGSdo1\n922Ppo59JFVCZwXwAqkrCNLe+oyq0S8F9pH0znze4TOk0PoN6dzAC5I+LmkNSfsBbxlm8euRjgKe\nyieFZ42mLUORNCBpR0njSN1tz7O6nUOp1+br8+u3AQfl92Nv6ncDVdbPJ/N8DmDo9bNernNFPnH/\n2QZqrUfAKZLWlPR2YB/g0vyZ+i7wjcoJe0lbSHpP3RmlE9iVk+trV5/QlnQSqXvqf0bEE6Oo10bJ\nATH27UvqVrlWq7/N9JNhpvkMsJjUd/0Yqb+40c9CvZOp2wK/kLQCuA74Vu6nhnRO4Av52zvHR8Tv\nSXuyZwF/IW1o3pv7yZ8nff3xn4HlwEHAXOC5IWr6Bmnv+1FSyPy0wZqLJP1U0ol1Xp4KXA48STpx\nOh+4aLjlDNXmPMpxpPdyOWnjOKfOfCrr5wjSe/ePwA+HaM4pwJuAJ0jrsXbckaybJbm+h4ELgaMq\n3zYCTgD+ANwg6QlgHrDdEPP6G/BUXv7vSEcgFV8mHV39oeozXe/9sDZSO28YJGkaqc92U9Jh7nci\n4puSZpJOaj2SRz05In7WtkJsTJN0A3BORFzQ7VrMXknaHRBTgakRsVDSRGABsB+pv3xFRJzRtoXb\nmCXpHcDdpCOCQ4CzgRkRsWzICc2spdr6LaaIWEr+VkdEPC3pLlafYGv5CURbTdIdvPTkpUiH80dF\nxA+6U1XDtif12a9D+rbV+x0OZp3X1iOIlyxI2hoYBHYkfevkcFI/7i3ApyPiyY4UYmZmDenISerc\nvXQ5cGz+ylqly2An0hGGu5rMzHpM248g8qUErgSuiogzC69PB+ZGxBsKr3XqCpZmZn0lIkbdjd+J\nI4jzgTurwyGfvK44ALij3sQR0bd/M2fO7HoNbl9zf/nT2cG/zv5f6Of37pXQvlZp60lqSbuTLpC2\nWNJtpE/6yaQfBO1E+urrfcBR7azDzMxGrt3fYrqOdO2XWv7Ng5lZj/MvqbtoYGCg2yW0Vb+3r5/1\n+3vX7+1rlY59zbUZL70Yp1nvSFct7eRnUy3tW7b+JokYIyepzcxsDHJAmJlZkQPCzMyKHBBmZlbk\ngDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzM\nrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JA\nmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWVFbA0LSNEnXSPqt\npMWSjsnDp0iaJ+luST+XNLmddZiZ2cgpIto3c2kqMDUiFkqaCCwA9gOOAB6LiNMknQBMiYgTC9NH\nO+sza5YkoJOfTeH/C9YoSUSERjufth5BRMTSiFiYHz8N3AVMI4XEBXm0C4D921mHmZmNXMfOQUja\nGtgJuAHYNCKWQQoRYJNO1WFmZo0Z34mF5O6ly4FjI+JpSbXHynWPnWfNmrXq8cDAAAMDA+0o0cxs\nzBocHGRwcLDl823rOQgASeOBK4GrIuLMPOwuYCAiluXzFPMj4rWFaX0OwnqSz0FYLxsT5yCy84E7\nK+GQXQEcnh8fBvy4A3WYmdkItPtbTLsDvwIWk3a3AjgZuAm4FNgSuB84MCKeKEzvIwjrST6CsF7W\nqiOItncxjYYDwnqVA8J62VjqYjIzszHIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIH\nhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZm\nRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPC\nzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggzMytqa0BIOk/SMkmLqobNlPSgpFvz397t\nrMHMzJrT7iOI2cBeheFnRMQu+e9nba7BzMya0NaAiIhrgeWFl9TO5ZqZ2eh16xzEJyQtlHSupMld\nqsHMzIbQjYA4G5gRETsBS4EzulCDmZkNY3ynFxgRf6l6+l1g7lDjz5o1a9XjgYEBBgYG2lKXmdlY\nNTg4yODgYMvnq4ho+UxfsgBpa2BuRLw+P58aEUvz408Bb46Ig+pMG+2uz6wZkoBOfjaF/y9YoyQR\nEaM+19vWIwhJFwMDwIaS/gzMBN4paSfgReA+4Kh21mBmZs1p+xHEaPgIwnqVjyCsl7XqCMK/pDYz\nsyIHhJmZFTkgzMysyAFhZmZFDggzMytyQJiZWZEDwszMihoKCElXNzLMzMz6x5C/pJa0NrAOsJGk\nKay+TPckYIs212ZmZl003KU2jgKOAzYHFrA6IJ4CzmpjXWZm1mUNXWpD0icj4psdqKd2ub7UhvUk\nX2rDelmrLrXR8LWYJL0V2Jqqo46I+N5oCxhmmQ4I60kOCOtlHb2aq6QLgVcDC4GVeXAAbQ0IMzPr\nnkYv9/33wOu8O29m9srR6O8g7gCmtrMQMzPrLY0eQWwE3CnpJuC5ysCI2LctVZmZWdc1GhCz2lmE\nmZn1Ht9RzqwJ/haT9bJOf4tpBav/N6wFrAn8NSImjbYAMzPrTQ0FRESsV3mstOu0H7Bru4oyM7Pu\na7qLSdJtEbFzi+upXYa7mKwnuYvJelmnu5gOqHo6jvS7iGdHu3AzM+tdjX6L6b1Vj18A7iN1M5mZ\nWZ/yt5jMmuAuJutlrepiavSGQdMkzZH0SP77oaRpo124mZn1rkYvtTEbuIJ0X4jNgbl5mJmZ9alG\n7wexMCJ2Gm5Yq7mLyXqVu5isl3W0iwl4TNIhktbIf4cAj4124WZm1rsaDYiPAAcCS4ElwAeAw9tU\nk5mZ9YBGv+b6ReCwiFgOIGkD4HRScJiZWR9q9AjiDZVwAIiIx4G2/orazMy6q9GAGCdpSuVJPoJo\n9OjDzMzGoEY38l8Drpd0WX7+j8CX21OSmZn1goZ/SS3pdcCe+ek1EXFn26pavUx/zdV6kr/mar2s\nVV9z9aU2zJrggLBe1unfQZiZ2StMWwNC0nmSlklaVDVsiqR5ku6W9HNJk9tZg5mZNafdRxCzgb1q\nhp0I/CIitgeuAU5qcw1mZtaEtgZERFwLLK8ZvB9wQX58AbB/O2swM7PmdOMcxCYRsQwgIpYCm3Sh\nBjMzG0YvnKT2VzPMzHpQN34NvUzSphGxTNJU4JGhRp41a9aqxwMDAwwMDLS3OjOzMWZwcJDBwcGW\nz7ftv4OQtDUwNyJen5+fCjweEadKOgGYEhEn1pnWv4OwnuTfQVgvGxM/lJN0MTAAbAgsA2YC/wlc\nBmwJ3A8cGBFP1JneAWE9yQFhvWxMBMRoOSCsVzkgrJf5l9RmZtZWDggzMytyQJiZWZEDwszMihwQ\nZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkV\nOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMysyAFhZmZFDggz\nMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCYmVmRA8LMzIrG\nd2vBku4DngReBJ6PiLd0qxYzM3u5rgUEKRgGImJ5F2swM7M6utnFpC4v38zMhtDNDXQA/yXpZklH\ndrEOMzMr6GYX0+4RsUTSxqSguCsirq0dadasWaseDwwMMDAw0LkKzczGgMHBQQYHB1s+X0VEy2c6\n4iKkmcCKiDijZnj0Qn1mtSSRDoI7tkT8f8EaJYmI0Gjn05UuJknrSJqYH68LvAe4oxu1mJlZWbe6\nmDYF5kiKXMP3I2Jel2oxM7OCnuhiqsddTNar3MVkvWxMdzGZmVnvc0CYmVmRA8LMzIocEGZmVuSA\nMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKHBBmZlbkgDAzsyIHhJmZFTkgzMys\nyAFhZmZFDggzMytyQJiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzMrMgBYWZmRQ4IMzMrckCY\nmVmRA8LMzIocEGZmVuSAMDOzIgeEmZkVOSDMzKzIAWFmZkUOCDMzK3JAmJlZkQPCzMyKuhYQkvaW\n9DtJv5d0QrfqMDOzsq4EhKRxwFnAXsAOwIclvaYbtXTT4OBgt0toq35vXz/r9/eu39vXKt06gngL\ncE9E3B8RzwOXAPt1qZau6fcPab+3r5/1+3vX7+1rlW4FxBbAA1XPH8zDzMysR4zvdgG94qGHHuLo\no4/u2PJ23HFHJkyY0LHlmZmNlCKi8wuVdgVmRcTe+fmJQETEqTXjdb44M7M+EBEa7Ty6FRBrAHcD\n7wKWADcBH46IuzpejJmZFXWliykiVkr6BDCPdB7kPIeDmVlv6coRhJmZ9b5u/Q7iPEnLJC2qGvYB\nSXdIWilpl2GmHyfpVklXtL/akRtN+yRNlnSZpLsk/VbS/+hM1Y0ZZdtOym1aJOn7ktbqTNWNq9O+\n0/L7sVDSDyVNqjNtz//4s9n2SZom6Zr8/i2WdExnK2/MaN6/PO5Y3LY0+vkc8balW19znU36kVy1\nxcD7gF82MP2xwJ2tLqqFRtO+M4GfRsRrgTcCvdb11lTbJE0HjgR2jog3kLo3P9SuIkeh1L55wA4R\nsRNwD3BS7URj6MefTbUPeAE4PiJ2AHYDPt5n7asYi9uWRts34m1LVwIiIq4FltcMuzsi7gGGPPMu\naRrwv4Bz21fh6DTbvpz8b4+I2XmaFyLiqbYWO0KjeO+eAv4bWFfSeGAd4OG2FdqkOu37RUS8mJ/e\nAEwrTDomfvzZbPsiYmlELMyPnyZtXHrut0ujeP/G8rZl2PY1u20Zixfr+zrwWaAfT55sAzwqaXY+\nzP2OpFd1u6hWiIjlwNeAPwMPAU9ExC+6W1VTPgJcVRjeLz/+rNe+VSRtDewE3NiBelptqPb1w7al\nXvua2raMqYCQtA+wLO/JiGGONsag8cAuwLciYhfgGeDE7pbUGpJmAJ8CpgObAxMlHdTdqkZG0ueB\n5yPi4m7X0g6NtE/SROBy4Nh8JDFmDNW+fti2DPP+NbVtGVMBAewO7CvpXuAHwDslfa/LNbXSg8AD\nEXFLfn456U3tB38PXBcRj0fESuBHwFu7XFPDJB1O6n6oF2oPAVtVPZ+Wh40JDbSP3DV4OXBhRPy4\nQ6W1RAPtG9Pblgba19S2pZsBMVRKF4dHxMkRsVVEzCCd4LwmIg5tV4Gj1Ez7lgEPSNouD3oXvXnC\nbMRtI/0wcldJa0sSqW29dgK+4iXtk7Q3qeth34h4rs40NwN/J2l6/nbWh4Ce/CYMzbUP4Hzgzog4\ns831jdaI2zeWty0Ntq+5bUtEdPwPuJh0gvI5Up/0EcD+pD7cv5F+XX1VHncz4MrCPPYAruhG/e1s\nH+nbBTcDC0l72ZO73Z4Wtu2zwG+BRcAFwJrdbk+D7bsHuB+4Nf+dXad9e5OC8B7gxG63pZXtI+1h\nr8yfy9vyeHt3uz2tfP+q5jHWti2Nfj5HvG3xD+XMzKxorJ2DMDOzDnFAmJlZkQPCzMyKHBBmZlbk\ngDAzsyIHhJmZFTkg+pyka7u03FMk7dnB5e0haW4b579fL129NLd3tyame5OkbwwzzmRJH22+uiHn\nPV3S4nbM21rPAdHnIuJtXVruzIi4ptOLbeO89yddxrtXDFDnUiVKt/QtiogFEXHcMPOeAnys+dKG\n5R9fjREOiD4naUX+dw9J86tuGHJh1ThvlnRdvuHIDZLWlTRB0vlKN/dZIGkgj3uYpDmS5km6V9In\nJH06XyHyN5LWz+PNlnRAfvwnSbPyfG6v/Nxf0kZ5PoslfVfSfZI2KLTh3Xnet0j6D0nr5OF757bc\nAhxQNX7d+Uo6WNKNud5zlIzL9S7K9R1bs/zdgH2B0/J020h6o6TrtfomLZMLdW8i6Ud5nNsk7ZqH\nH59rW1RZVt6zvlPpKpt3SPqZpAn5tWOUbvCyUNLFSvfWOBo4Lteze67/HEk3AKfm9/Q3eZ1fK2nb\nqs/B3Px4ptINaOZL+oPSbYAB/gWYked9ak2bhqqzuE6UjloWSroN+HjVvMYp3ezmxvz6keVPsXVN\nt3867r+2/zT/qfzvHqTryG9Guo7Lb0h7oGsCfwR2yeNNBNYAjgfOzcO2J/2Ufy3gMOD3pPs5bAQ8\nCRyZxzsDOCY/ng0ckB//CfhYfvxR4Dv58TeBE/LjvUiXctigpv4NSTcielV+/jng/wITSJcamJGH\n/wf58gj15gu8hnR9pDXya98CDiFdtGxe1TInFdbjqvbk57cDb8uPTwG+Xpjmkqr1IWC9vKzbgbWB\ndYE7SJdAmE66X8brq9pzUH78EPmyJJXagJmkG/hU13dF1fOJwLj8+F3A5VWfgyuq5nEt6UqfGwKP\n5vd+OrCozudpqDpr18kZVcN3z49Pq8ybdAOpk/PjtUiXgZje7f8z/lv95yOIV5abImJJpP+RC4Gt\nSRv/hyPiVkg3g4l0tdW3ARflYXcD9wGVC33Nj4hnIuJRUuhcmYcvzvMsmZP/XVA1zttIG1Ei4ufU\n3Agl2xV4HXBd3gM9lLSReg1wb0Tcm8e7qGqaevN9F2kDfXOe157ADOBeYBtJZ0raC1hRpw3Aqpuv\nTI508xZI15V6e2HUPYFzch0REStybXMi4tmI+CvpmjiVaf8UEZX++er1dDtwsaSDSWFXz2VVj9cH\nLlfq7/86aR2W/CTSzWMeA5YBmw4x/4qX1VlnnbwjH0VMjojr8vALq+bzHuDQ/F7cSArxbRtYvnXI\n+G4XYB1VfaXHlax+/xu59n31ONXziarnL1L/M1UZZ+UQ45TqEGnv/uCXDJTeWGf8oeYr4IKI+PzL\nRkjz2ws4CjgQ+N8NznsoI+1rr31/1s6P9wHeQerm+rykHetM/9eqx18iXZH0gNwlNb+BZQ71/jVS\n50iv8CvgkxHxXw0s07rARxD9b7iN6N3AVElvgnRDGKWTnL8GDs7DtgO2zOO20nXAB/My3kPa6611\nA7C7pFfn8dbJ/em/A6ZL2iaP9+EG5ns18AFJG+fXpkjaStKGpG6nOcAXgJ0LdawAJgFEulXjckm7\n59f+ifL9uK8mn+zN/e2TSOt1f6XLnq9Lupf3r/P49d6rrSLil6QbvEwidR+tqqeOSay+H8URQ4xX\nsoLUHVbPy+rM6+Tx2nUSEU+S1lXlhPohVZP9HPiY0n0mkLSt+uQOiv3CAdH/6u3FBkCk+yd/EDhL\n0kLSDdAnAGcDa0haRLqBymF53BHNf5hxTgHenZfxfmApNd07uRvrcOAHkm4nnTvZPtJ1748Cfqp0\nknrZcPONiLtI5y/m5XnNA6aSbg06mLs6LqR8p61LgM/mk77bkM7FnJ7X2RuBLxamOY5045lFwC3A\nayPiNuDfSf3t15POx9xebz3ljedFud4FwJl5YzwXeF/lJHVh2q8C/yppAY3/P698Jh4ndektqj1J\nXa/O7HDK6+QjwNmSbq2Z9lzSPQluzV1h38a9Gj3Fl/u2rlG6sc7KiFiZv+FzdqTbIfbkfM1eaZzW\n1k1bAZdKGkfq127V1xzbNV+zVxQfQZiZWZHPQZiZWZEDwszMihwQZmZW5IAwM7MiB4SZmRU5IMzM\nrOj/A5X7wQkyjHA2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1007342d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "total incoming edges\n",
      "300\n"
     ]
    }
   ],
   "source": [
    "# First define the input to constraint node connectivity\n",
    "\n",
    "N = 60 # Should be a multiple of 12\n",
    "print \"N = %i\" % N\n",
    "\n",
    "N_c = int(5*N/12)\n",
    "print \"N_c = %i\" % N_c\n",
    "\n",
    "z = 5 # Should be 5\n",
    "print \"z = %i\" % z\n",
    "\n",
    "z_c_avg = int(N*z/N_c)\n",
    "print \"z_c_avg = %i\" % z_c_avg\n",
    "\n",
    "r = 8\n",
    "print \"r = %i\" % r\n",
    "\n",
    "M_typical = 2**(z_c_avg-r)\n",
    "print \"M_typical = %i\" % M_typical\n",
    "\n",
    "input_to_constraint_node_connectivity = matching_assignment(N, N_c, z, z_c_avg)\n",
    "    \n",
    "constraint_node_to_input_connectivity = []\n",
    "for i in range(N_c):\n",
    "    temp = []\n",
    "    for j in range(N):\n",
    "        if i in input_to_constraint_node_connectivity[j]:\n",
    "            temp.append(j)\n",
    "    constraint_node_to_input_connectivity.append(temp)\n",
    "    \n",
    "plt.figure()\n",
    "plt.title(\"z histogram: should all be 5\")\n",
    "plt.xlabel(\"outgoing edges\")\n",
    "plt.ylabel(\"count\")\n",
    "plt.hist([len(q) for q in input_to_constraint_node_connectivity])\n",
    "plt.show()\n",
    "\n",
    "print \"\\ntotal outgoing edges\"\n",
    "print np.sum([len(q) for q in input_to_constraint_node_connectivity])\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"z_c histogram: should all be 12\")\n",
    "plt.xlabel(\"incoming edges to constraint node\")\n",
    "plt.ylabel(\"count\")\n",
    "plt.hist([len(q) for q in constraint_node_to_input_connectivity])\n",
    "plt.show()\n",
    "\n",
    "print \"\\ntotal incoming edges\"\n",
    "print np.sum([len(q) for q in constraint_node_to_input_connectivity])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated number of stable states: 7.17465e-43\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/adammarblestone/anaconda/envs/py27/lib/python2.7/site-packages/ipykernel/__main__.py:7: RuntimeWarning: divide by zero encountered in log10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEPCAYAAABhkeIdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFOW1x/HvQUWjoqKokEEUBCIa2cUFhBFFlNwouF28\nSVxA4y4mahSNMiGJQbka3DWJghKV6DUqBGRTRqORTZZBg4BxhajEuMQlRoRz/3hrsMUZZunqql5+\nn+fpZ7rf6u46Uw9wOPVu5u6IiIjEpUnaAYiISHFRYhERkVgpsYiISKyUWEREJFZKLCIiEislFhER\niVXBJhYzO8rMXjKzlWZ2WdrxiIhIYIU4j8XMmgArgcOBvwMLgKHu/lKqgYmISMFWLL2AVe7+uruv\nAyYBx6Yck4iIULiJpQx4M+P16qhNRERSVqiJRURE8tSWaQfQSGuANhmvW0dtX2FmhdeBJCKSMne3\nbD5fqBXLAqC9me1pZk2BocDkmt74q185LVo4d9zhbNjguJfmY9SoUanHkC8PXQtdB12L2h9xKMjE\n4u7rgfOBmcCLwCR3X17Tey+/HCor4Xe/g4ED4Y03EgxURKQEFWRiAXD36e7+LXfv4O5jNvfe/faD\n556D8nLo0QN++1uIKTGLiMgmCjaxNNSWW8IVV8CcOXDnnXDUUaVVvZSXl6cdQt7QtQh0Hb6kaxGv\ngpwgWV9m5jX9fuvWwXXXwbhxMGYMDBsGllVXlYhIcTAzPMvO+5JMLNWqquC002D33cPtsdatk4tN\nRCQfxZFYSuZWWE06d4Z58+CQQ6BbN7j7bvW9iIhkq6QrlkxVVXDqqdCqFfzmN6peRKQ0qWKJUefO\nMH8+HHRQqF7Gj1f1IiLSGKpYarBkSeh7KSsL1UuZViETkRKhiiVHunYN1csBB4Tq5Z57VL2IiNSX\nKpY6LF4cqpc2bcL8l29+M57YRETykSqWBHTrBgsWhJ9du8LEiapeREQ2RxVLAyxaFKqXvfYK1Uur\nVrF9tYhIXlDFkrDu3WHhQujSJTzuu0/Vi4jIplSxNNLzz4d5L+3bwx13QMuWOTmNiEiiVLGkqEeP\nkFz22y9UL/ffr+pFRARUscRi4cJQvXTsGKqX3XfP+SlFRHJCFUue6NkzdOx36hRm8D/wgKoXESld\nqlhiNn9+GDnWqRPcfjvstluipxcRyYoqljzUq1eoXjp0CNXLgw+mHZGISLJUseTQvHmhevn2t+HW\nW1W9iEj+U8WS5w48MCwJ07ZtqF4eeijtiEREck8VS0Lmzg3VS+fOoXrZdde0IxIR+TpVLAXkoINC\n9bLnniG5PPxw2hGJiOSGKpYU/OUvcPrpYWHLW26BFi3SjkhEJCjqisXMRpnZajNbFD2Oyjg20sxW\nmdlyMzsyzTgb45BDwmZiZWWw//7wxz+mHZGISHzytmIxs1HAR+5+wybtnYD7gQOA1sBsoENNpUm+\nViyZnn02VC89eoTqZZdd0o5IREpZUVcskZp+uWOBSe7+hbu/BqwCeiUaVYx69w7VS6tWoXp55JG0\nIxIRyU6+J5bzzWyJmf3OzHaM2sqANzPesyZqK1jbbgs33BAmU/7kJ/C978E//5l2VCIijbNlmic3\ns1lA5pKNBjhwJXAbMNrd3cx+AVwPnNHQc1RUVGx8Xl5eTnl5eRYR51afPrB0KVxxRRg5dtttcOyx\naUclIsWssrKSysrKWL8zb/tYMpnZnsAUd+9sZpcD7u7XRsemA6PcfV4Nn8v7PpbaPP00DBsGBx8M\nN94IO++cdkQiUgqKuo/FzDK3zjoOeCF6PhkYamZNzawt0B6Yn3R8uda3b6hemjcPfS9TpqQdkYhI\n/eRtxWJm9wJdgQ3Aa8BZ7v5OdGwkMBxYB4xw95m1fEfBViyZnnoqVC+9e4fqpXnztCMSkWIVR8WS\nt4klDsWSWAA++QQuvzyMGrvjDviv/0o7IhEpRkosdSimxFKtsjJUL337wrhxsNNOaUckIsWkqPtY\npGbl5VBVBdtvH/pepk1LOyIRka9SxVLAnnwShg8PyebXv1b1IiLZU8VS4vr3D9XLN74RqpfHH087\nIhERVSxF44knQvVy+OFhFv+OO9b9GRGRTalikY0OPxyWLYOttgrVy4wZaUckIqVKFUsRmj0bzjgD\nBgyA66+HHXZIOyIRKRSJVCxmtquZnWtmk8xsrpk9Fz0/z8x2y+bkkhtHHBH6XsxC9TKzxumjIiK5\nsdmKxczuAtoB0wnLpvydsFBkK8JS9UcDL7t7gxeHTEKpViyZZs6EM8+EgQPhf/9X1YuIbF7OJ0ia\nWWd3r6ojiDrfkxYlluDDD+GSS0KSueuuUNGIiNREM+/roMTyVTNmhOpl0CAYOxaaNUs7IhHJN0n1\nsXQzs7FmNs/M3jGzt6K+lrFm1i2bk0uyBg4MI8fWrQt9L088kXZEIlKM6roVNg14D5hCzX0sxwA7\nuft3ch9qw6liqd3jj8MPfxgWs7zuOlUvIhIk0ceym7uvrSOIXd39H9kEkStKLJv3wQfw4x/DnDmh\n76V//7QjEpG0Jd7HYmbbEHZv/E82J02KEkv9TJsWqpdjj4Vrrw0LXIpIacp5H4sFQ8zsITNbA7wC\nvG5mq6O2IWaWVQCSvkGDQt/LJ59A585haX4Rkcaq61bY08DThD6Wxe7+edTeFOhG6GPp4+79Eoi1\nwVSxNNzUqXDWWTBkCIwZA9ttl3ZEIpKkJPpYmlYnk2zekxYllsZ5/30YMQKefRbuvhv65eV/G0Qk\nF3J+KyyjQplYw8knZr5Hikfz5nDvvWGHyv/5H7jwwnCbTESkPuq7uvF+mS/MbAugR/zhSD757ndD\n38v770OXLvDnP6cdkYgUgro670ea2UdAZzP7V/T4CFgLPJZIhJKqnXeGiRPDHi9Dh8JFF8Gnn6Yd\nlYjks3oNNzazX7n7yATiiZX6WOL13nvhtti8eTB+PPTpk3ZEIhK3ROexmFlzoAOwTXWbuz+dzclz\nTYklNx59FM49F/77v+GXv4Rtt007IhGJS2I7SJrZGYRhxzOAn0U/K7I5cfS9J5jZC2a23sy6b3Js\npJmtMrPlZnZkRnt3M6sys5VmNi7bGKThBg8OfS/vvANdu4bRYyIi1erbeT8COAB43d0PI8xh+SCG\n8y8DhgBPZTaaWSfgJKATYc+X2zImYt4ODHf3jkBHMxsYQxzSQLvsAvffH+a6nHACXHwx/PvfaUcl\nIvmgvonlM3f/DMDMtnb3l4BvZXtyd1/h7qsIC1tmOhaY5O5fuPtrwCqgl5m1BJq5+4LoffcCg7ON\nQxrvuONC9bJmTahennsu7YhEJG31TSyrzWwn4FFglpk9Bryeu7AoA97MeL0maisDVmfGFbVJilq0\ngEmT4JprQqK59FJVLyKlbMv6vMndh0RPK8xsDrAjYbviOpnZLGD3zCbAgSvdfUoDYm2UioqKjc/L\ny8spLy/P9SlL1vHHQ9++cP750K0bTJgABx2UdlQisjmVlZVUxrxAYENGhW1BSBAbk5G7vxFLECFZ\nXezui6LXl4ev92uj19OBUYQqaY67d4rahwL93P2cWr5Xo8JS8tBDcMEFcMopMHo0bLNN3Z8RkfQl\nOSrsAuAdYBYwNXr8KZsT13SajOeTgaFm1tTM2gLtgfnu/jbwoZn1ijrzT0ETNfPSiSdCVRW8+mqo\nXubNSzsiEUlKfSdIvgwc6O7/jPXkZoOBm4EWhFFmS9z96OjYSGA4sA4Y4e4zo/YewATCfJpp7j5i\nM9+viiUPPPhgmFh52mlQUaHqRSSfJTZBMrpVNcDdv8jmZElTYskfa9eGSZXLl4e+lwMOSDsiEalJ\nEsvm/zh6uh9hePFUYOPuke5+QzYnzzUllvziDn/4Q1iSf9iwUL1svXXaUYlIpiT6WJpFjzcI/StN\nM9qaZXNiKT1mYSHLqipYsQK6d4cFC+r+nIgUlgbteV/rl5jd7O4XxBBPrFSx5C93eOAB+NGP4Iwz\n4OqrVb2I5IPERoXVQ++YvkdKhFnYRGzpUnjxRejZE55/Pu2oRCQOcSUWkUZp2RIeeQQuvxwGDYKr\nroLPtSepSEFTYpHUmcH3vgdLloQKpmdPWLQo7ahEpLHiSixZ3Y8TAWjVCh57LKw1dtRRod9F1YtI\n4anvzPtDoyVdMtsy90+5MdaopGSZwQ9+EKqXxYvDfJclS9KOSkQaor4TJD8FFgAnuvvaqG2Ru3ff\n/CfTpVFhhc0dJk6ESy4JkyuvuAKaNk07KpHiluSosBXAWOApMzuk+vzZnFikLmZhEcvFi2HhQujV\nK/TBiEh+q29icXf/E3AMcIuZnU9Y+l4k58rKYMoUuOgiGDAAfvYzWLcu7ahEpDb1TSwGEO32eCjQ\nF+icq6BENmUWFrFcvDislHzggWEGv4jknzoTS7QVcTcz2xrA3T9x95OAdjmPTmQTZWUwdWrYTOzw\nw+HnP1f1IpJv6lOx3GFm2wC3ZTbGtcmXSEOZhUUsFy2CZ58Nu1QuW5Z2VCJSbbOJxcz6AguBPwPP\nR69F8sIee8Djj4cRY/37wy9/CV8U1MYOIsWproqlen96jQCTvGQGw4eH6uXpp0P18sILaUclUto2\nm1jc/SmgF6HDvqe7P51IVCINtMceMH06nH02HHYYXHONqheRtNQ5QdLMtgfOBG53988SiSommiBZ\nmt54I1QxH3wA99wD++6bdkQihSORCZLu/jHwP4WWVKR0tWkDM2fCmWdCv34wZoyqF5Ek1XdJl18D\nWwF/AD6pbnf3vF6DVhWLvP56qF4++gjGj1f1IlKXnO95n3GiOTU0u7v3z+bkuabEIhDWHLvzzrDX\nyyWXwMUXw5Zbph2VSH5KLLEUKiUWyfTaa6F6+eQTmDAB9tkn7YhE8k8+bU0skvf22gtmzQoLW/bp\nA2PHwvr1aUclUnxSTSxmdoKZvWBm6zP3dzGzPc3sUzNbFD1uyzjW3cyqzGylmY1LJ3IpVE2ahAmV\nCxbAtGlw6KGwYkXaUYkUl7QrlmXAEOCpGo697O7do8e5Ge23A8PdvSPQ0cwGJhGoFJe2beGJJ8KW\nyH36wPXXq3oRiUuDE4uZ/Sauk7v7imjF5Jru532tzcxaAs3cfUHUdC8wOK54pLQ0aQLnnRdWS54y\nBfr2hZUr045KpPA1pmLpGXsUNdsrug02x8z6RG1lwOqM96yO2kQarV07ePJJOPlkOOQQuOEGVS8i\n2WjMoMu1DXmzmc0Cds9sIqw/dqW7T6nlY38H2rj7+1Hfy6Nm1qgZCBUVFRufl5eXU15e3pivkSLX\npElYiv/oo+H00+GPfwzzXjp0SDsykdyqrKyksrIy1u/Mi+HG0TyZi2ubcFl9nJBw5rh7p6h9KNDP\n3c+p5XMabiwNtmED3HILjB4NP/0pXHhhSDwipaDYhhtv/EXMrIWZNYmetwPaA6+4+9vAh2bWy8wM\nOAV4LJVopWg1aRKSydy58PDDYVmYl19OOyqRwpH2cOPBZvYmcBDwJzN7PDrUF6gys0XAg8BZ7v5B\ndOw84C5gJbDK3acnHbeUhvbt4amn4Pjjw3L8N94YqhkR2by8uBWWK7oVJnFZtSr0vWyxBdx9N+y9\nd9oRieRGIrfCzGygmd1uZpOjx+1mdlQ2JxUpNB06hOpl8GA48EC4+WZVLyK12WzFEs1s70iYL1I9\nzLc1oW9jlbuPyHmEWVDFIrmwYkWoXpo2DdVLu3ZpRyQSn5wvQmlmK6MZ7pu2G7DS3fN6MKYSi+TK\n+vUwblzY66WiAs45RyPHpDgkcSvsMzM7oIb2AwBt/CUla4stwvL7zzwDv/89HH44vPpq2lGJ5Ie6\nEstpwC1m9lczmxk9lgM3RcdEStq3vhWSy3e+A716wW23qe9FpL4bfbXky6VT1kTzSfKeboVJkl56\nCU47DbbbDu66KyzTL1JoEpsg6e5vu/vz0aMgkopI0vbZJ1QvAwfCAQfAHXeE3StFSk2j57GY2SJ3\n7173O9OjikXS8te/hpFjzZqF6mXPPdOOSKR+Ul3SJd+Tikia9t0Xnn0WjjgCevaEO+9U9SKlQzPv\nRXLsxRdD30vz5vC730GbNmlHJFK7VCsWM1uWzYlFSsV++8Fzz8Fhh0GPHvDb36p6keJW1wTJ42o7\nBNzh7rvmJKqYqGKRfPPCC6F6adEiJJg99kg7IpGvSmLm/TrgPsLGXJs6wd2bZXPyXFNikXy0bh1c\nd92XM/eHDQPL6q+xSHySSCzPA6e6+ws1HHvT3fP6/1tKLJLPli0L1ctuu4XqpXXrtCMSSaaP5SLg\nX7UcG5LNiUVK3f77h83EeveG7t3Dgpb6f5AUA40KE8kDVVVw6qnQqlWoXsrK6v6MSC7kvGIxs1Fm\ntvtmjrc0s4psAhAR6NwZ5s8PO1V26wYTJqh6kcJVVx/Ld4BLgK2ARcBbhBFhLYEehBWO/9fdH6/1\nS1KkikUK0dKloe+lrCxMrFT1IknKeed9xolaA72B6oUpXgeedffVtX8qfUosUqg+/xyuuSasljx2\nLJxyikaOSTISSSxm1srd3zKzloW2AKUSixS6JUtC30ubNqF6+eY3045Iil1SM+9Hm1kL4GfZnEhE\nGq5rV1iwIIwa69oVJk5U34vkv7r6WE4F2gCnA3cDb7j7vQnFljVVLFJMFi0KfS9t24Yl+Vu1Sjsi\nKUZJVCyVwCvAX6OfT2VzMhFpvO7dYeHCMIKsa1e47z5VL5KfNptY3P114FDgGKBP9Do2ZnadmS03\nsyVm9rCZ7ZBxbKSZrYqOH5nR3t3MqsxspZmNizMekXzXtCn8/OcwbVpYDmbIEHi7oHo+pRTU2cfi\n7me7+wZ3PzcH558J7OfuXYFVwEgAM9sXOAnoBBwN3Ga2cUzM7cBwd+8IdDSzgTmISySv9egRqpdv\nfxu6dIH771f1Ivmjrj6W8dS8ACWAu/vw2AIxGwwc7+4/MLPLo++/Njr2OFBBGOb8pLvvG7UPBfq5\n+zm1fKf6WKToLVwY+l46doTbb4fda53SLFK3JPpY/gRM3eRRBRwOxF0pDAOmRc/LgDczjq2J2sqA\nzLkzq6M2kZLVsyc8/zzss0+oXiZNUvUi6dpycwfd/eHq52bWDrgC6AuMAe6qzwnMbBaQ+X8oI1RB\nV7r7lOg9VwLr3P2BBkVfDxUVFRufl5eXU15eHvcpRFK39dZhQuXgwaF6eeihUL3stlvakUm+q6ys\npLKyMtbvrM8EyX2AnwLdgLHA7939i9gCMDsNOBPo7+7/ido2vRU2HRhFuBU2x907Re26FSayic8+\ng4qKsN7YTTfBSSelHZEUkiT2Y3mIsCbY9cCDwPrM4+7+XlYnNzsq+u6+7v7PjPZ9CRuMHUi41TUL\n6ODubmZzgQuBBYRbcze5+/Ravl+JRUrWvHmhevn2t+HWW1W9SP0kkVhe48vO++qf1Sd0d2+X1cnN\nVgFNgeqkMrd69JmZjQSGA+uAEe4+M2rvAUwAtgGmufuIzXy/EouUtM8+g1Gj4J574Oab4cQT045I\n8l1ii1AWKiUWkWDu3FC9dOkCt9wCu+6adkSSr5LYj6XOiqQ+7xGRdB10ECxeHBaz7NwZHn647s+I\nNFZdt8ImAdsBU4CFfHU/lp7AscC/3P3k3IfacKpYRL7uuedC9dKtW6heWrRIOyLJJ0ktm783MBTo\nQ9iPxQmjs54BJrn7K9kEkEtKLCI1+/e/4aqrwoz9W26B445LOyLJF+pjqYMSi8jmPfssnH56mGR5\n882wyy5pRyRpS2Q/FjP7RfRzdDYnEpH807t32EysZUvYf3949NG0I5JiUJ9bYccCRwIz3H1yIlHF\nRBWLSP0980yoXnr1ChMrVb2UpiRGhY0C+gMnA/3N7OpsTiYi+atPH1i6NAxF7twZHnss7YikUNWn\nYrkRuA04x90vSiSqmKhiEWmcP/85VC8HHww33gg775x2RJKUpPa8v9vdVwDjszmRiBSOQw8N1cvO\nO4e+lylT0o5ICkl9Ekv1IhDH5zIQEckv220XqpUHHoAf/QhOPRXefz/tqKQQ1CexLDCzWwkTJEWk\nxPTtG6qXHXcM1cvUqWlHJPmurpn3o4CdgR8A9wLvuXvBDDtWH4tIvCorYdiwkGzGjYOddko7Iolb\nzvtY3P1n0dODo9cFk1REJH7l5VBVBdtvH6qXadPq/IiUoPqMCuvq7kvMrIu7L00orlioYhHJnTlz\nQvVy2GFwww2qXopFUqPCXjGzvTdNKmbWOZsTi0hhO+wwWLYMttkmzHuZXuN2e1KK6upjOQkYB6wF\ntgJOc/cF0bFF7t49kSgbSRWLSDKefBKGD4f+/UP1suOOaUckjZVExXIF0MPduwKnAxPNbEj1+bM5\nsYgUj/79Q9/LVluFvpcZM9KOSNJUV8WyzN33z3jdCvgTcA+helHFIiJfMXs2nHEGDBgA118PO+yQ\ndkTSEElULB9F+7EA4O5vAeWEDb72y+bEIlKcjjgiVC9NmoTqZdastCOSpNVVsXQBPnX3VZu0bwWc\n5O735Ti+rKhiEUnXrFmhejnqKBg7VtVLIUiiYqnaNKkAuPu66qRiZuprEZEaDRgQRo5t2BBGjs2e\nnXZEkoS6KpanCH0qj7n7yk2OdQSGAIPcvV9Oo2wkVSwi+WPGDDjzTBg0KFQvzZqlHZHUJImKZQDw\nLnCrmb1lZivMbKWZ/R24BXg7eo+IyGYNHBiqly++CNXLE0+kHZHkSr33vDezJkCL6OW77r4h65Ob\nXQd8F/gP8DfgdHf/l5ntCSwHXoreOtfdz40+0x2YAGwDTNvcHjGqWETy0+OPww9/CN/9Llx3XVgi\nRvJDUjPvMbPdga5Aa0IyyjqpRGYC+0XzZFYBIzOOvezu3aPHuRnttwPD3b0j0NHMBsYUi4gk5Oij\nQ/Xy2WehepkzJ+2IJE51bU3c1czmApXAddHjKTObG1UOWXH32RlJai4hcW08fQ3xtASaVc/+J6y4\nPDjbOEQkeTvtBHffDbfcAqecAuedBx9/nHZUEoe6KpYJwAh37+TuR0SPfYCLiH9HyWHA4xmv9zKz\nRWY2x8z6RG1lwOqM96yO2kSkQA0aFOa9fPJJqF4qK9OOSLK1ZR3Ht3P3eZs2uvtcM9uuPicws1nA\n7plNgANXuvuU6D1XAuvc/f7oPX8H2rj7+1Fl9KiZ7Vuf822qoqJi4/Py8nLKy8sb8zUikkPNm8OE\nCWETse9/H4YMgTFjwi6WkluVlZVUxpzN6xpufBOwN+GW05tR8x7AKcCr7n5+1gGYnQacCfR39//U\n8p45wMWEhDPH3TtF7UOBfu5+Ti2fU+e9SIF5/3246CJ45hkYPz5sKibJiaPzvj77sRxNWMKl+pbT\nGmCyu2e9xY+ZHQVcD/R1939mtLcg7Fa5wczaAU8B+7v7B1Gfz4XAAmAqcJO717hgtxKLSOGaMgXO\nOQeOPx6uuUbVS1ISSSy5ZGargKZAdVKZ6+7nmtlxwGjgc2ADcHV1IjOzHnx1uPGIzXy/EotIAXvv\nvVC9/OUvoXo59NC0Iyp+qSYWM/uNu/8wm5PnmhKLSHGYPDlULyeeGKqXbbdNO6LilfN5LGa2cy2P\nXYBB2ZxYRKS+jjkmzHt5913o0iX0v0j+qqvzfj3wOl+dU+LR6zJ3b5rb8LKjikWk+Dz6KJx7Lgwd\nCr/4haqXuCUx8/4VoNzd22Y82rl7W+CdbE4sItIYgweH6uXtt6Fr19D/IvmlrsQyDmhey7HrYo5F\nRKRedtkF7r8frr0WTjgBLr4Y/v3vtKOSaqmOCss13QoTKX7vvgsXXACLFoVJlgcfnHZEhS2xUWHR\n8N9NfQgsc/e12QSQS0osIqXj4Yfh/PPDzP3Ro+Eb30g7osKUZGKZChwMVK9BWg48D7QFRrv7xGyC\nyBUlFpHS8o9/hOSydGmoXg46KO2ICk9iy+YT1hTr5O7Hu/vxwL6E0WEHApdlE4CISFx23RX+8Icw\nWmzIELjssrA0vySrvollD3fPHAW2Nmp7D1gXf1giIo13wglhxeRXX4Xu3WH+/LQjKi31TSyVZvYn\nMzvVzE4FJkdt2wEf5C48EZHG2XVXePBBqKgIEywvv1zVS1Lq28diwHFA9b4ozwIP53sHhvpYRARg\n7dowqXL58tD3csABaUeUvxJdKyzanrgXoW9lfj6PBqumxCIi1dxD/8uIETB8OIwaBVtvnXZU+SfJ\nPe9PAuYDJwAnAfPM7IRsTiwikiSzsAxMVRW89BL06AELF6YdVXGq762wpcCA6irFzHYFZrt7lxzH\nlxVVLCJSE3eYNCksyX/GGXD11apeqiU53LjJJre+/tmAz4qI5BUzOPnkMN/lxRehZ094/vm0oyoe\n9U0O081shpmdFm0lPBXIegdJEZE0tWwJjzwSRowNGgRXXQWff552VIWvIZ33xwO9o5d/dvdHchZV\nTHQrTETq66234Oyzw9yXCRPC/JdSVPBbE+eaEouINIQ73HdfWC357LPhyiuhaV7vOhW/nCcWM/uI\nMLz4a4cAd/cdsjl5rimxiEhjvPUWnHUWvPFGqF66dk07ouSoYqmDEouINJY7TJwIl1wSJldecUVp\nVC9JjgoTESkpZnDKKbBkSZjv0qtXGEUmdVNiERHZjG9+E6ZMCXNeBgwIe72s09K7m6VbYSIi9bRm\nDZx5Jrz9duh76dw57YjiV/C3wsxstJktNbMlZjbbzFpnHBtpZqvMbLmZHZnR3t3MqsxspZmNSydy\nESlFZWUwdWrYCvmII+DnP1f1UpNUKxYz297dP46eXwB0dvczzWxf4D7gAKA1MBvo4O5uZvOA8919\ngZlNA2509xm1fL8qFhHJidWrQ/Wydm2oXvbfP+2I4lHwFUt1UolsR1gqBuAYYJK7f+HurwGrgF5m\n1hJo5u4LovfdCwxOKl4RkWqtW8O0aXDeedC/P/zyl/DFF2lHlR9S77w3s1+Y2RvAacCvouYy4M2M\nt62J2sqA1Rntq6M2EZHEmcGwYbBoETz9NBx0ELzwQtpRpW/LXJ/AzGYBu2c2ESZdXunuU9z9p8BP\nzewyYBxwepznr6io2Pi8vLyc8vLyOL9eRIQ99oDp0+Guu+Cww+DHP4ZLL4Utc/4vbPYqKyuprKyM\n9TvzZlS4g5BbAAAHdUlEQVSYme0BTHP3/c3scsLM/mujY9OBUcDrwBx37xS1DwX6ufs5tXyn+lhE\nJFFvvBGW4v/gg9D3su++aUfUMAXfx2Jm7TNeDgaWRM8nA0PNrKmZtQXaE3atfBv40Mx6RdslnwI8\nlmjQIiKb0aYNzJgRkku/fjBmTOn1vaQ9Kuz/gI7AeuAV4JyMzcRGAsOBdcAId58ZtfcAJgDbECqc\nEZv5flUsIpKa118P2yB/9BGMH18Y1YvWCquDEouIpM0d7rwz7PVy6aVh5eQttkg7qtopsdRBiUVE\n8sVrr4Xq5ZNPQt/LPvukHVHNCr6PRUSkVOy1F8yaBaeeCoceCmPHwvr1aUeVG6pYREQS9uqroXr5\n7LPQ9/Ktb6Ud0ZdUsYiIFKC2bWH2bPj+96FPH7j++uKqXlSxiIik6JVXwuz9detC9dKxY7rxqGIR\nESlw7drBk0/CySdD797w618XfvWiikVEJE/87W+helm/PlQvHTokH4MqFhGRIrL33jBnDpx0Ehxy\nCIwbBxs2pB1Vw6liERHJQy+/HKoX91C9tG9f92fioIpFRKRItW8PlZVwwglw8MFhY7FCoYpFRCTP\nrV0Lu+2WzLm0pEsdlFhERBpGt8JERCTvKLGIiEislFhERCRWSiwiIhIrJRYREYmVEouIiMRKiUVE\nRGKlxCIiIrFSYhERkVgpsYiISKxSTSxmNtrMlprZEjObbWato/Y9zexTM1sUPW7L+Ex3M6sys5Vm\nNi696EVEpCZpVyzXuXsXd+8KPAZUZBx72d27R49zM9pvB4a7e0ego5kNTDDeglVZWZl2CHlD1yLQ\ndfiSrkW8Uk0s7v5xxsvtgHczXn9tETQzawk0c/cFUdO9wODcRVg89BfnS7oWga7Dl3Qt4rVl2gGY\n2S+AU4BPgQMzDu1lZouAD4Gr3P0ZoAzI3JVgddQmIiJ5IucVi5nNivpEqh/Lop/fBXD3n7p7G2A8\nUN1n8hbQxt27AxcD95vZ9rmOVUREspc3+7GY2R7ANHffv4ZjcwgJ5u/AHHfvFLUPBfq5+zm1fGd+\n/HIiIgUk2/1YUr0VZmbt3f3l6OVgYEnU3gJ4z903mFk7oD3wirt/YGYfmlkvYAHhFtpNtX1/thdH\nREQaLu0+ljFm1hFYD7wCVFcefYHRZvY5sAE4y90/iI6dB0wAtiFUONOTDVlERDYnb26FiYhIcUh7\nHktOmNlRZvZSNInysrTjyTUzu8vM3jGzqoy25mY208xWmNkMM9sx49hIM1tlZsvN7Mh0os4NM2tt\nZk+a2YvRQJELo/aSux5mtrWZzTOzxdH1uCZqL7lrAWBmTaIJ15Oj1yV5HQDM7LVocvpiM5sftcV3\nPdy9qB6EZPkysCewFaHfZp+048rx79wH6ApUZbRdC/wken4ZMCZ6vi+wmHAbdK/oWlnav0OM16Il\n0DV6vj2wAtinhK/HttHPLYC5QO8SvhY/An4PTI5el+R1iH7HV4Dmm7TFdj2KsWLpBaxy99fdfR0w\nCTg25ZhyysMcn/c3aT4WuCd6fg9fTiQ9Bpjk7l+4+2vAKsI1Kwru/ra7L4mefwwsB1pTutfj0+jp\n1oT/dL1PCV6LaLmoQcDvMppL7jpkML5+xyq261GMiaUMeDPjdalOotzN3d+B8I8tsFvUvun1WUOR\nXh8z24tQyc0Fdi/F6xHd/lkMvA1UuvtfKc1r8WvgUiCzU7kUr0M1B2aZ2QIzOyNqi+16pD0qTJJT\nUqM0ogm1/weMcPePa5jTVBLXw903AN3MbAdghpmV8/XfvaivhZl9B3jH3ZdEv39tivo6bKK3u79l\nZrsCM81sBTH+uSjGimUN0CbjdeuordS8Y2a7w8Y11tZG7WuAPTLeV3TXx8y2JCSVie7+WNRcstcD\nwN3/BUwDelJ616I3cIyZvQI8APQ3s4nA2yV2HTZy97ein/8AHiXc2ortz0UxJpYFQPto6f2mwFBg\ncsoxJcH46sKdk4HTouenElaPrm4famZNzawtYfLp/KSCTMjdwF/d/caMtpK7HmbWonpkj5l9AxhA\n6IQtqWvh7le4ext3b0f49+BJd/8BMIUSug7VzGzb6iWyzGw74EhgGXH+uUh7dEKORjwcRRgNtAq4\nPO14Evh97ycsd/Mf4A3gdKA5MDu6DjOBnTLeP5IwsmM5cGTa8cd8LXoTJtwuIfwjuij687BzqV0P\nYP/o918MLAUuidpL7lpk/H79+HJUWEleB6Btxt+PZdX/RsZ5PTRBUkREYlWMt8JERCRFSiwiIhIr\nJRYREYmVEouIiMRKiUVERGKlxCIiIrFSYhFJiJltMLOxGa8vNrOr04xJJBeUWESS8x/gODPbOe1A\nRHJJiUUkOV8AvwF+nHYgIrmkxCKSHAduBb5nZs3SDkYkV5RYRBLkYfOxe4ARaccikitKLCLJuxEY\nDmybdiAiuaDEIpIcA3D394EHgTM2/3aRwqTEIpKczKXErwd2obR2LZQSoWXzRUQkVqpYREQkVkos\nIiISKyUWERGJlRKLiIjESolFRERipcQiIiKxUmIREZFYKbGIiEis/h+AQV6wZbHdTwAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f6016d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Estimating the number of stable states: From equation 2 of the paper, requiring r*z/z_c <= 1 for exponential scaling\n",
    "z_hat = float(z)/float(z_c_avg)\n",
    "num_states = 2**(N*(1.0 - r*z_hat))\n",
    "print \"Estimated number of stable states: %g\" % num_states\n",
    "\n",
    "plt.figure()\n",
    "plt.plot([np.log10(2**(Nprime*(1.0 - r*z_hat))) for Nprime in range(1000)])\n",
    "plt.xlabel(\"N\")\n",
    "plt.ylabel(\"log10(2**(N*(1.0 - r*z_hat)))\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1])]\n"
     ]
    }
   ],
   "source": [
    "def generate_random_parity_states(input_size, num_parity_states_to_generate):\n",
    "    parity_states = []\n",
    "    for i in range(num_parity_states_to_generate):\n",
    "        parity = 1\n",
    "        state = []\n",
    "        while parity == 1:\n",
    "            state = np.random.randint(0,2,size = input_size)\n",
    "            parity = sum(state) % 2\n",
    "            if state in parity_states:\n",
    "                parity = 1 # A hack to make sure all the permitted parity states in a node are distinct\n",
    "            if parity == 0:\n",
    "                parity_states.append(state)\n",
    "\n",
    "    return parity_states\n",
    "\n",
    "print generate_random_parity_states(12, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def M_for_constraint_node(constraint_node_index):\n",
    "    z_c_for_this_node = len(constraint_node_to_input_connectivity[constraint_node_index])\n",
    "    return 2**(z_c_for_this_node - r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{   'bias_term': 6,\n",
      "    'inputs': [3, 13, 14, 20, 26, 29, 32, 36, 43, 44, 50, 52],\n",
      "    'num_inputs': 12,\n",
      "    'offs_preffered': [13, 20, 29, 36, 43, 52],\n",
      "    'ons_preferred': [3, 14, 26, 32, 44, 50],\n",
      "    'preferred_parity_config': array([1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0])}\n"
     ]
    }
   ],
   "source": [
    "def generate_constraint_node_neurons(constraint_node_index):\n",
    "    neurons = []\n",
    "    z_c_for_this_node = len(constraint_node_to_input_connectivity[constraint_node_index])\n",
    "    M_for_this_node = M_for_constraint_node(constraint_node_index)\n",
    "    for i in range(M_for_this_node):\n",
    "        neuron = {}\n",
    "        neuron[\"inputs\"] = constraint_node_to_input_connectivity[constraint_node_index]\n",
    "        neuron[\"num_inputs\"] = len(neuron[\"inputs\"])\n",
    "        neuron[\"preferred_parity_config\"] = generate_random_parity_states(neuron[\"num_inputs\"], 1)[0]\n",
    "        neuron[\"ons_preferred\"] = [neuron[\"inputs\"][j] for j in range(neuron[\"num_inputs\"]) if neuron[\"preferred_parity_config\"][j] == 1] \n",
    "        neuron[\"offs_preffered\"] = [neuron[\"inputs\"][j] for j in range(neuron[\"num_inputs\"]) if neuron[\"preferred_parity_config\"][j] == 0]\n",
    "        neuron[\"bias_term\"] = neuron[\"num_inputs\"] - np.sum(neuron[\"preferred_parity_config\"])\n",
    "        neurons.append(neuron)\n",
    "\n",
    "    return neurons\n",
    "\n",
    "pp = pprint.PrettyPrinter(indent=4)\n",
    "pp.pprint(generate_constraint_node_neurons(N_c-1)[0]) # Show what the first neuron in the last constraint node looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_full_index(constraint_node_index, index_within_node):\n",
    "    return N + sum([M_for_constraint_node(k) for k in range(constraint_node_index)]) + index_within_node"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total number of neurons 460\n"
     ]
    }
   ],
   "source": [
    "# Generating the full weight matrix and biases\n",
    "\n",
    "last_constraint_node_index = N_c-1\n",
    "\n",
    "last_neuron_in_last_constraint_node_index = M_for_constraint_node(last_constraint_node_index) - 1\n",
    "\n",
    "total_number_neurons = get_full_index(last_constraint_node_index, last_neuron_in_last_constraint_node_index) + 1\n",
    "total_number_neurons = int(total_number_neurons) # Not sure why it is making me do this\n",
    "print \"total number of neurons %i\" % total_number_neurons\n",
    " \n",
    "full_weight_matrix = np.array([[0 for j in range(total_number_neurons)] for i in range(total_number_neurons)], dtype = np.int32)\n",
    "\n",
    "biases = np.array([0 for j in range(total_number_neurons)], dtype = np.int32)\n",
    "\n",
    "constraint_node_preferred_states = []\n",
    "\n",
    "for node in range(N_c): # Index over constraint nodes\n",
    "    neurons_in_node = generate_constraint_node_neurons(node) # Generate neurons in that constraint node\n",
    "    z_c_for_this_node = len(constraint_node_to_input_connectivity[node])\n",
    "    M_for_this_node = M_for_constraint_node(node)\n",
    "    this_node_preferred_states = []\n",
    "    for m in range(M_for_this_node):\n",
    "        neuron_full_index = get_full_index(node, m)\n",
    "        neuron = neurons_in_node[m]\n",
    "        for g in neuron[\"ons_preferred\"]: # Constraint neuron to input neuron and vice versa\n",
    "            full_weight_matrix[neuron_full_index][g] = 1\n",
    "            full_weight_matrix[g][neuron_full_index] = 1\n",
    "        for g in neuron[\"offs_preffered\"]: # Constraint neuron to input neuron and vice versa\n",
    "            full_weight_matrix[neuron_full_index][g] = -1\n",
    "            full_weight_matrix[g][neuron_full_index] = -1\n",
    "        for p in range(m, M_for_this_node): # Inhibitory internal connectivity within the same node\n",
    "            p_full_index = get_full_index(node, p)\n",
    "            if p != m:\n",
    "                full_weight_matrix[neuron_full_index][p_full_index] = -1*(neuron[\"num_inputs\"] - 1)\n",
    "                full_weight_matrix[p_full_index][neuron_full_index] =  -1*(neuron[\"num_inputs\"] - 1)\n",
    "        biases[neuron_full_index] = neuron[\"bias_term\"]\n",
    "        this_node_preferred_states.append(zip(neuron[\"inputs\"], np.ndarray.tolist(neuron[\"preferred_parity_config\"])))\n",
    "    constraint_node_preferred_states.append(this_node_preferred_states)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counting paths...\n"
     ]
    }
   ],
   "source": [
    "# Trying to actually count stable states\n",
    "\n",
    "def intersection(list1, list2):\n",
    "    return list(set(list1) & set(list2))\n",
    "\n",
    "pairs = []\n",
    "\n",
    "for k in range(N_c-1):\n",
    "    states_first = constraint_node_preferred_states[k] # Zip of the input neurons and their preferred parity states for each neuron in the node\n",
    "    states_second = constraint_node_preferred_states[k+1]\n",
    "    \n",
    "    input_identities_first = [q[0] for q in states_first[0]] # The input neurons\n",
    "    input_identities_second = [q[0] for q in states_second[0]] # The input neurons\n",
    "    \n",
    "    input_overlap = intersection(input_identities_first, input_identities_second)\n",
    "    \n",
    "    pairs_from_k = []\n",
    "    \n",
    "    for f in range(len(states_first)):\n",
    "        for s in range(len(states_second)):\n",
    "            sect = intersection(states_first[f], states_second[s])\n",
    "            if len(sect) == len(input_overlap): # For any shared input neurons, the neurons must assign the same 0/1 preference\n",
    "                pairs_from_k.append(((k, f), (k+1, s)))\n",
    "    \n",
    "    pairs.append(pairs_from_k)\n",
    "\n",
    "last = len(pairs) - 1\n",
    "\n",
    "def build_paths(k, edge):\n",
    "    if k == last:\n",
    "        return 1\n",
    "    else:\n",
    "        s = 0\n",
    "        for edge_outgoing in pairs[k+1]:\n",
    "            if edge[1] == edge_outgoing[0]:\n",
    "                s += build_paths(k+1, edge_outgoing)\n",
    "        return s\n",
    "\n",
    "count_paths = True\n",
    "if count_paths:\n",
    "    print \"Counting paths...\"\n",
    "    total_full_paths = 0\n",
    "    for edge in pairs[0]:\n",
    "        total_full_paths += build_paths(0, edge)    \n",
    "    print total_full_paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (18, 18)\n",
    "plt.figure()\n",
    "plt.title(\"Full weight matrix\")\n",
    "cmap = cm.get_cmap('gist_heat', z_c_avg + 1)\n",
    "plt.imshow(full_weight_matrix, cmap = cmap, interpolation='none')\n",
    "plt.colorbar(ticks = np.arange(-1*z_c_avg, 2, 1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (21, 7)\n",
    "plt.figure()\n",
    "plt.title(\"Input to constraint node mapping\")\n",
    "adj = full_weight_matrix[:N][:]\n",
    "plt.ylabel(\"Input neurons\")\n",
    "plt.xlabel(\"Constraint neurons\")\n",
    "cmap = cm.get_cmap('gist_heat', 3)\n",
    "plt.imshow(adj, interpolation='none', cmap = cmap, aspect = 0.4*float(total_number_neurons)/float(N))\n",
    "plt.colorbar(ticks = [-1,0,1])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (21, 8)\n",
    "plt.figure()\n",
    "plt.title(\"Neuron biases\")\n",
    "plt.xlabel(\"Neuron indices\")\n",
    "plt.xlim(0,total_number_neurons)\n",
    "plt.ylabel(\"Bias term\")\n",
    "plt.plot(biases)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Defining some Hopfield aspects\n",
    "\n",
    "def random_state(num_neurons):\n",
    "    return np.random.randint(0, 2, size = num_neurons)\n",
    "\n",
    "def time_step(state, weight_matrix, biases): # Asynchronous mode\n",
    "\n",
    "    x = state[:]\n",
    "    \n",
    "    k = np.random.randint(0, total_number_neurons, size=1)[0] # Pick a random single neuron to update\n",
    "    \n",
    "    a = 0\n",
    "    \n",
    "    for j in range(total_number_neurons):\n",
    "        a += weight_matrix[k][j] * x[j]\n",
    "        \n",
    "    a += biases[k]\n",
    "        \n",
    "    ns = x\n",
    "        \n",
    "    if a > 0:\n",
    "        ns[k] = 1\n",
    "    elif a < 0:\n",
    "        ns[k] = 0\n",
    "    elif a == 0: # If the activation is exactly zero\n",
    "        ns[k] = random.randint(0,1)\n",
    "                        \n",
    "    return np.array(ns, dtype = np.int32)\n",
    "\n",
    "def energy(s, weights, bias):\n",
    "        \n",
    "    E = 0.5*np.dot(np.matmul(weights, s), s) + np.dot(s, bias)\n",
    "        \n",
    "    return -1*E\n",
    "\n",
    "def small_perturbation(perturb_size, state):\n",
    "    new_state = np.zeros(len(state))\n",
    "    prob = float(perturb_size)/float(len(state))\n",
    "    for i in range(len(state)):\n",
    "        if np.random.rand() < prob:\n",
    "            new_state[i] = (state[i] + 1) % 2\n",
    "        else:\n",
    "            new_state[i] = state[i]\n",
    "    return np.array(new_state, dtype = np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Test out some Hopfield dynamics\n",
    "initial_random_state = random_state(total_number_neurons)\n",
    "\n",
    "num_timesteps = global_num_timesteps\n",
    "state = np.array(initial_random_state, dtype = np.int32)\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "    \n",
    "    energies.append(e)\n",
    "    \n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Initial settling\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)\n",
    "\n",
    "print \"\\n___Recovery from Perturbation___\"\n",
    "perturb_percent = 2.0\n",
    "print \"Percentage perturbed: %f\" % perturb_percent\n",
    "perturbed_state = small_perturbation(perturb_percent * total_number_neurons / 100.0, state)\n",
    "num_timesteps = global_num_timesteps\n",
    "state = perturbed_state\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "        \n",
    "    energies.append(e)\n",
    "    \n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Recovery from perturbation\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initial_nonrandom_state = np.ndarray.tolist(random_state(N)) # Visible nodes random\n",
    "initial_nonrandom_state.extend([0 for j in range(total_number_neurons - N)])\n",
    "for node in range(N_c): # Set just one neuron in each constraint node to be 1\n",
    "    M_for_this_node = M_for_constraint_node(node)\n",
    "    m = np.random.randint(0,M_for_this_node)\n",
    "    neuron_full_index = get_full_index(node, m)\n",
    "    initial_nonrandom_state[neuron_full_index] = 1\n",
    "    \n",
    "num_timesteps = global_num_timesteps\n",
    "state = np.array(initial_nonrandom_state, dtype = np.int32)\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "    \n",
    "    energies.append(e)\n",
    "\n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Initial settling\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)\n",
    "\n",
    "print \"\\n___Recovery from Perturbation___\"\n",
    "perturb_percent = 2.0\n",
    "print \"Percentage perturbed: %f\" % perturb_percent\n",
    "perturbed_state = small_perturbation(perturb_percent * total_number_neurons / 100.0, state)\n",
    "num_timesteps = global_num_timesteps\n",
    "state = perturbed_state\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "        \n",
    "    energies.append(e)\n",
    "    \n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Recovery from perturbation\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initial_nonrandom_state = np.ndarray.tolist(np.zeros(N)) # Visible nodes zero\n",
    "initial_nonrandom_state.extend([0 for j in range(total_number_neurons - N)])\n",
    "for node in range(N_c): # Set just one neuron in each constraint node to be 1\n",
    "    M_for_this_node = M_for_constraint_node(node)\n",
    "    m = np.random.randint(0,M_for_this_node)\n",
    "    neuron_full_index = get_full_index(node, m)\n",
    "    initial_nonrandom_state[neuron_full_index] = 1\n",
    "    \n",
    "num_timesteps = global_num_timesteps\n",
    "state = np.array(initial_nonrandom_state, dtype = np.int32)\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    energies.append(e)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "    \n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Initial settling\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)\n",
    "\n",
    "print \"\\n___Recovery from Perturbation___\"\n",
    "perturb_percent = 2.0\n",
    "print \"Percentage perturbed: %f\" % perturb_percent\n",
    "perturbed_state = small_perturbation(perturb_percent * total_number_neurons / 100.0, state)\n",
    "num_timesteps = global_num_timesteps\n",
    "state = perturbed_state\n",
    "energies = []\n",
    "e = energy(state, full_weight_matrix, biases)\n",
    "print \"energy %f\" % e, state[:2*N]\n",
    "\n",
    "constraint_node_satisfactions = []\n",
    "sats = []\n",
    "for p in range(N_c):\n",
    "    M = M_for_constraint_node(p)\n",
    "    ind1 = get_full_index(p, 0)\n",
    "    ind2 = get_full_index(p, M-1)\n",
    "    sats.append(np.sum(state[ind1:ind2+1]))\n",
    "constraint_node_satisfactions.append(sats)\n",
    "\n",
    "for i in range(num_timesteps):\n",
    "    new_state = time_step(state, full_weight_matrix, biases)\n",
    "    e = energy(new_state, full_weight_matrix, biases)\n",
    "    energies.append(e)\n",
    "    state = new_state\n",
    "    \n",
    "    sats = []\n",
    "    for p in range(N_c):\n",
    "        M = M_for_constraint_node(p)\n",
    "        ind1 = get_full_index(p, 0)\n",
    "        ind2 = get_full_index(p, M-1)\n",
    "        sats.append(np.sum(state[ind1:ind2+1]))\n",
    "    constraint_node_satisfactions.append(sats)\n",
    "    \n",
    "print \"energy %f\" % e, new_state[:2*N]\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"Hopfield energy function\")\n",
    "plt.title(\"Recovery from perturbation\")\n",
    "plt.plot(energies, linewidth = 2)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"step\")\n",
    "plt.ylabel(\"# neurons active in constraint node\")\n",
    "plt.title(\"constraint node satisfaction\")\n",
    "for p in range(N_c):\n",
    "    plt.plot([c[p] for c in constraint_node_satisfactions], label = \"node %i\" %p, linewidth = 2, alpha = 0.5)\n",
    "plt.legend(loc=9, bbox_to_anchor=(0.5, -0.1), ncol=4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}