block_gibbs.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9749d3c9078f66cd447fa3e7953afc05fec6ad8e5d86d54a3acfbebea376fb65"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib\n",
      "%matplotlib inline\n",
      "from matplotlib import pyplot as plt\n",
      "from itertools import izip, izip_longest\n",
      "from block_gibbs import tree_iter, neighbours\n",
      "from numpy import random"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ImportError",
       "evalue": "No module named compat",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-3-6d96d1e44aff>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmagic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mu'matplotlib inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mitertools\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mizip\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mizip_longest\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/usr/lib/pymodules/python2.7/matplotlib/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m    155\u001b[0m \u001b[0;31m# definitions, so it is safe to import from it here.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    156\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mis_string_like\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 157\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msubprocess\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    158\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    159\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mImportError\u001b[0m: No module named compat"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size = 30\n",
      "theta = 0.25\n",
      "num_iter = 1000\n",
      "vis_step = 100"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\psi_{0, 0} = \\exp(\\theta x_{0, 1} x_{0, 0} + \\theta x_{1, 0} x_{0, 0}) $$\n",
      "$$p(x_{0, 0}) \\propto \\exp(\\theta x_{1, 0} x_{0, 0} )$$\n",
      "$$\\psi_{1, 0} = \\exp(\\theta x_{0, 0} x_{1, 0} + \\theta x_{1, 1} x_{1, 0} + \\theta x_{2, 0} x_{1, 0}) $$\n",
      "$$p(x_{1, 0} | x_{0, 0}) \\propto \\exp(\\theta x_{0, 0} x_{1, 0})$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rng = random.RandomState(1)\n",
      "vars = rng.binomial(1, 0.5, (size, size)) * 2. - 1.\n",
      "probs = np.zeros((size, size))\n",
      "\n",
      "def do_ancestral_sampling(vars, iterator, theta):\n",
      "    cur_iter = iterator()\n",
      "    next_iter = iterator()\n",
      "    next_iter.next()\n",
      "    for (i, j), nind in izip_longest(cur_iter, next_iter):\n",
      "        neigh = neighbours(i, j, size)\n",
      "        if nind is not None:\n",
      "            neigh = neigh.difference(nind)\n",
      "        energy = 0.0\n",
      "        for ii, jj in neigh:\n",
      "            energy += vars[ii, jj]\n",
      "        prob_pos = np.exp(theta * energy)\n",
      "        prob_neg = np.exp(-theta * energy)\n",
      "        prob_pos, prob_neg = (prob_pos / (prob_neg + prob_pos),\n",
      "                              prob_neg / (prob_neg + prob_pos))\n",
      "        probs[i, j] = prob_pos\n",
      "        vars[i, j] = rng.binomial(1, prob_pos) * 2. - 1.\n",
      "        \n",
      "        \n",
      "ret = []\n",
      "plt.figure(1, figsize=(15, 6))\n",
      "for i in xrange(num_iter):\n",
      "    do_ancestral_sampling(vars, tree_iter(0, size), theta)\n",
      "    do_ancestral_sampling(vars, tree_iter(1, size), theta)\n",
      "    \n",
      "    if i % vis_step == vis_step - 1:\n",
      "        plt.subplot(250 + (i / vis_step))\n",
      "#        ret += [vars.copy()]\n",
      "        #plt.figure(1)\n",
      "        plt.imshow(vars.copy(), interpolation='none')\n",
      "plt.savefig(\"sample_block.eps\")\n",
      "plt.show()\n",
      "#        plt.figure(1)\n",
      "#        plt.subplot(121)\n",
      "#        plt.imshow(probs, interpolation='none')\n",
      "#        plt.subplot(122)\n",
      "#        plt.imshow(vars, interpolation='none')\n",
      "#        plt.show()\n",
      "    \n",
      "#for image in ret:\n",
      "#    plt.imshow(image, interpolation='none')\n",
      "#    plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAFlCAYAAACJNVsDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3W+MZFd55/HfBDvqASO32ERmsNx9ySqsULKRjbRsVnTk\n0S6LzBsS3iRCimQlI8QLAoi8CEyk8R/mDYkE4h1vsJEDEVEEwju8IMFY9MStFUQktjFgAqxcozG2\nx0jJKFiaTJyk98WtdlffOlX11DnPOfecme9Hak1PTdW55977u6funarnXAkAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAgCR3SfqBpB9J+sjIfQGsyC1aQ2bRGjKL1pBZNOVVkn4sqZN0\no6QnJL15zA4BBuQWrSGzaA2ZRWvILEbzc5Gve6v60E4kvSzpLyT95uwTbr31v+5L4oeflJ9d+VqZ\n2+3x15mftn925YvM8pP7Z1e+VmZW2h57nflp+2dXvlZmdmPjP4+9zvy0/bOrBW5Y9A8r3Crp4szf\nn5X032ef8JOfPCXp3unfdiWdDDZ0r+6P7MKhxa1L97/Sh/QlWPsaWuay1/atr/86ez8OljAvZZ0O\nHbY/bC92++/sbGlv79SdUS9ebGVuL6hP7a4WZ8pq2f5c1r51m4X2nWWfW1kym8ts+8Pt4T1meLW/\nsbmp05cvN53ZZdZtf/1jf90lrGu99tcfG9drf32H7XscAwful4pnNpRa33Wy7ZPYZYbaz9P/w2Wk\n93WRvv1Yq/q1rPX487P7i2f2X/7l/8nrnDb3Plmtb3/5OcR6ZtvqW09jOefMfU7uadk4G3shtm97\n2u70z8n0p4tcHK4Pk+mPdOHCzTkWYMrt7rQXu+oT2+XoCa4JEx0kVrrhypUciyCzcDXRYWYzWeP8\nYKLD1AKLTJQ5tZzTwtVE9sTGXoj9RNJtM3+/Tf3/IAycnP65KwKL1Tod5GR7e0sXL57zXoAptyeV\n//+jcG3odDiybRw/rkevXvVeBJmFq05H343P+y9ijfODXZFarNYpc2o5p4WrTvbExtaIfVvSL0+X\n8/OSfkfSkrPmbvE/Ocjbev4l5G29xBJyt+/GnNsuc0dyt996Zltv3xGZpf1C7btZ4/ygy9yV1tsv\nsYy87edt3U1V57Stt5+39fxLyNv6vNhPxP5N0h9I+mv1s808IOnpxU/vFv5L6LubxwKfEt+jj61s\n3ace7FC+73zbfUz3zD22bFvM8ux/bPuhvobWadjWlna0F9+1RUy5ne2fdVuHWNbd2n7oecsz1SX1\na5XFrafxPoa92k85DhOZMrtqvaxjau7tH5J7nNrXsbnHLGNQiDUHludZj7nc2ycDU2Z9amlW6dZa\nXkr7IbHjRu59HpvZRc8b8hxHwnVM7kyZjT2G19dJCq+7dTyztO/J41zerlvr2eseT+u1ni72QkyS\nvjr9AVpCbtEaMovWkFm0hsxiFLFfTQQAAAAAROJCDAAAAAAKS/lqYjaFai/WlvK959jvWlu3hWcN\ni3U9Lc8LfUe49pqHMzr7yu8p29Xz+9GhtqyGyzyrM3PPyb1Pxqh19KxTqHVMOlD7MZWDd02CZ/1g\nyvOuV5btn1IzU0v9aWw/xng/iX3/9qhDrlHK+7CX0Pu3pXZtjPeI0LFZi9zH4Xz7i7c/n4gBAAAA\nQGFciAEAAABAYVyIAQAAAEBhXIgBAAAAQGE5K+n27xk0n1KsOSyQq+VmpFbhmxLGFSeHhNqqdZID\nSw52dra0t3dKypvRkH05rX9scXXNNyYPGfY3300c7X1IEXvMbW5u6PLl09IImfVa+1rG0NwZ8rkp\n6npKF9Tb9+X9UgWZtewTay5y7998N9pN78eQ9zFtORcLsW4f202Tx8ms5dzAsk88bzLvLfZc0rpO\nsWN7qH3LxCVS/DGQsv0XTNYRzCyfiAEAAABAYVyIAQAAAEBhXIgBAAAAQGFciAEAAABAYTfkbDz2\nLuS13JV9WKhXc4GlZZtZi2o9t/8Y+y3V7P5LKcCOvat8yvYvPTFHyBmdnXvMUpS9iCXH1vW27M9a\nJqzwlLL9xxB77ITUMPmQ9Xkp+6n13Fq2jzUXseceIdaxJbRMz0kOahF7bFrXyXPfebNNJDIv9ti0\njgex53qx/Qi1b81/KD+h7WOZFMZzPVMmxkvFJ2IAAAAAUBgXYgAAAABQWOpXEyeS/lnSv0t6WdJb\nUzsEZDYRmUVbJiKzaMtEZBZtmYjMYgSpF2L7kk5K+sf0rgBFkFm0hsyiNWQWrSGzGIXHZB3udzcf\no1jTUqhn7Ze1wDJ2MpBQW9Y7jFvEFh2HVFqIbM6s57ZIWWaI9Q7ylqJXTymTF+QuyPXM3jAHG9rU\nabfW5yzNbOsTNbQi99iVkvXYcSnjaFYks54T0aSM7Z4TK+R+j7FOouA5GdkYk5hFWJpZy3up51hs\nndgilrUtzwlarMfmcJm1vMdZJ2wa5n1Z+lNrxPYlfV3StyW9N7EtoAQyi9aQWbSGzKI1ZBajSP1E\n7G2Snpf0i5IekfQDSY8d/OPuzBO7xAXhejCZ/kgXLtycayFLMyvN57bL1RM0b6KDxEo3XLmSazEr\nM0tqYTXRYWYzIrNwNFGB1K59brCdu0do1kT2xKZeiD0//fOnkr6svrjxleCeTGwc15tOB2/G29tb\nunjxXI6FLM2sRG5h1+nw9HHj+HE9evVqjsWszCyphVWno5c85/MshszCUacCqV373KDeO55hbJ3s\niU35auKrJb12+vtrJL1D0lMJ7QG5kVm0hsyiNWQWrSGzGE3KJ2K3qP9fg4N2/lzS12afEFtcZ3ld\nLYV7IbGTboQmaYgtbKzFGBOvLLEys0OeheAprAWiIZ7ZSOnHUMrkNF6vsxqx0HztzIbUXPxv5Tnm\ne+Y4xLI9vCf+sEwCUYhLZq2s45tln1smO1rUVko/VvVr0es8j4ncE1ONsU5rMGU2Z99iJ4GQbGOL\ndTyIPU6sPNdpDDmOk5QLsWck3Z60dKAsMovWkFm0hsyiNWQWo0mdNREAAAAAsCYuxAAAAACgMC7E\nAAAAAKCw1OnrlxoW/dVSbBeSu0DUUrhuLab07Kv3JAcVFYy7SOl/7H6yFofnniBhjIlKYo+TMzo7\n95jncWJpa1Mbkk67LbMlNU+eFDLGe9HwuE4p8LZs71omGsot94QPKZMdWfbB9bKfrNvRclyMNnXS\nQOykP9Zj35rt4dhizY+1H8PnpZxnWPtWeqK3lOPQNt4s3tZ8IgYAAAAAhXEhBgAAAACFcSEGAAAA\nAIVlrRGr7Ga+o8p9k7xYnjfVXdSepa2a6gdn+2Lta+zNEr3rAzxvNhjal2PcVHfYj5SbrtaUs1aw\nHZeLPb5ScuzZj9aNkUXrOBg7ho4xxsW+F13Pdea5z3FrfV87qzNzj6XUrFv2sXdNfOwcFjnO0/lE\nDAAAAAAK40IMAAAAAArjQgwAAAAACuNCDAAAAAAKi79T4Wr7w5I2z5vNWguTPZdZQ5HkIrabIMYX\nGeYuBB/2bWdnS3t7p6S8GQ3ZV8R28izYD0kp4o8Vm5fcN5r2Pg4tRbuW/buxuanTly9LI2S29DQ/\nY0wsFDvRjecY7TkOWo+J/GOvpAoy6zlGxLLmx3OM896/w36kbMMxJoCx9f9+qYLMhljGRu/38+E2\nC02mMcakJ7nPb1Leh4bLPKOzc89JuUn1/GsXZ5ZPxAAAAACgMC7EAAAAAKAwy4XYg5IuSXpq5rHX\nSXpE0g8lfU3Spn/XgGhkFq0hs2gNmUVryCyqY7kQ+6ykuwaPfVR9cN8k6dHp34FakFm0hsyiNWQW\nrSGzqM4Nhuc8JqkbPPYuSXdOf39I0q4C4Y29Q7qlEHaMYnHr3cRDck9WkHt7hNofYwIJo+jMDo1R\nrBySsl0tx2Hu/ISyElvw7ltAG38cDrfZpjYknY5qS4mZXbX/vI9Vy8QZ3obrMMaxOcb7TsXcxlnJ\nNkZY30dTJt2wtDXG++31IpQDR66ZjeV9njTcZqHjJOUcLncePSfmGON9wTJ2LetVbI3YLeo/3tX0\nz1si2wFKIbNoDZlFa8gsWkNmMSrLJ2Kr7E9/Av+wO/O3zmFRuLZNpj/ShQs351zQwsz2dl/5bSKS\ni2UmOsjslSsew+lC5sz2ie1y9gUNm+ggsdktzezuzO+dSCyWmahQasksXExkT2zsmcMlSa+X9IKk\nE5JeDD3pmE5GNo/rU6eDoW17e0sXL57zbNyU2d7JmR6d9+wDrjmdDjJ7/PiGrl591LPxqMwCy3Q6\negLpPMKZM3vSd7m4pnXKmFoyC3ed7ImN/WriOUl3T3+/W9LDke0ApZBZtIbMojVkFq0hsxiV5c7k\nX1BfyPgL6v/n4B5J/0fSX0raUv/p229Lujx4neku5KGiQstdwccoZvUuAiy9DrGTIyx6rafhttjZ\n2dLe3inJltGh2MxK0r6c9oulENZzEpda5J7ExbrNQv04o7Om18bY3NzQ5cun+0WvLymzY5f21zKO\nWOWePMnzGPY8nsLvOff3i1lf05kN8X5Pjs17Lfms4XhtLbOxGap48rNRJvOKPVeyZtZzEjNL/+8/\neGqA5auJ71nw+NsNrwXGQGbRGjKL1pBZtIbMojqxX00EAAAAAETiQgwAAAAACuNCDAAAAAAKiyl2\ntNo/ZiiGC/EsDLQW7lkKAVOKKWPXyVpom7tgPLQdLf2IXe/EyTpSuE3W4SllopWWWI4xa/5Lb5/E\nyTpSZJ34wLq9h0LH/hjF/9b+jzEBVCzfCTymTZY1N87aCt7j91Hu9nNnu6XxvuIJZlJknawjJDZT\n3sfJcJKrlMlecufYmr3S+2nZOMsnYgAAAABQGBdiAAAAAFAYF2IAAAAAUFjWGrHYb2B6fncz5bvK\nw++yWr/HGrvMlO/OWr6LW+v3yUPGrBG7Z2aRKfV4ludZ2/KsdayF5/ffa9g+NdWIWXJmZdneY9RX\ned8ANXa8zF2fG7tMe63E9OllzdWI1XvzYJvcN+TNfYzlrufx3L+1ZDa32ExZ3yPHOF+w1KBZ+zHG\n2BsS2rbz+2BxXSOfiAEAAABAYVyIAQAAAEBhXIgBAAAAQGFciAEAAABAYTfkbLyGG2SGihGt/RoW\nFeYuJg4VMXq2l1KsGTtBRa0TRSwzu91S9rll3UP7PGWZlgylFLha+pZSABw6NmOPw9gbEYfUUhRs\nFdoHFtax0fK8lIloLBmyrqPneF/LZDje48YYLDdTHq5nykRJuaW8l1p4jkHWrFgmi0gZ2y3LHGNf\npvB8j/TM1BjHiTWzZ3Vm7rHhdqxl7A1JvdbhEzEAAAAAKIwLMQAAAAAozHIh9qCkS5KemnnsPknP\nSnp8+nOXe8+AeGQWrSGzaA2ZRWvILKpjuRD7rOaDuS/pk5LumP78lXO/gBRkFq0hs2gNmUVryCyq\nY5ms4zFJXeDxlZXSnsV2sYXIKcv0LGSMLUS2FjvGFmJaJwixbrPY4urh/t3SjvZMSwyKzqx0dF/F\n3tneyrrfLPlZ9FrLseNZCG7tl7XAdXgMW7fFGZ2NXqalfWdJmR3uqxombkjZZp7b23OCB2u2PYvI\nPfdlOP/R7Sdl1rKPh8dwaLtaj/PcY7k1B8O+jZGplO1j2W/eE4M5SspsrOF6pqyjZTvGZtGb56RC\n1rZqntRjkZQasQ9IelLSA5I2fboDZEVm0Royi9aQWbSGzGI0sdPXf1p65ZL7rKRPSDo1fNLuzO+d\npO3IheH6MJn+SNLNFy54N2/KbG/3ld/2ta1jwf9AA6TZ1F654n43EHNm9+dGW2CRiQ5HWnfmzO7O\n/N7l6g2uCRNlTGzkuUGf2i5bp9C6iaypjT1zeHHm989I+kroSScHf8/+5R40rdPhsLa1va1zFy96\nNm/KbO/kK78V+EoamtbpILXHj2/o6tVHPRs3Z/bY3GgLLNLp6Ankec/GzZk96blUXNM6ZUxs5LkB\nsFwna2pjv5p4Yub3d+voDDRAjcgsWkNm0Royi9aQWYzK8onYFyTdKekXJF2UdK/6/xa4Xf2HXM9I\nel/ohbkLAS1yTwaSIlRUOGS54/hYLP0I9d8yccOOtqT46TqiM2vpa4jnZDKxk1ik8CyqDUkplq0h\n7wXGsqTMrprQJGWiAktBtHcheOxrU4r/Y7dZ7k/Na3gfXSAps5b1Gm7b3JNvpUySEduPlGPH0t8x\nxk/PzDpPMON2brCO0pPC5B6Tck8il2KMiTksx9iyZ1guxN4TeOxBw+uAsZBZtIbMojVkFq0hs6hO\nyqyJAAAAAIAIXIgBAAAAQGFciAEAAABAYe43vhmTtUDUWrxaukjae3kxxdBS/ruyWyaG2NJO/FQd\niWbXP7SOKRMfWNq6Vu4WP8uaKVvRa97jMravG9rU6Wy9SuM5sYs0PwlPyiQH3uN26bas7VsmOAmx\nblvLJBCt8TzWLdss5f3QariMWibTiH1f894+NYvdV7E5Lr08K+tkIN7vOxax509jnnfxiRgAAAAA\nFMaFGAAAAAAUxoUYAAAAABSW8wuc+/dEfic+pIbaESvr93pjawZqZqlTsHwXfWtnR6f29vqnl7Wv\nFTnK/Z1+600QrXkfbu/YG9x6t59b6bq6zc0NXb58WroOMxviPR575myMG40OpdTG+Wb7fqnCzKbw\nrEkNCe27Wo+xlPzUe941TmYta1q6HizEc/+G2vM+l69hPLaKz/vizPKJGAAAAAAUxoUYAAAAABTG\nhRgAAAAAFMaFGAAAAAAUlvWGzpbias/iee+i5tzFt7m3z7D/3gWcsca4yd86ShZdlyh0zr29LUW1\nnsW41rZib87a4g2dvTJbc9F0bD9qmcRpKGWfWW+o2pIxJruwHPue7VuXMcYNecfY/q2pYSKOkNh+\neU7IZT1XtS4z9tgc4z3MNqnNYnwiBgAAAACFcSEGAAAAAIWtuhC7TdI3JH1P0nclfXD6+OskPSLp\nh5K+JmkzVweBNZFZtIjcojVkFq0hs6jOqguxlyV9WNKvSPp1Se+X9GZJH1Uf2jdJenT6d6AGZBYt\nIrdoDZlFa8gsqrNqso4Xpj+S9JKkpyXdKuldku6cPv6QpF0ZgjtGgbFnUWruIsDc2yd28oJFz7MY\noQDYNbPe+zy2qHaMSVU8eR4nKW1ZjrGRCtmTchtT2Jx7QoCUgu7Qa2uddCO3itfbday1rKf1Pcwq\n9/vaGZ2de8xyrKa8L8dO8JU7Z96TqUVyzaxV7klhhrzzY+mv92RzlknGKh4b17JOjVgn6Q5J35J0\ni6RL08cvTf8O1KYTmUV7OpFbtKUTmUVbOpFZVMA6ff1Nkr4k6UOSfjb4t/3pz5zdmd+76Q+wyGT6\nI0k3X7iQ2lxUZqWjud3XRMdILhaY6DCzN1y54tFkVG7350ZbYJGJDlPrInKs3Z35vRO5xWITkVm0\nZCJ7Yi0XYjeqD+znJD08feySpNer/4j3hKQXQy88aewEIB0d1ra2t3Xu4sXYpqIzKx3N7d8w0GKJ\nToeZ3Th+XI9evZrSXHRujzHawqzT0RPI8ymNJYy1J1OWi+tKJzKLlnSyJ3bVVxOPSXpA0vclfWrm\n8XOS7p7+frcOwwyMjcyiReQWrSGzaA2ZRXVWfSL2Nkm/K+k7kh6fPnZa0scl/aWkU+o/ffttz055\nFsZ7TrZgKR70llIAaSkQ9S5+roBrZq37PHfRqLWo2ZIX7+L24bqnZCr3hDghlUx6kpRbr7Ep99hr\nFTthyhjHibX9GgrLQ9swYY+7jrWWfZKS89zHuTUHllzlnjDIM4sp/bJMnlRTZr3GoND2T3nfjN0H\nnjmztmUdo2u4Dig1Zq+6ENvT4k/N3u7cF8ADmUWLyC1aQ2bRGjKL6qwzayIAAAAAwAEXYgAAAABQ\nGBdiAAAAAFCY9T5i1Ukpoost3BujqN9SzLrIsChyjIkQQiz7bkdb6r/OXV7JovqUyVishevD51kL\nhWPzfkZn5x7zLk7GUTH7Kvd4lpJti5Tj9KzOzD3W0iROIZVMOmOWs7+eY7h1nIod97zHvGF7NUwS\nkyLc/3qzbnkvtU5YYd13lkzVMj54vgekTAZimWSs1DbjEzEAAAAAKIwLMQAAAAAojAsxAAAAACiM\nCzEAAAAAKGz0yTqsBd2eBa3WYj7PZdYwUUZKUbm1/7GFwcOiyC3tjDRVx1G5izVTJmOx9m24n0L7\nMuWYiN1GtUxycD3IPVlK7glgQnIfmymTHHiNg5L92Bw+Zh+f6yjiD6lhggHrODVmsf+yfnhPmmNZ\nJ88JQkLLGz8VPet7pOXYTDn2a30vjd0+Vp45G3Mb8okYAAAAABTGhRgAAAAAFMaFGAAAAAAUlrVG\nzPL9zdw3dbV+HzX0/VDLDd+sfbXe7NGTpTbI+h3e3PVsw20x5g2dZ3nX9uWuGbD0d4y6hZRlWr5f\nbzXGDVZrZRnzpPhxL+U797n3U+xxnfsm1SG11B6NIfbYt+6na3E75j6v8ByPLWqqa7RsW89zhpQa\n8lix+Qn1NXTee1Zn5h6LratL4Xl+nIpPxAAAAACgMC7EAAAAAKCwVRdit0n6hqTvSfqupA9OH79P\n0rOSHp/+3JWpf8C6yCxaRG7RGjKL1pBZVGdVjdjLkj4s6QlJN0n6O0mPSNqX9MnpD1ATMosWkVu0\nhsyiNWQW1Vl1IfbC9EeSXpL0tKRbp393ufuZtWDcInchtfcNkYdSbm4dep6l0DPlRr65i4IjJWd2\ndl1D26LS9ZYULo7NXfSau/1hRq3b3zrJxLD/I918PSm3q/ZB7puFek96YNlPKROLxI7luSfmyJ09\n55vjZj8/iF13636yjF0p+8Rz0pnYibW8b6obOx57jhFjZdaybS1ji+ekSFalJ4cr9bxYlrz7Tji1\nOLXr1Ih1ku6Q9M3p3z8g6UlJD0jaXKMdoJROZBbt6URu0ZZOZBZt6URmUQHrhdhNkr4o6UPq/xfh\n05LeKOl2Sc9L+kSW3gHxyCxaRG7RGjKL1pBZVMNyH7EbJX1J0uclPTx97MWZf/+MpK+EX7o783s3\n/QEWmUx/pAsXbk5pKCGz0u6RrwhMRG6xyEQHiZVuuHIltbno3O7PjbVA2ESHmXUQndndmd87kVos\nNlEtmT08N+j8+oNr0kTW1K66EDum/mPa70v61MzjJ9T/r4EkvVvSU+GXnzR1Auh1Ohjetre3dPHi\nuZhGEjMrnZz5XvDfMNxiiU6Hb8gbx4/r0atXY5tKyu0xxloYdTp6Enk+vqmkzJ6MXy6uM51qyaxL\n6SOuC52sqV11IfY2Sb8r6Tvqp/SUpD+W9B71H+HuS3pG0vui+uksdyG1VWxRZGyBsXc/rAWKluJb\n64QAjpIzu2qSE89JYbwLaC39SCkUtvR3pMku3JZpmfAkpf0FknI7PM5iJyHwNMakP577JGWbxb42\nJXvx2zF66oOkzFr6GzuxRUjsZBEpE1pZ2vN+D7BM0uV57KccJyO8VyRl1pLH2Am+UiaGsPB8D0iZ\neKXmyc5ipa7TqguxPYXryL6atFQgHzKLFpFbtIbMojVkFtVZZ9ZEAAAAAIADLsQAAAAAoDAuxAAA\nAACgsJzVgfv3rCggT2EtPLQWFeaeVCLUj1UF9ilthdpL2RZndWbusdii2lBh47BvWzs7OrW3J+XN\naMj+sZntZi2eH4N1gg1LzlKKb2PbtxZqxx4n1u1jYSnG3dzc0OXLp/tFl7U/dvlzzQXYnjkY433C\nc9uG25c0Qma9zg9S9m/NubWwjNsp65h7gof44+l+aYTMyikvY0yeNIbYcxQp//YY9sN7ecO8Lxtn\n+UQMAAAAAArjQgwAAAAACuNCDAAAAAAK40IMAAAAAApbdUPnJDmL7bzv3B4qKhwW26UUFOYuCg71\nP1ZoYg7rpBVekyHsaEv9vRfLm93PKRnOXZAbymOIZzY827dmJXf7Fpai+A1t6rTbEtczzJXnurck\nZdIZi9zHUuuTR6wjZiy0TvgQ2o6WiQOuhQkTPDMUu82s733XyzhlGZfGmMBjjHOU3GNc7HuA9wSA\n6+ATMQAAAAAojAsxAAAAACiMCzEAAAAAKGz0Gzrn/o5q7pqB6+k7/TFib7Y55g2dY27aGHuz7Nw3\nEQyx3FB7ndfGis1GqA+h13nehNxiY3NTpy9flq7DGzpb5b5xp+fNm62sefRcT9/2x7k5rldmrWPS\n9XIT3aHY9ybJNobG1uMtYtuf42TWck6b+5zTooYbJC9aZu6bhFtZ6qpTzpXm25I05g2dJ5kLnKVJ\nw62XWELr7Y9hkrn1vMfEJGvr+ZeQt3USW2P7rWfqWtgDpU2aX0Lu9vMvI2/r0v41mVvev5e3zjn/\nOgpdiLW9hLytl1hC6+2PYdJw6+1nNm/rJLbG9lvP1LWwB0qbNL+E3O3nX0be1sssobRJ4+2TqTFb\nn0eNGAAAAAAUxoUYAAAAABSWs9hxV9KdGdvHte+8pJOFl7krcot4ZBatIbNoDZlFa8bILAAAAAAA\nAAAAAAAAAAAAAAAAAACgcXdJ+oGkH0n6SKZlTCR9R9Ljkv7Wob0HJV2S9NTMY6+T9IikH0r6mqRN\n5/bvk/Ss+nV4XP12i3GbpG9I+p6k70r64PRxz/4vWsZ98lmHGuTO7URkdlbu3JLZdBOR2VlkNh2Z\nXd3+fSKzNSGzq9u/T2S2Gq+S9GNJnaQbJT0h6c0ZlvOM+p3i5Tck3aGjwfpTSX80/f0jkj7u3P69\nkv4woc0Dr5d0+/T3myT9g/pt7tn/RcvwWoexlcgtmT0qd27JbDoyexSZTUNmbe2T2XqQWVv7ZHZN\nOe8j9lb1oZ1IelnSX0j6zUzL8pyG/zFJ/zR47F2SHpr+/pCk33JuX/JZhxfUDw6S9JKkpyXdKt/+\nL1qGlPd2CKWUyi2ZPZQ7t2TWB5k9RGbTkFlb+xKZrQWZtbUvkdm15LwQu1XSxZm/P6vDFfS0L+nr\nkr4t6b0Z2pekW9R//Krpn7dkWMYHJD0p6QGlfVR8oFP/PxXfUr7+Hyzjm9O/e6/DGErklswu1ilv\nbg/aJ7PrIbOLdSKz6yKzdmS2DmTWjsyuIeeF2H7Gtme9Tf3Ge6ek96v/qDSnffmv26clvVH9R6TP\nS/pEYntsXA9RAAAgAElEQVQ3SfqSpA9J+tng37z6f5OkL06X8ZL812EsJXJLZsNy55bMxiOzYWQ2\nDpm1IbP1ILM2ZHZNOS/EfqK+EO7Aber/B8Hb89M/fyrpy+o/PvZ2Sf13SSXphKQXndt/UYdh+ozS\n1uFG9YH9nKSHp4959/9gGZ+fWYbnOoypRG7J7LzcuSWzacjsPDIbj8zakNl6kFkbMrumnBdi35b0\ny+o/7vt5Sb8j6ZzzMl4t6bXT318j6R06WjTo5Zyku6e/363DneXlxMzv71b8OhxT/zHq9yV9auZx\nz/4vWobXOowtd27J7LzcuSWzacjsPDKbhszakNl6kFkbMluZd6qfheTHkk5naP+N6gvtnlA/9aTH\nMr4g6TlJ/6r++8C/p34Wm6/LZ7rPYfu/L+nP1E9Z+qT6QMV+33VH0n+o3x6z02569j+0jHfKbx1q\nkDO3ZHZe7tyS2TRkdh6ZTUdml7dPZutDZpe3T2YBAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAa8Bdkn4g6UeSPjJyXwArcovWkFm0\nhsyiNWQWTXmVpB9L6iTdKOkJSW+efcJ/+9Vf3ZfEDz8pP7vytTK3W1u/NvY689P2z658rczsr21t\njb3O/LT9sytfjLP85P7Zla+Vmf2lX7p97HXmp+2fXS1wbNE/rPA/JN2r/n8QJOmj0z8/PvOc/Xun\nv+xKOhm5oEXu170zf1tvCffq/hXtDS1uP9RWSKj9g9cetL68D8stX6d+Cda+Wg3bX9SPxa87NHzd\n1s6OTu3tSfEZDTHltn/Krkpts6H18xly2H7sOliPidh9vtxh+0Mpx5yl/Vibmxu6fPm0NEJmDxMb\nZjnmFj3v0OIlpBwn1nFk2Lf1c7C8/VSHrecx237Ke8W8+6WRx1lpnH1i2Y7L+rWq/Vih85vY94VV\n27Vv3dIPW3vzr+uX4Ps+NGZmpRx7fnb7HLTue5zPOljC8n4cWHfsXdx62BjnB8v70LfvOSZNWwpm\n9uci27xV0sWZvz87fQyoGblFa8gsWkNm0Royi9HcEPm6fcuTdqd/TqY/XeTCcH2YTH8k6eYLF3Is\nwpTbPrkTSbvkFitMdJDaK1dih9OlzGPtZPpnJzKLZSY6HGmzWHucJbFYbqJ6MitxVotVJrInNvbM\n4SeSbpv5+23q/wfhiJMzHeoiF2STt/XWe9/K9ulmWtra3ta5ixcXPzmOKbd9cieSOnU6792HqS5T\nu7Rftv3ulbaOH9/Q1auPOrX7CvNYO1HrI1Xb7edt3bP9btCa+xi39jibS76Wy7RfYil5W/daQqd6\nMiu1nlsylb/9TvbExn7H9gZJ/yDpf0l6TtLfSnqPpKdnnrPv+y32XN+VtcvxHfZZ+WrEFj8nRez3\nxy0y1YiZcquI/XAs8B9q9+hjUZ0M2Q9sho/pHtNrY/eJdZnDdQ+td0r/Vy1v0TJDLP0ItWXpa6Ya\nsaix1rJOIdYxKE/d4SHP4ym0TEv7oX3ueUynqLxGLGqczf3+GjLcjrHHvpT/PSDEc1yN5b3elrFr\nWb1NJFNm7xkssvS29uZTb93L/b5sNcY4YnmvWJbZ2E/E/k3SH0j6a/WzzTygo4EFakRu0Royi9aQ\nWbSGzGI0KUUNX53+AC0ht2gNmUVryCxaQ2YxithZEwEAAAAAkbgQAwAAAIDCssy3XEKoIM+7cHW+\n2C6+QNqzoDVcvFp+MhNLYXMtxe3rWFXsac2ZZZIA6/YJ5cdq2A9rFkPPiy2E9SxqPqOzc4+lTAJh\nWSfLcza0qdOmXuQXyqin3OON5biwTqZhnYDBc51iJ/rwXKeUiSdKsBxTpd9PUsZZ6zEXu5+sY/RQ\n7nOlUFuxE+QseqwWNR0/ucSOgyn7LXbypJCU7FnkmMSJT8QAAAAAoDAuxAAAAACgMC7EAAAAAKAw\nLsQAAAAAoLAqJ+uILRZMmVwg9x3qPQtQrQXFsRMrjDHxR0usObNkKqXwMzbboWXmKEBd1ZZ1mcPn\n5c5/iOWY2NSGVM10He2I3U8pk7ZYWPtlXWZs32Jfdy1MLFDzxA1DseNx7vyEeGYj92Qg4fX2G9vX\nYZlAotaJdKzLtOzPlL6OcX4ZO3mVta+p5098IgYAAAAAhXEhBgAAAACFcSEGAAAAAIVxIQYAAAAA\nhcVVsNnsx5bkhQrkPAvvrcscSimw9Lyzd4il/9Y+WAttYydWCL1u2P+dnS3t7Z3qmyhrZW5T7txu\nKXpNybqlbym5y91+rJR9YmE5JjY3N3T58ul+0WXNZbbWCXdyj+O5J10KZSo0aUgNE2XYJxuRVEFm\nY6Uc+2NM7uI58YFF7sx6H9OWiS1az2xIzRMN1bDMlElhPDMau8xlmeUTMQAAAAAojAsxAAAAACgs\n9T5iE0n/LOnfJb0s6a2pHQIym4jMoi0TkVm0ZSIyi7ZMRGYxgtQLsX1JJyX9Y3pXgCLILFpDZtEa\nMovWkFmMIvVCTFpSMDksyrMW0MYW1qXc3doy8URKUW3snb1TWO4Cn1KsGdqfsftu+Lot7WgvqiUT\n951hzVlom8Wy7jtLcbjncWKVejf6VWLbsvZrmNkNbep01BJNlmbWq+g6dv9aX5d7n3seXyHWvloK\ny2uY0COztcZZz2zEFtTnnuxFKp9R63FoaWud18YaY7KnGaNlNnayl1A/vCdrshwnKcuMPW+0HkuW\n/WTtv+d574HUGrF9SV+X9G1J701sCyiBzKI1ZBatIbNoDZnFKFI/EXubpOcl/aKkRyT9QNJjh/+8\nO/PUTtJ24uJwLZtMfyTp5gsXci1mRWbnU9vl6gmaN9FhZm+4ciXXYlZmltTCaqLDzGa09jgLLDIR\nmUVbJrJnNvVC7Pnpnz+V9GX1xY0zwT05eHrej+TRtk6Hg9vW9rbOXbyYYzErMjufWmCRToeZ3Th+\nXI9evZpjMSszS2ph1enoSeT5PIthnIWbTmQWbelkz2zKhdirJb1K0s8kvUbSO6TlX5RMuVFw7PdP\nre2f1ZmVz8v93Wjv+glLf1O+l1zrzX2XWDuzuW8OmPsGtyEp9Yqex4BnNjy3o/Vmp4WybcrsqvW3\njr3W/Ru7vVPeAyxy38g6pa+xx461TmeMsWSBtcdZKX6fW7ePZdxLGd9yZ8+r/rqE3Od1hjitKyqz\nseNZ7qzkZp3vYYxxKnddY475HlIuxG5R/78GB+38uaSvJfcIyIfMojVkFq0hs2gNmcVoUi7EnpF0\nu1dHgALILFpDZtEaMovWkFmMJnXWRAAAAADAmrgQAwAAAIDCct5leF+VFiTGFgtaCyxzFyN6Fnpa\niy5DLIWqoW1h6f/Ozpb29k71XSxrP2dqY29ynpI9y41kPfuRkv/YyRBSbsboVUi9ubmhy5dP94sp\nKyqzofXMPTZa1dKPoZRsl76hs/2Ynj69rP1jg/X3PM5DapkYKXa/p7xXe4rNcey5QNj90giZtZzT\nWnLm+X4uze8D7wkrxpiELfZcydJWqL3c67RsnOUTMQAAAAAojAsxAAAAACiMCzEAAAAAKIwLMQAA\nAAAoLOtkHZ5l06ULnVN4FtWG1vOMzpqe58m6Tl6FvFs7Ozq1t9cvuqyVuY2dUCKF52QU3v3wLE4O\n8SyqDRU6n9UZl/Y3Njd1+vJlqYLMWoq3XTvgOCnBIpYc1DKhUojnRCKe/R9rsg6vrWHNXi2TXVhy\nYJ1sIXeOLe91uSfICRtnsg6vNfUep4btpUyMYpkgZPieKflOnBFqz/LevU4/SmOyDgAAAACoCBdi\nAAAAAFAYF2IAAAAAUBgXYgAAAABQ2A1jd8BqjMk5hkV/1j6EigpDLAX1KQWKw/Zj79Iu2ScIGfbD\nusxhAeeOtiTtzT2vBqFtETuZhndhqaU978lGck8MYT2eLEJFwZ7t16J0wXJou4ZYi89jJ//wncQi\n7yQE1gkZsFzs+6GV93hpOa8IHU+5J+aw5t3zuGjtGBiuu+cEMDWPXcP1tL5nWvthOSe0GmMytVR8\nIgYAAAAAhXEhBgAAAACFWS7EHpR0SdJTM4+9TtIjkn4o6WuSNv27BkQjs2gNmUVryCxaQ2ZRHcuF\n2Gcl3TV47KPqg/smSY9O/w7UgsyiNWQWrSGzaA2ZRXUsk3U8JqkbPPYuSXdOf39I0q4C4fUs4h8W\nMsYWcy/qh6VgNqWw1LOAMGXdh1KKgmP7by3sTxCdWYuUYuIxCnktRa+xd7uX5reHd6Fw7iL1odj+\nb2pD0unYxWbNbEjuombrfvPsh2dbKZMbWY4T7wkZck4ksED2zFrWyTqRVMzyUsW2l9KP2POW0uPs\nIrGTpBklZdarb6F9kvIebJH7HKIWuSffyvG+GVsjdov6j3c1/fOWpF4A+ZFZtIbMojVkFq0hsxiV\nx/T1+9OfwD/szvytc1gUrm2T6Y904cLNORe0MLOS5lLb5ewJGjfRQWavXMl6NxAyCxcTHSQ2OzIL\nJxMVSu2KzB7+U1egM2jXRPbExp45XJL0ekkvSDoh6cXQk47pZGTzuD51Ohjetre3dPHiOc/GTZmV\nRGqxhk4HmT1+fENXrz7q2TiZhbtOR08iz/s2T2aRQaeMqV0js9fePSeRRyd7YmMvxM5JulvSn0z/\nfDj0pJw3UbN+vz53PZX1+9eeN6yzsnyf2dp/z3050veNTZm1qOW79CGx3/23Hie56/vG2Laex4kz\nt8yGjHGTy5SxxZINa1ux/bDm3/K83OOg/cbzrseca2Zjb35sZdkHKe+RnuOG53t1LTexr6T2yJzZ\nYdZi369SbkJu2edj3BS79dpwqxzvm5YasS9I+r+S/ouki5J+T9LHJf1v9dN9/s/p34FakFm0hsyi\nNWQWrSGzqI7lE7H3LHj87Z4dARyRWbSGzKI1ZBatIbOoTuysiQAAAACASFyIAQAAAEBhWedbjr2h\ns2fRn/dNM3O273mz0EWPWVhvkGnZn9Z9OdxmW9rRnumV/mbXNXfxv3dbsUXY1oloPI9N63FSQ0F3\ngZuQXxdSxjiv1y3qh2fOcmc29j0m9Lo6S+J7pSc9sW5Xa84s/a91UgLJ930tdj+1llmLlAkfQmOX\nJWe5J0qysr4HxN6g3pozy43PrVKPEz4RAwAAAIDCuBADAAAAgMK4EAMAAACAwrgQAwAAAIDCsk7W\nMZz0wVrQlnuyi1pZ+xq6K3tomw23h7V9a1GkRc2FyIvErKu1+N9SIJpSNJp7go3YPlizZ+mrNVOx\n40GLmV213VKOacs+sbafu+g7xbC/Y7x3tJi9WDVMwuO5vXOf38T2Y4wJp1Jci8dA7D7xNMYyQ6zv\nFcNt5j2ZXez7R45xi0/EAAAAAKAwLsQAAAAAoDAuxAAAAACgMC7EAAAAAKCwrJN1eN652iI0iYW1\ncM/SN8vkC9a2FrVnETsxhLUI1rp9YiejuJ5ZCkRrKGJPETomQmLXM6VAN1aoreHxtKFNnXZboi/r\nmJG7fc+JgKys6+k5kYulsNzaVi3b0dtw/T0n24odD1ImrPCePMbCc9yz9r+GiSfGYll3z3HVMgmb\ndZm1TMbi2Y9rJXt8IgYAAAAAhXEhBgAAAACFWS7EHpR0SdJTM4/dJ+lZSY9Pf+5y7xkQj8yiNWQW\nrSGzaA2ZRXUsF2Kf1Xww9yV9UtId05+/cu4XkILMojVkFq0hs2gNmUV1LJN1PCapCzxevjJ1BWvh\nXuyd5q2Fk7HFiCnth8QWMp7R2bnHrAWow/6OVCDaTGZrmZgjZVKb4fOs6xSb95R8WljbGj5vUxtS\n/HQdWTObezITK+uEAJ6TD4WkvNYiNJlG7DKt+264zQpM6OGaWeskP0O5J6yoZYwOic1ZSv6H42/N\n2ycgKbPDfOSecCrE8zgZwxiT2uSW+n6SUiP2AUlPSnpA0mZSL4AyyCxaQ2bRGjKL1pBZjCZ2+vpP\nS69cXp+V9AlJp4ZP2p35vVP4vyGAA5PpjyTdfOGCd/OmzErkFuuY6CC1V6643w3EnFlSC6uJDsfZ\nDBhn4W6iOjJLamE1kT2zsWcOL878/hlJXwk96WRk47g+dToc1ra2t3Xu4kXP5k2Zlcgt1tHpILXH\nj2/o6tVHPRs3Z5bUwqrT0dPH877NM87CXac6MktqYdXJntnYC7ETkp6f/v5uHZ2Bpiq5b9CYu87C\n+j1/Sz9SbgwaW0NkqWUI2dGWpL2Vz1tDFZm1bJ/crN9nDu1zy03TresU+/16z2Oicm6Z9d5mln1u\nbd9ys2zJrz6jBM+b3Xvuk/A2c62XKz7OWrdP7rrA2Pa966gtx4V1PBjjGLMts+3MptQietaferLm\nOLavNayjlKcflguxL0i6U9IvSLoo6V71/y1wu/rZZp6R9D73ngHxyCxaQ2bRGjKL1pBZVMdyIfae\nwGMPencEcERm0Royi9aQWbSGzKI6KbMmAgAAAAAicCEGAAAAAIW5z7c8K3cxvkXKzeNy3xgxlvWG\nfsPiyZTifOs+GU7mcI1OojDHup6Wgv2U7eOZx9A+Dx1PsTeY9JSyzTwnnmhJ7kkOrO1bb1o/3Aex\nNwS38p5EwfJelDI5Qg3HYW65j8OUTNU8ecxQymQRNUw4dT2z5HGMmzfHjuNWNR9fqevJJ2IAAAAA\nUBgXYgAAAABQGBdiAAAAAFAYF2IAAAAAUFjWyTpiizotRW61FO6lFPZbiuBT7lY+LJ4cTqRh7cMi\nlskcrocCcil+PUP7xDM/ISmFpZZlphTU1zD5zRh98Nb6OsROshQ6DmO3RUrBu+ekErFtWdd7rKTE\nvIenjJexPCcCsrZnHaNznwflHkdamwRpuD1iJwfyXu8zOruyLet4FvvalKyMcT7vme3UyVH4RAwA\nAAAACuNCDAAAAAAK40IMAAAAAArjQgwAAAAACss6WYdn0eLweZbJKax9WMRyt/JhkeSiZcYW83kW\nWOYu1pTanyTAYoz967ldYydCWGRYaBuaMMGzODmlsPd6yKcUv41ix+PY1y163hj7ybLNcmfbukzP\niZ5aEhq7vCehilVT8f8qKX21nBelHPstiR0Pck9eZT0XDvW11snyUs75a3/f5xMxAAAAACiMCzEA\nAAAAKGzVhdhtkr4h6XuSvivpg9PHXyfpEUk/lPQ1SZu5OgisicyiReQWrSGzaA2ZRXVWXYi9LOnD\nkn5F0q9Ler+kN0v6qPrQvknSo9O/AzUgs2gRuUVryCxaQ2ZRnVWTdbww/ZGklyQ9LelWSe+SdOf0\n8Yck7SoQ3GHRX6iIroai5kWGhYApBYq5C6c9J2BIaWuMIs6BpMxKR9fBWoybmzU/sdvfM58pk+bk\nPk4qnsAgKbfDcW64D6zbPzbvodelPC+U4+FrvfdbrQXdsRMfhCaSCo8P0eudPNZ6CI031nHQkinr\nJBPW/bTqWJXGeR+NXWbKeZdl24baSjhSkzJber/kfl/OvcyQ2MldUrb9GO8n61inRqyTdIekb0m6\nRdKl6eOXpn8HatOJzKI9ncgt2tKJzKItncgsKmCdvv4mSV+S9CFJPxv82/70Z87uzMPd+n3DdWcy\n/ZEuXLg5tbGozPZ2Z564rWOkFwtNdJDZK1dc7gYSOdYe6iRte/QE16iJDjLrJHKs3Z35vRNnCVhk\nIufEkllkNpE9s5YzhxvVB/Zzkh6ePnZJ0uvVf8R7QtKLoReedL5fEa51nQ4Gtu3tLV28eC62oejM\n9k6+8pv1q1S4XnU6yOzx4xu6evXRlMYSxtqjSC0W63T0BPJ8SmMJY+3JlOXiOtLJMbFkFgV0smd2\n1YXYMUkPSPq+pE/NPH5O0t2S/mT658PzL40/ia31u/opPGu4PG9Yl3Ljxdjvv2e+iWNSZocqqR9y\nrSPwvrFm6eM15caOnsfhsB8b2tTp+OZcc+v5Hwix36/3rBORbHVvnnWTISk1o5aa6ZRjc7h9rDc6\nTjh6kzIb814xxk1dY7f/otcO93Fr/9nnOd5bjifnukbXcTa32Lou63FiHRs9j0Nr3ofjV8oyQ4Y1\ntN61cfPbdnFmV12IvU3S70r6jqTHp4+dlvRxSX8p6ZT6T99+O6KfQA5kFi0it2gNmUVryCyqs+pC\nbE+LJ/R4u3NfAA9kFi0it2gNmUVryCyqs86siQAAAAAAB1yIAQAAAEBhLvMtW3lPEtC64brnLjQP\niS00twoVjNeeg9n+eU8IMMaNOj2Lq2u9MaJnzmL30aY2pJTpOjKybgvrunvuY+sEGMP+phSfx2Yj\nZb0t7Vsn2LCoZaIhq9gbA8fmZ9HzPMVODpTSf8/30jEmTmtp8hLP7Z/SliXHuY8T637zzLF1zIid\n1MNyviPZJqtadiTxiRgAAAAAFMaFGAAAAAAUxoUYAAAAABTGhRgAAAAAFFZ0so7YYlyrlCJAS7Fd\nyl28Q2ILea08J4awFijGPEea39Zb2tGerWvuZtd1ePd1aZwJN2opyh6jeNsiNmfXilUT/7RUAC/l\n76/nsRM7QUioDynvCZb2a5rAw7LdYicCCuXHsm29x3bLeJOyn2InmBnjPSxWaBvWMorHntOGsmid\ncCrmuFnUluc4a81UaN1zXxvEToJkPee3rfvi5fGJGAAAAAAUxoUYAAAAABTGhRgAAAAAFMaFGAAA\nAAAUVnSyjhSritEXPRYSKsCz3EHb2n7uiRVyG6P/w+2/oy1ptOk6Do2x37wnlKhh0pYUlmN/jIkJ\nhvtpQ5s6nXWJiw3HtNh9nnsyE8/9FOpryoRKsetuzWPsWGLdZsNjM/dkUKliJi+xbgtr/i3Pq2U8\njh33Ut7DLO1bzp1Cr0tpv5bpOmLXPff7qPUc12r4WuvkGtbxODbHnsf+GOc2B/hEDAAAAAAK40IM\nAAAAAApbdSF2m6RvSPqepO9K+uD08fskPSvp8enPXZn6B6yLzKJF5BatIbNoDZlFdVbViL0s6cOS\nnpB0k6S/k/SIpH1Jn5z+ADUhs2gRuUVryCxaQ2ZRnVUXYi9MfyTpJUlPS7p1+veilcCliuZmndHZ\nI3+3FgFa75BuKaJNKRQevjalcLIhyZn1nNxiKPdkCJ5qmXTGayKKyiXl1qtAP7Yg2lroPBxTQ20t\nas9y7MS+LiR2koxFy7RMPBFiHbdHeI9Myuxwu8Vus9zjgfe2js2jZ7ZDYsd76+us+2n42tA6Jqx1\n9nPa0pNEhVjHjNznIymTetQ64VQO69SIdZLukPTN6d8/IOlJSQ9I2vTtFuCiE5lFezqRW7SlE5lF\nWzqRWVTAOn39TZK+KOlD6v8X4dPSK5f9ZyV9QtKp4Yt2Z37vpj/AYpPpj3Thws2pjUVltrc783sn\nkotFJjpIrHTDlSseTUaOtYf/89h59ALXrIkOM+uE8wNkNVEdmSW1sJvImlrLhdiNkr4k6fOSHp4+\n9uLMv39G0ldCLzxp6gJwoNPBwLa9vaWLF8/FNhSd2d7J2OXiOtPp8K144/hxPXr1akpzCWNt3feM\nQj06HT19PJ/WHOcHyK5THZkltbDrZE3tqq8mHlP/Me33JX1q5vETM7+/W9JT63QPyIjMokXkFq0h\ns2gNmUV1Vv036o6kv5H0HemV7778saT3SLp9+tgzkt4n6dLgtfvHBsWB1rtgx97h2vvO2zXwvNu3\n93pbiiJj98nOzpb29k5J6xfQpmRWkvbvmVmktf/W5w0nckkp7K0lx8OMehcre66n5XgKLc+yzzc3\nN3T58ul+MetLGmuHPa4hG7mLplMmUbC8NlRAbp2IySJlH3lu22lLo2fWwrrNxphoKPa92vMcyNp+\nimH2xsjxWJnVYF0t7/GLnjeUOwcpYidnim3fm6W/Ke8ntslR7u+7ErDqq4l7Cn9q9lVD34AxkFm0\niNyiNWQWrSGzqM46syYCAAAAABxwIQYAAAAAhVmnr48S+t7kUOx3uWNvDpiyzBSxN6u03ujO8zu7\nnvUHsTcW3NKO9tx6sZ7ZbWndFtbnWY6JkNzfobbeZLGGG8mmHL+x27/2m5wP+2f5/nvKPq/hppnW\nfRmb2VrqM0Is61TDsZoqdtwb43j13Lax41RIynEeYtknLb2frCP2ZsqxYrdF7u3qWYtrXUbKMW1Z\nZu73k2X4RAwAAAAACuNCDAAAAAAK40IMAAAAAArjQgwAAAAACvOrKpw3d/O7FMMiyZRi09w8JwPx\nLrS1SClSXzVpgNXWzo5O7e313Skr6kajIZaiZuu29rxppnd+PG/waWHdZqXHiMQbOqeIGmtrmHCj\nhFAea1j33MehVcLNcVO4ZbbmcwFPuSeP8Ry389+QevHNcTPaP+Z0fhOSMglb7AQqtRhjMjJLH0Kv\ni+/r4szyiRgAAAAAFMaFGAAAAAAUxoUYAAAAABTGhRgAAAAAFHZDyYVZJ7GwFNaFihitchf8ntWZ\nucc8J/CwFnDGtn9GZ+ceiy1QtN6FfNj+jrYk7UUtM9Vsnz2LcUNCufOcmCMklB/PY8KaxdjMWreZ\ndT0t7dde/D/clsP+Wgu1YwvGrfs3t9xF31aW7Zh7+9SyT3JLORewSBkbLa+1ji2hZdbK+32tFrH7\nIHYSszG24xhjVw0TiYzZh0KfiE0abr3EEmi/NvuNb7O8rZdYQt7287Y+jknmE9Lcx0Te1vMvIW/r\nZZZQ3qTp9nMfcwdLabf1Mksob9J0+3lbp/11cSFWxRJovz6TptvP23qJJeRtP2/r45g0voS8redf\nQt7WyyyhvEnT7edtvcxS8rZeZgnlTZpuP2/rtL8uasQAAAAAoLCsNWJvecsJSdJzz92kN7zhhE7o\nLfPP0Ym5x0LPW+am557TiTe8IdhWSKj90HdxD9o76L9V6Dvry9Z9WfvWbeHZf+t+OrCs/XW3xYE3\nvek/aW+cEjG95S0n9Nxzr12YWXM7kdssxUF/D46JFNZMDbfRsteF+rqq/VjLsrds+1j7P/Ta1/68\nztHoT58AAAQHSURBVJ9fs5NOTrzlLXrtzDoN+2vNcWg9D7bjwTERas9zHE9hPeZSjutFblqy/aX5\nPHqOs6H2zPvk7/9+YZs59eNs+ntTsG3HcWTZ+cKyY27Va4d9XcTj/GCZdc+f1pXjvW6kyJrOaYOv\nM4zHs8/Jsc1ml7lsrI0du2Z5jOW1th97nCzLbM7qz11Jd2ZsH9e+85JOFl7mrsgt4pFZtIbMojVk\nFq0ZI7MAAAAAAAAAAAAAAAAAAAAAAAAAgMbdJekHkn4k6SOZljGR9B1Jj0v6W4f2HpR0SdJTM4+9\nTtIjkn4o6WuSNp3bv0/Ss+rX4XH12y3GbZK+Iel7kr4r6YPTxz37v2gZ98lnHWqQO7cTkdlZuXNL\nZtNNRGZnkdl0ZHZ1+/eJzNaEzK5u/z6R2Wq8StKPJXWSbpT0hKQ3Z1jOM+p3ipffkHSHjgbrTyX9\n0fT3j0j6uHP790r6w4Q2D7xe0u3T32+S9A/qt7ln/xctw2sdxlYit2T2qNy5JbPpyOxRZDYNmbW1\nT2brQWZt7ZPZNeW8ofNb1Yd2IullSX8h6TczLctzGv7HJP3T4LF3SXpo+vtDkn7LuX3JZx1eUD84\nSNJLkp6WdKt8+79oGVLe2yGUUiq3ZPZQ7tySWR9k9hCZTUNmbe1LZLYWZNbWvkRm15LzQuxWSRdn\n/v6sDlfQ076kr0v6tqT3Zmhfkm5R//Grpn/ekmEZH5D0pKQHlPZR8YFO/f9UfEv5+n+wjG9O/+69\nDmMokVsyu1invLk9aJ/MrofMLtaJzK6LzNqR2TqQWTsyu4acF2L7q5/i4m3qN947Jb1f/UelOe3L\nf90+LemN6j8ifV7SJxLbu0nSlyR9SNLPBv/m1f+bJH1xuoyX5L8OYymRWzIblju3ZDYemQ0js3HI\nrA2ZrQeZtSGza8p5IfYT9YVwB25T/z8I3p6f/vlTSV9W//Gxt0vqv0sqSSckvejc/os6DNNnlLYO\nN6oP7OckPTx9zLv/B8v4/MwyPNdhTCVyS2bn5c4tmU1DZueR2Xhk1obM1oPM2pDZNeW8EPu2pF9W\n/3Hfz0v6HUnnnJfxakmvnf7+Gknv0NGiQS/nJN09/f1uHe4sLydmfn+34tfhmPqPUb8v6VMzj3v2\nf9EyvNZhbLlzS2bn5c4tmU1DZueR2TRk1obM1oPM2pDZyrxT/SwkP5Z0OkP7b1RfaPeE+qknPZbx\nBUnPSfpX9d8H/j31s9h8XT7TfQ7b/31Jf6Z+ytIn1Qcq9vuuO5L+Q/32mJ1207P/oWW8U37rUIOc\nuSWz83LnlsymIbPzyGw6Mru8fTJbHzK7vH0yCwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAKAJ/x8B61paJBrbwgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3cfbe10>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import block_gibbs\n",
      "gibbs = reload(block_gibbs)\n",
      "a = [i for i in gibbs.tree_iter(0, 4)()]\n",
      "print a\n",
      "b = [i for i in gibbs.tree_iter(1, 4)()]\n",
      "print b\n",
      "c = [i for i in gibbs.rev_tree_iter(0, 4)()]\n",
      "print c, c == [x for x in reversed(a)]\n",
      "d = [i for i in gibbs.rev_tree_iter(1, 4)()]\n",
      "print d, d == [x for x in reversed(b)]\n",
      "print gibbs.prec(2, 0, 4)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[(0, 0), (1, 0), (1, 1), (1, 2), (2, 0), (3, 0), (3, 1), (3, 2)]\n",
        "[(3, 3), (2, 3), (2, 2), (2, 1), (1, 3), (0, 3), (0, 2), (0, 1)]\n",
        "[(3, 2), (3, 1), (3, 0), (2, 0), (1, 2), (1, 1), (1, 0), (0, 0)] True\n",
        "[(0, 1), (0, 2), (0, 3), (1, 3), (2, 1), (2, 2), (2, 3), (3, 3)] True\n",
        "[(3, 0)]\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rng.binomial(1, 0.5, 100) * 2 - 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "array([ 1,  1,  1,  1, -1,  1, -1, -1,  1, -1, -1, -1,  1,  1, -1, -1, -1,\n",
        "        1, -1,  1,  1,  1,  1,  1,  1, -1, -1,  1, -1,  1, -1,  1,  1, -1,\n",
        "        1, -1,  1,  1,  1, -1, -1, -1,  1,  1, -1,  1, -1, -1,  1,  1,  1,\n",
        "       -1, -1, -1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1,  1, -1,  1,  1,\n",
        "        1, -1, -1, -1,  1,  1, -1,  1,  1,  1,  1,  1, -1, -1,  1, -1,  1,\n",
        "        1,  1,  1, -1, -1, -1, -1,  1,  1, -1, -1,  1, -1,  1,  1])"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = 3\n",
      "np.exp(0.45 * x) / (np.exp(0.45 * x) + np.exp(-0.45 * x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 37,
       "text": [
        "0.9370266439430035"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}