{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# BAC: Block Arithmetic Coding\n",
    "\n",
    "I developed BAC (Block Arithmetic Coding) back in the 1980's and 1990's.  For several reasons, I've become interested in this problem again.  \n",
    "\n",
    "This notebook discusses BAC and presents code to calculate its coding efficiency. BAC was presented in Ref [1] below.\n",
    "\n",
    "Note: see companion notebook for encoder and decoder code.  \n",
    "\n",
    "-- Charles Boncelet, <boncelet@udel.edu>, 3/29/17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#python 3.6\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#utility function\n",
    "def entropy(p):\n",
    "    \"\"\"binary entropy function with parameter p\"\"\"\n",
    "    return -p*np.log2(p)-(1-p)*np.log2(1-p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## BAC (Block Arithmetic Coding) Basics\n",
    "\n",
    "BAC (Block Arithmetic Coding) is a variable to fixed (VtoF) encoding.  Here we describe the coder for IID binary inputs, but BAC works perfectly well for non-binary, non-stationary, and/or Markov inputs.\n",
    "\n",
    "For each input bit $x$, BAC recursively divides the set of available output codewords into two non-empty groups, one for $x=1$ and one for $x=0$.  The size (number of codewords) in each group is proportional to the probability of that group, i.e., $k$ codewords are divided into $k_0$ and $k_1$ with $k_0+k_1=k$ and $k_1 \\approx pk$ and $k_0 \\approx (1-p)k$.  The recursion continues until the group contains only one codeword.  That codeword represents that input phrase.\n",
    "\n",
    "The expected length of a BAC input phrase is\n",
    "\n",
    "$$ E(L(k)) = 1 + p E(L([pk]) + (1-p)E(L(k-[pk]))$$\n",
    "\n",
    "where $k_1=[pk]$ is a rounding of $pk$ with the restrictions that $1\\le k_1 \\le k-1$.  This also guarantees $1 \\le k_0 = k-[pk] \\le k-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## BAC Example\n",
    "\n",
    "Let $p=\\Pr(x=1)=0.3$ and let $k=16$ (i.e., 4 bit codewords).  Index the output codewords from 0000 to 1111, or, more simply, from 0 to 15.\n",
    "\n",
    "Generate some random bits with probability $p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#uncomment to generate bits\n",
    "#p = 0.3\n",
    "#bits = 1*(np.random.rand(10)<p)\n",
    "#bits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The sequence I got was [1, 0, 0, 0, 0, 0, 0, 0, 0, 0].  Index the output codewords from 0 to 15.\n",
    "\n",
    "1. The first bit is a 1.  Let $k_1 = [0.3\\times 16] = [4.8] = 5$.  Thus, $x=1$ corresponds to the last five codewords, $\\{11, 12, 13, 14, 15\\}$.\n",
    "\n",
    "2. The second bit is a 0.  Let $k_0 = [0.7\\times 5] = [3.5] = 3$.  (Both $k_0=3$ and $k_0=4$ are the same distance from $k_0=3.5$; we arbitrarily choose $k_0=3$.)  Thus, the set of codewords is now $\\{11, 12, 13\\}$.\n",
    "\n",
    "3. The third bit is a 0.  Let $k_0 = [0.7\\times 3] = 2$.  The set of codewords is $\\{11, 12\\}$.\n",
    "\n",
    "4. The fourth bit is a 0.  Since the set only has two codewords, the division is trivial: codeword 11 is chosen (if the bit were a 1, codeword 12 would be chosen).\n",
    "\n",
    "Thus the input sequence 1000 is encoded by codeword 11 = 1011.  \n",
    "\n",
    "The full codebook is below:\n",
    "\n",
    "|Index|Input Phrase|Output Code|\n",
    "|-----|-------------|----------|\n",
    "| 0  | 0000000 | 0000 |\n",
    "| 1  | 0000001 | 0001 |\n",
    "| 2  | 000001  | 0010 |\n",
    "| 3  | 00001   | 0011 |\n",
    "| 4  | 00010   | 0100 |\n",
    "| 5  | 00011   | 0101 |\n",
    "| 6  | 0010    | 0110 |\n",
    "| 7  | 0011    | 0111 |\n",
    "| 8  | 0100    | 1000 |\n",
    "| 9  | 0101    | 1001 |\n",
    "| 10 | 011     | 1010 |\n",
    "| 11 | 1000    | 1011 |\n",
    "| 12 | 1001    | 1100 |\n",
    "| 13 | 101     | 1101 |\n",
    "| 14 | 110     | 1110 |\n",
    "| 15 | 111     | 1111 |\n",
    "\n",
    "Note, usually we don't generate the codebook.  We generate the parsing for the specific input on the fly.\n",
    "\n",
    "The expected input length can be calculated by the standard method:\n",
    "\n",
    "$$ E(L) = 7(1-p)^7+7(1-p)^6p + 6(1-p)^5p^2 + \\cdots + 3p^3 = 4.41 \\text{ bits per input phrase}$$ \n",
    "\n",
    "Since each codeword uses 4 bits, the normalized rate is $4.41\\times h(p)/4 = 0.972$, about 3% worse than entropy. The BAC code is slightly worse than the entropy (as it must be)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sum(p) =  1.0\n",
      "EL =  4.412719 bits per input phrase\n",
      "Coding Efficiency = 0.972 and entropy = 0.881\n"
     ]
    }
   ],
   "source": [
    "p = 0.3\n",
    "lengths = np.array([7,7,6,5,5,5,4,4,4,4,3,4,4,3,3,3])\n",
    "ones = np.array([0,1,1,1,1,2,1,2,1,2,2,1,2,2,2,3])\n",
    "pvector = p**ones * (1-p)**(lengths-ones)\n",
    "EL = pvector @ lengths\n",
    "Eff = EL*entropy(p)/4\n",
    "#Check: the p's must sum to 1\n",
    "print('Sum(p) = ', np.sum(pvector))\n",
    "print('EL = ', EL, 'bits per input phrase')\n",
    "print('Coding Efficiency = %.3f' % Eff, 'and entropy = %.3f' % entropy(0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Recursive Algorithm for Expected Input Phrase Length\n",
    "\n",
    "The standard algorithm requires us to actually compute the codebook.  For large codebooks, this is expensive (slow).  Fortunately, the recursive definition of the BAC code yields a recursive solution to the expected input length:\n",
    "\n",
    "$$ EL(k) = 1 + pEL(k_1) + (1-p)EL(k_0)$$\n",
    "\n",
    "with the boundary conditions $EL(0)=EL(1)=0$.  \n",
    "\n",
    "E.g. $EL(16) = 1+0.7 EL(11)+0.3EL(5)$, and $EL(11) = 1+0.7EL(8)+0.3EL(3)$, etc.\n",
    "\n",
    "The code below also gives 4.4127."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## BAC code for calculating Expected Input Phrase Length\n",
    "\n",
    "The BAC class below works with binary inputs and calculates the expected input phrase length of a BAC encoding with $k$ codeworks, i.e., each input phrase is encoded with $\\log_2(k)$ bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#reset if recursion problems arise, default = 1000\n",
    "#import sys\n",
    "#sys.setrecursionlimit = 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "class BAC:\n",
    "    def __init__(self,p,loopfirst=0):\n",
    "        self.p = p\n",
    "        self.els = {0:0,1:0}\n",
    "        #initialize dictionary with loop for small k\n",
    "        if loopfirst > 0:\n",
    "            self.calc_loop(loopfirst)\n",
    "        \n",
    "    def calc_loop(self,k):\n",
    "        \"\"\"use loop to calculate codelengths, fast for small k\"\"\"\n",
    "        for k in range(1,k+1):\n",
    "            k1 = round(self.p*k)\n",
    "            if k1 == 0:\n",
    "                k1 = 1\n",
    "            if k1 == k:\n",
    "                k1 = k-1\n",
    "            k0 = k-k1\n",
    "            self.els[k] = 1 + self.p*self.els[k1] + (1-self.p)*self.els[k0]\n",
    "    \n",
    "    def calc(self,k):\n",
    "        \"\"\"calculate expected length for k codewords\"\"\"\n",
    "        if k not in self.els:\n",
    "            k1 = round(self.p*k)\n",
    "            if k1 == 0:\n",
    "                k1 = 1\n",
    "            if k1 == k-1:\n",
    "                k1 = k-1\n",
    "            k0 = k-k1\n",
    "            #call smaller value first, minimizes recursion depth\n",
    "            if k0 < k1:\n",
    "                e0 = self.calc(k0)\n",
    "                e1 = self.calc(k1)\n",
    "            else:\n",
    "                e1 = self.calc(k1)\n",
    "                e0 = self.calc(k0)\n",
    "            self.els[k] = 1 + self.p*e1 + (1-self.p)*e0\n",
    "        return self.els[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Examples\n",
    "\n",
    "First, compute the example above with $p=0.3$ and $k=16$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.412718999999999"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "code = BAC(0.3)\n",
    "code.calc(16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Example with $p=0.95$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.547954674515136"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "code = BAC(p=0.95, loopfirst=2**16)\n",
    "code.calc(k=64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(109.7022288182992, 165.5529325173578, 221.41910059802296)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "code.calc(2**32), code.calc(2**48), code.calc(2**64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Notice the answer above: we did an exact calculation of the expected input phrase length for a codebook with $2^{64} = 18446744073709551616 = 1.8\\times 10^{19}$ codewords.  We couldn't do that calculation if we had to compute the codebook first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(18446744073709551616, 19.265919722494797)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**64, 64*np.log10(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Calculate the efficiency of the 64 bit code.  It's over 99%. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99083994779102502"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "code.calc(2**64)*entropy(p=0.95)/64"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 48 bit code\n",
    "\n",
    "Example is given in companion notebook on BAC Encoder and Decoder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "165.5529325173578"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "code = BAC(p=0.95, loopfirst=2**16)\n",
    "code.calc(2**48)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Plot the Input Phrase Length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAGLCAYAAAALRKIcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VNUWxuHfovdepfcmVUDEBoKKBfWiYq/YrwWxgb1c\nRbFiveJV7AUBFQERUFAQC9J7C72EXpJA6r5/nIkOMZAzyUySSb73eebJzJmZMyuhLfbZe3/mnENE\nREQkvyuS1wWIiIiI+KGmRURERKKCmhYRERGJCmpaREREJCqoaREREZGooKZFREREooKaFhEREYkK\nalpE8hEz62FmLsPtkJnFmNlIM2uVxft/D7znXR+f1dzM3jSz5WYWb2YHzWylmY0wsy7h+65yxsze\nz+RnEnz7XzbPW8nMHjezHmEuWUQipFheFyAimfoMmBi4XxpoB9wAXGhmbZ1z6zO+wcyOBboCa4D+\nZnancy4+s5Ob2QDgLeBQ4LPmAylAc+BC4EYza+OcWxrebytHbgXiMjm+OpvnqwQ8Frg/PZvnEJFc\npKZFJH+a65z7OPiAma0ChgP9gJczec8A4ABwJfAr0B8YmfFFZtYbGAEsBc50zm3J8PwQ4I4wfA/h\nNto5tzMvCzCz8s65A3lZg0hhpstDItEjvblIyviEmZXAa1ZGO+d+A+bhNTGZeQ4w4JKMDQuAcy7F\nOffy0UZZzOzWwKWZ8zJ5roiZbTKz+UHHupvZd2a2LXC5a7OZTTSzbkf7hkMVuNzjzKyFmT0TqCPR\nzBaY2dlBr+sBrA08fCzoUtO6wPMNA48fN7NLzGyOmR0EXgs6Rzsz+8rMdgW+p6Vmdr+ZFc1QU/rl\nrepm9mHg9fFm9oOZdQp6XQ0zSzKzT47wvb1hZmlm1jBMPy6RqKORFpH8qYyZVQvcLw0cCzwN7ATG\nZPL684FqwAeBx+8Dw82shXNuRfqLzKwR0AmYkcNLP5/jjfZcDYzL8FwvoA7wYuAzWwBTgG14I0Wx\nQE3gJKA98JvPz6xiZpkd3++cy9jIfQAkAy8AJYCBwNdm1tw5tw5YBtwd+B6+AsYG3pfx8tMFwJ14\nl9L+C+wPfE+dgZ8Cn/FG4Hvri9cQtgeuyKTOScBu4HGgFnA78JOZneCcW+yc225m44B+ZlbJObc3\n/Y1mVgq4HJgaqF+kcHLO6aabbvnkBvQA3BFuS4CWR3jfd3gjBxZ4XA1vROa5DK/rGzjXq2Go9Uu8\nOTGVMxz/CO8f8xqBx3cGPrNrNj/n/aP8TBxwUdBrHw8cG5/+swgc7xI4PjToWMPAsccz+cz055KB\nVpk8/wveHKB2QccMGBV4X69M6h+boabjgDRgUtCxMwKvvS3D510RON4/r3+P6qZbXt50eUgkfxoB\nnB649QUewGtEJppZg+AXmlk9vH/sPnTOOQDnzf2YAFxtZsEjqhUCX/eHocYPgJLAJUG1lAP+hfcP\n8fbA4X2Br+cHRgyy60L+/pkE337O5LXD038WAM652XijKM1C/MwJzrllwQfMrAbQHRjnnFsY9BkO\nbzQMvJ9BRsMy1DQHbwSqd+DnRuDxWv55aW8AsAv4OsT6RQoUNS0i+dMq59zUwG28c24YcB7QCO8S\nRLBr8f4s/2JmTdNvwI94lyHODnpterNSPgw1TgK2410iSnchUBb4MOjY58BU4EFgt5n9aGYPZGy+\nfPg56GcSfNueyWtjMjm2C6ga4meuzORYo8DXJZk8twxv9KTxEZ7LaClQFGgAfzU+/wM6mVkHADNr\njDcC95H752UwkUJFTYtIlHDO/Y43anFa+jHzJnlcF3j4PbAq6PZq4Pj1QadZHPjaMQz1pACfAicE\nmiTwGpg9BM1zcc4lOudOB44HhgKpwJPAcjPLbEQiHFKPcDzTSTFHkZDTQrLhPbxLT+mjLdfj1Z2t\n/WhEChI1LSLRpRiHj5L0xPuf/yvAxZncxgLnmFlNAOfcWryVRSeaWcsw1JM+8ffqwGWqHsAXzrnE\njC90zv3hnHsq0MA0BeKB/4ShhuxyWb8kU+mrjtpk8lxLvL9XMxvpyWxjwNZ4DdZf++4457YB3wJX\nmFlZvJG0351zmY3siBQqalpEooSZnY536WVO0OEBeP/oPeOcG53xhjfaUozDL+E8EPj6uZnVyuRz\niprZQDNrnVVNzrn5wEK85dZX4f2d8kHwa4JWQQXbBOwAqmT1GRGUvlIopBoCl6NmAX0DG/oBf416\nDQk8/CqTt95vQcufAsudewM/OOcyrlp6B6iMt2KpDhplEQG05Fkkv+pkZlcG7pfE+1/9TXirWR4G\nbxt6vI3mZjjndhzhPDPw5p1cDzwP4JybYmY34S3jXWFmwTviNsWbl9IEb5m1Hx/gLW9+AFjpvH1i\ngj1sZmfgrehZi3epoy/eqMQwn58BcJGZZbYj7nbn3OQQzgOAc26Xma0GLjWzNXhLseOdc9/6ePtd\neEueZ5hZ+pLnc4EzgU+dcz9k8p4GwPeBZc218ZY8HwTuy+S13+ONvlyJ11x9HtI3J1JAqWkRyZ8u\nC9zAm9i5C5iMt2R3duD4FUAp/t5j5B+cc2lm9jVwk5l1d87NChx/18xm4u1f0gtvJKYIsBHvH+NL\nnf99XD7BmxxcgcybkK/x/pHuj7c/y0G8OTc3AllmJAV56wjHf8H72WTHFXh7tTwDlMFrFLJsWpxz\nf5pZd+AJ4Da8EbAYvMbtxSO8rQ/wUuA9pfH2p7kveAVS0PnTzMuPehIYlclIjEihlL6ng4iIRICZ\nvQ9c45wLaRKwmd2P1wx2d879GonaRKKN5rSIiOQzgb11bgYWqWER+ZsuD4mI5BOBmIUT8GIZGvP3\nJUIRQU2LiEh+cipeMvdO4EnnnCbgigTRnBYRERGJCprTIiIiIlEhKi8PVatWzTVs2DCvyxAREZEw\nmDNnzk7nXPWsXheVTUvDhg35888/87oMERERCQMzW5/1q3R5SERERKKEmhYRERGJCmpaREREJCqo\naREREZGooKZFREREooKaFhEREYkKalpEREQkKqhpERERkaigpkVERESigpoWERERiQpqWkRERCQq\nqGkRERGRo0uKh6Xf5HUValpERETkKFZPhTe7wahrYNeaPC1FTYuIiIj8U9wOGHMDfHwhFC0J106A\nqk3ytKRiefrpIiIikr84B/M+hskPe5eFTh0MJw+CYiXzujI1LSIiIhKwczWMHwjrZkD97tD3Faje\nIq+r+ouaFhERkcIuJQl+eQV+fgGKlYK+w6Hj1VAkf80iUdMiIiJSmG34Db69C3Yshzb9oM+zUL5m\nXleVKTUtIiIihdHBvTD1cZgzEirWg8tHQfMz87qqo1LTIiIiUpg4B0u/hu8egPgdcMLt0GMIlCyX\n15VlSU2LiIhIYbF3I0y8F1ZOglrt4PIv4JiOeV2Vb2paRERECrq0VPhjBPzwFODgjP/A8bdCUX9t\nQEpqGn+u30O3xlUjW2cW1LSIiIgUZFsXeBNtt8yDpqfDOS9C5Qa+37548z4eGLOQ5dsOMHXQqTSq\nVjaCxR6dmhYREZGCKCkepj8Lv74BZarARe95q4PMfL39UHIqr0xdxTszYqhcpgSvXdaRhlXLRLjo\no1PTIiIiUtCsmgoT7oa9G6DTNXD6E1C6su+3z1qzkwfHLmLdrgQu6VyPB89uRcUyxSNYsD9qWkRE\nRAqKuO0waQgsHg1Vm8G1E6Hhib7fvu9gMkMnLuPz2RupX6UMn95wPN2bVotgwaFR0yIiIhLtgvOC\nkhOylRc0afFWHvlmCbviErn5lMYM7N2c0iWKRrDo0KlpERERiWY7V8H4u7OdFxS7/xCPfrOY75fE\n0rp2BUZe24Vj61SMYMHZp6ZFREQkGv2VF/Q8FC8dcl5QWprjiz838szEZSSlpPFAn5bccHIjihfN\nX3lDwdS0iIiIRJv1v3rLmHeuyFZe0Nqd8QwZu5DfYnbTrXEVhvZrl6dLmf1S0yIiIhIt/pEX9CU0\nP8P325NT03hnRgyvTF1FyWJFeLZfWy7pUg/zuQw6r6lpERERye/CkBe0aJO3SdzSrfvp06YWT5zf\nhpoVSkWw6PBT0yIiIpKf5TAv6GBSKi9PXcn/ZsRQrVxJ/ntlJ/ocWzuCBUeOmhYREZH8KC0Vfn8b\nfvwPXl7Q03D8Lb7zggB+Wb2TIWMXsWF3Apd1rcfgs1pRsXTebxKXXWpaRERE8psc5gXtS0jm6YlL\nGfXnJhpWLcOnNx5P9yb5Z5O47FLTIiIikl8kxcP0ofDrm9nKC3LOMXHRNh4bt4Q9CUnc2qMJd/Vq\nRqni+WuTuOxS0yIiIpIf5DAvaNu+QzzyzWKmLI3l2DoVeP+6/LtJXHapaREREclLwXlB1ZrDdd9B\ng+6+356W5vhs9gaenbic5LQ0Hjy7Jdef2Ihi+XiTuOxS0yIiIpIXnIN5H8HkR7y8oB5D4KS7Q8oL\nWrMjjiFjFvHHut10b1KVof3a0qBq/t8kLrvUtIiIiOS2navg24Gwfma28oKSU9MY8XMMw39YRali\nRRh2YTsu7lw3ajaJyy41LSIiIrklJRFmvgIzXgjkBb0KHa/ynRcEsGDjXh4Ys5Dl2w5wTtvaPHZe\na2qUj65N4rIrV5sWM1sHHABSgRTnXGczqwJ8ATQE1gH9nXN7crMuERGRiMthXtCh5FRenLyCd2eu\npXr5koy46jjOaFMrggXnP3kx0tLTObcz6PFg4Afn3LNmNjjw+IE8qEtERCT8Du6FqY/BnPehYv2Q\n84IAFm/ex91fzGfV9jgu61qfIWe3pEKp6N0kLrvyw+Wh84EegfsfANNR0yIiItEuDHlBKalpvDV9\nDcN/WEWVsiX44PqunNq8egSLzt9yu2lxwFQzSwXeds6NAGo657YGnt8G+B8rExERyY+C84Jqtw85\nLwhg7c54Bo2az7wNezm3XW3+c8GxVCpTIkIFR4fcblpOcs5tNrMawBQzWx78pHPOmZnL7I1mdhNw\nE0D9+vUjX6mIiEioMuYFnfkMdL05pLwg5xwf/76BZyYso3hR49XLOnJe+2MiV3MUydWmxTm3OfB1\nu5l9BXQFYs2stnNuq5nVBrYf4b0jgBEAnTt3zrSxERERyTNbF8C4O2Hr/GzlBQHE7j/EfaMX8vPK\nHZzcrBrPX9SeWhULx8ogP3KtaTGzskAR59yBwP0zgCeBccA1wLOBr9/kVk0iIiI5dlheUNWQ84LS\nfbtgCw9/vZjElFSePL8NV3VrUOD3XQlVbo601AS+CvwCFAM+dc5NMrPZwCgzGwCsB/rnYk0iIiLZ\nl8O8IPASmR/5ZjHjFmyhfb1KvNy/PY2r+5+sW5jkWtPinIsB2mdyfBfQK7fqEBERybEc5gWlm7Fq\nB/d9uZCdcYkMOr05t/VoUiAzg8IlPyx5FhERiQ5hyAsCOJiUytDvlvHhr+tpUr0s71x9Im3rFqxE\n5khQ0yIiIuLHP/KChkP15iGfZt6GPdwzagExO+O5/sRG3N+nBaWKF41AwQWPmhYREZGjCUNeEHgh\nh6/9sIo3pq+hZvmSfHrD8XRvWi1CRRdMalpERESOJDgv6NgL4cyhIeUFpVu9/QB3f7GARZv30a9T\nHR7r24aKpQvfNvw5paZFREQko4x5QVeMhmanh3yatDTH+7PW8dyk5ZQpUZS3rujEWW1rh7/eQkJN\ni4iISLrM8oJ6PgglyoZ8qs17D3LflwuYtWYXp7WswbMXtqVGeW0UlxNqWkRERMDba2XCvbDqe6jd\nAS4fBcd0CPk0zjm+mreZx75ZQqpzDO3Xlku71NNGcWEQctNiZpWAw2YfOed2h60iERGR3BSGvKB0\nu+OTeHDsIiYt2UbnBpV5qX8H6lctE/6aCylfvyJm1gD4L9ADCI6YNLzkZq3VEhGR6BOcF9TsDC8v\nqFL2Qnl/XB7L/aMXse9gEg/0aclNpzSmaBGNroST3zZyJFAJGABswWtUREREolNSPEx7Bn57K5AX\nNBLa/CvkvCCA+MQU/jNhKZ/9sZGWtcrz4fVdaX1MhQgULX6blq5AN+fc4kgWIyIiEnGrpsD4QbBv\nAxx3LfR+POS8oHSz1+3mnlEL2LgngZtPbcyg05tTspguPkSK36ZlLRDaHsUiIiL5Sdx2mDQYFo+B\nai2ynRcEkJiSystTVvH2z2uoU6k0X9x0Al0bVQlzwZKR36blLmComd3mnFsdyYJERETC6q+8oIch\n+SD0eBBOGhhyXlC6JVv2ce+XC1m2dT+XdK7HI31bU66kFuPmhiP+lM3sAIfPXSkFrDCzRCAl+LXO\nOV28ExGR/Cc4L6jBiXDuK9nKCwJISknj9WmreXPaaiqVKc47V3fm9Nah744r2Xe01vD2XKtCREQk\nnDLmBZ33GnS4MuS8oHSLNu3jvtELWL7tABd0OIbH+rahctkSWb9RwuqITYtz7oPcLERERCQs1s8K\n5AWthGMvgj5DoVyNbJ0qMSWVV39YxX9/iqFq2RL87+rO9NboSp7xu09LKlDbObc9w/GqwHbnnKZK\ni4hI3grOC6pUH64YA816Z/t08zfu5b4vF7BqexwXHVeXR85pTcUyCjnMS35nDh1p4XpJIClMtYiI\niITOOVjylbcyKH4HdL8DegzJVl4QwKHkVF6espJ3ZsRQs0IpRl7XhZ4tsjdSI+F11KbFzAYF7jrg\nFjOLC3q6KHAysDxCtYmIiBxdmPKC0s1Zv5v7Ri8kZkc8l3apx4PntKJCKY2u5BdZjbTcEfhqwA1A\natBzScA64JbwlyUiInIUh+UFAWcOha43ZSsvCOBgUiovTF7Be7+s5ZiKpfloQFdOblY9jAVLOBz1\nV9c51wjAzKYB/Zxze3KlKhERkSM5LC/oTDjnhWznBQH8HrOLB8YsZN2uBK7sVp/BZ7XSviv5lK9f\nFedcz0gXIiIiclQZ84Iufh9aX5CtvCCAhKQUhk1awfuz1lGvSmk+vfF4ujepFt6aJaz8rh567whP\nOeAQsBr4wjm3JVyFiYiI/CWMeUEAs9bs5IExC9m4+yDXdm/IfWe2oKxGV/I9v79C1fEm3aYB6aGJ\nx+LNdZkD9AOeNLOTnXPzw16liIgUTv/IC5oEDU7I/ukSUxg6cRmf/L6BhlXLMOpmZQZFE79Ny09A\nHDDAOZcAYGZlgHeAhcDZwIfAi0CvCNQpIiKFSVqalxc05REvL6jnQ3DiXdnOCwKYsWoHg8csYsu+\ng9xwUiPuOaMFpUtom7Fo4rdpuRvond6wADjnEszsaeAH59xzZvYcMDUSRYqISCGyYyWMHwjrf/Hy\ngvoOh2rNsn26/YeSeWbCMj6fvZHG1csy+pYTOK6BRleikd+mpQJQG1iW4XgtoFzg/v4QziciInK4\nw/KCyuQ4Lwhg2ortPDh2EbH7D3HzqY25u3dzShXX6Eq08ttkfAW8a2b3A7MDx7oAw4CxgcddgZXh\nLU9ERAqFMOYFAexLSOapCUsZPWcTzWqU463bTqRDvUphLFjygt+m5RbgJeDjoPekAO8B9wYeLwNu\nDGt1IiJSsB3cA1Meg7kfhCUvCGDq0lge/GoRu+KT+HfPJtzZqxkli2l0pSDwu09LAt42/vcATQKH\n1zjn4oNeo1VDIiLij3OwZCx8NxgSduU4LwhgT3wST3y7hK/nb6FlrfK8e00X2tatGMaiJa+FNAcl\n0KQsjFAtIiJSGOzdABPugVWTvbygK0dD7fY5OuWkxdt4+OvF7E1I4q5ezfh3z6aUKJb9uTCSP/nd\nXK4UcBfecuYawGG/E5xz7cJfmoiIFCipKfD7f2Ha04DlOC8IYFdcIo+NW8L4hVtpXbsCH1zfhTbH\naHSloPL7O+VN4F/Al8AsvJ1wRURE/NkyH76908sNanYmnPMiVKqX7dM55xi3YAtPfruU/YeSuef0\n5tzSownFi2p0pSDz27RcAFzsnNM+LCIi4l9iHEwfCr+9CWWq5TgvCGDx5n08+e1S/li3m3Z1K/Lp\nRd1oUat8+GqWfMtv05IAbIxkISIiUsCsnOzNXdm3AY67LpAXlP1lx7viEnlh8ko+n72BymVKMLRf\nW/p3rkfRItlvgCS6+G1ahgGDzOwW55wuDYmIyJGFOS8oOTWNj35dzytTV5KQlMp13RtxV+9mVCxd\nPIxFSzTw27ScjheY2MfMlgLJwU86584Ld2EiIhJlIpQX9OS3S1m1PY6Tm1Xjsb6taVpDl4IKK79N\ny068XXFFRET+acdKb0fbDbOgwUnQ95Uc5QWt3xXPfyYsY8rSWBpULcM7V3emd6saWA7mwkj087u5\n3HWRLkRERKJQSiLMfBlmvBjIC3odOl6Z7Ym28YkpvDFtNf+bsZZiRY37+7RgwEmNtKOtACFuLmdm\nnfF2xB3vnIs3s7JAonMuJSLViYhI/rXuFy+NOQx5Qc45vp6/mWe/W07s/kT6dazDA2e1pGaFUmEu\nWqKZ383lagLf4IUiOqAZEIOXR3QIb+M5EREpDMKcF7Rw014eH7eEuRv20r5uRd668jg61a8cxoKl\noPA70vIyEAtUBTYEHf8SeC3cRYmISD4U5rygHQcSef775Xw5ZxNVy5Zk2EXtuKhTXYpoCbMcgd+m\npRfQyzm3J8MkqDVA/bBXJSIi+UsY84KSUtJ4f9ZaXv1hNYkpqdx0cmNuP60p5UtpCbMcnd+mpTSQ\nlMnx6niXh0REpCDKmBfU51kvL6hI9ibGTlu+nafGLyVmZzyntazBw+e0onH1cuGtWQosv03Lz8C1\nwIOBx87MigIPAD9EoC4REclrW+Z5y5jDkBcUsyOOp8YvZdqKHTSuVpaR13WhZ4vsTdqVwstv03I/\n8JOZdQFKAi8CbYCKwIkRqk1ERPJCcF5Q2eo5ygs6cCiZ135czchf1lKyWFEeOrsV13RvSIliCjaU\n0Pndp2WpmbUFbgUSgVJ4k3DfcM5tDeUDAyM0fwKbnXPnmlkV4AugIbAO6O+c2xPKOUVEJEzClBeU\nluYYM3cTz01awc64RC4+ri739WlBjfJawizZ53ufFufcNuCx4GNmVtXMbnPOvRnCZ94FLAMqBB4P\nBn5wzj1rZoMDjx8I4XwiIpJTB2K9vKAlY3OcFzRvwx4e/3YpCzbupWP9Srx7TWfa18t+UKJIupA2\nl8tEXbwlz76aFjOrC5wDPA0MChw+H+gRuP8BMB01LSIiuSMtDeZ9CFMezXFeUOz+Qzz33XLGzttM\njfIlefmS9pzfvo6WMEvY5LRpCdUrePNjgtOuagZdYtoG1MzlmkRECqcdK+DbgTnOC0pMSeXdmWt5\n48fVJKc6buvRhH/3bErZkrn9T4wUdLn2O8rMzgW2O+fmmFmPzF7jnHNm5o7w/puAmwDq19fWMCIi\n2ZaSCDNegpkv5SgvyDnHd4u38dyk5azflcDprWvy8DmtaFA1e5vNiWQlN9vgE4HzzOxsvIm8Fczs\nYyDWzGo757aaWW1ge2Zvds6NAEYAdO7cOdPGRkREshCcF9T2YjhzKJSrHvJpZq/bzTMTlzFvw16a\n1yzHRwO6cnKz0M8jEoqjNi1mNuhozwPH+P0g59wQYEjgvD2Ae51zV5rZ88A1wLOBr9/4PaeIiPh0\ncI83b2XuhznKC1q9PY7nJi1nytJYalbwtt6/sFNdimreiuSCrEZa7vBxjg1Zv+SongVGmdkAYD3Q\nP4fnExGRdGHKC9p+4BDDp67i89kbKV28KPed2YLrT2xE6RLZ2xlXJDuO2rQ45xpF4kOdc9PxVgnh\nnNuFl20kIiLhtGe9t+fK6inZzguKT0zhnRkxjPg5hqSUNK7q1oA7TmtK1XKhry4SySlN7RYRKWjC\nkBeUkprGF39u5OUpq9gZl8g5bWtz35ktaFhNk2wl76hpEREpSILzgpr3gbNfCCkvyDnHlKWxPDtp\nOTE74unSsDIjrj6OTvUrR7BoEX/UtIiIFAT/yAv6AFqfH9Iy5rkb9jB04jJmr9tDk+pleefqzvRu\nVQPLRuaQSCSoaRERiXYrJ8OEQbBvY7bygtbujOf575czcdE2qpUryTP/akv/znUpVlShhpK/qGkR\nEYlWB2Jh0gOw5Kts5QXtjEvktR9W8cnvGyhRrAh3927ODSc30k62km/5/p1pZqWAc4EmwNvOub1m\n1gTY45zbHakCRUQkgxzmBR1MSuXdmTH896cYDianclnXetzZq5kSmCXf89W0mFlTYApeZlAl4Etg\nL3Br4PENkSpQRESC5CAvKDXNMWbOJl6csoLY/Ymc0bom9/dpSdMa5SJctEh4+B1peQWvabkVr1lJ\nNw4YGe6iREQkg/S8oBkvehvDhZAX5Jxj+oodDP1uGStj4+hYvxKvX96JLg2r5ELhIuHjt2npDnRz\nzqVmmEW+gRC28hcRkWxY94u3jHnXqpDzghZu2svQicv5NWYXDauW4a0rOtHn2FpaESRRKZTZVsUz\nOVYf2BemWkREJFjGvKArx0BTf3lBG3cn8Pz3Kxi3YAtVy5bgyfPbcFnX+hTXiiCJYn6blsnAIGBA\n4LEzswrAE8CESBQmIlJoOQeLx8CkIYG8oDuhx2BfeUF74pN4fdpqPvx1HUWLGHec1pSbTmlM+VKZ\n/b9TJLr4bVoGAdPMbAVQCvgCaArEooBDEZHwCc4LOqajN7pSu12WbzuUnMr7s9bxxrTVxCem0L9z\nPe4+vTk1K2hFkBQcvpoW59wWM+sAXAZ0AooAI4BPnHMHI1ifiEjhkJoCv78F054h1LygGat28NBX\ni9mwO4FeLWvwwFktaV6zfORrFsllvue0BJqT9wI3EREJly3zYNydsG1hSHlBO+MS+c/4pXw9fwuN\nqpXlkxuO58Sm1XKhYJG84Xeflv7AXufc5MDjR4GbgCXAtc65rZErUUSkgEqM80ZWfn8LytbwnRfk\nnGPUnxt5ZuJyEpJSuLNXM27r0YRSxf2nOItEI78jLY8DAwHMrBPwIPAo0Ad4Ebg8EsWJiBRYK7/3\n5q7s2widr4dej/nKC1q9PY4Hv1rEH2t307VRFZ75V1ttDieFht+mpQGwInD/X8DXzrlhZjYZ+D4i\nlYmIFETBeUHVW8L130P9blm+7VByKm9OX8Nb01dTpkQxhl3YjouOq0uRItpvRQoPv03LIbwt/AF6\n8fe8ln0U84TUAAAgAElEQVRBx0VE5EgOyws6BD0fDuQFlcjyrbPW7OThrxYTszOeCzocw8PntqZa\nOX85QyIFid+mZQbwopnNBDoDFwWONwc2RqIwEZECY8cKb0fbDb9Cw5Ph3FegWtMs37Y7PomnJyxj\nzNxNNKhaho8GdOXkZv52whUpiPw2LbcDb+E1K7c457YEjp+FLg+JiGQuOC+oZDk4/w3ocIWvibZj\n527mPxOWcuBQCv/u2YQ7TmumibZS6Pndp2UT0DeT4wPDXpGISEFwWF5QfzjzGV95QWt3xvPQV4uY\ntWYXxzWozDP/akuLWroKLwKhZQ+JiEhWDssLauA7LygpJY23f1rDa9NWU7JYEZ7+17Fc1qW+JtqK\nBPG7T0sJ4CG8HXHrkyE80TmnMUsRKdz+ygsaDAm7vUm2pw6GEmWyfOsfa3fz4FeLWL09jnPb1ebR\nc1tTQ9vvi/yD35GWp4BLgKHAy8B9QEPgUuCRiFQmIhIt/pEXNNZXXtDehCSe/W45n8/eSJ1KpRl5\nXRd6tqiRCwWLRCe/TUt/vAm4k8zsBeAb59waM1sGnA68HbEKRUTyq+C8ICsCfZ6DrjdmmRfknGPc\ngi08NX4pexKSufmUxtzVuxllSuiKvcjR+P0TUhNYGrgfB6Rv2zgJeC7cRYmI5HuH5QWdBee8ABXr\nZvm29bviefjrxcxYtZP29Srx4fVtaX1MhVwoWCT6+W1aNgDHBL6uBs4E5gAnAEp5FpHCI2NeUP8P\nodV5WS5jTk5N450ZMQyfuoriRYvwxHltuLJbA4pqoq2Ib36blq/wdsL9DRgOfGZmNwJ1gOcjVJuI\nSP5yWF7QAOj9GJSqmOXb5qzfw4NjF7Ei9gB92tTi8fPaUKuiJtqKhMrvPi1Dgu6PNrONwInASufc\n+EgVJyKSL2QzL2jfwWSe/345n/y+gVoVSvHO1Z05vXXNXChYpGDKsmkxs+LAx8CDzrk1AM6534Hf\nI1ybiEjeSkuDuR/AlMcg5RCc9jB0zzovyDnHhEVbeeLbpeyKS+S67o0YdEZzypXURFuRnMjyT5Bz\nLtnMzgCGZPVaEZECI5t5QRt3J/DoN4uZtmIHx9apwHvXdKFt3awvIYlI1vy2/WOBfsALEaxFRCTv\nJR+CmS95mUEly8H5b0KHy7OcaAvw1bxNPPTVYgAeObc115zQgGJFi0S6YpFCI5TVQw+b2cnAn0B8\n8JPOuZfCXZiISK5bNxO+HRhyXtDBpFQeH7eEL/7cSNeGVXj50g7UqVQ6FwoWKVz8Ni3XAnuAdoFb\nMAeoaRGR6JWw28sLmvdRIC9oLDTt5eutq7cf4N+fzGPl9gPc3rMpA3s30+iKSIT4XT3UKNKFiIjk\nuhzkBQGMnrOJR75eTJkSRfnw+q6c3CzrURkRyT5NZReRwmnPukBe0FQ4ppPvvCCAhKQUHv1mCaPn\nbKJb4yoMv7QjNRVwKBJxvpsWM7sEb4O5GsBhY5/OufPCXJeISGSkpsBvb8L0oSHlBaVbGXuAf38y\nl9U74rizVzPu6tVMu9qK5BJfTYuZPQ8MBKYBW/DmsYiIRJfNc71lzCHmBYG398qXf27i0XGLKVey\nOB8POJ4Tm1aLcMEiEszvSMvVwGXOudGRLEZEJCIS42Da0/D7f0PKC0oXn5jCI18vZuy8zXRvUpVX\nLu1AjfK6HCSS2/w2LUWA+ZEsREQkIlZMgon3hpwXlG75tv38+5O5rN0Zz929m3P7aU11OUgkj/ht\nWkYAVwKPR64UEZEwOrANvnsAln4dUl5QOuccn8/eyOPjllChdHE+uaEbJzSpGsGCRSQrR2xazOzV\noIdFgCvM7HRgIZAc/Frn3J2RKU9EJETZzAsKFpeYwkNfLeKb+Vs4uVk1XurfgerlS0awaBHx42gj\nLW0zPE6/PNQyw3FNyhWR/GH7cm+i7cbfQsoLCrZ0y35u/3Qu63bFc+8ZzbmtR1OK6HKQSL5wxKbF\nOdczNwsREcm25EMw40WY+XLIeUHpnHN88vsGnhy/lMplivPZjd04vrEuB4nkJ1nOaTGzBsDpQAlg\nunNuacSrEhHxa+0MGD8Qdq2Gdpd4eUFlQ1uKfOBQMkPGLmL8wq2c0rw6L/dvT9Vyuhwkkt8ctWkx\ns1OAiUD6ntYpZnaNc+6ziFcmInI0CbthyiMw7+OQ84KCLd68j9s/ncvGPQe5v08LbjmliS4HieRT\nWY20PAX8CNwMJALPAMOAkJsWMysF/AyUDHzuaOfcY2ZWBfgCaAisA/o75/aEen4RKSScg0Wj4fsh\ngbyggXDqA77zgv4+jeOj39bzn/HLqFquBJ/f1I0uDatEqGgRCYesmpa2wCnOua0AZnYPcKOZVc5G\nY5EInOacizOz4sBMM/sO6Af84Jx71swGA4OBB0I8t4gUBnvWwfhBsOYHLy/oqq+gVsY1A1nbfyiZ\nwWMWMnHRNk5rWYMXLm5PlbL+VxeJSN7IqmmpBGxPf+CcizezhMDxkJoW55wD4gIPiwduDjgf6BE4\n/gEwHTUtIhIsNQV+ewOmDfUygkLMCwq2cNNebv90Hpv3HuTBs1tyw0mNdTlIJEr42VyunZntDnps\nwLFmVjn9gHNurp8PM7OiwBygKfCGc+53M6uZPpIDbANq+itdRAqFzXMCeUGLoMXZcPbzvvOCgjnn\neH/WOp6ZuIzq5Uoy6uYTOK5B5azfKCL5hp+m5Xu8RiXYN0H3HeDrvzvOuVSgg5lVAr4ys2MzPO/M\nLNN9X8zsJuAmgPr16/v5OBGJZv/IC/oIWvUNaRlzun0Hk7l/9AK+XxJL71be5aBKZXQ5SCTaZNW0\nNIrEhzrn9prZNKAPEGtmtZ1zW82sNkGXozK8ZwRenACdO3fWhnYiBdmK72DCvbB/M3QZAL0eDSkv\nKNj8jXu5/dO5bNt3iIfPacWAkxph2Wh8RCTvHbVpcc6tD9cHmVl1IDnQsJTG2/vlOWAccA3wbODr\nN0c+i4gUaAe2wXf3w9JvoHqrQF7Q8dk6lXOOd2eu5blJy6lRvhRf3nICHevrcpBINPMbmBgOtYEP\nAvNaigCjnHPjzexXYJSZDQDWA/1zsSYRyQ/S0mDOSJj6RCAv6BHofmdIeUHB9iYkce+XC5m6LJYz\nWtfk+YvaU7FM8TAXLSK5LdeaFufcQqBjJsd3AaHvCCUiBcP2ZYG8oN+9vKC+w6Fqk2yf7veYXQz8\nYj474xJ5rG9rru3eUJeDRAqI3BxpERH5W/IhmPECzHwl23lBwVJS03jtx9W89uMq6lcpw9hbT6Rt\n3ezNgxGR/ElNi4jkvrUzvNGV3WuynRcUbNOeBAZ+Pp8/1+/hwk51eeL8NpQrqb/eRAoaX3+qzew9\n4C7n3IEMx8sCrznnro9EcSJSwATnBVVu6O1o2+S0HJ1y4qKtDB6zkDQHr1zSgQs61glPrSKS7/j9\nr8g1eNvrH8hwvDRwNaCmRUSOLD0vaNJgOLgHTrobTrk/5LygYAeTUnly/BI++2Mj7etV4tVLO9Cg\natkwFi0i+U1WKc9V8DaWM6CymaUEPV0UOAeIjVx5IhL1dq+FCYNgzY9Q5zi4+uts5QUFW7Z1P3d8\nNo81O+K4tUcTBp3enOJFi4SpYBHJr7IaadmJt+OtA5Zm8rwDHgt3USJSAKQmw29v/p0XdNbz3kZx\n2cgLSuec48Nf1/P0xGVULF2cj64/npOaZX8ujIhEl6yalp54oyw/AhcCwRlEScB659yWCNUmItFq\n8xwYdxfE5iwvKNju+CTuH72Aqcu207NFdV64uD1Vy5UMU8EiEg2y2hH3JwAzawRsCCQ1i4hkLvEA\n/Pgf+GMElKsJl3zs5QXl0KzVO7l71Hz2xCfz6Lmtue5E7b0iUhj5nYhbFah6pL8k/KY8i0gBtuI7\nmHAP7N+S47ygdMmpabwydSVvTl9Do2plee/aLrQ5RnuviBRWfpuWP/HmrwR3LcGjLtm/SC0i0S04\nL6hGa7j4fajXNcen3bg7gTs+m8f8jXu5pHM9HjuvNWVKaO8VkcLM798AGdOei+Ntyf8QMCSsFYlI\ndPgrL+hxSEnMcV5QsHELtvDQ2EVg8PrlHTm33TE5r1dEop6vpuUIac+rzWwf3uqh78JalYjkb8F5\nQY1OgXNfyVFeULr4xBQeH7eEL+dsolP9Sgy/tCP1qmR/LxcRKVhyOta6FugQjkJEJAoclhdUHi54\nC9pflu28oGCLN+/jzs/msXZXPHec1pS7ejWjmPZeEZEgfrfxr5LxEFAbeBxYEeaaRCQ/WvszfDsw\nkBd0KZz5dI7ygtKlpTne+2Utz01aTtWyJfn0hm6c0KRqGAoWkYLG70hL+iZzwQzYCFwS1opEJH9J\n2A2TH4H54csLSrczLpF7v1zA9BU7OL11TYZd2I7KZXM+J0ZECia/TUvPDI/TgB3AaudcSiavF5Fo\n5xws+hImDQlbXlCwn1fuYNCoBew/lMxT57fhym4NtPeKiByV34m4P0W6EBHJR/6RF/QN1Do2LKdO\nSknjxckrePvnGJrVKMfHN3SlZa0KYTm3iBRsvifimllt4FagdeDQMuAtbeMvUoCkJsOvb8D0Z8OW\nFxRs3c547vx8Hgs37eOK4+vz8DmtKV1C2zyJiD9+J+KeDnyDN4fl98Dhi4F7zOwC59zkCNUnIrnl\nsLygcwJ5QXXCdvqxczfxyNeLKVrE+O+VnehzbO2wnVtECge/Iy2vAv8D7grOHzKz4cBwoFUEahOR\n3JCeF/T721C+VtjygtIdOJTMo98s4at5m+nasAovX9qBOpVKh+38IlJ4+G1aGgKvZxKY+AZwY1gr\nEpHcs3wiTLw3kBd0A/R6JMd5QcHmb9zLXZ/PY+PuBAb2bsbtPZtq7xURybZQsofaAiszHG8LzAtr\nRSISefu3enlBy8aFNS8oXez+Q7w4eQVfztlE7Qql+OLmE+jSMON2TyIiofHbtLwJvGxmzYDfAse6\n4U3MHWxmndJfqMRnkXwsY15Qr0fhhDvCkhcEkJCUwjs/r+W/P60hJS2NASc24o5ezahYunhYzi8i\nhZvfpuWTwNdnjvIceBvQaSmASH4Uobwg8Ha1HTtvM89/v5zY/YmcdWwtBp/VkgZVy4bl/CIikP2U\nZxGJFv/IC/ovtL80LHlBALPW7OTpCctYsmU/7etW5PXLO+lSkIhERE5SnkUkv4tQXhDAmh1xDJ24\nnKnLYqlTqTTDL+1A33bHUKSIdrUVkcgIZXO5usApQA3gsOn/zrmXwlyXiOREBPOC9sQnMfyHVXz8\n23pKFS/KfWe2YMBJjShVXFeGRSSy/G4udwXwHpCClzkUvPTZAWpaRPKDv/KCBsOhfWHNC0pMSeXD\nWet57cdVxCWmcGnX+tzduznVy5cMQ+EiIlnzO9LyJPAi8IhzLjWC9YhIdh2WF9QZ+g4PS16Qc47v\nFm/j2e+Ws2F3Aqc2r85D57Siec3yYShaRMQ/v01LTeB/alhE8qHD8oKKhTUvaN6GPTw9YRl/rt9D\ni5rl+eD6rpzavHoYihYRCZ3fpmUicDwQE8FaRCRUwXlBLc+Fs4aFJS9o054Ehk1awbgFW6hWrgRD\n+7Xl4uPqajdbEclTfpuWKcBzZtYGWAQkBz/pnBsb7sJE5CgilBd04FAyb05fw7sz12LA7T2bckuP\nJpQr6XvOvohIxPj9m+jtwNcHM3lOG8qJ5KYI5AWlpKbx+eyNvDxlJbvik/hXxzrcd2YLjlGwoYjk\nI373adGYsEhe+0de0AdQr0uOTumcY/rKHTwzYRmrtsfRtWEVRl7XinZ1K4WpaBGR8NGYr0h+l5YG\nc96DqU9AapKXF9T9Tiiaszyf5dv28/SEZcxYtZOGVcvw3yuP48w2NbEw7ZQrIhJuR21azKyfn5No\nTotIhByWF3QqnPtyjvOCth84xEuTVzLqz42UL1WcR85tzVXdGlCimAZURSR/y2qkZbSPc2hOi0i4\nRSAv6GBSKv+bEcNbP60hKSWNa7s34s5eTalUJjwJzyIikXbUpkVzWUTyQHBeUPvL4IynoWzVbJ8u\nNc3xzfzNPP/9CrbuO8SZbWoy+KxWNKqmBGYRiS6a0yKSXyTshskPw/xPoHIjuOpraNIz26c7lJzK\n6DmbeHfmWtbujKdtnYq8ckkHjm+c/QZIRCQvqWkRyWv/yAsaBKfeD8Wzt9x4V1wiH/66no9+W8/u\n+CTa1a3I65d35OxjayuBWUSimpoWkbyUMS/ovFehZptsnSpmRxzvzlzL6DmbSExJo3erGtx4cmO6\nNqqiFUEiUiCoaRHJCxnzgs5+ATpfn628oDnrd/P2TzFMWRZL8SJF6NepDjec3IimNRRoKCIFi5oW\nkdyWMS/o7OehwjEhnSI1zTFl6TZG/BzD3A17qVi6OP/u0ZSruzegRvlSESpcRCRvqWkRyS2H5QXV\nhks+gVbnhnSKg0mpjJ67iXdnxLBuVwL1qpTmifPacHHnupQpoT/OIlKwHfFvOTNbi7cHS5acc43D\nVpFIQRScF9T1RjjtEShVwffbd6ZPrv11HXsSkmlfrxJv9mnJmW1qUVSTa0WkkDjaf81eD7pfDhgE\n/AH8Gjh2AtAVeDEypYkUAIflBbUJOS9ozY44/jdjLWPmbiIpJY3erWpy0ymN6dKwsibXikihc8Sm\nxTn3VzNiZu8Dzznnngl+jZkNAXwtdTCzesCHQE28EZwRzrnhZlYF+AJoCKwD+jvn9oT0XYjkN//I\nC3oMut/hKy/IOcef6/cw4ucYpi6LpXjRIlzYqQ4DTmpM0xrlcqF4EZH8ye9F8H5Ap0yOfwkM8XmO\nFOAe59xcMysPzDGzKcC1wA/OuWfNbDAwGHjA5zlF8p/YpTB+oJcX1LiHlxdUJesrqKlpjslLtvH2\nzzHM37iXSmWKc0fPplx1QkOqly8Z8bJFRPI7v01LPNADWJ3heA8gwc8JnHNbga2B+wfMbBlQBzg/\ncB6AD4DpqGmRaJR8CH5+Hn55BUpWgH+9De0uyTIvKCEphdFzNvG/GWvZsDuB+lXK8OT5bbjoOE2u\nFREJ5vdvxJeBN8ysM/Bb4Fg34Brg8VA/1MwaAh2B34GagYYGYBve5SOR6HJYXtDlcMZ/sswL2nEg\nkY9+XceHv61nb0IyHepVYshZLTlDk2tFRDLlq2lxzg0zs3XAXUD/wOFlwDXOuVGhfKCZlQPGAAOd\nc/uDJxM655yZZbpiycxuAm4CqF+/figfKRI5GfOCrv7GuyR0FKu3x/HuzBjGzN1Mcurfk2s7N9Dk\nWhGRo/E99hxoTkJqUDIys+J4DcsnzrmxgcOxZlbbObfVzGoD24/w+SOAEQCdO3f2tRRbJGKcg4Wj\n4PshXl7QyffAKfcdNS9ozvo9vDV9NVOXbadEsSJc2KkuN5zciCbVNblWRMQP302LmZUCzgWaAG87\n5/aaWRNgj3Nut4/3G/AusMw591LQU+PwLjM9G/j6TQj1i+S+3TEwfhDETIO6XaDv8KPmBa2MPcCw\nSSuYuiyWymWKc2evZlx9QgOqldPkWhGRUPhqWsysKTAVb7+WSnirhvYCtwYe3+DjNCcCVwGLzGx+\n4NiDeM3KKDMbAKzn78tPIvlLajL8+nogL6h4lnlBW/Ye5OUpKxkzdxNlSxTj3jOac/1JjTS5VkQk\nm/z+7fkKMBmvSdkbdHwcMNLPCZxzM4EjXbDv5bMOkbyxaQ58eyfELs4yL2hvQhJvTl/D+7PWgYPr\nT2zEbT2bUqVsidytWUSkgPHbtHQHujnnUjNMFNwAhJb0JhJNEg/AD0/BHyOyzAs6mJTKyFlreWv6\nGuISU+jXsS53n96MupXL5HLRIiIFUyjj1Jlt5Vkf2BemWkTyl+UTYOJ9WeYFpaSm8eWcTbwydSWx\n+xPp1bIG9/VpQcta/rOFREQka36blsl42UMDAo+dmVUAngAmRKIwkTyzf0sgL+jbo+YFOef4fsk2\nhn2/gpgd8XSqX4nXLutE10ZV8qBoEZGCz2/TMgiYZmYrgFJ4WUFNgVg0cVYKirQ0+PNd+OHJLPOC\nfovZxbPfLWf+xr00rVGOt686jjNa19Q+KyIiEeR3c7ktZtYBuAwvg6gI3p4pnzjnDkawPpHcEbsU\nvr0LNv1x1LygpVv2M+z75UxfsYNaFUrx3IVtubBTXYoVLZLrJYuIFDZ+lzyfAsxyzr0HvBd0vJiZ\nneKc+zlSBYpEVPLBQF7Q8KPmBW3cncBLU1by9fzNlC9ZjMFnteTa7g0pVTzz5c4iIhJ+fi8PTQMy\n2622YuA5/c0t0SfmJy+NeXfMEfOCdsUl8vq01Xzy2wbM4KZTGnPbqU2pWCazeekiIhJJfpsWAzLb\nOr8qXgK0SPSI3+XlBS341LsElEleUHxiCu/OXMuIn2NISEqhf+d63NW7GbUrHnmbfhERiayjNi1m\nNi5w1wEfm1li0NNFgWOBWRGqTSS8nIOFX8D3Dx4xLyg5NY3P/9jA8B9WszMukTPb1OS+M1vQtEb5\nPCxcREQg65GWXYGvBuwBgifdJgEzgXciUJdIeO2OgfF3Q8x0qNMZznv1sLygtDTHhEVbeXHyCtbt\nSqBrwyq8fdVxHNegct7VLCIihzlq0+Kcuw7AzNYBzzvnEnKjKJGw8ZEXNHPVTp6btJxFm/fRslZ5\n3ru2Mz1b1NDyZRGRfMbvnJbRePuyLAw+aGbtgBTn3NJwFyaSY5v+9JYxxy6GVn3hrGGH5QUt3ryP\n5yYtZ8aqndSpVJoXL27PBR3rULSImhURkfzIb9MyAniVDE0L0Bq4HTgpnEWJ5Mih/fDjU/DHO5nm\nBa3fFc8Lk1fy7YItVC5TnIfPacWV3Rpo+bKISD7nt2lpB8zO5PhsoG34yhHJoeUTYMK9cGBrpnlB\no2Zv5KGvF1GsSBFu79mUm05tTIVSWr4sIhIN/DYtqUBmMxIr403SFclb+7d44YbLx3t5QZd8BHU7\n//W0c46Xp6zk1R9Xc3Kzarx4cXtqVCiVhwWLiEio/DYtPwEPmdnFzrlU8HbDBR4CtBuu5J30vKCp\nT0BacqZ5QUkpaQweu5Cxczdz8XF1eaZfW4pr230Rkajjt2m5H29582ozmxk4dhJQDjglEoWJZCl2\nSSAvaPYR84L2H0rmto/nMnP1Tgb2bsZdvZppVZCISJTyG5i4IrBS6HagQ+DwJ8CbzrktkSpOJFPB\neUGlKh4xL2jrvoNcN3I2q7fH8fxF7bi4c708KlhERMLB70gLzrmteJeDRPJOzHRvk7ij5AWBl8Z8\n/fuziUtMYeR1XTi5WfXcr1VERMLKd9NiZm2Bm4HGwADn3FYzuwBY75ybF6kCRQBfeUHpZqzawa0f\nz6VcyWJ8ecsJtKpdIdPXiYhIdPHVtJjZGcA44DugF5Ae1tIEuBa4IBLFifyVFzRpCCTuzzQvKNiX\nf25kyNhFNK1RjpHXdVHAoYhIAeJ3pOUpYJBz7k0zOxB0fDpwT9irEoHD84LqdoG+r0LN1pm+1DnH\n8B9W8crUVZzYtCpvXXmc9l8RESlg/DYtxwITMzm+G6gSvnJE8PKCZr0GPz0XlBc0AIpkvkw5OTWN\nB8cu4ss5m7iwU12G9mtLiWJa0iwiUtD4bVp2A3WAdRmOdwI2hbMgKeQ2/Qnj7oTtSzLNC8rowKFk\nbvtkLjNW7eTOXs24u7eWNIuIFFR+m5ZPgefNrD/ggGJmdirwAjAyUsVJIZIxL+jST6HlOUd9y7Z9\nh7h25B+s2h7HsAvb0b+LljSLiBRkfpuWh4H3gfV42/YvDXz9FHg6IpVJ4bFsvLcF/4Gt0PUmOO3h\nw/KCMrN8236uGzmb/QeTee/aLpzaXEuaRUQKOr+byyUDV5jZo0BHoAgwzzm3KpLFSQGXRV7Qkfyy\neie3fDSH0iWKMuqWE2hzTMVcKFZERPKa731aAJxza8wsNnA/LjIlSYGXlgp/vvd3XlDvx+GE2w/L\nCzqSMXM28cCYhTSuXpaR13WlTiUtaRYRKSxC2VxuIDAIb0IuZrYFeAl4xTnnIlOeFDg+8oIy45zj\n9R9X8+KUlXRv4i1prlhaS5pFRAoTv5vLDQNuAp4Hfg0cPgF4FKiNF6gocmTJB+GnYTDr1UBe0Aho\n1/8feUGZvjU1jUe+XsznszfSr2Mdnr2wnZY0i4gUQn5HWm4AbnDOjQ469qOZrQDeRk2LHE3MdPh2\nIOxZe9S8oMzEJabw70/m8tPKHdxxWlMGnd5cS5pFRAqpUOa0LDzCMf2XVzIXvwsmPwQLPssyLygz\nsfsPcd3I2ayIPcDQfm25rGv9iJUqIiL5n9+m5UPg38BdGY7fCnwU1ook+v0jL+heOOXeI+YFZWZl\n7AGufe8P9h5M5n/XdKZnixoRLFhERKKB36alJHC5mZ0J/BY4djxwDPCJmb2a/kLn3J3hLVGiyq41\nXl7Q2p+gblfoO/yIeUFHMmvNTm7+aA6lihdl1M0ncGwdLWkWERH/TUtLYG7gfoPA122BW6ug12kV\nUWEVnBdUtASc8yIcd/0R84KO5Ot5m7lv9AIaVi3LyOu6ULdymQgVLCIi0cbv5nI9I12IRLGNs71l\nzD7zgjLjnOPN6Wt4/vsVdGtchbev7EzFMlrSLCIif/O75Lmhc27dEZ7r7pybFdaqJDpkIy8oMymp\naTzyzRI++2MD53c4hmEXtaNksaIRKFhERKKZ38tDC8zs3865j9MPmFkR4HHgAbw5L1KYBOcFHX+z\nlxdUsnzIp4lPTOH2T+cybcUObuvRhHvPaEGRIlrSLP9v787Dq6qu/4+/l0EEAqjIDIJAGAWVQQVs\nrYgiiNZWq0IdQBx+OFStU6v2q7SiBQeslToVFRwKWi1WAUHUOhWHiEZmCjJDSAjzPCTr98c5wSTc\nkBDvkJt8Xs9znpt77jlnr82QrOyzz14iIgcqbdJyF/CMmZ0LDAXqAa8CTYFD/9VaklfBekENOsKl\nr0DTrmW6VPaWXQwZl868NVt48JcduezU5iWfJCIilVZp57Q8a2YfEyQqc4CjgOnAue6+IYbxSXlx\nQFNe62gAABoaSURBVL2gP0KPG0tVLyiSRVlbGfxiOht37OH5QSfTq50eaRYRkYM7lEc7MoFlQAOg\nOjBVCUslkTUXXjgHptwRVGG+4XP4ya1lTlg+W5TDhU/PYPe+PF67rocSFhEpk6eeeooWLVpQrVo1\nunbtyqeffnrQ47du3cqtt95K8+bNqV69Oj179iQ9Pb3QMcOGDcPMCm0NGzaMZTfkEJQqaTGz04HZ\nBLeDjgeGAI+Y2UQzK9167JJ89u4MRlaePR02LIEL/w5XTCxVgcPi/OPLFQx68SsaH1mdt27sSaem\nWoNFRA7da6+9xi233MI999zDt99+S8+ePenXrx8rVqwo9pxrrrmGadOmMW7cOGbPnk2fPn0466yz\nWL16daHj2rZtS2Zm5v5t9uzZse6OlFJpR1reJ1gV9zR3X+zuLwOdCea26G+zIlryETzVAz4bBSdc\nCjd9XeoCh5Hk5jkPTp7HPRNn85O0urxxfQ+twSJSAeXk5GBmPP7445x88slUq1aN1q1bM3Xq1Ki2\nM2rUKAYPHsy1115L+/btefLJJ2nUqBFPP/10xON37tzJm2++yYgRIzjjjDNIS0tj2LBhpKWlHXBO\nlSpVaNiw4f6tXr16UY1dyq60SctZ7v4Hd9+Xv8PdlwKnA0/FJDJJjO3rYeJQeOmCIEG58m34xVNQ\no06ZL7ljzz6GvjKTv3+6lEE9mvP8oG7UqqY1WEQqooyMDADGjBnDyJEjmTVrFieeeCKXXXYZO3bs\nKHTsQw89RM2aNQ+6Rbrls2fPHmbOnEmfPn0K7e/Tpw8zZkRegWPfvn3k5uZSrVq1QvurV6/OZ599\nVmjfkiVLaNy4MS1atGDAgAEsWbLkkP8cJDZKOxH3k2L25wHDoxqRJIY7fDcBpt1T5npBkazdvIur\nx6UzP3MLw87vwODTWkQpYBEpjzIyMkhJSWHixIm0adMGgJEjR5KWlsbChQvp3Lnz/mOHDh3KJZdc\nctDrNWnS5IB9OTk55Obm0qBBg0L7GzRowPvvvx/xOrVq1aJHjx4MHz6cjh070rBhQ8aPH8/nn39O\nWlra/uNOPfVUxo4dS7t27cjOzmb48OH07NmTuXPncswxmg2RaAdNWsxsBsETQpvC938GHsmfgGtm\ndYFv3F3ld5NZFOoFRTJn9WauHpfOtl379ISQSCWRkZHB+eefvz9hAahdu3bEY+vUqUOdOmUfxT1U\nL7/8MkOGDKFp06akpKTQpUsXBg4cyMyZM/cf069fv0Ln9OjRgxYtWjBu3Dhuu+22uMUqkZV0e6g7\nULXA+xsJHnfOl0IwObdEZvaCmWWb2ZwC++qY2XQzWxS+Hl3awCUKcvfCp6Pg6Z6w5tugXtCQaVFJ\nWKbPy+LiZz4nxYw3ru+phEWkksjIyCg0mgIwY8YMqlWrRtu2bQvtL+vtobp165KSkkJWVlah/VlZ\nWQd90qdVq1Z8/PHHbNu2jZUrV/LVV1+xd+9eWrYs/uGC1NRUjj/+eBYtWlSa7kuMlXZxuXyRZmGW\ntkjiWGA0wYTefL8HPnD3EWb2+/D97w4xJimLA+oFPQK1G/3oy7o7Yz5dykPvzueEJkfy9yu7Ub92\ntZJPFJGkt2vXLhYuXEheXl6h/Y899hgDBgygRo3Ck+/LenuoatWqdO3alenTp3PxxRfv3z99+nQu\nuuiiEuNMTU0lNTWVjRs3Mm3aNB5++OGD9mnBggX06qUSfOXBoSYtZebun5jZcUV2XwCcEX49DvgI\nJS2xVbBeUO3GMGA8tDs3Kpfem5vHfWENoX4dGzLqkpOoXlU1hEQqizlz5uDujB8/nt69e1O/fn2G\nDx/O4sWLef311w84/sfcHrrtttu44oorOOWUUzjttNN45plnWLNmDUOHDt1/zOjRoxk9ejQLFiwA\nYNq0aeTl5dGuXTsWL17MnXfeSbt27bjqqqv2n3PHHXdw/vnn06xZM7Kzs3nggQfYvn07gwYNKlOc\nEl0lJS3OgSMppR1ZKY0G7p4Zfr2WYOE6iZUo1QuKZPPOvdz46jd8tjiH689oxZ2qISRS6WRkZNC6\ndWvuv/9+Bg4cyIYNGzjzzDP59NNPo75A26WXXsr69esZPnw4mZmZdOzYkSlTptC8+Q/lQHJycli4\ncOH+95s3b+buu+9m1apV1KlTh4suuogHH3yQww//4WnGVatWMXDgQHJycqhXrx7du3fniy++KHRd\nSRxzLz4HMbM8guX6d4e7+gEfA/nPrR1B8Dh0qX6dDkdaJrl7x/D9Jnc/qsDnG9094rwWM7sOuA6g\nWbNmXZcvX16aJgWK1AvqFEy0LWO9oEhWrN/BkHHpLMvZzkMXduKSbsdG7doikjxuuukmsrOzI46q\niByMmc10924lHVfSSMu4Iu9fiXDMSxH2lVaWmTVy90wzawRkF3eguz8HPAfQrVu3aI72VFyF6gXt\n+9H1giKZuXwD1740k9w85+WrT6VHKz0SKFJZZWRk0Ldv30SHIRXYQZMWd7/qYJ9HwdvAIGBE+Prv\nGLdXeWTNhbdvhtVfQ6szof8oqBPdNVL+nbGaO9+YReMjq/HC4JNpWa9mVK8vIsnD3Zk1axZ33XVX\nokORCixuE3HNbDzBpNu6ZrYKuJ8gWXndzK4GlgMHn0YuJdu7Ez5+GGb8FaodFdQL6nRxmZffj8Td\neeKDRfzl/UWc0qIOz17elaNTq5Z8oohUWGbGli1bEh2GVHDxfHpoYDEf9Y5XDBXe9/8JFonbuBRO\nuhz6PPCjlt+PZNfeXH7/5izeyljDRV2a8tCFHTmiip4QEhGR2Itb0iIxtD0Hpt0LsyZAnVZBvaCW\nP4t6M+u37ea6l2cyc/lG7jynLTec0QqL4giOiIjIwShpSWZF6wWdfmdQM+jw6C/mtjh7K1eNTSd7\ny25G/7oz553QOOptiIiIHIySlmQVo3pBkXy2KIfrX53JEVUOY8J13encTNUWREQk/pS0JJvcvcEk\n248fhpSqwVNBXa+Cw0oqI1U2479awR/emkNavZo8P7gbTY+uUfJJIiIiMaCkJZmsTId3bobsedD+\n59Dv4ajUC4okN88ZOXUBz32yhJ+1qcfoX3emVrXore8iIiJyqJS0JINdW+CDP0H6mKjXC4pkx559\n3DIhg+nzsriyR3PuO68DVVJiM5IjIiJSWkpayrv574T1gtbCqUPhzHujVi8okrWbd3HNS+nMW7OF\nYed3YPBp0V2QTkREpKyUtJRXm1fDu3f9UC/o0lejWi8okjmrN3P1uHS27drHmEHdOLOd6leKiEj5\noaSlvMnLhfTng9tBMaoXFMn0eVncPP5bjqpxOP8c2pMOjWvHtD0REZFDpaSlPFk7B965Jab1gopy\nd57/bCkPTplPpyZHMubKbtSvHf11XkRERH4sJS3lwd6d8PFImPFkzOoFRbJ99z6GT57P+K9W0Pf4\nhjx+6UlUr6ol+UVEpHxS0pJocagXVFRenvNWxmpGTl1A1pbdDP1ZK+46py2HHaYl+UVEpPxS0pIo\ncaoXVNQ3Kzbyx3fm8d3KTZzY9EieuqwrXZtrhVsRESn/lLTEmzt8Nz5IWHZvjWm9oIIyN+9k5LsL\neCtjDfVrHcFjF5/ILzs30eiKiIgkDSUt8bT+e5h0Kyz9BI49NagXVL99TJvctTeX5z5ZwtMffU+u\nOzf2asUNZ6SReoT+6kVEJLnoJ1c87NsT1Av65JG41AuC4KmgSbMyGfHuAlZv2sm5nRpyd7/2HFtH\ntYNERCQ5KWmJtZVfBY8xZ8+DDhdA35ExqxeUb87qzfzxnbmkL9tI+0a1eeySE+ne8piYtikiIhJr\nSlpiZdfmsF7Q83GpFwSQvXUXj05byD9nrqJOjar8+cJOXNLtWFI0b0VERCoAJS2xEOd6Qbv35fLi\nf5cx+sPF7N6XyzU/acFveremtqoyi4hIBaKkJZo2rw6SlYWTg3pBA16FJrGrF+TuvDcvi4emzGf5\n+h2c1b4+9/bvQIu6qTFrU0REJFGUtERDXi6kj4EPHgjqBZ39J+h+Q0zrBS1Yu4UHJs3jv4vX07p+\nTV4acgqnt6kXs/ZEREQSTUnLj7V2DrxzM6yeCa16w3mj4OjjYtbchu17GDV9If/4cgW1qx/Ony44\nnl+f0owqKbF7EklERKQ8UNJSVgfUCxoDnX4Vs3pBe3PzePnz5fzl/f+xfU8uV/Y4jlvPas1RNarG\npD0REZHyRklLWXz/YVgvaBl0vhzOjm29oP8szGb4pHl8v247P21dl/vO60DrBrGb2CsiIlIeKWk5\nFAXrBR2TBoMmQYufxqy5xdnbGD55Hh8tXEeLuqk8P6gbZ7arj8W4+rOIiEh5pKSlNA6oF3QX/PT2\nmNUL2rxjL098sIiXPl9G9cNT+EP/9lzZ4ziqVtG8FRERqbyUtJSkUL2g7mG9oHYxaWpfbh4T0lfy\n2HsL2bRzLwNObsbtfdpQt+YRMWlPREQkmShpKU5+vaCPH4Yq1eC8x6HL4JjVC5qxOIc/TZrHgrVb\n6d6yDveddzwdGteOSVsiIiLJSElLJIXqBf0C+o2EWg2j3syaTTuZMjuTybMz+XbFJpoeXZ2nL+tC\n344NNW9FRESkCCUtBRWqF9QEBk6Atv2i2kTm5p1Mmb2WybPW8M2KTQAc37g2f+jfnsu7N6fa4SlR\nbU9ERKSiUNKSb8cGeLonbMuC7tdDr3vhiJpRufTazbuYMjuTKbMz+Xr5RgDaN6rNnee05dxOjbTs\nvoiISCkoaclXow50uRLanBOVekFZW35IVNKXBYlKu4a1uKNPG87t1IiW9aKTEImIiFQWSloK6nXP\njzo9e8su3p2zlsmzMklfvgH3IFG57ewgUUmrr0RFRESkrJS0/EjZW3cxdc5aJs3KJH1ZkKi0aVCT\nW3u3of8JDUmrr5VrRUREokFJSxms27qbqXOCp36+XBokKmn1a3JL79b079RIS+yLiIjEgJKWUsrZ\ntpup4a2fL5euJ8+hVb1UfnNma847oRFtlKiIiIjElJKWg1i/bTdT5waJyhdLgkSlZb1UbuqVRv8T\nGtOmQU2tpyIiIhInSlqKWL9tN9PmZjF59ho+/z5IVFrUTeXGXmmc26kR7RrWUqIiIiKSAEpaQtt2\n7+P6V2Yy4/v15OY5xx1Tg+vPaEX/To1p30iJioiISKIpaQmlVk2hymHG/zu9Jf1PaESHRrWVqIiI\niJQjSlpCZsaLV52S6DBERESkGLEpWSwiIiISZUpaREREJCkoaREREZGkoKRFREREkoKSFhEREUkK\nSlpEREQkKZSLpMXM+prZQjNbbGa/T3Q8IiIiUv4kPGkxsxTgb0A/oAMw0Mw6JDYqERERKW8SnrQA\npwCL3X2Ju+8BJgAXJDgmERERKWfKQ9LSBFhZ4P2qcF8hZnadmX1tZl+vW7cubsGJiIhI+VAekpZS\ncffn3L2bu3erV69eosMRERGROCsPSctq4NgC75uG+0RERET2Kw9JSzrQ2sxamFlVYADwdoJjEhER\nkXLG3D3RMWBm5wJ/AVKAF9z9wRKOXwcsj0dscVQXyEl0EHGmPlcO6nPlUNn6XNn6C7Htc3N3L3Hu\nR7lIWgTM7Gt375boOOJJfa4c1OfKobL1ubL1F8pHn8vD7SERERGREilpERERkaSgpKX8eC7RASSA\n+lw5qM+VQ2Xrc2XrL5SDPmtOi4iIiCQFjbSIiIhIUlDSkmBmdqyZ/cfM5pnZXDO7JdExxYOZpZjZ\nt2Y2KdGxxIuZHWVmb5jZAjObb2Y9Eh1TLJnZb8N/03PMbLyZVUt0TNFmZi+YWbaZzSmwr46ZTTez\nReHr0YmMMdqK6fMj4b/rWWY20cyOSmSM0RapzwU+u93M3MzqJiK2WCmuz2b2m/Dveq6ZPRzvuJS0\nJN4+4HZ37wB0B26sJFWubwHmJzqIOHsCmOru7YATqcD9N7MmwM1AN3fvSLAG04DERhUTY4G+Rfb9\nHvjA3VsDH4TvK5KxHNjn6UBHdz8B+B9wd7yDirGxHNhnzOxYoA+wIt4BxcFYivTZzHoRFDQ+0d2P\nBx6Nd1BKWhLM3TPd/Zvw660EP8gOKBhZkZhZU6A/MCbRscSLmR0JnA48D+Due9x9U2KjirkqQHUz\nqwLUANYkOJ6oc/dPgA1Fdl8AjAu/Hgf8Iq5BxVikPrv7e+6+L3z7BUE5lgqjmL9ngMeBu4AKNzm0\nmD5fD4xw993hMdnxjktJSzliZscBnYEvExtJzP2F4D96XqIDiaMWwDrgxfC22BgzS010ULHi7qsJ\nfgtbAWQCm939vcRGFTcN3D0z/Hot0CCRwSTAEODdRAcRa2Z2AbDa3b9LdCxx1Ab4qZl9aWYfm9nJ\n8Q5ASUs5YWY1gTeBW919S6LjiRUzOw/IdveZiY4lzqoAXYCn3b0zsJ2Kd9tgv3AexwUEyVpjINXM\nLk9sVPHnweOZFe638OKY2b0Et7xfTXQssWRmNYB7gPsSHUucVQHqEExluBN43cwsngEoaSkHzOxw\ngoTlVXf/V6LjibHTgJ+b2TJgAnCmmb2S2JDiYhWwyt3zR9HeIEhiKqqzgKXuvs7d9wL/AnomOKZ4\nyTKzRgDha9yH0BPBzAYD5wGXecVfS6MVQUL+Xfi9rCnwjZk1TGhUsbcK+JcHviIYLY/rBGQlLQkW\nZqnPA/PdfVSi44k1d7/b3Zu6+3EEEzM/dPcK/xu4u68FVppZ23BXb2BeAkOKtRVAdzOrEf4b700F\nnnhcxNvAoPDrQcC/ExhLXJhZX4Jbvj939x2JjifW3H22u9d39+PC72WrgC7h//OK7C2gF4CZtQGq\nEueikUpaEu804AqCEYeMcDs30UFJTPwGeNXMZgEnAQ8lOJ6YCUeU3gC+AWYTfK9J+Gqa0WZm44HP\ngbZmtsrMrgZGAGeb2SKCEacRiYwx2orp82igFjA9/B72TEKDjLJi+lyhFdPnF4CW4WPQE4BB8R5V\n04q4IiIikhQ00iIiIiJJQUmLiIiIJAUlLSIiIpIUlLSIiIhIUlDSIiIiIklBSYuIFGJmd4QLZsW6\nnWGRquZGuQ03s1/Fso1EtCVSWSlpEUlSZtbAzJ4ws+/NbLeZrTazd8vDOj9mlmJmvzOz+Wa2w8w2\nmtnXZnZzgcMeBX4W41AaAe/8mAuY2dgwIcnfcsxskpm1K64tMzsuPLbbj2lbRAqrkugAROTQhcU1\n/wtsBe4GviP4JaQ38AzQLFGxhe4HbgBuAr4CUgmKgTbPP8DdtwHbYhlEFFcofZ9gEUgIaik9AkwE\n2segLREphkZaRJLTU+FrN3d/3d0Xuvt8dx8NnJB/kJk1M7OJZrY13P5lZk0LXsjM7jKztWa2zcxe\nAmoWbczMrjKzeWa2y8z+Z2a/NbODff/4OfCMu09w9yXhsucvufsDBa65//ZQgZGJotuyAsd3MLPJ\nYT+yzWx8SbVeCt6yKdDGRWY2PRwBmmdmZx/sGqHd7r423L4BHgfamVn1SG0BS8PX9HD/R+Exnczs\nAzPbEv55f2dmvUrRvoigpEUk6ZhZHaAv8LdwtKIQd98UHncYQd2bBgT1QnoRjBK8lV+Z1cwuAYYT\njIx0ARYCtxVp71qCkgP3EYws3A78jmAkpThrgTPMrEEpu7WS4PZK/tYGWA58FMbQCPgEmAOcQrA8\nfk3g3yUkT5E8CPwVOBFIByZYUGW9VMysFnApMNvddxZz2Cnha9+wPxeG7/8BZIafnwQMA3YdYvwi\nlZZuD4kknzTAKLkAYW+CUZdW7r4MwMx+DSwOP3sfuBUY5+7Phuc8GP7mn1bgOv8H3OXub4Tvl5rZ\nCIKkZXQxbd9GUHso08zmE9QwmQJMjFSrxN1zCRKd/GRrDMEP96HhIdcD37n77/LPMbMrgQ1AN4Jb\nUKX1uLvnzz25B7iSIIH47CDn9DWz/AQxlSDJOtjcoXXh6/oit42aA4+6+4Lw/eJDiFuk0tNIi0jy\nsVIe1x5Yk5+wALj7EmAN0KHAMZ8XOW//ezOrBxwLPBveztgW/vAeAbQqrmF3nwd0BE4lSECOAV4H\nJpdiZGQkQbL1C3fPH4XoCpxeJIaV4WfFxlGMWQW+XhO+1i/hnE8IEpuTCEZJPgDeM7NjD7HtUcAY\nM/vQzO6NMJlXRA5CSYtI8lkEOAUmgZZBaSul5n+PGMoPP7RPIkhIjj9oA+557p7u7o+7+y+BwUA/\n4PTizjGzQWFb57t7VpE4JheJ4SSgNTCplH3Jt7dAjPl/DiV9L9zh7ovDLR24BqgNXHcoDbv7MIKE\n8S2gJzDLzIYcyjVEKjMlLSJJxt03ANOAmyLNxTCzo8Iv5wONwyeN8j9rSTCvZV6BY7oXucT+92Hi\nsIbgFtPiotshhp7fZsT5I2bWE3gauNzdvyvy8TcESdLyCHFsPcQ4osHDrUYxn+8JX1MOONF9kbv/\n1d37A88TJEAiUgqa0yKSnG4keOT5azP7P4JbHkYw2fZugkee3w/3v2pmt4TnPUmQAHwYvn8CeMnM\n0gkmvf6K4JbOhgJt3Q88aWabCOalHE4wabeJu/85UnBm9kYY3wyCuSotgD8DWeG+osc3JHiE+Cng\nywJPBeW6+zrgb8C1wGtmNpJgzkhL4BLg9jgkLkcUiOlogke5Uyl+DZhsYCdwTvgE1C6CROZR4J/A\nMoIJ0j8BvoxZ1CIVjEZaRJJQODelCzCdYA7ILIJE5JfAb8NjHLiA4Af8f8JtLcFcEQ+PeY3gCZYH\ngW+BTgTzLgq2NQYYQrBOyXfApwS3RZZSvGlAf+Bt4H/AywRPA/UOR4qKakcwr+R2ggm4+Vt6GMMa\n4DQgD5gKzCVIZHaHW6ydVSCmL4GTgYvd/aNIB7v7PuBmglGUNQRPceUSJDxjCZ7Smkgwf+i2SNcQ\nkQNZhIn8IiIiIuWORlpEREQkKShpERERkaSgpEVERESSgpIWERERSQpKWkRERCQpKGkRERGRpKCk\nRURERJKCkhYRERFJCkpaREREJCn8f4WE0sspLOCUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111f8c828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "logks = np.arange(1,17).astype('int')\n",
    "codelengths = [code.calc(2**logk) for logk in logks]\n",
    "plt.figure(figsize=(9,6))\n",
    "plt.plot(logks,codelengths)\n",
    "plt.plot(logks,logks/entropy(code.p))\n",
    "plt.ylabel('Expected Input Phrase Length',fontsize=14)\n",
    "plt.xlabel('Code Size in Bits',fontsize=14)\n",
    "plt.annotate('$p=0.95$', xy=(12,10),fontsize=14)\n",
    "_ = plt.title('BAC vs Entropy',fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Comments and Relationship to Tunstall Coding\n",
    "\n",
    "The most famous VtoF code is the Tunstall code.  The Tunstall code builds a parsing tree by continually breaking the highest probability leaf into two leaves ($m$ leaves for an m-ary alphabet).  The analagy isn't perfect but BAC is more top-down while Tunstall is bottom-up.  BAC is similar to Shannon coding; Tunstall is similar to Huffman coding.\n",
    "\n",
    "As suggested by the plot above, $\\log(k)/h(p) - EL(k) \\le C$ for some value of $C$ and for all $k$.  I proved this result in the [cite Trans. Info. Th. paper]. While Tunstall is \"optimal\", Tunstall originally proved the weaker result that $$\\frac{h(p)EL(k)}{\\log(k)} \\rightarrow 1$$ (Showing the ratio approaches 1 does not imply the difference between $EL(k)$ and $\n",
    "\\log(k)/h(p)$ is bounded.)\n",
    "\n",
    "BAC has many advantages over Tunstall coding:\n",
    "\n",
    "1. BAC accomodates much larger alphabets.  The encoding is top-down and the codebook never needs to be stored.\n",
    "\n",
    "2. BAC gracefully allows non-identically distributed bits. \n",
    "\n",
    "3. BAC allows dependency between bits as long as $p(x_n) = f(x_1, x_2, \\ldots, x_{n-1})$.  E.g., the input sequence can be Markov or hidden Markov."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "## References\n",
    "\n",
    "1. Boncelet, Charles G. \"Block arithmetic coding for source compression.\" IEEE transactions on information theory 39.5 (1993): 1546-1554.\n",
    "\n",
    "2. Boncelet, C.G., 1993, January. Block arithmetic coding for Markov sources. In Information Theory, 1993. Proceedings. 1993 IEEE International Symposium on (pp. 118-118). IEEE."
   ]
  },
  {
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}