blackbody.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import bisect\n",
    "import itertools as it\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import numpy.linalg as npla\n",
    "import random\n",
    "import sys\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the transition probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "weights = np.array([0, 0.01, 0.0001])\n",
    "weights[0] = 1 - weights.sum()\n",
    "# weights now has the values [0.9899, 0.01, 0.0001] which are the probabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a class that randomly chooses one of a set of elements with specified weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# adapted from http://stackoverflow.com/a/4322940/56541\n",
    "class random_chooser(object):\n",
    "    '''Implements a weighted random choice, one where the probability\n",
    "    distribution of the choices is not uniform. The probability of choosing\n",
    "    a particular element is proportional to the weight associated with\n",
    "    that element; weights are provided along with the elements, in the\n",
    "    class constructor.'''\n",
    "    def __init__(self, choices):\n",
    "        '''Constructs an object to randomly choose from the given options.\n",
    "\n",
    "        The argument can be a dict in which the keys are the choices and\n",
    "        the values are the associated weights, or it can be a list-like\n",
    "        object of length N, in which case the choices will be 0,1,...,N-1\n",
    "        and the weights will be the elements of the list. In either case\n",
    "        the weights must be numeric.'''\n",
    "        if len(choices) == 1:\n",
    "            if isinstance(choices, dict):\n",
    "                self.choose = lambda: list(choices.values())[0]\n",
    "            else:\n",
    "                self.choose = lambda: 0\n",
    "        if isinstance(choices, dict):\n",
    "            weight_iter = choices.values()\n",
    "            value_iter = choices.keys()\n",
    "        else:\n",
    "            weight_iter = iter(choices)\n",
    "            value_iter = range(len(choices))\n",
    "        self.cumulative_weights = list(it.accumulate(weight_iter))\n",
    "        self.values = list(value_iter)\n",
    "    def choose(self):\n",
    "        '''Returns an element randomly chosen according to the probabilities\n",
    "        provided in the constructor.'''\n",
    "        x = random.random() * self.cumulative_weights[-1]\n",
    "        i = bisect.bisect(self.cumulative_weights, x)\n",
    "        return self.values[i]\n",
    "    def __len__(self):\n",
    "        return len(self.values)\n",
    "\n",
    "def test_random_chooser():\n",
    "    choices = random_chooser({0: 10, 1: 45, 2: 25, 3: 20})\n",
    "    results = [0] * len(choices)\n",
    "    for n in range(100000):\n",
    "        results[choices.choose()] += 1\n",
    "    print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test should return something close to `[10000, 45000, 25000, 20000]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10004, 45247, 24969, 19780]\n"
     ]
    }
   ],
   "source": [
    "test_random_chooser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Monte Carlo simulation starts with a state representing the frequency of a photon. Instead of a vector as in Chad's code, my implementation represents the frequency by an integer. For the three-frequency simple case, 0 is the main frequency, 1 is the medium frequency, and 2 is the low frequency.\n",
    "\n",
    "At each step, the code makes the following choices:\n",
    "\n",
    "1. A photon in state $\\lvert i\\rangle$ can be absorbed with probability $P_i$, or not absorbed with probability $1 - P_i$.\n",
    "2. If it was absorbed, it will be reemitted into state $\\lvert j\\rangle$ with probability $P_{ji}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def monte_carlo(state, absorb_weights, emit_weights):\n",
    "    '''Create a generator to run the Monte Carlo simulation.\n",
    "    \n",
    "    absorb_weights is an iterable of the probabilities P_i\n",
    "    emit_weights is a matrix (an iterable of lists) of the\n",
    "    probabilities P_ji, such that emit_weights[i][j] = P_ji'''\n",
    "    absorb_chooser = [random_chooser({True: w, False: 1-w}) for w in absorb_weights]\n",
    "    emit_chooser = [random_chooser(w) for w in emit_weights]\n",
    "    while True:\n",
    "        if absorb_chooser[state].choose():\n",
    "            state = emit_chooser[state].choose()\n",
    "        yield state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_with_monte_carlo(initial_freq, steps, absorb_weights, emit_weights, iterations=10000):\n",
    "    total = np.zeros(3)\n",
    "    for n in range(iterations):\n",
    "        for s in it.islice(monte_carlo(initial_freq, absorb_weights, emit_weights), steps):\n",
    "            pass\n",
    "        total[s] += 1\n",
    "    return total/total.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transition matrix "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transition matrix method calculates the transition probabilities instead of simulating individual absorptions and emissions. In this case, the state is a NumPy array in which element `i` is the probability of a photon having frequency `f_i`. For the three-frequency simple case, a photon at the main frequency is represented by the state vector `[1,0,0]`, a photon at the medium frequency is represented by `[0,1,0]`, and one at the low frequency by `[0,0,1]`.\n",
    "\n",
    "Each absorption/emission step is represented by the application of a transition matrix, which is computed as\n",
    "\\begin{equation}\n",
    "T = \\underbrace{\\begin{pmatrix}\n",
    "P_{00} & P_{01} & \\cdots \\\\\n",
    "P_{10} & P_{11} & \\cdots \\\\\n",
    "\\vdots & \\vdots & \\ddots\n",
    "\\end{pmatrix}}_\\text{emit}\n",
    "\\underbrace{\\begin{pmatrix}\n",
    "P_0 & 0 & \\cdots \\\\\n",
    "0 & P_1 & \\cdots \\\\\n",
    "\\vdots & \\vdots & \\ddots\n",
    "\\end{pmatrix}}_\\text{absorb}\n",
    "+\n",
    "\\underbrace{\\left[I - \\begin{pmatrix}\n",
    "P_0 & 0 & \\cdots \\\\\n",
    "0 & P_1 & \\cdots \\\\\n",
    "\\vdots & \\vdots & \\ddots\n",
    "\\end{pmatrix}\\right]}_\\text{not absorbed}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def transition_matrix(absorb_weights, emit_weights):\n",
    "    '''Create the transition matrix that represents one step.\n",
    "    \n",
    "    absorb_weights is an iterable of the probabilities P_i\n",
    "    emit_weights is a matrix of the probabilities P_ji,\n",
    "    such that emit_weights[i][j] = P_ji'''\n",
    "    emit_matrix = emit_weights.T\n",
    "    absorb_matrix = np.diag(absorb_weights)\n",
    "    identity = np.identity(len(absorb_weights))\n",
    "    return emit_matrix.dot(absorb_matrix) + (identity - absorb_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_with_transition_matrix(initial_freq, steps, absorb_weights, emit_weights):\n",
    "    initial_state = np.zeros(len(absorb_weights))\n",
    "    initial_state[initial_freq] = 1\n",
    "    full_transition_matrix = npla.matrix_power(transition_matrix(absorb_weights, emit_weights), steps)\n",
    "    return full_transition_matrix.dot(np.array([1,0,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Test data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reproduces the first plot in the blog post"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "weights = np.array([0, 0.01, 0.0001])\n",
    "weights[0] = 1 - weights.sum()\n",
    "# weights now has the values [0.9899, 0.01, 0.0001] which are the probabilities\n",
    "\n",
    "absorb_weights = np.array([\n",
    " # probability of absorbing a photon in state 0\n",
    " 1,\n",
    " # probability of absorbing a photon in state 1\n",
    " weights[1],\n",
    " # probability of absorbing a photon in state 2\n",
    " weights[2]\n",
    "])\n",
    "emit_weights = np.array([\n",
    " # probabilities of emitting a photon from state 0 back to...\n",
    " weights[[\n",
    "  0, # state 0 (P_00)\n",
    "  1, # state 1 (P_10)\n",
    "  2  # state 2 (P_20)\n",
    " ]],\n",
    " # probabilities of emitting a photon from state 1 back to...\n",
    " weights[[\n",
    "  1, # state 0 (P_01)\n",
    "  0, # state 1 (P_11)\n",
    "  2  # state 2 (P_21)\n",
    " ]],\n",
    " # probabilities of emitting a photon from state 2 back to...\n",
    " weights[[\n",
    "  2, # state 0 (P_02)\n",
    "  1, # state 1 (P_12)\n",
    "  0  # state 2 (P_22)\n",
    " ]]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.3658  0.6271  0.0071]\n",
      "[ 0.36497205  0.62866848  0.00635947]\n"
     ]
    }
   ],
   "source": [
    "print(compute_with_monte_carlo(0, 100, absorb_weights, emit_weights))\n",
    "print(compute_with_transition_matrix(0, 100, absorb_weights, emit_weights))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "counts = [int(pow(10,n/3.)) for n in range(30)]\n",
    "results = np.array([compute_with_transition_matrix(0, n, absorb_weights, emit_weights) for n in counts])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x105f6ca90>,\n",
       " <matplotlib.lines.Line2D at 0x105ff7b00>,\n",
       " <matplotlib.lines.Line2D at 0x105ff7c88>]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOX9/vH3LYIgiztaKYgKtkDcrUrrktaNxLXaal1r\nFay1gkWrqPVnMm2tu1WgiOJSlbpgq3WDWGsNUpeKdhFcgW8VUIviUqGKIHx+f5yJTJIzyUxyZs6Z\nmc/ruuaKeeaZMzfJeJ6c82wyM5xzzlWudeIO4JxzLl7eEDjnXIXzhsA55yqcNwTOOVfhvCFwzrkK\n5w2Bc85VOG8InHOuwnlD4JxzFW7dQh5cUk9gEvAZ0Ghmdxby/ZxzzuWv0FcERwLTzOw04LACv5dz\nzrkOyLshkHSLpCWS5rQoHyHpVUnzJI1LF/cDFqX/e3UnszrnnCuAjlwR3AqMyCyQ1AWYmC4fChwr\naQiwGOjfifdyzjlXYHmfnM1sFvBhi+Ldgflm9oaZrQLuBg4H7gOOkjQJeLCzYZ1zzkUvqs7izFtA\nEFwJ7GFmnwCntPVCSb78qXPOdYCZKYrjRHW7ppMnc0s/ahrMTJkPGPHo2ufrmtUFWwdsfbC+cMLf\nwuud8x7Ys2CL4f+tAfsP2Gyw+2Hn/2Y7dssca/OQyvZcLvValmd+n8ux6cvB7MJ89gXq04/deJOD\nuIZ6plDPK9TzMfX8qVW9pscGzG913A1Y8MXzma8ZROvfySAebVWvrbphGQbzKOvyK+rZkHr6Us+X\n2ZS3QjPsyivUcxr1nE09ddRzJfVcz3CWhB57b1ZQz1+pZxr1XEs957EjL+b8s8jzd5KEz0U+9bLV\nba+sHH8WuWZISo4WGSIT1RXBW6ztCyD934tzf3k98OzbsHhC6+cWjYdR28KUQVCdLhu5ABZOMMOA\nT4FPpaUfrX1NdcbrX37BjBoA6Zr94OevAV8OMvbbA9ig9XsO3FZiGPBy+j0IXl9VCzvXSkdXw7IV\nsGi82dzpWf5RjTmWZ6sXrg9jOIxB/Duj7BC2YgbHAZcA1wNzrM4+1yTV0psbCa7YAg+ygFXc1Oq4\nq5jCg4zkMAY1q/sxrX8nHzOeB9mW7XOsa1S1yvBfxvM5n1id/bepWJP1drN6Tf7Lm1ZnU1oWa6q2\n4Rsc1Kr+Up4D/h+wZfp4W7I+WzCwVU3YjHWV0kDgTav7YnOORm2uWvowho35qgZrTz5mvC2xbL9r\nyP33mK1ey/Jcj9fR14TVzbUsqhzZ6rUsT0KGpORolFRN85Ncp8k6sDGNpIHAQ2a2ffr7dYHXgP2A\nt4HngGPN7JUcjmXBX/cLJ2Q7qQYn4AGjoVcPWP5pWN2gzvDrggajycgF8OyY7MetaYAZ6RNJffoB\nMHoxTFgJ9AZmAo1w+mpYfU7z44+aD8+c1UZjkBdJ9WZW32adHfUsR7JHqyemMdNesupW9YMT2mi6\n0YOVfMrHTMh2Qvui7gcMYWNeyaluPsdtp642Vy1f5joOYxBPAN8kaDQWM6bd+k2y1NdgNXBCSKPR\nwLuMYDXBH0XPAbOZifExJ3MoW2fkmM9izmqnMSiIXD4XlZIjCRmSkkOSRXVlkHdDIOkuggv3TYB3\ngYvN7FZJNcC1QBfgZjO7NMfjRfePyaHBaF2/qfFoJGhk1zYeEgMI/q3VMO57cPn6rY9S22A2vSaa\n/Ko2s8bQ51IS8H3+zA3sT7dWFabSYPOs4DkK7YtG4zO+xHq801YD06x+Po1Mk4xGQyn1Ixj0sDsN\nnMYINgbg38DW6foR/ozzEefvI2k5kpAhKTlibQiilu4sThHMPG4s/vvn1nhIRzfCtH1bH+HomWbT\nqguaMaWNgMnAMB5hEqsZm8tfwa65nBuNYWrkaFr/rh/iZQ5lZ6uzlUUJ7FyIjFtDdVE1BAVdYiJX\ncV5ipU/6OZxAl60IL99isMRAM96IMlcTpbQvcDvwAHCyPWefanO9wdTcbsu4tdI/o/Z/TisJ/113\nYwvgTaV0CzDF6uwN+KKBGUM3urOSFTn0JzjXYek/mBsl1UV1zEQ0BKUhs9O6yaj/g288C7wgcTtw\niRlLo3g3pdSV4ErpZGCk1a09seR8QnMd09QZHnLVRXCz6IfAC0rpGWYymy9zQou622pzUezGoEqq\n7Q9jekP3ZbBiEYyfa+EZvG6ycuSbOXJmFuuDYLxmPVAdd5b2sw6rhZoZ8N3G4Ouw2qDcNgebCLYU\n7CKwnkHdEQ1B3RENTXVDj9uXWgbRwFAaGUQDOzOSemZTzyPU0zfuf3clPtK/kxnp38kM+tLs90c9\n61PPyQznI+qxVo9BzOhshmFQOwIavguNI6BhGGT9DA2D2pEwz8CaHiNhXthrvG6ycuSbmeC2UH1w\n+o7o8x77/3AR/mPifoANArsL/vw+jHk34/dqMHJeWGNAX2rZhXnNTiJ78TkH8xvqgz4cfyT3wVAa\nQxuCoTR25rjtnhygp8F2Bt8yOHEkvG7NP3BmYD+Bdw0eNnjEYLrBjJ/Ae2F1xwblj2Y+fgJLy7Vu\nUnJkq1tD239MRHnu9FtDETJjPnCsdMPTMG1482enDILa0bS8pdM0LyDT/nRhKtvYw8Fv2yVYtv6E\n1azpzGH7w5gpNP9cTIFB58HdSKuB7gTzdxYDi9cPvm9lOSwhGGgA6RmTy+AyYNOWdT8Ohn5fk1m2\nDK4gGCFYdnWTkiNb3V7Qo2VZoXhDUBhZRpX0av2L7Rb+PzDdivchcJ0Q1p/QwIfswc5K6XjgzoyJ\najnbBDYMK18SzNcZAXzQ9GchwOtSA80ndQLwFizG7OFmZdIYYPuWdd+GtzF7tEXdscAO5Vg3KTmy\n1V0eTJYtikSsCCqpPj0kqkxkG2G0vPUvNttflCuL9yFwHWdLbDqLOYupNDCNmUylgf/jBLZhf+B8\n4AGltGXOB5S2QrqxP+wS9vR7sBSz9zMbAYBFMH4UzM8sGwkLFtJ6trfXTVaOfDNLqpZUH/ZcRyVi\nHoFFvG5G3MJnOY9bBXNONJt+T7O6h+h6PuIUDsiYJObzAsqCUuoGXAScDvwUuCPr1YE0ALgQ+C4w\n+WB4cUv4ZebtoZGw4FkY09bIkwEwuhf0WA6fLoQJXrftuknJkW9mKMMJZeXWEEDYRLXjP4XjNwYO\nNAtuHSmlw4HruY+f8QlH+7yA8qSUdgZ+CyzqO4FXt1nO93uKrv8zVq1j3PbUSnoA3wOmAFdhthQ6\ndnJwlcMbghIk0QX4PfAxcDL1+hrwCFBjdfZ8rOFcwSmlbhv+kUcPfJXqezJuBp4PbL4OT41dw7cx\ney+2gK7keENQoiR6Ao0MmDWLU/Y5FjjN6uyhuHO54hjeXe8981nr0TrD12PpMytsszgyudIV5bkz\nEaOG0h0fjZaAxaQKyYz/aa8rTmLXG1/kpe9MtWn3eiNQQXqKrvmUOxcmMctQRxqgkq4IUuoOPMaH\nA+dz3b8PBo4yY1bcuVxxfK+rlt/9OT1blu/eh1Wzz2aY1dm8OHK50hTluTMRw0crgVJah6DD8B02\neuNU4ATgXonBsQZzxSFdcMLnfHZiVz7PLP5eN1bNH8DjwDNK6adKKRFX6a6yeENQPJcS7Ix2ktXZ\nGjP+RLCL1iNS61mFrkxIIrj1edIhsP3jXfjl8PVYun93/jt8PZY+uQ6XfDDHaoA9gFrgaaVUFWtm\nV3H81lCBNFuauBubsT192JadrM7eb1ZPXA6PHAyT3oae3XLYAtOVCknAr4BDgP0xW9Jm9WDzoZHp\n10xkCv+gB2f48tYuTNmNGiLGjWkKIXQ3rId5k4Wc0WoLRe1wMBx4N1zVa21ptFtguhgEjcA1BJ16\nBzTNDcjppSl9mVe5j0Xs2GKiYWzbZbrkKMTGNIloCMrtiiDr/rghWx023zc5U3RbYLoik9YBJgK7\nAQdh9mHeh8jjM/TFa/LYIMfrJitHRzY3Krvho2Unr4XkeofXDVugziVS5qYiy+Gzk8GOhl4Et4M+\n7tBBs32GtmR/pfQ48PeMxzwmMSJkT+bQDXKy7N9ccXWTkiPfzIXgDUEh5LWQXB4L1LnEqZJqh8N1\nmWsCnQefToYT/9LRRgCyf4aW8CTBssW7AEcCvwQ2Y1s+5yA2alb3MAZxP5copfWBNQTLUK9hC+pa\nLX1+GIP4Az9XSs1vEWzOz8u2blJyZKs7ldbL1heINwSFMIC/8Wf2Y/+Mn++DLODjsNUEQ7fAXAAL\nQ1cedMkStm/AFdCjNuj0/UOHD5xtu8wPudrq7FHgi6WMldJGrOAJaNEQAPSkH8E6RuukH6IXA0Lf\nsw9bA2NCysqzblJyZKtbxKXovSGImFLagn34IY9xAVPZr72F5MzmTpeqCDat6dUD+m8Pg6eZTfEO\nwRLQO8uGMJ3dVMSW2HRtLpjK6HY/Q3X2oabqP8COrQ70Li9YnX0ns0hT1QAh/Q9LeM7qWvRhlXHd\npOTIWreIS9H7PIIIZUwau9H+alfZPKuxl6za5llNW/f6zOZON5teYzatGq7+MZy3Z7Eyu85ZRvgt\nnCg2FbElNj3Xz1D6CqLZmvZZr0K9brJy5Ju5ABIxaogyGT6qlMYAxwN7WZ2t6tAxRDdgEbCvGa9G\nmc9Fb5x0zjpw+aXQpamsvX0DCiU98qTdKwivm7wcedX14aPJlZ4N+gQw3Opsfnv12zyWuBTobsbY\nSMK5wpC2AJ7/JUx5Gvb0fQNcMZXdhLJSbwjSi8k9B1xrdXZLp48nBgLPAwPM+KSzx3MFIHUF/gL8\nGbNU3HFc5fFF55LnUuB14NYoDmbGG8DfgGOiOJ4riKuAj4BfxB3Euc7yhqCTlNKBwHcINpmJ8vLq\neuBHER7PRUU6gWCBuBMxWxN3HOc6yxuCTlBKmxJcBZxsdfZBxIefAfSV2DXi47rOkHYEfg0cidlH\nccdxLgreEHRQeqXIm4A7rc4ej/r4ZqwGbsSvCpJD2hi4DxiN2Zy44zgXFW8IOm4ksBVwUQHf42bg\nKIkNC/geLhdSF+B3wAOY3R13HOei5DOL8/DFCoHrsxF92IkVnGWz7bNCvZ8ZSyQeBU6E4k0ucaHq\nCGYLj4s7iHNR8+GjOcqyQmDB14eXqAYmAcPMiPeXVamkw4DfALu1t7mMc8VSdstQK9jKL9kzi/sw\nJqYVAmemv+6T8d+uwJqWlt4MNuoHO62C86/yRsAlQMbM4sgkoiEws/q4M7Qrrz0GomOGSUwm6DT2\nhqAIwpaWHgVnVEmv+YxhF7f0H8yNkuqiOqZ3FudqFeF9AcVZIfB24CCJzYvwXhUvbGnpKTBoAIyO\nK5NzheQNQa6G8B/+1GK5hyKtEGjGRwRr259a6PdyhVta2rmkSsStoaRTSl9lFw7mD/yYqRyT66qG\nEbseuE/i8vQcA1cgq7P8gRTF0tLOJZE3BO1I7zFwM1BnL9pvCfYbKDozXpBYAtQAD8eRoSJI3UbC\nFmfBe9fBZk3FI2HBQh/C68qUNwTtO4Ngr9fr4w7C2vWHvCEonF/WwCvnwg21MNqXlnaVwOcRtEEp\nbQW8AHzD6uy12POI9Qk2rdk1vUKpi5J0AMHaUTthtjTuOM61xZehLoL0WkI3AFcnoREASO9NcAdw\nWtxZyo60GcFtv5O8EXCVxm8NZXcisDnBuvNJMhkef1q6Zjfo2Q2WrYBF483m+m2LjpJEcCVwO2Z/\niTuOc8XmDUEIpbQ5cCVQ09G9hwunahuo7Q6PHLC2bNS2UhXeGHTYmUBf4OK4gzgXB781FG4CcKvV\n2d/jDtJa/zFwRYvx7FMGwQCf7NQR0g4EDcCxWNIafeeKo6ANgaStJd0k6d5Cvk+UlNIRwE5AQveh\n7R2+1AW9fLJTvqT1gbuBszFbEHcc5+JS0IbAzP5tZiML+R5RUkobAhOBkVZnCZ08tGxFePnyhOZN\ntGuAv2N2R9xBnItTTg2BpFskLZE0p0X5CEmvSponqRzWab8SeNDq7Mm4g2S3aDyMmt+8bOQCWOiT\nnfIhHQkcQDBPxLmKlmtn8a0E981vbypQsGPTRGB/4C1gtqQHgd2AXYArzeztaONG74vNZnqzBevz\nFd7hRCJb0y96ZnOnS1XA0RfC4F3gHzNh4QTvKG5f09LSG8MG/WGX1XDBlWYfx53LubjlPKFM0kDg\nITPbPv39cKDOzEakvz8fwMwuy3jNxsCvgP2Am8zs8pDjxjahLK7NZqIg0QVYCgw145248yRdlqWl\n5z8DZ/mMYVeKkrIxTT+CWa5NFgN7ZFYwsw+A09s7UHpjmibF26Amvs1mOs2M1RKzCDasuSfuPEmX\nbWnpWpL/u3YOCrMhTZPONASRrU0R28Y0MW02E6GZwL54Q9AuX1ralbqmDWmavk/KxjRvAf0zvu9P\ncFWQN0n16dauuFbGutlMFJoaAteOFYQv3e1LS7tSI6m6xV2UTutMQ/A8MFjSQEndgGOABztyIDOr\nj2W/4u2Yz2M0H45ZpM1mIvJP4MvS2uWSXQhJp0D3s+GjzGJfWtqVIjNrjPouSk6dxZLuIvjLcxPg\nXeBiM7tVUg1wLdAFuNnMLs07QEydxUppC2AOD3MxH3FYTJvNdJrEdOBmM/4Qd5bEks4ATt0Tfr4x\nnO5LS7tyEOW5MxHLUBPM4i1eJzGglO4C3rA6u6BY71kIEuOAfmaMiTtLIkk7AY8BX8dsXtxxnOus\njE7jurJqCIp9RaCUDiRYYnqY1dkn7dVPMok9gRvM2DHuLIkj9SK4hfkLzH4XdxznolR2VwTFbAiU\nUg9gDjDG6kr/toBEV+B9YKAZH8SdJ1GkWwEw+0HMSZyLXNltTFPkUUM/A/5RDo0AgBmrgGeBvePO\nkijSCcBwgiWmnSsbhRg1VFFXBEppCDAL2MHqkr/8Ra4kfgZsYsbZcWdJBGkw8DSwP2b/ijuOc4VQ\ndlcExZDeevJ6IFVOjUCazydoIq1HsLR0nTcCzuWmYhoC4PtAT2BS3EEKYDbwFYkN4g6SAJcDbxI0\n+s65HCRiq8r0/a6CDR9VSpsSnCBqrc5CZ5iWMjM+k3gO2At4JO48sZEOA44Adibue57OFUgh1hyq\niD4CpXQLsMzq7KxCvk+cJOqAnmacF3eWYmq5tLTBhZebXR13LucKLSmrj5YEpbQPwQYkw+LOUmAz\nCa56KkaWpaVPr5Je8RnDzuWurPsIlFI3YDLwE6sr+w1I/gYMk+gdd5Biyba09IBgaWnnXI4S0RAU\ncB7BucD/AfcV4NiJYsanwN+Br8edpVh8aWlXiZK2+mhkCrH6qFLaFhgLnGl1FdNxWFHDSC3LrU1f\nWtqVs0KsPpqIhiBq6TkDk4ArrM7eiDlOMT1JsGNZ+ZM2PQ22GQNLMot9aWnn8ldWo4a+2Ih+A/qz\nHgN4g+PsLXsoimOXAomeBCfGvmaU9GJ6bQr2v/gT8EwVzBoAo31paVdpfNG5sOOU8Eb0UZJ4BviZ\nGX+JO0tBSCJYOXYL4AjM1sScyLlY+BITYbJtRN+n4kaQlHs/wZkEi8kd742Ac9FIREMQyaih0t+I\nPirl2xBIBwIXAodhtizuOM7FwUcNtaX0N6KPylPAblKWhrFUSV8B7gCOwezfccdxLi4+aqgtwUb0\nzU/6pbURfSTM+Bh4Bdg97iyRkTYCHgQuxOzJuOM4V27KorNYKW0OzOURLubD0t2IPioSVwH/NeMX\ncWfpNGldYDrwEmZj447jXFL4qKGWx0jpNuBdq7NzI4pV0iQOBc4yY/+4s3SadB3wFeAQzD6PO45z\nSeGLzmVQSvsC3wKGxJ0lQWYBd0p0M2Nl3GHy0bSaaG/o3hM2PRY2OBC290bAucIp6T4CpdSVYAbx\nWKuz5XHnSQozPgLmA7vFnSUfTauJzoCDpsG+t8Kw+2F1VQWtn+RcHBLREHRi+OhPgEXAH6JNVBZK\nbhhp2Gqi18NWvpqoc2v58NEMSqk/MI7KWlQuHyXXEPhqos61z4ePNvdrYKLV2fy4gyTULODrUun0\nA3XJcsL31USdK6ySbAiUUg2wExW2I1c+zFhKsIn7LnFnyYl04Cj4ypnwdmaxrybqXOGV3PBRpdQD\nmEtwS2hG4ZKVPomJwJtmXBl3ljZJRwA3At+ugg18NVHn2lfpw0fHAf/wRiAnM4HvQ4IbAul44Gqg\nBrMX5galfuJ3rohKqiFQSoMJVp/cOe4sJeJJYIpEFzNWxx2mFek04GJgP8xeijuOc5WqJPoItLlq\nNVgNPMbfeJD3mcT2cWcqBWYsAd4Bdow7SyvSOcAFQLU3As7FK/FXBCEbzmzEg1ynzUUlriPUATMJ\ntq/8e1wBMmcLL4MVx8D7JweT3fbBbFFcuZxzgcR3FmuwGjiBg1o9MZUGm2c1BQ1XBqSrL4eFJ8M7\nr8CyFbBovNncojWgTbOFMyeKjYOVr8MP7je7s1g5nCs3ZddZnJ4l1xg6qcw3nOkwqaoW9joaJvcF\n+galo7aVqihWYxA2W/hy6FYLJwLeEDiXp/QqDNVRHjMRfQRtziy2LI1V5W040wH9x8Dkgc3LpgyC\nAUVbssFnCzsXrYqbWayUxM70YQbvNXuiAjec6ZjeWXYp61Wck7DUdUPoF/aUzxZ2LjkScWuoDcex\nHat5glOZyhmVvuFM/patCC9fXviTsLQbcPMh8OHp8MZkGNj0lM8Wdi5ZEttZrJQ2INhy8Sirs2eK\nn6z0BX0Ew68Lbgc1OfM/0HhqwfoIpPWBFHAScA7wuyqo8dnCzkWrInYoU0rXAr2szkbGEKtsBI3B\ngNHB7SB1gx8OhW/tbEb0G8BL3wSmALOBszB7N/L3cM4BFdAQKKUdgceAoVZnS+NJVp4kxgLHAHub\nsaqjx8mcG7ACVh8PK4+BKuAMzB6KLLBzLlRZNwRKaR2CpRHusDq7Ib5k5UlCwEPAXDPO78gxwuYG\nnAMfz4WRj5rdG1VW51x2UTYESRw1dBLQDbgp7iDlyAwDfgAcL3FgR44xEMa2nBtwNfTpAqdEENE5\nV2SJagiU0kbAZcAZVmfJWyStTJjxHsGErt9KbJHTiyQh7YF0wy5ZJrP43ADnSlPSho9eAtxndfZ8\n3EHKnRmNElOAOyQOMmNNyzWBFsH4ufA8QaNxCsGV2i0vwlOEbIPpcwOcK02JaQiU0q7AkcCQuLNU\nkF8AfwHOq5JebHnf/1zY6wmwb8LvgdOBv2Jm86V/jYJ+mXV9boBzpavgncWSDgcOBvoAN5vZYy2e\nN+rpAjwDTLY6u7WggSpM6F/5GWP4JfqLNc/XssnbD/PRTi1ffyj8+SGzA8KO63MDnItPSS06Z2YP\nAA9I2hC4imBYaHN/5EWGIrbjtkLnqSRho3tGwbbVUvdGeB8YbrDnCtZb7+esbtUIAPSArmHl6ZO+\nn/idKwM5dxZLukXSEklzWpSPkPSqpHmSxrVxiIuAiaHPHMEwXqc3kxiRax7XvrCVP6fAoG/ANOBy\nYAvgru58tsOTDFkYdoxs9/2lqlqppkE6ujH4WlWbLUc+dZ1zxZfPFcGtBPeAb28qkNSF4OS+P/AW\nMFvSgwSbjuxCsFfuOwQjgWaY2T+zHv0Q+jOV0fhfmdGQttsWtgt7aj48jdk+mWULtMVtx7DBRffw\n3y8uNb9Hr1X/pOuzrQ8dtnRF+PLW+dRdW7//mGDBvPb3T8i3vnOutZwbAjObJWlgi+Ldgflm9gaA\npLuBw83sMuCOdNkYYD+gj6RBZiGTxJ5If/2AIZKqsy5J7YA27vtLQ4HvAEcBfQ1Whr1+GfyvZdl/\n2Hn36ZyprzGBnqzgf3TnVUZ3Xc79x0h8AKxOPz6H3X7S/MQOwfcn/kJinYy6q2H7+vC6R42TeJkg\n46rg617fhK9fBTduu7Zue41G4RoZ55KkEPsQNOlsH0E/IHOrwcXAHpkVzGw8ML7No3wz/fUtXrH3\nvRFoS9h9/7Gw+6PSsoOCW31/AH4MPDMTDhrVom720T29uy/nYJ7n4Bblj/UBvgJ0WfvYZLPwdFts\nBfyIIEe67lahVyXw1a8BjQR9EN2Cx7694JIWtyunDILz7paYCywDlq/9ulctTN66df1DzqLFlWW+\njYZzSZP+A7mx6XtJdVEdu7MNQXRDjnyPgZyE3ff/NWz0PXj1INgLszVN5XNhepVEbU6je7ItWf36\nHDPOzCyRXm6AkO1DeWm2WfOWRPpXlrr/mGlGTfO68xoJmZ8AS14DzgV6Ab3TX3vBuln6GnY9QOJ9\n4O21j732gcnbNK83ZRDU+u1IV/E62xC8BfTP+L4/wVVBfn7NfFZxk/3Phx+2ZzPYOKx8DazMbASa\n5D66Z9F4GLVt87+YRy6AhSGNc6HqZmuM3ltqxlMtS6UFhwBbt64/+0/ACcCWBFetW0LX6vBjDxoS\ntJXMTs+4Th/bbyO5ZCrELaLONgTPA4PTfQdvE6xqeWy+B7GPbHAnc1QG6fh+sGPYU52d1Ws2d7pU\nRfAXcq8eweY1CyeEnfwKVTe/RqPN+uPNWAosBV4EkOZ/l4zNcdZa+TkwFvhaui/kObhxGexbA7/J\n2F3NbyO5ZGi6RRTlrSHMLKcHcBfByf4zgn6BH6TLa4DXgPnABbkeL+O4lu9rKu4B6xvcZPDqGXDG\nSJhnYE2PU2H+MKiNPWcEDxhWCzUz4LuNwddhbf67cq0f1Bs5L+PHZnDq/Kb6YOuAfQXsRPjRm83r\nNT1qZsT98/GHP5oeUZ47E7EMNcGOVo3mo4Vak74K3Evwl+3pmC3zWb0d03yTnuxXJtLRjTAtpK9i\nzJswfi+zDtz+dC4iGbeG6qxc9yNwGaQTgWuAC4CbifuXVSGkmgaYEdLBfeZCmNgbmAlMBh4zY433\nJ7g4lNQSE659LecF/A9ueBIOAfYC9sPsxbgzVpZsfQ/PjgFmEfSDXQZMkm6aCd+ohhsyOq29P8GV\nlkRcEVDBt4bC5gWcDyu/Bk8dBYdjtizOfJWqvdtI6Z3edocz/wAT+7U+Qm2D2fSa1uXOdY7fGipD\nNVLDjJBx9rXQMN3MTyQJl70/4eiZZtOqi53HVY5y36qyovSG7mHlvttXqcg290HdipvDuY7zhiBm\nyyD0ROK7fZWKReNh1PzmZaOXwGlDJC6RWC+eXM7lLhENgaT69H2virMNPHhBsEDbF3y3r9IR9Bs8\ncxbUNsBbbEKIAAAO8ElEQVTRM4OvT5wC+w0BhgIvSHwt7pyufEiqllQf6TG9jyBG0pbAM1fD3Y/D\nDj4voLykO5SPAa4DbgFSZuFXgM7lK8pzpzcEcZF6Ak8C92P2y7jjuMKR2ByYBAyBn06Blw7yOQeu\ns3weQakLNvS5i2C28CUxp3EFZsYSie/Alb+CrlfCjC5rn/U5By5+3kcQj2uAnsAPfbZwZTDD4C87\nw6Vdmj8zZVAwX8G53BSijyARVwRmVh93hqIJdmw7APg6ZqE7iLly1Tt0qDD06VXcHK6UWQFWH03E\nFUHFkA4FzgdqMfso7jiu2LLNORi4q8QBxc3i3FreEBSLtCvByJEjSO/x7CpN2JyDkQugx6+AmyRu\nkdgolmiuovmooQLJXEhuFejHMHT/oE/gvrizufhkW8NIojdwKfBt4Mdm/DHmqC7hfPhowoUtJPcT\nWPpn+L7PD3BtkdgHuAn4BxzxAHx2kg81dWHKbvhouge8bFYfDdtg/lrYtBZ8o3TXJjOelNgR7rgd\nhtzefJSRDzV1ydyzOBLlNmrIF5JznWHGp9KdvZvPN4BgqOkRPyXLHxP5bJDjdZOVI5+6hRg1VFb7\nbiblMQIaQja8tRrwPW/9kdMj2IPZQj5GF60CewfsEbBfgH0bbKvwPZlHzgvbw9nrJitHvpnXvg6L\n6vOWiCuCcmNwywWw36UZV1y+kJzLT7ahpi/8GTgN2CX9OAWYCIduBpd2bV53yiA4/TqJ3YA1gAVf\nh5/cfPe1prqjrpYY0rx8z1HlWzcpObLVrS3arWRvCKImqQGOugum10I3X0jOdUy27TIXTjBjEbAI\neKDpGWnh08Dw1sfp1pXg//N10g/B+lluUfbsBWzZvKxXz/Ktm5Qc2er2KtqtZG8IoncCUHUs7Has\nme8p4DokGFJaRfBXYfh2mc198HF4+fxXzLg4s0R6fWegf+u6r88145zmdV8bBny5HOsmJUf2usuL\nd/6I/15odPe5Yn/A1gbvGewUexZ/VNQj/D7zqfNzvyddeXWTkiPfzGtfh0X1+UnEPALKYfN6aV2g\nEfgjZlfFnMZVoGyT1bxuW1dSyciRX13fvD65pIuAbwIHYLYm7jjOufJWdhPKSp60O8FksV29EXDO\nlRpfdK6zpF7AVOBMzBbHHcc55/Llt4Y6S7oR6IbZyXFHcc5VDr81lBTSEcD+wE5xR3HOuY7yhqCj\npC8Bk4GjMMsyhts555LP+wg6QloHuBW4EbOn4o7jnHOd4VcEeWjabGYb2LYHbDYTfjM77lDOOddJ\niWgISmE/grDNZkbBNVXSal9DyDlXLIXYj8BHDeWoRmqYAQe1LK+FhulmNXFkcs5VrijPnd5HkCPf\nbMY5V668IcjRaugSVr4cfIVR51xJ84YgF9I6p8BGZ8F7mcW+2Yxzrhx4H0EupDOA43eGS78EP/bN\nZpxzcYvy3OkNQXukAcDfgb0xeyXuOM45B95ZXDySCGYPX+uNgHOuXHlD0LbjCLaQuyLuIM45Vyh+\naygbaTNgDnAoZj6B2DmXKN5HUAzSncBbmJ0bdxTnnGvJl6EuNOkQYHdgh7ijOOdcoRW0IZD0VeAs\nYFPgcTObXMj3i4TUB5gEfB+zT+KO45xzhVaUW0MKlm2+zcxODHkuWbeGpOuBLpidFncU55zLpujD\nRyXdImmJpDktykdIelXSPEnjsrz2UOBhIPkTr6R9gEOB8+KO4pxzxZLTFYGkvYHlwO1mtn26rAvw\nGsFWjW8Bs4Fjgd2AXYArzeztjGM8bGaHhBw71iuCpj0G+sD6W8GuG8LVF5pdHFce55zLRdE7i81s\nlqSBLYp3B+ab2RvpUHcDh5vZZcAd6bJ9gSOB9YBHoggcpSx7DBxbJT3rS0c45ypFZzqL+wGLMr5f\nDOyRWcHMZgIz2ztQemOaJkXboKY/jMlsBACmwKBaGE0p3MpyzlWMQmxI06QzDUFkvcxmVh/VsfLh\neww450pF+g/kxqbvJdVFdezOLDHxFtA/4/v+BFcFJWMZrAgr9z0GnHOVpDMNwfPAYEkDJXUDjgEe\n7MiBJNWnL3uKqg/87kJYk1nmeww455JMUnWL2+mdP2aOo4buAvYFNgHeBS42s1sl1QDXEuzedbOZ\nXZp3gDhHDUm/vw1W3gMb+R4DzrlS4msNRfPGBwPXAdtj5reCnHMlpezWGkpf5hRttBDS+sBE4DRv\nBJxzpaQQo4cq84pAugwYgNlxRX1f55yLSNldERSVVAWcCmwfdxTnnEuCRDQERbs1FCx+dwNwMWb/\nKeh7OedcAfitoc6/2UhgJPB1zNa0V90555LKRw117I36AnOBAzD7V8HfzznnCsgbgo690e3Au5j9\ntODv5ZxzBVZ2ncUF7yOQvkVwT21oQY7vnHNF4n0EHXuD9YAXgfMwe6Bg7+Occ0VU9B3KStw44BVv\nBJxzLlx5XxFIg4FngF0wW1iQ93DOuRh4H0GOBwWuBy71RsA5Vy68jyC/Ax8PnAvshtnnkR/fOedi\n5MNH2z/oRsDLwBGY/S3SYzvnXAKU3a2hqFRJtf1hzBDY4TNYNRM2mRt3KOecS7iyaQiqpNrhcF3m\nZvSj4LoqCd9oxjnnskvE8NEotqrsD2MyGwGAKTBoAIzuVDjnnEuQQmxVmYgrAjOr7+wxekP3sPJe\n0KOzx3bOuaRIj65slFQX1TETcUUQhZVZypeD70DmnHNtKJuG4CTQWPgws2wkLFgIE+LK5JxzpaA8\nho9KBwI37ANn94LTekGP5fDpQpjgHcXOuXLk8wiaH2A9YA5wNmYPRxbMOecSzBeda+5c4GVvBJxz\nrmMSMWqow2sNSdsAY4Fdok/lnHPJ42sNtXgh8DDwV8wujTyYc84lmC8xETgc2Ab4dtxBnHOulJVm\nQyD1BK4DfoBZtikEzjnnclCqncUXEdwS+kvcQZxzrtSV3hWBNBQYCewQdxTnnCsHpXVFEHQQ/wb4\nOWbvxB3HOefKQUlcETTtMzAABvSGfk/ANS/EHco558pE4huCLPsMXFMlrfblI5xzrvMSf2vI9xlw\nzrnCSsQVQVszi32fAeecW6sQM4sT0RC0tTHNcvgsS7nvM+CcqzgVuTHNCHhlHKzILPN9BpxzLjrJ\nXmtI2gR4eQz8Yj4c7PsMOOdcoHL2I5BuAFZi5h3DzjmXoTIWnZN2Bw4DhsQdxTnnylky+wikLsAk\nYBxmH8UdxznnylkyGwI4DfgEuCPuIM45V+6S10cg9QXmAvthNie2YM45l2Dl3Vks3QJ8hNnZ8aVy\nzrlkK9/OYukbwIHA0LijOOdcpUhOH4G0LkEH8U8x+zjuOM45VykK3hBI6ilptqSD26n6Y2ApcE+h\nMznnnFurGFcE59HeyV36EsH2k2cSU6dFeiGnWCUhAyQjRxIyQDJyJCEDJCNHEjJAcnJEJaeGQNIt\nkpZImtOifISkVyXNkzQu5HUHAC8D77V1/NHwwh3wBGav5BM+YtUxvneT6rgDpFXHHYBkZIBk5KiO\nO0BaddwBSEYGSE6OSOR6RXArMCKzQMGkr4np8qHAsZKGSDpR0q8lbQnsC+wJHAeMUrDVZCsT4Euz\nYNcqqbatELm2wtnqtSzvSKuehAydzZFrWSEzhJUXMkO2uv6zyL8sqhxJ+FnkmiEpOQp1JZJTQ2Bm\ns4APWxTvDsw3szfMbBVwN3C4md1hZmPN7G0zu8jMxgJ3AjdaG7d9boRtcthspjqXvG3Ua1me6/GS\nlqGzOXItK2SGsPJCZshWN9eyqHJkq9eyvJAZstXNtSyqHNnqtSxPQoak5OhIhnblPI9A0kDgITPb\nPv39d4CDzGxU+vsTgD0szwXiJMU7kcE550pUEuYRRHICj+of4pxzrmM6M2roLaB/xvf9gcWdi+Oc\nc67YOtMQPA8MljRQUjfgGODBaGI555wrllyHj94FPA1sJ2mRpB+Y2efAmcCjBENE77F4h38655zr\ngNgXnXPOORev5Kw15JxzLhaJawjSaxPdJulGScfFmGNrSTdJujfGDIenfw53p2dpx5Xjq5Kul3Sv\npNNjzJHrulWFzFAtaVb657FvTBkk6RJJ4yWdFEeGdI690j+HKZKeiinDAEn3S7o5bHWDImUYKuke\nSZMkHRXD+7c6V+V7Hk1cQwAcCUwzs9MI9iyOhZn928xGxvX+6QwPpH8OpxN0xseV41Uz+1E6wzfi\nykEu61YV3hpgGbAe8Y2SOwLoB6yMMQNm9tf05+Jh4LcxxagCfm9mpwI7x5RhBDDBzM4Ait4wZzlX\n5XUeLUpDkOdaRf2ARen/Xh1jjoLoYIaLCJbziC2HpEMJ/oefHkcG5bhuVaFzALPMrBY4H0jFlGE7\n4Ckz+ynwo6gydCBHk+MIVg+II8PfgFMlPQ40xJThDuB7kq4ANonh/cPkdx41s4I/gL0JWus5GWVd\ngPnAQKAr8E9gCHACcHC6zl1x5ch4/t4YfxYCLgf2i/N30uJ1D8f0s/gl8GuCUWp/JD3QIcbPRbco\nPxt5/iyOB76brnNPnJ8LYADB8jFx/T9yDrB3uk4sv48Wz/8xxve/N+O/8zqPFmWHMjObpWCJikxf\nrFUEIOlu4HBgPDAxfR840nkJ+eSQtAT4FbCTpHFmdnmxMwD7A/sBfSQNMrMbosiQbw4F+0gfSXA7\n5JE4MpjZRenvvw+8Z+lPeLFzSPoqcBCwITAhjgzAdcAESXsDjVFl6ECOV4BTgFtizPAQUJ++D/7v\nODJI+gS4EOgJXBHD+4edq+4jj/NonFtVZl66QHCvcw8z+4TgwxV3jg8I7s3HmWE0EZ5sOpFjJjAz\nzgxN35jZbXHmMLPLgPtjzvApUMz+q6y/EzOrjzODmb0EfDfmDG8CP4zx/Vudq/I9j8bZWZyUCQxJ\nyJGEDJCMHEnIAMnIkYQMkIwcnqGA7x9nQ5CUtYqSkCMJGZKSIwkZkpIjCRmSksMzFPD942wIkrJW\nURJyJCFDUnIkIUNSciQhQ1JyeIZCvn+Uvf1t9IDfBbwNfEZwj+sH6fIa4DWCnvALKiFHEjIkJUcS\nMiQlRxIyJCWHZyj++/taQ845V+GSOLPYOedcEXlD4JxzFc4bAuecq3DeEDjnXIXzhsA55yqcNwTO\nOVfhvCFwzrkK5w2Bc85VuP8PJfdwqj4JoSUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105b842e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(counts, results, 'o-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the second plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "absorb_weights = np.array([\n",
    " # probability of absorbing a photon in state 0\n",
    " 1,\n",
    " # probability of absorbing a photon in state 1\n",
    " weights[1],\n",
    " # probability of absorbing a photon in state 2\n",
    " weights[2]\n",
    "])\n",
    "emit_weights = np.array([\n",
    " # probabilities of emitting a photon from state 0 back to...\n",
    " weights[[\n",
    "  0, # state 0 (P_00)\n",
    "  1, # state 1 (P_10)\n",
    "  2  # state 2 (P_20)\n",
    " ]],\n",
    " # probabilities of emitting a photon from state 1 back to...\n",
    " weights[[\n",
    "  0, # state 0 (P_01)\n",
    "  1, # state 1 (P_11)\n",
    "  2  # state 2 (P_21)\n",
    " ]],\n",
    " # probabilities of emitting a photon from state 2 back to...\n",
    " weights[[\n",
    "  0, # state 0 (P_02)\n",
    "  1, # state 1 (P_12)\n",
    "  2  # state 2 (P_22)\n",
    " ]]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "counts = [int(pow(10,n/3.)) for n in range(30)]\n",
    "results = np.array([compute_with_transition_matrix(0, n, absorb_weights, emit_weights) for n in counts])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10617ee80>,\n",
       " <matplotlib.lines.Line2D at 0x105d02208>,\n",
       " <matplotlib.lines.Line2D at 0x105d02390>]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXEW5x/Hvm31fWYQwMYSwhQHZkStIlCXJICASNgHZ\nAoKQIIgGBJ2MVxEUxCyCLIIQFBDkXgEhiGgQRLgEBVk1CVsCEsMSSSAhC+/94/SQnpnuSffM6a7q\n7t/nefqBqak555dO0pVzqk695u6IiEjt6hI6gIiIhKWBQESkxmkgEBGpcRoIRERqnAYCEZEap4FA\nRKTGaSAQEalxGghERGpct1Ie3Mz6AlcAHwBz3P2XpTyfiIgUr9RXBF8AfuXupwIHl/hcIiLSAUUP\nBGZ2nZktNrOnW7WPM7MXzGyemU3JNA8DFmb+f20ns4qISAl05IrgemBcdoOZdQVmZtpHA0eb2bbA\nIqCuE+cSEZESK/rD2d0fAt5p1bw7MN/dX3b31cAtwCHAHcBhZnYFcGdnw4qISPrSmizOvgUEyZXA\nHu7+PnBSez9oZtr+VESkA9zd0jhOWrdrOvlh7pnX+NnubtkvGHffuu83tukL3hN8KBzzaO5+5ywG\nvx/8ObhgJfgH4C+B/xl2Wd7esdtm2e5AGL4ADn8wybXdgbn6ZX5zmgppz/46388Uc/xC+q2vrTPH\nLvd70dn3TO9Fbb8XhWaIJUerDKlJ64rgNdbNBZD5/0WF//hU4NHXYdGMtt9bOB1O2QKuGQVjMm0T\nF8CrMwDcWQW8bfbWf9b9zJisn3/+b+6MBzD73hj47uPAJsAw2PAGoF/bc47exYyzgceAv7mzwqy+\nAfacBsdsAWO2SPqdsoVZPe7P3JPjFzUnzy+2dXu+futT6M/l6ldoWykz5GovZYZ8fQttSytHvn6t\n20uZIV/fQtvSypGvX+v2GDLEkmOOmY2h5Ydc57l70S9gBPB01tfdgAWZ9h7Ak8C2BR7LYfy9sF1D\n/j7bNSR9Dp+Tr2/SZ+I8cF/3Onl++8cdN3td38asnzvh7+A/AZ8L/h74E/CVV1seu/k1/t6OvId5\n3oupaR2r0nPEkCGWHDFkiCVHDBliyZF8fKdzrKKvCMzsZmAfYKiZLQS+7e7Xm9mZwH1AV+Bn7v58\n4YPRPePb//4z9wC5/tXdoo9ZPdAwCfr1huUr4NUZef61npHvauPx89yT85nRG9gJuDH3Mfr1bi9X\nkeakeKzOmBM6AHFkgDhyzAkdIGNO6ADEkQHiyZEKy4ws4QIkk8VNJE8ezyn/+esbYPh6Bw+z8bPh\n3rFtj/C1t+CyrwG/dmd5GSKLSA3LujXU6CnNFUQxEKT1iymldXME14xa13rKi7DlLfCN7YFPkyyR\nvRH4I9SPhbrJ0L8XLFsJC6e3f3UiIlK4ND87NRAUob2rBzM2Ao4Gjoffbwa/7QaXD17306fMh7+c\npcFARNJQdQMBAW8NlYLZ4Q/DbZ9q+52G2eubDxERaU8pbg1Fse2Du0+tlkEgYWtyt398pBm9yptF\nRKqJu89x96lpHjOKgaD6LFuZu73XYOAlM843Y3DuPiIi5aWBoCQWTk/mBLJNXAD3nwAcAGwDLDDj\nMjPqzOobzMbPNjtiTvLf+obyZxaRWlXSwjSFMrOpVNEcQQHPNBxvRh3wVfjDs8mmrZf2X3eEdp9Y\nFpEaVooni6OYLK6UVUOlYHbwA3DnZ9t+RxPLIpJfmp+dujUUXK+uudt3+awZvzDjJDM+nv0d6z+k\n0eoGLbEtBi+1ukFLrP+QxnxHL6aviNSmKG4N1bZ8E8vPPkryGPsBwMVmvAs8QO+DBrPF6kOZsHzd\n793t/S6w/kPwZW83ZR/B+g9pZMTqC5iwvPv6+n7Uf9CHZ9LDurPKV7O0y8xc/USkuujWUGC5n1ie\nuAAenZz1sFoXoB7Yl2F1F3PKoh5tDnR7vzVsO+4Ruq1cRddVH9Bt5SqeeLyBw1b0bNP35j7vcvT7\nhwHLgHeBd5k58Ct8bM3XmfBe9qCxmpe7f0+Dhkh80vzsjOKKoNomi4vh/sw91n/Ibgy6LftDdZYv\ne3vdRPFU2xz4FLA3D9A954G6OcBaVg7sw9qeg1jbvSdd/pa77wAfwMv7/JZeS52e70KPZcbWb/Zg\n/1b9Jizvzh29zrdzN11DnyX/pOuaRcAiLh90MiPWfLPQKw0RSY8mi6uQbWwNbMY0DmbdFcFdvMhg\nbmEvhpLcGuoF/A64n+v7T+fEZUPaHOjagW/6wqUbtjh23aAlTPzPBm37DnqTRe/UAwMyr/5svtHd\nHL+kb5u+dwx06o7+gH5v9GDQK2vo/xo8+u8e7JfjF3N9/7f85XfbnE9XDyLpq7orgpo2gMktBgGA\ngxjJbE4FLgJmAM95YzJi26VDRnF7v5b3/W/ru5qlXWa2OfbSLjNz97WZ7iwGFjc3W92qFUDbgeBt\n3vKnrtwwsx33psAwNt/oXljSp03f4cuH2vkD3+P9DV7kw25z6ffGA9z24e6M8NNa3XLS1YNIRHRF\nEJhtZ3M4gn3afONXPOjP+picP1PEv7AL7ZtzYvm2vqt5pUebOYK8VxrXbPAf7M7pbPzUXgx+cTRD\n523IK3d1YdzatsFyXMGISOGqbtO5mh4IdrSnOJQd2nzjJmb7PC/rcwRpDxpm9GTzDd7k+LfalgP9\nXVdnx61/Tf/Xb6H30ge80ZcWk0Gk1lXdQECV7T5aCGuynsBlLOBQnmctn8uq+XwnC1jEZF/s0T5Z\nXPCgke/q4cahKxl50hI2fWJT6v4C72+wmKdWvsU7y0ZzyMp1z1a0s3JJpBapME2VsCYbDtwGvA6c\nwBV8igFMoge9WcUK3mVGzINAMdZ39WBGf/q+MYZRs49i7VlHMeHdtg85zuq3lOOWD/NGfz/n8Qu8\ngti0tzV+3Dmzr9H9PWf1K8bM11d4VfSNJUel9Y0lR7GZIeXPzmoqwFwJL6Yylqm8wVTOZWoyEFf7\ni36DG9ls4BJGDlrKZgOX0G9wY85+IwctZSre5rVnd+fCnms5Z9O3OL3+L3x5xx8zZXADmwy4pN/w\nXqt33QLf5+P4rlvg/Yb3Wp3r+Jv0ovHIHqxx8ObXkT1Ys0kvKr5vLDkqrW8sOYrN/NHflxQ/O3VF\nUCbWZF2AbwGnAkd7o/8pcKTo5F/uOvgt3n7hWEY8uD9D//lJBr28NUPnDe736INdGhbDre+s63rk\nYLhnQI9Vy8fW/xzzZdjapXRZ85/df/7cJY+9T+/Wh96/Jyv23Y3TzDGALg73z+XK361q2/eAnqwY\nvxOnZ7fd+9fwfWPJUWl9Y8mRr+8ne/KfR1f6oNbtzbR8tALYxtbAACbTg16s5UN2ZgBb8x6wqzf6\nv0Lni1Le5a7M8Pc2mg2Hz25uNqPXNr26L7915ZoWezXd+g58admqHl/4zTvHDV39frcNVq/oOnTV\nyi4zP8h9yh1W0fuA5/vcAND8T6I31rS5AwVA/Wp67zW/z8+z2xZG0DeWHJXWN5Yc+fr2pkvb5dyl\nktalRUdfVOGtITaigZ2Z1+L2xid5h004KHS22F+F3kZy6H5kV2txOd38OqarrXU4yuGzDvUOG+/e\ns9uqXH1369ltFXivrFfv3drv2yf7FUPfWHJUWt9YcrTXt92/K+Cp/b0L/hc/xV9MLC9GMTvnve5R\n3Bs6W8W/oL/D2Q6vTIL3c/0F2su6zGv9cwO61f3ziO7dWvQ7vHs3H9Ct7p+V3jeWHJXWN5YcxWZu\nfqX52RnFraGq22uoR566xD3a3geUlurNGupgcn/otQxWLoTpz7jfg9kmwGTgFOABYMIfYMPjYNYs\n+GjLjePg7Xf8w7NaH/fdNQO+em/XQTfu1vP1oX1Zw3t04x9rNn172doPv1rpfWPJUWl9Y8lRbGbt\nNVQhbFubw1E5nhYO8JBYJak3a9gTpl3Dui03TodXPwMvHAG7Ab8ALsf9xeyfGQ6T+kHv5bDiVZjx\njOdeepvs9Do8X9W4iu4bS45K6xtLjmIzJz9TZQ+UVdNAYE3Wj/k8xTwGMp6hH32jAh4SC2282ex7\nYWzr9i/DvKtgT9zfCpFLJEZaNRSpzBLRWYziQX7H7dxUnQ+JlUp/ct9Sewde1yAgUjoaCNJ1Ecn9\n6iN9sa8C9MFfhGWQs1rbclhR7iwitUQ1i1NiTXYCMAE4zBt9VeA4lces64lgU1oNBhNhwavJVtwi\nUiKaI0iBNdmnSfYO2scb/YXQeSqOWQ/gJmDIZ+HKXjCxkMlfkVqmyeKIWJNtAfwZOM4b/f7QeSqO\nWR/g1yRXAkfjnvP2kIi0lOZnp24NdYI12SDgbuA7GgQ6wGwgcB+wBDhcg4BIGBoIOsiarBtwK3C/\nN/oVofNUHLMNgT8ATwEn4L4mcCKRmhXFqqFKebK4xUZyfdmcUSxha84JnavimG0G3E9yS+hbhL4/\nKVJB9GRxQLaxNbAZ01oUmr+LF1nIJD0f0L7sbSPWQJevwFb7wWW4/zB0NpFKpTmCEAYwucUgAHAQ\nIxnApECJKkLzthH3wthfwT53wN53gdXDs6GziUhCA0GhtJFch9TB5Oy9gwCmwUbD0QAqEgsNBIVa\nRe7SJqv01Gt78m0b0Q8NoCKx0EBQqK15iftbbYFwJwt4V0+9tmc5uQdQbRshEo8oVg3FzppsNLsx\ngd8wiZs4TBvJFe4IePMbsOIHWVcA2jZCJC5aNbQe1mQ9gMeAK73Rrw6dp6KYfQ648jA4bwUcq20j\nRNKjbajLqwlYBFwTOkhFMdsSuA445NfufyEpKiMiEdJA0A5rsr2BE4AdvVEPPRXMrC9wB9BIMgiI\nSMR0aygPa7IBJNsfTPZGvyt0nophZsAvSTaRO0lPDYuUhm4Nlcd0kn2ENAgU56vAVsBeGgREKkNJ\nBwIz2xy4ABjo7oeX8lxpsiY7DPgUsFPoLBUl2QPlPGAP3LU8VKRClPQ5And/yd0nlvIcabMm2xS4\ngqS+wPLQeSpGspHcL4FjcX85cBoRKUJBA4GZXWdmi83s6Vbt48zsBTObZ2ZTShOxfKzJjGSly0+9\n0R8NnadimPUEbgdm4KrLIFJpCr01dD3JA0A3NjeYWVdgJrAf8BrwuJndCewK7Az80N1fTzduyZ1B\nUnz+u6GDxC57R9ENYeT+sPDzcHHoXCJSvIIGAnd/yMxGtGreHZjvmdsAZnYLcIi7XwzMyrQNAS4C\ndjSzKe5+SUq5U/NRjYE+DGEAO/IeZ/hffXXoXDFr3lE0ezO5L8PqC2H8M6AHxUQqTGcmi4cBC7O+\nXgTskd3B3d8GTlvfgTKFaZqVrUBNzhoDd/IN29he09YR+eXaUfQqGNmQ7Ciq902kBEpRkKZZZwaC\n1JYGuvvUtI5VlFw1Bg5mFDfpA6092lFUpPwy/0Ce0/y1mTWmdezODASvAXVZX9eRXBUULVipStUY\n6JBlkLPIvHYUFSm9UlwZdGb56FxgSzMbYWY9gCOBOztyIHefGqResWoMdMh/wV/OgxbF5rWjqEh5\nuPuctO+iFHRFYGY3A/sAQ81sIfBtd7/ezM4E7gO6Aj9z9+fTDFdy2/Ay97OC/bOuAFRjoH1mG34L\nvnwBXNgAY7SjqEjli2KvIZIdPst6a8iabAQwl99wPsv4gmoMFCDZR+gOYB7u3wgdR6QWZd0aakxr\nr6EoBoJybzqXeXBsNvBHb3StfS+U2UnAZJItJHLfVhORstCmc513LLAxcFnoIBXDbCRwCfAZDQIi\n1SWKgaCcq4asyTYCLgUO9EY9OFYQs24kDwlehPszoeOI1LJSrBqquVtD1mS/AP7ljX5uuc5Z8cwu\nAD4DHID7h6HjiIhuDXWYNVkD8Elg+9BZKobZrsBZwM4aBESqU80MBNZk/YErgZO90d8PnacimPUB\nbgIm496hhwVFJH5RDARlmiP4HvAHb/Tfl/Ac1eYHwBO43xI6iIgkNEfQ0XM02Z4k69+380Z/u5Tn\nqhpm44CrgE/gvjR0HBFpSXMERbAm6wFcC5ylQaB9zTUGBkP/4bBLN5j6XQ0CIlWv6gcCkhq6C4Db\nQgeJWa4aA6fAyfVmf9fWESLVraQ1iwtlZlMz973SPW6TjSbZI/8r3hj4HljkctUYuAZGDU/ePxGJ\nhJmNaVXDpdOiuCIoRT0Ca7IuJLeEGr1RK17WRzUGRCpDc12CNOsRRHFFUCKnAx8CPw0dpBIsJ/eW\n3KoxIFL9orgiSMtH9Yd7M5CB7MIKJvlcPQRViIPgxSmw8pKsKwPVGBCpDVWzfDRP/eH5LOIsbSu9\nHmZbAX8+Cc5/Aw5TjQGR+KW5fLR6BoItbTbHMrbNN25its/z8Z09ftVKNpR7GJiF+09CxxGRwlTd\ncwSpPFms+sMd9Q1gOcn2GyISuVI8WRzFQJDKqqFVuQuqq/5wO8w+AZyNNpQTqRhaNdSejbmf37O2\nRZvqD+dn1hO4ETgX94Wh44hIOFUxR2BN1h2YyyPM5kV2UP3hAphdBIwGDiX0HwIRKVrVzRGkYDLw\nb/6L8/w+faitl9mewEkkG8rp/RKpcRU/EFiTDQfOB/bUNhIFMOsL3ACcgfvi0HFEJLxqmCOYBkz3\nRp8XOkiFuBj4P9x/HTqIiMQhiiuCji4ftSY7mOQ+91EliFV9zPYFPg/sEDqKiHSMCtNk/1yT9QWe\nJSk9+UD6yapDc42BQdB3OOzaB77X6P7d0LlEpHM0WZxoBB7WIJBfnhoDx9eb/VVbR4hIs4qcI7Am\n2x44Efha6CwxU40BESlExQ0EmToDPwW+5Y1a9dIe1RgQkUJU3EBAsv69C3B16CCxU40BESlERQ0E\n1mQbAhcBp3mj9sZZn4PhpSm03INJNQZEpLWKWDX0UcGZDdmR1azgVc7Q1hHrYbYD8MAJcN6/YYJq\nDIhUl5qqR6CCMx1g1huYC/wA9xtCxxGR9KU5EMR/a2gAk1sMAgAHM4oBWvnSjh8CT5PsLioi0q4o\nniNo98liFZwpjtnngM8BO2pDOZHqU5uFaT4k96WPCs60ZfYx4BrgCNyXho4jIumrzcI0O9KD2bzd\nok0FZ9oy6wL8HLgW94cCpxGRChLFFUE+1mQHsQ1DeZiTuInTVHCmXZOBgcB3QgcRkcoS7aoha7I+\nJJvKneqNfn/5k1WQpPbw74E9cH8xdBwRKb1aWTX0TeAxDQLrkSwV/SVwjgYBEemIKG8NWZNtDZyG\n9s0vxKXAU8BNoYOISGWKbiCwJjPgJ8B3vdFfD50nRs01BjaDYUNgiyfh+PtC3+MTkYoV3UAAHAls\nCMwMHSRGeWoMXFRv9p62jhCRjohqjsCabCBwGXC6N/qa0HlipBoDIpK2qAYCoAmY7Y3+SOggsVKN\nARFJWzS3hqzJdgS+SFKMXvLoCn1ytavGgIh0VMmvCMzsEDO72sxuMbP9c/ZJqo5dAVzgjf5mqTNV\nLLNNJsIWZ0CLSXTVGBCRzijbA2VmNgi41N0ntmp3duQZtqMbW7KdCs7kYdYNeAB4oB7mDodJqjEg\nUruC1CMws+uAA4F/u/v2We3jgB8DXYFr3f2SPD9/KXCTuz/Zqt2ZCtzNQl7lNG0dkYfZ94GdgfG4\nBkuRWhfqyeLrgXGtgnQlWeY5juTe/tFmtq2ZHWdml5vZppa4BLi39SDQwueoU42BPJKtpY8BjtUg\nICJpK3iy2N0fMrMRrZp3B+a7+8sAZnYLcIi7XwzMyrRNBvYFBpjZKHe/qs3B/5j579tsa2ZjctYl\nqFXJe/4z4FDcl4QNIyKhlKIOQbPOrhoaBizM+noRsEd2B3efDkxv9yifyfz3NZ73tzQIfMSsJ/Ar\n4BJcS2pFallzHYLmr2OqR5DeTLNqDORyKcngennoICJSvTp7RfAaUJf1dR3JB1dxLmc+q7nW39NE\n8UfMjgTGA7uq5KSINCvFLaKilo9m5gjual41ZMmSxn+QzAG8DvwfcLS7P1/EMVOb+a4aZlsDDwMH\n4P630HFEJD5BVg2Z2c3AI8BWZrbQzE509zXAmcB9wHPArcUMApKDWR/gduACDQIiUg5RVCgj2WNo\nTq2uFmreVro/9NoIRo2Bf06AfXVLSERay7o11Fj2B8pKpdZvDeXaVvpUePERmKSnhUUkn1opVVkT\ncm0rfTWM1LbSIlIuUew+amZTqdFbQ9pWWkSKUYpVQ1EMBO4+NXSGUJbBylzt2lZaRHJpfrAspgfK\npJPGwrPnwarsNm0rLSLlpMnikMy2Bf70Ffj2y3CwtpUWkUIF2Ya6VGp2+ahZb5IH8Kbhfm3oOCJS\nGbR8tJqYXQX0B47R8wIiUqw0PzujmCyuOWZHAZ8FdtEgICKhaSAoN7NRJBPBY3F/N3QcEZEoBoKa\neY4gqS9wK9CE+19DxxGRyhN899FSqKk5ArPpJMV8JuiWkIh0huYIKpHZocBBwE4aBEQkJhoIyiGp\n43AVcBDuS8OGERFpSQNBiWRtLd27Dj6xNdx2qvtjoXOJiLSmOYISyLW19Ckw/y9wlp4YFpE0VN02\n1GY2NTMTXhVybS19DYzS1tIi0llmNiaz0jI1UdwaqrbdR7W1tIiUinYfrRDaWlpEKokGghL4LDx5\nHqzObtPW0iISK00Wp81sE+DJr8H3n4ex2lpaREqh6rahrpqBwMyAu4G/4v6t0HFEpHpV3aqhKnIy\nsAnw36GDiIgUKopVQ1Wx6ZzZ5sD3gc/gvmp93UVEOkKbzsXKrAvwB+C3uP8wdBwRqX66NRSfs4Cu\nwI9CBxERKZauCDrLbDTwJ2AP3BeEjiMitUFXBLEw6w7cCFyoQUBEKpUGgs75JvAmyRbTIiIVKYpV\nQxXJbBfgDFRoRkQqnAaCIjTXGBgAfYbDTh+DK7/m/lroXCIinaGBoEB5agwcWm82R1tHiEgli2KO\noBLqEajGgIjEQPUIAlKNARGJgeoRBKQaAyJSrTQQFGgbuPd8WJvdphoDIlIN9GRxIZIHx+b+CH73\ne6hXjQERCU31CMrN7HxgH2C8nhkQkRhoICgns62BPwO74P5K6DgiIqC9hson2V76auA7GgREpFpp\nIGjfKUBP4Cehg4iIlIpuDeVjNgx4EhiD+7Oh44iIZNOtoVJLitBfAfxEg4CIVLsoniyO0ARgS+CI\n0EFEREqtpAOBmW1DUsZxA+ABd/9pKc+XCrMhwDTgMNw/CB1HRKTUyjJHYMnqmxvc/bgc34trjsDs\nOmA57pNDRxERyafscwRmdp2ZLTazp1u1jzOzF8xsnplNyfOzBwF3A/E/gWu2H7AvcEHoKCIi5VLQ\nFYGZ7Q0sB2509+0zbV2BfwD7Aa8BjwNHA7sCOwM/dPfXs45xt7t/Lsexg14RNBebGQh9Pg679IVL\nvu3+nVB5REQKkeZnZ0FzBO7+kJmNaNW8OzDf3V/OhLoFOMTdLwZmZdr2Ab5Ashb/t2kETlOeYjPH\n1ZvN1R5CIlIrOjNZPAxYmPX1ImCP7A7u/iDw4PoO1KrIwpzMftsll6/YTENSbEYDgYhEI1O8a0wp\njt2ZgSC1WeZQhWlUbEZEKkVzQZrmr2MpTPMaUJf1dR3JVUHFULEZEZHODQRzgS3NbISZ9QCOBO7s\nyIFC1SweCL/8JnyY3aZiMyISs1LULC501dDNJPvxDwX+DXzb3a83s/HAj4GuwM/c/ftFBwi5asjs\njp/Dil/BEBWbEZFKonoE6Zz4IOAyYAfcc94iEhGJVdmXj5Za5jKnbKuFMOtLcvvnZA0CIlJJSrF6\nqDavCMx+AGyK+7FlPa+ISEqq7oqgrMy2B04Atg+cREQkClEMBGW7NZRsfncV8C3cF5f0XCIiJaBb\nQ50/2akkVwN74f7henqLiERLq4Y6dqKNgGeA/XD/e8nPJyJSQhoIOnaiWcAbuH+95OcSESmxqpss\nLvkcgdm+wKeB0SU5vohImWiOoGMn6AX8HTgX9w5tgSEiEpuyVyircFOAZzUIiIjkVt1XBGZbAY8A\nO+G+cH3dRUQqheYICjwocAVwkQYBEakWmiMo7sDHAF8HdsV9TerHFxEJSMtH13/QwcBzwOdxfyzV\nY4uIRKDqbg2lpd6soQ4mbws7fACrH4Shz4QOJSISuaoZCOrNGvaEadnF6E+BafVmqNCMiEh+USwf\nTaNUZR1Mzh4EAK6BUcNhUqfCiYhEpBSlKqO4InD3qZ09Rn/olau9H/Tu7LFFRGKRWV05x8wa0zpm\nFFcEaViVp305rChrEBGRClM1A8GXwM6Gd7LbJsKCV5OSlCIikkd1LB81OwC46tNwTj84tR/0Xg4r\nXoUZmigWkWqk5whaHqAn8DRwDu53pxZMRCRi2nSupa8Dz2kQEBHpmChWDXV4ryGzkcDZwM7ppxIR\niY/2Gmr1g8DdwMO4fz/1YCIiEdMWE4lDgJHAoaGDiIhUssocCMz6AtOAE3HP9wiBiIgUoFIniy8k\nuSX0h9BBREQqXeVdEZiNBiYCO4SOIiJSDSrriiCZIP4J8B3c/xU6johINaiIK4LmOgPDYXh/GPZH\n+NEToUOJiFSJ6AeCPHUGflRvtlbbR4iIdF70t4ZUZ0BEpLSiuCJo78li1RkQEVmnFE8WRzEQtFeY\nZjl8kKdddQZEpObUZGGacfD8FFiZ3aY6AyIi6Yl7ryGzocBzk+G/58OBqjMgIpKonXoEZlcBq3DX\nxLCISJba2HTObHfgYGDb0FFERKpZnHMEZl2BK4ApuC8NHUdEpJrFORDAqcD7wKzQQUREql18cwRm\nGwHPAPvi/nSwYCIiEavuyWKz64CluJ8TLpWISNyqd7LY7FPAAcDo0FFERGpFPHMEZt1IJojPxf3d\n0HFERGpFyQcCM+trZo+b2YHr6XoG8CZwa6kziYjIOuW4IvgG6/twN9uEpPzkmQSatMhs5BRUDBkg\njhwxZIA4csSQAeLIEUMGiCdHWgoaCMzsOjNbbGZPt2ofZ2YvmNk8M5uS4+f2B54DlrR3/EnwxCz4\nI+7PFxM+ZWMCnrvZmNABMsaEDkAcGSCOHGNCB8gYEzoAcWSAeHKkotArguuBcdkNljz0NTPTPho4\n2sy2NbPjzOxyM9sU2Af4JPBF4BRLSk22MQM2eQh2qTdraC9EoaNwvn6t2zsyqseQobM5Cm0rZYZc\n7aXMkK/CeDgjAAAEY0lEQVSv3ovi29LKEcN7UWiGWHKU6kqkoIHA3R8C3mnVvDsw391fdvfVwC3A\nIe4+y93PdvfX3f1Cdz8b+CVwtbdz2+dqGFlAsZkxheRtp1/r9kKPF1uGzuYotK2UGXK1lzJDvr6F\ntqWVI1+/1u2lzJCvb6FtaeXI1691ewwZYsnRkQzrVfBzBGY2ArjL3bfPfD0BGOvup2S+PhbYw4vc\nIM7Mwj7IICJSoWJ4jiCVD/C0fiEiItIxnVk19BpQl/V1HbCoc3FERKTcOjMQzAW2NLMRZtYDOBK4\nM51YIiJSLoUuH70ZeATYyswWmtmJ7r4GOBO4j2SJ6K0edvmniIh0QPBN50REJKx49hoSEZEgohsI\nMnsT3WBmV5vZFwPm2NzMrjWz2wJmOCTzPtySeUo7VI5tzOxKM7vNzE4LmKPQfatKmWGMmT2UeT/2\nCZTBzOx7ZjbdzL4UIkMmx16Z9+EaM/tzoAzDzex/zOxnuXY3KFOG0WZ2q5ldYWaHBTh/m8+qYj9H\noxsIgC8Av3L3U0lqFgfh7i+5+8RQ589k+E3mfTiNZDI+VI4X3P30TIZPhcpBIftWld6HwDKgJ+FW\nyX0eGAasCpgBd3848+fibuDngWLUA7e7+8nAToEyjANmuPtXgLIPzHk+q4r6HC3LQFDkXkXDgIWZ\n/18bMEdJdDDDhSTbeQTLYWYHkfyFvydEBitw36pS5wAecvcG4DygKVCGrYA/u/u5wOlpZehAjmZf\nJNk9IESGx4CTzewBYHagDLOAo8zsB8DQAOfPpbjPUXcv+QvYm2S0fjqrrSswHxgBdAeeBLYFjgUO\nzPS5OVSOrO/fFvC9MOASYN+Qvyetfu7uQO/Fd4HLSVap/S+ZhQ4B/1z0SPPPRpHvxTHA4Zk+t4b8\ncwEMJ9k+JtTfka8Be2f6BPn9aPX9/w14/tuy/r+oz9GyVChz94cs2aIi20d7FQGY2S3AIcB0YGbm\nPnCqzyUUk8PMFgMXATua2RR3v6TcGYD9gH2BAWY2yt2vSiNDsTksqSP9BZLbIb8NkcHdL8x8fTyw\nxDN/wsudw8y2AcYCg4AZITIA04AZZrY3MCetDB3I8TxwEnBdwAx3AVMz98FfCpHBzN4Hvgn0BX4Q\n4Py5PqvuoIjP0ZClKrMvXSC517mHu79P8ocrdI63Se7Nh8wwiRQ/bDqR40HgwZAZmr9w9xtC5nD3\ni4H/CZxhBVDO+au8vyfuPjVkBnd/Fjg8cIZXgC8HPH+bz6piP0dDThbH8gBDDDliyABx5IghA8SR\nI4YMEEcOZSjh+UMOBLHsVRRDjhgyxJIjhgyx5IghQyw5lKGE5w85EMSyV1EMOWLIEEuOGDLEkiOG\nDLHkUIZSnj/N2f52ZsBvBl4HPiC5x3Vipn088A+SmfDzayFHDBliyRFDhlhyxJAhlhzKUP7za68h\nEZEaF+OTxSIiUkYaCEREapwGAhGRGqeBQESkxmkgEBGpcRoIRERqnAYCEZEap4FARKTG/T/PnRh9\nwEtkygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105f95e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(counts, results, 'o-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}