julia/Tutorials/coding_2_entropic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Entropic Coding and Compression\n",
    "===============================\n",
    "\n",
    "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes.\n",
    "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n",
    "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n",
    "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n",
    "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n",
    "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n",
    "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
    "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n",
    "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n",
    "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n",
    "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n",
    "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n",
    "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n",
    "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n",
    "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n",
    "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n",
    "$\\newcommand{\\CC}{\\mathbb{C}}$\n",
    "$\\newcommand{\\RR}{\\mathbb{R}}$\n",
    "$\\newcommand{\\EE}{\\mathbb{E}}$\n",
    "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n",
    "$\\newcommand{\\Ww}{\\mathcal{W}}$\n",
    "$\\newcommand{\\Vv}{\\mathcal{V}}$\n",
    "$\\newcommand{\\Nn}{\\mathcal{N}}$\n",
    "$\\newcommand{\\NN}{\\mathcal{N}}$\n",
    "$\\newcommand{\\Hh}{\\mathcal{H}}$\n",
    "$\\newcommand{\\Bb}{\\mathcal{B}}$\n",
    "$\\newcommand{\\Ee}{\\mathcal{E}}$\n",
    "$\\newcommand{\\Cc}{\\mathcal{C}}$\n",
    "$\\newcommand{\\Gg}{\\mathcal{G}}$\n",
    "$\\newcommand{\\Ss}{\\mathcal{S}}$\n",
    "$\\newcommand{\\Pp}{\\mathcal{P}}$\n",
    "$\\newcommand{\\Ff}{\\mathcal{F}}$\n",
    "$\\newcommand{\\Xx}{\\mathcal{X}}$\n",
    "$\\newcommand{\\Mm}{\\mathcal{M}}$\n",
    "$\\newcommand{\\Ii}{\\mathcal{I}}$\n",
    "$\\newcommand{\\Dd}{\\mathcal{D}}$\n",
    "$\\newcommand{\\Ll}{\\mathcal{L}}$\n",
    "$\\newcommand{\\Tt}{\\mathcal{T}}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\al}{\\alpha}$\n",
    "$\\newcommand{\\la}{\\lambda}$\n",
    "$\\newcommand{\\ga}{\\gamma}$\n",
    "$\\newcommand{\\Ga}{\\Gamma}$\n",
    "$\\newcommand{\\La}{\\Lambda}$\n",
    "$\\newcommand{\\si}{\\sigma}$\n",
    "$\\newcommand{\\Si}{\\Sigma}$\n",
    "$\\newcommand{\\be}{\\beta}$\n",
    "$\\newcommand{\\de}{\\delta}$\n",
    "$\\newcommand{\\De}{\\Delta}$\n",
    "$\\newcommand{\\phi}{\\varphi}$\n",
    "$\\newcommand{\\th}{\\theta}$\n",
    "$\\newcommand{\\om}{\\omega}$\n",
    "$\\newcommand{\\Om}{\\Omega}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This numerical tour studies source coding using entropic coders (Huffman and arithmetic)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#from __future__ import division\n",
    "\n",
    "#import numpy as np\n",
    "#import scipy as scp\n",
    "#import pylab as pyl\n",
    "#import matplotlib.pyplot as plt\n",
    "\n",
    "#from nt_toolbox.general import *\n",
    "#from nt_toolbox.signal import *\n",
    "\n",
    "#%matplotlib inline\n",
    "#%load_ext autoreload\n",
    "#%autoreload 2\n",
    "\n",
    "using PyPlot\n",
    "using NtToolBox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Source Coding and Entropy\n",
    "-------------------------\n",
    "Entropic coding converts a vector $x$ of integers into a binary stream\n",
    "$y$. Entropic coding exploits the\n",
    "redundancies in the statistical distribution of the entries of $x$ to\n",
    "reduce as much as possible the size of $y$. The lower bound for the\n",
    "number of bits $p$ of $y$ is the Shannon bound :\n",
    "\n",
    "$$p=-\\sum_ih(i)\\log_2(h(i))$$\n",
    "\n",
    "where $h(i)$ is the probability of apparition of symbol $i$ in $x$.\n",
    "\n",
    "Fist we generate a simple binary signal $x$ so that $0$ has a probability $p$\n",
    "to appear in $x$.\n",
    "\n",
    "Probability of 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Size.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "512"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 512"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Signal, should be with token 1,2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "512-element Array{Int64,1}:\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 1\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " ⋮\n",
       " 2\n",
       " 1\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 2\n",
       " 1\n",
       " 2\n",
       " 2\n",
       " 2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = (rand(n) .> p) .+ 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can check the probabilities by computing the empirical histogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the entropy of the distribution represented as a vector $h$ of proability that should sum to 1.\n",
    "We take a max to avoid problems with null probabilties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Empirical p = 0.11\""
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = [sum(x .== 1), sum(x .== 2)]\n",
    "h = h/sum(h)\n",
    "\n",
    "@sprintf(\"Empirical p = %.2f\", h[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Entropy = 0.50\""
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e = - sum(h.*log2(collect(max(e,1e-20) for e in h)))\n",
    "@sprintf(\"Entropy = %.2f\", e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huffman Coding\n",
    "--------------\n",
    "A Hufman code $C$ associates to each symbol $i$ in $\\{1,...,m\\}$ a binary code $C_i$\n",
    "whose length is as close as possible to the optimal bound\n",
    "$-\\log_2\\left(h(i)\\right)$, where $h(i)$ is the probability of apparition of the\n",
    "symbol $i$.\n",
    "\n",
    "We select a set of proabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.1 \n",
       " 0.15\n",
       " 0.4 \n",
       " 0.15\n",
       " 0.2 "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = [.1, .15, .4, .15, .2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tree $T$ contains the codes and is generated by an iterative algorithm.\n",
    "The initial \"tree\" is a collection of empty trees, pointing to the symbols numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = length(h)\n",
    "#T = [0] * m # create an empty tree\n",
    "T = zeros(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build iteratively the Huffman tree\n",
    "by grouping together the two erees that have the smallest probabilities.\n",
    "The merged tree has a probability which is the sum of the two selected\n",
    "probabilities.\n",
    "\n",
    "Initial probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "MethodError: Cannot `convert` an object of type Tuple{Float64,String} to an object of type Float64\nThis may have arisen from a call to the constructor Float64(...),\nsince type constructors fall back to convert methods.",
     "output_type": "error",
     "traceback": [
      "MethodError: Cannot `convert` an object of type Tuple{Float64,String} to an object of type Float64\nThis may have arisen from a call to the constructor Float64(...),\nsince type constructors fall back to convert methods.",
      "",
      " in setindex!(::Array{Float64,1}, ::Tuple{Float64,String}, ::Int64) at .\\array.jl:415",
      " in macro expansion; at .\\In[19]:4 [inlined]",
      " in anonymous at .\\<missing>:?"
     ]
    }
   ],
   "source": [
    "#we use the symbols i = 1, 2, 3, 4, 5 (as strings) with the associated probabilities h(i)\n",
    "\n",
    "for i in 1:m\n",
    "    T[i] = (h[i],string(i))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterative merging of the leading probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "while len(T) > 1: \n",
    "    T.sort() #sort according to the first values of the tuples (the probabilities)\n",
    "    t = tuple(T[:2])\n",
    "    q = T[0][0] + T[1][0]\n",
    "    T = T[2:] + [(q,t)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We trim the computed tree by removing the probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trim(T):\n",
    "    T0 = T[1]\n",
    "    if type(T0) == str:\n",
    "        return T0\n",
    "    else:\n",
    "        return (trim(T0[0]),trim(T0[1]))\n",
    "\n",
    "T = trim(T[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We display T using the ete3 package (install it in the terminal with \"pip install ete3\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAGhCAYAAACap/TMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAN1wAADdcBQiibeAAAIABJREFUeJzt3Xt0VfWd9/HPOSfkRi5AEiCAEEIIlyoOknZGRx8Cjs+A\nWgpLxopWBjtTp+rYzkJW1+DQJXHNxRE7q4w42nm8jLjExRIHqzJoESEdtS4LOFM7FLmUy1i5JOR6\ncs7JSXL288eXzc6BBAKEnGT7fq2110myT/bZu9VPvn73b/9+AcdxHAEAfCWY6hMAAPQ+wh0AfIhw\nBwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfIhw\nBwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfCgt\n1SeA3vfGG9LmzdKHH0oHD9rPSkula6+V5s6V5s1L7fkBuPwCjuM4qT4J9I7du6WbbpIiEclxpJwc\nKStLSiSkcFhqaZHS0qTcXGnrVqm8PNVnDOByIdx9YulS6emnpeJiafRoC/CMDNvX1iZFoxbwTU1S\nba3U3Cw9+KC0alVqzxvA5UFbxgeWL5eee066+mppzBipqMjCPS1N6ujwgr2xUWposGq+rk566ikp\nO1uqqvKOdfLkSR0+fFg5OTkqv4TSvreOA+DiEO4D3IsvSmvWWLBPmCCNGyeNGGEtGUmKxaxKb2y0\nUE9Pl0IhKRCw1s0TT0iTJkl33mnv37JlixYtWqQ5c+Zo8+bNpz/n4MGD2rBhgz755BMdOnRIjuOo\ntLRUCxcu1IIFC846r+6OA6BvEO4D3D33WLCXlEhlZbaNGmVBHo9bG6a+3kI9eGpsVCIhtbdbu6a1\nVVq82Av3goKCpFdJamhoUGlp6envCwoK1NLSoo8++kjr1q3TN77xDW3cuFGBQCDpPWceB0DfYSjk\nAHb33VJhofXZi4ulsWOter/ySmniRAv5wkJp6FApP1/Ky7N2TU6ObYMH22tGhv2RkLoP5XHjxmnN\nmjU6ceKEamtrdfLkSW3cuFHDhw/XT3/6Uz3zzDNJ7yfcgdQi3AewjRu98B42zF4LC6WCAgvx7Gzb\nsrJsy8zsesvOll57zY7ZVSjn5+dr//79euCBB1RUVCRJys7O1vz58/UXf/EXkqStW7cmnRvhDqQW\nbZkB6uhRa6vk51twZ2TYDdS2NrtpKlnrxXGsvx4MWq89Lc3bBg2yLSPDevLHj3cdyoFAQGlpXf+j\ncu2110qSamtrk35OuAOpReU+QO3c6YV1IGB99FjMQrqmxkbDNDfbSJl43EbNJBIW9pKFfefQDwal\nHTuknJwcrVq1SpWVlT06j3feeUeSzhoRc6HHAdC7er1yr66u5l/oPvE9DR68Wo5jFXo0amFeV2ff\nu8MgIxHb19pqVX17u4V8IpF8NMeR9u+3r5ctW9ajM4jFYnr99dclSUuWLDlrf0+PA6D39Xrlvn37\ndj3yyCNyHIftMm5vvrla7e0JxePeOPaGBntAqabGXk+etEq+udmeTnVDPh73gt6t6AMBG2lzIe6/\n/34dPnxYd955p6677rre/kcJwCWg5z5AzZghOU4w6QGlzEzbF4tZL12yAI/FLNzdgI/FLODjcdvv\nVvMVFT3//OXLl+uFF17QV7/6Vf3kJz/p/QsEcEkI9wGquNjGrjc1eWPZO7di0tOtGu/osGo9FrM/\nAuFwcqvG3bKy7OGnnlixYoUee+wxTZ8+XW+//bZy3CemAPQbhPsAtmCB9Prr0pAhNuIlGLQq3A33\nYNAqcvdhpUjEq+DdgI/F7HXRovN/XiKR0AMPPKBnnnlGM2fO1BtvvKG8vLzLf6EALhjhPoCtXWsj\nXerqrGqXrNWSnZ38RGp7u/28c3vGnSWypcX2P//8uT8rEonoW9/6ljZu3Kj58+frlVdeUabbBwLQ\n7xDuA9zatdK991rIJxJWobvhHgrZezo6vOo9GrWq3W3RxGLSyy+f+zOOHTumr3/969qxY4eWLl2q\nVatWKRhkFC3QnxHuA9xdd0l799rUvR0dFtZuuLvVvONYuLvVezTqVe0PPyzdfvu5P+PVV1/Vjh07\nNGzYMB05ckTf/OY3u3zfn//5n+uP//iPe/kKAVwMwt0HqqossFevtmGPeXleuLuzP3au3sNhC/qH\nHpJWrOj559TV1WnDhg3d7p81a1YvXA2A3tDri3VUVVXJcRytXLmyNw+LHti3T7rxRhsW2dpqwyHT\n05Mr96wsuwG7bZs0fnyqzxjA5ULl7iMTJ0pHjkibNtkaqh98IB04YPvKyqTrr7c1VG++ObXnCeDy\nI9x96JZbbAPw5cWQBwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0A\nfIhwBwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0A\nfIhwBwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0A\nfIhwBwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0A\nfIhwBwAfItwBwIcIdwDwIcIdAHyIcAcAH0pL9QkAwED2xhvS5s3Shx9KBw/az0pLpWuvlebOlebN\nS815Ee4AcBF275ZuukmKRCTHkXJypBEjpERCOnpUeukl6ZVXpNxcaetWqby8b8+PcAeAC7R0qfT0\n01JxsVXpublSRobta2uTolEpHJaamqTaWmnaNOnBB6VVq/ruHAl3ALgAy5dLzz0nXX21NGaMVFRk\n4Z6WJnV0eMHe2Cg1NEhZWVJdnfTUU1J2tlRV5R3r5MmTOnz4sHJyclR+CaV9V8ch3AGgh158UVqz\nxoJ9wgRp3DhrxeTk2P5YTGputmDPypLS06VQSAoErHXzxBPSpEnSnXfa+7ds2aJFixZpzpw52rx5\nc9JnRSIR/eQnP9HHH3+sQ4cOadq0aZo1a5buuOOOs86rq+MQ7gDQQ/fcY8FeUiKVldk2apQFeTxu\nbZj6egv14KmxiImE1N5u7ZrWVmnxYi/cCwoKkl5dtbW1+vrXv66PPvpImZmZGjdunJ599ln967/+\nqz744AOtXr1awaA32LGr4zAUEgB64O67pcJC67MXF0tjx1r1fuWV0sSJFvKFhdLQoVJ+vpSXZ+2a\nnBzbBg+214wM+yMhdR/uDz30kD766CP95V/+pRoaGrRnzx6dOHFCFRUVWrNmjV5++eWk9xPuAHCR\nNm70wnvYMHstLJQKCizEs7Nty8qyLTOz6y07W3rtNTtmV6Hc1NSkDRs2qLS0VD/+8Y+VcepObUFB\nwelQf/7555POjXAHgItw9Ki1VfLzLbgzMuwGalub3TQNh6314jjWXw8GrdeeluZtgwbZlpFhvfnj\nx7sO5U2bNikSiWjWrFkKhUJJ51FeXq6SkhJVV1erpqbm9M8JdwC4CDt3emEdCFgfPRazG6c1NTYa\nprnZRsrE4zZqJpGwsJcs7DuHfjAo7dgh5eTkaNWqVaqsrDz9WZ9//rkkadq0aV2ey7Rp0+Q4zun3\nSV0fhxuqwABWXV2d9C80LpfvafDg1XIcq9CjUQvzujr73h0GGYnYvtZWq+rb2y3kE4nkozmOtH+/\nfb1s2bKkfcePH5ck5ebmdnkm7s+PHTuW9PMzj0PlDgxg27dv1yOPPCLHcdgu4/bmm6vV3p5QPO6N\nY29osAeUamrs9eRJq+Sbm6WWFi/k43Ev6N2KPhCwkTZdOXHihCSrxrsyePBgSUpqy3SFyh0AzmPG\nDMlxgkkPKGVm2r5YzHrpkgV4LGbh7gZ8LGYBH4/bfrear6jo+rOys7MlSW1tbV3uj8fjkryQ7w7h\nDgDnUVxsY9ebmryx7J1bMenpVo13dFi1HovZH4FwOLlV425ZWfbwU1dGnNoRiUS63B8Oh5Pe1x3C\nHQB6YMEC6fXXpSFDbMRLMGhVuBvuwaBV5O7DSpGIV8G7AR+L2euiRd1/jhvaZ/bUXUePHk16X3cI\ndwDogbVrbaRLXZ1V7ZK1WrKzk59IbW+3n3duz4TD3tft7dIZw9ST/P7v/74k6ec///lZ+6LRqH75\ny1+qqKhIJSUl5zxfbqgCQA+tXWvhfvKkdzPV3WprvRurdXXWumlstDZOc7MFfCwmnfFw6VlmzJih\nadOmqbq6Wr/85S+T9j311FOKx+NavHixBrmN/m5QuQNAD911l7R3r03d29FhYe1W7m417zjWmnGr\n92jUq9offli6/fbzf86KFSt0xx13aN68eVqyZIkmTZqkX/ziF3ruuedUWFio++6777zHOGe4NzZK\nu3ZJP/uZDbiXpK99TfqjP5KmT7feEwB8mVRVWWCvXm0VeV6eF+7u7I8dHV7vPRy2oH/oIWnFip59\nxp/8yZ9o0KBBuvfee/XYY4+d/vnVV1+tDRs2aMKECec9RrfhvmiR9NZb3gD9zEw76Q8/lP7pn+xi\n5s2z1UYA4Mvk8cel73xHuvFGK4JbW204ZHp6cuWelWVzz2zbJo0ff2GfMX/+fM2fP1+HDx/WoUOH\ndNVVV2nYsGE9/v2zwn37dunWW+1EJ070VhgJBLwxnOGw/cX693+3u8ebNkn/5/9c2IkDwEA2caJ0\n5Ijl3+bN0gcfSAcO2L6yMun6620N1ZtvvrTPGTdunMaNG3fBv5cU7m6wX3GFzVfsrjAyaJA3ftMN\n9oYG6zXV1NgFvPOOXQwAfJnccott/c3pcP/Nb+wvzJQpNkfx2LHS8OE2/3AgYP+J0XnpKHdWtGDQ\nXv/oj6Rf/zqVlwIAcJ0O91mz7Cms8ePtPzdKS+377Gxrx7jVemZm8o0Dd5WRaFSaOVO6995UXg4A\nQDoV7itWWMvlyittwddx4yzgS0qsJdPUZGM3QyGdnhXNvRPsDvUZMkQ6cULavr1SM2duS/FlAcCX\nW1CSfvITb0WRggLrtbvLSBUV2VAfd5moM1cb6bzqSHq69MtfdjMbDgCgzwTDYXuSatgwL7zd0TEt\nLVaZO07yyiLuiiLu0B93y8qSIpFstbWd+8kpAMDllbZzpxfO7hNWbW3WiklLs0CPxawF405V2Xkp\nKXcLhdzwd/TFF6NSe1UA8CWX9qtfeRPedB7uWF9v36eleTOfuUtIdZ50vvMKI4GAJAV0/PjIVFwL\nAOCUtMJCqa0toba24OkHlBoaLNTdp64SCfu689SV7jJSHR3J6wU6jpSd3ZLq6wKAL7W0GTOkQCB4\nenrKpibrnUvWjun8OK0b/p3nJu68VmB7uyQlNGrUF6m8JgD40ksrL7equ7nZgr2hwavWo1FvhRF3\njuLOE9B3Dnh3GSnHkYYNq0v1dQHAl1qaZA8f7dplY9UzM8+9wkgsZj93J58/c4WR8eMPpvqaAOBL\nL02Snn5auuoqu4nqzv/eeYWRzg8vuQ8udV5dxA34jg7p1lvfSuX1AAB0KtzLymzagP/3/yzI3Ruo\nncNd8uYo7jwBvbsIbFOT9MADUl5eoxwnlZcEADi9zN7q1Tb9wBdf9Gz5qIYGb/mohgZbvONHP0rl\npQAAXElT/n78sXTNNTZD5NChNmrmzBVGOs8rE41a1X7VVdIvfpGqSwD8ZedO2ySposL+nQQu1FmL\ndezaZauM/M3fWKhnZHg3Wd3lo2Ix72Gmxx+Xli5NxakD/rF+vbRypbRnj/37Zg8E2r9zsZhNxf3o\no9LChSk9TQwgXS6z94MfSPfdJy1fLm3ZIh086N1obWuz6YD/7/+VHnvM+vIALo7jSDfcIP3qVza/\nU0WFN7eTW0hFItYWXbJEevJJqbo61WeNgaDbNVRzc6U1a7zv9+2z14kTL/cpAV8Ozc3W/iwstHtW\nhYXeymeOY63PSMTubeXn272uTz6x/fX1NlMr0J1uw/1MhDrQexzHgr2szLYrrvBWPguFvAcG3QcL\ns7Is1NPSbJDDsGH2HqA7PQ53AL3nhhusUi8rk8rLbQW0kSNt2u2ODhtm3Nhobc9Bg+yelzuooaPD\nW/ls9uxUXwn6K8Id6GPr11uPffp0q9jHj5cmTbKvMzK8WVndqT8SCe8Zk9ZW2woKpB07pOLiKZo8\neXeqLwn9UPD8bwHQm1autLZKYaG1YkaOtGCfMsWWuSwosKlAuloBrfMqaIMGSe+9NyvVl4N+inAH\n+tiePbZ8ZW6uF94ZGdZDTyQ6L3yTvNJZerq9z92ys6WamqJUXw76KdoyQB/audPGsWdkWOUdClnL\nJRy2G6WhkN1IbWs7e9WzUOjMVc+ktLQOHTtWnOrLQj9EuOO0qqqqVJ+C7+3cOUOBwK2nb47G43bz\ntL7evg+FLNjd9YvdhwW7WvXMtiDLWqJLhDskWbBv27ZNlZWVqT4VX3Mc5/SIF3cce2OjtVza2qwa\nTyS8hXE6r5fQeXlLd9UzoDuEO06rrKzUypUrU30avrZrl/Tuu96Tp01N1jsPBJJXPovHk1c+6xzy\n7spnHR2S47DyGbpGuAN96JprkoP9XCufuZW9u26CuyiOu+qZBXxII0ceTfVloR8i3IE+NmWKzRWT\nn29DGoNBq8LPXPnMrd47L4wTjXoL1EciUlHRiVRfDvopwh3oY48+apOAuSufBQLWaum8OE4i4VXv\nnRfG6bzyWSIhzZ79XqovB/0U4Q70sYULbXbHTz6xG6juJGHZ2TZEMnjq6ZOODq96dwO+udmbmuC6\n66TJkz/jxiq6RLgDKVBdbVV7TY03ciYry8LdXdbSXZQ+Hrdwj0Qs2JuabN/WrRKjV9Edwh1Ikfp6\nm4YgGrUpBzqHe+eVz+Jxrz3T2GjBXleX6rNHf0e4AymSk2PBPXOmTQI2aJC1Zjova9nRYQEfiVio\nX3edVezA+TC3DJBi1dXS2rXS6NHSsWPSiRP2evy4LUpfUyOVlEjr1hHs6Dkqd6AfuO022yS70bpj\nh31dUWFTAwMXinAH+pnp0wl0XDraMgDgQ4Q7APgQ4Q4APkS4A4APEe4A4EOEOwD4EOEOAD5EuAOA\nDxHuAOBDhDsA+BDhDgA+RLgDgA8R7gDgQ4Q7APgQU/4CwCV44w1p82bpww+lgwftZ6Wl0rXXSnPn\nSvPmpea8CHcAuAi7d0s33WRLIDqOLZs4YoQth3j0qPTSS9Irr0i5ubaCVnl5354f4Q4AF2jpUunp\np6XiYqvSc3NtcXNJamuzxczDYampSaqtlaZNkx58UFq1qu/OkXAHgAuwfLn03HPS1VdLY8ZIRUUW\n7mlptqC5G+yNjVJDg5SVJdXVSU89ZQugV1V5xzp58qQOHz6snJwclV9Cad/VcbihCgA99OKL0po1\n0lVXSZMmSVOnSldeKX3lK/Z9aak0dqwtdj5ypDR8uFRYKBUU2B+BJ56whc5dW7Zs0YwZM/T973+/\n289MJBL63ve+p5deeqnb93R1HMIdAHronnuksjKppMRep0yxcJ8yRZowQbriCmvVFBVZqA8bJg0d\nKuXn25abKy1e7B2voKAg6bUr69at05NPPqlt27Z1+56ujkNbBgB64O67LbCLi20bO9YCfcIEKRiU\n6uul9HR7b3u79d5bW6VYzFo10ajddG1utj8SL7zQfbj/5je/0apVq/TrX/9aO3bsOO+5dXUcKncA\n6IGNGy3chw71KnK35ZKba/307GzrsWdlSZmZXW/Z2dJrr9kxuwv3xsZG7d27V+np6Zo6daokKRAI\ndHtuVO4AcBGOHrVKPD/fgjsjw26gtrXZTVPJqnXHkQIBq+RDIXuPuw0aZFtGht1sPX68+3D/gz/4\nA73//vuSpG3btmn27NnnPD8qdwC4CDt3emEdCNhY9ljMQrqmxkbDNDdb6yUet1EziYSFvWRh3zn0\ng0Fpxw4pJydHq1atUmVl5SWdX1fHoXIHBrDKykpVVlaqqvP4OlwG39PgwavlOFahR6MW5nV19r07\nDDISsX2trVbVt7dbyCcSyUdzHGn/fvt62bJlvXKGZx6HcAcGsJkzZ8pxy0NcNm+9JS1cmFA8Hjw9\njr2hwarw1lZrtyQSVrVHIlJLixfy8bgX9G5FHwjYaJvLiXAHgPOYMUNynGDSA0qZmbYvFrNwlyzA\nYzELdzfgYzEL+Hjc9rvVfEXF5T1nwh0AzqO42IY5NjXZVl+f3IpJT7dqvKPDG/4YDtvWuVXjbllZ\nNg/N5US4A0APLFggvf66NGSIjXgJBq0Kd8M9GLSK3B3f7rZnWlq8gHfHvC9adPnPl3AHgB5Yu9Z6\n7HV1VrVL1mrJzvbCXbLAj8eT2zPhsPd1e7v0/POX/3wJdwDoobVrpXvvtZBPJKxCd8M9FLL3dHR4\n1Xs0alW726KJxaSXX+6bcyXcAaCH7rpL2rvXpu7t6LCwdsPdreYdx8Ldrd6jUa9qf/hh6fbbz/85\nBw4c0CuvvCJJ+u1vfytJ+q//+i/97d/+rSSpoqJCc+bMOecxCHcAuABVVRbYq1fbWPe8PC/cAwEL\n987VezhsQf/QQ9KKFT37jH379umHP/xh0s927dqlXbt2SZLuv/9+wh0Aetvjj0vf+Y504402LNId\n656enly5Z2XZ3DPbtknjx/f8+HPmzLnk5xcIdwC4CBMnSkeOSJs22RqqH3wgHThg+8rKpOuvtzVU\nb745NedHuAPAJbjlFtv6GyYOAwAfItwBwIcIdwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAH\nAB8i3AHAhwh3APAhwh0AfIhwBwAfItwBwIeY8tfHdu60TZIqKqRrrknt+QDoO4S7z6xfL61cKe3Z\nI2Vm2rJfkq0OE4tJU6ZIjz4qLVyY0tMEcJkR7j7hONINN0i/+pU0bJhV6hkZFu7uQr6RiFRbKy1Z\nIj35pFRdneqzBnC5EO4+0NwsDR0qFRZK06fba26urenoOLa+YyQiNTVJ+flSfb30ySe2v75eyslJ\n9RUA6G2E+wDnOBbsZWW2XXGFNHy4BXYoZIv0usHe0GAL9g4aZCu119RYlR+Pp/oqAPQ2wn2Au+EG\nq9TLyqTyclthfeRIafBga8e0tNjq7NnZFurBoLVqHMf2R6PSzJnS7NmpvhIAvYlwH8DWr7ce+/Tp\nVrGPHy9NmmRfZ2RI4bC1XdLTLdATCQv0tjZr1bS2SgUF0o4dUnHxFE2evDvVlwSglzDOfQBbudLa\nKoWF1ooZOdKCfcoUacwYC+4hQ6S8PGvTDB5sFby7ZWV5bZr33puV6ssB0IsI9wFszx6pqMhunrrh\nnZFhPfREwnruoZD119PTk7eMDG/LzpZqaopSfTkAehFtmQFq504bx56RYZV3KGQtl3DYbpSGQnYj\nta3Ngt5xrDUTDNo+99UN/7S0Dh07VpzqywLQSy5LuFdXV6uqqupyHBqn7Nw5Q4HAradvjsbjdvO0\nvt6+D4Us2FtabIx7PC61t1vQJxLecQIBdwvqiy9Gpe6CAPSqXg/3Rx55RI7jyHGc3j40OrH/ja1a\nd8exNzZay6WtzarxRMKCPRy219ZWL+Q7OryKnv+rAP+5LJX7ypUrL8dh0cmuXdK773pPnjY1We88\nELCfpad7Fb0b8JFIcsi3tXlB7zgJjRr1RaovC0Avoec+QF1zTXKwNzRY7z2RsLHr7vDH9navsg+H\nrU0TjXqtGq+SD2nkyKOpviwAvYRwH8CmTLG5YvLzbUhjMGhVeCRi4R4MWti71XtLS3LAR6Ne8BcV\nnUj15QDoRYT7APboozYJWH29Ve2BgLVasrMt3EMhC3e3eo9GkwM+ErGfJRLS7NnvpfpyAPQiwn0A\nW7jQZnf85BO7gepOEpadbUMkg6eeYujo8Kp3N+Cbm72pCa67Tpo8+TNurAI+QrgPcNXVVrXX1Hgj\nZ7KyLNxDIXtPImEVfTxu4R6JWLA3Ndm+rVslRq4C/kK4+0B9vU1DEI3alAOdw90dB9/ebuHutmca\nGy3Y6+pSffYALgfC3Qdyciy4Z860ScAGDbLWTFpa8gyQ7e1WtScS1orZujXVZw7gcmFuGR+prpbW\nrpVGj5aOHZNOnLDX48elkyetdVNSIq1bR7ADfkfl7jO33WabZDdad+ywrysqbGpgAF8OhLuPTZ9O\noANfVrRlAMCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfIhwBwAfItwBwIcIdwDwIcIdAHyI\ncAcAHyLcAcCHCHcA8CGm/AWAS/DGG9LmzdKHH0oHD9rPSkula6+V5s6V5s1LzXkR7gBwEXbvlm66\nyZaudBxb7nLECFvG8uhR6aWXpFdekXJzbeWz8vK+PT/CHQAu0NKl0tNPS8XFVqXn5tqi9JLU1maL\n0IfDUlOTVFsrTZsmPfigtGpV350j4Q4AF2D5cum556Srr5bGjJGKiizc09JsIXo32BsbpYYGKStL\nqquTnnrKFq6vqvKOdfLkSR0+fFg5OTkqv4TSvqvjEO4A0EMvviitWWPBPmGCNG6ctWJycmx/LCY1\nN1uwZ2VJ6elSKCQFAta6eeIJadIk6c477f1btmzRokWLNGfOHG3evPn05xw/flwbNmzQxx9/rIMH\nD6qlpUWjRo3S7Nmz9e1vf1v5+flJ59XVcQh3AOihe+6xYC8pkcrKbBs1yoI8Hrc2TH29hXrw1FjE\nREJqb7d2TWurtHixF+4FBQVJr65JkyapsbFRkpSXl6f8/Hz993//t9566y09+eSTev/99zVq1KjT\n7+/qOAyFBIAeuPtuqbDQ+uzFxdLYsVa9X3mlNHGihXxhoTR0qJSfL+XlWbsmJ8e2wYPtNSPD/khI\n3Yd7YWGh/v7v/15HjhxRY2Ojjhw5oi+++EK33HKLDh48qL/+679Oen9Xx6FyB4Ae2LjRAn3oUGnY\nMHstLJQKCqwyb221nnpWlm2ZmV1v2dnSa69JL7zQfbjv3r1b6enpST8bPny4Hn/8cW3atEkff/xx\n0j7CHQAuwtGj1lbJz7fgzsiwG6htbXbTVLKAdxzrrweD1mtPS/O2QYNsy8iwnvzx492H+5nB7mpp\naZEkjR54qpMZAAARh0lEQVQ9OunntGUA4CLs3OmFdSBgffRYzEK6psZGwzQ320iZeNxGzSQSFvaS\nhX3n0A8GpR07pJycHK1atUqVlZXn/HzHcfTpp5/qr/7qr5SZmamVK1cm7e/qOFTuwABWXV193mBA\nb/ieBg9eLcexCj0atTCvq7Pv3WGQkYjta221qr693UI+kUg+muNI+/fb18uWLev2Ux3H0c0336xd\nu3bpxIkTuuOOO7Ru3TqNGzfurPeeeRwqd2AA2759ux555BE5jsN2Gbc331yt9vaE4nFvHHtDgz2g\nVFNjrydPWiXf3Cy1tHghH497Qe9W9IGAjbQ5H8dxtGPHjtMjZ/bs2aOtW7fKcf+T4BwIdwA4jxkz\nJMcJJj2gVFdngV5TY9vJkzYM8syAj8Us4ONxC3i3mq+oOP/nBoNB1dTUKBaL6bPPPtMVV1yhP/uz\nP1NV5yehuvvdXrhuAPC14mIbu97U5I1l7xzsbuVeV2f7mpos4MPh5FaNu2Vl2cNPF6K8vFzr169X\nZmamHnvsMbW3t5/z/YQ7APTAggVeZe5W7W5bpnPA19dby6apycK9pcUL+FjMXm+77eLOISsrS2Vl\nZWptbdV+t2nfDW6oAkAPrF1rI13q6uwGqmStluzs5CdS29vt57GYBXtLixfyLS22//nnL+4cHMfR\n//7v/yoQCCQ9odoVwh0AemjtWuneey3kEwnvwSV3DhnJbpq6Uw1Eo1a1h8O2xWLSyy+f+zM++OAD\nHT9+XPPnz1cwmNxc+Yd/+Ac1NjaqoqJCeXl55zwO4Q4APXTXXdLevTZ1b0eHhbUb7m417zgW7m71\nHo16VfvDD0u3337uzzhw4ID+9E//VGVlZbrhhhs0depUtba26mc/+5l+/vOfKz09Xc8+++x5z5Vw\nB4ALUFVlgb16td00zcvzwt2d/bFz9R4OW9A/9JC0YsX5j19RUaGFCxdq8+bNeuGFF5L2zZo1S//8\nz/+sK6+88rzHIdwB4AI9/rj0ne9IN95oN1hbW21qgfT05Mo9K8vmntm2TRo/vmfHnjp1ql599VW1\nt7fr888/1+HDh5WZmany8nINHTq0x+dIuAPARZg4UTpyRNq0ydZQ/eAD6cAB21dWJl1/va2hevPN\nF3f8tLQ0lZSUqKSk5OJ+/+I+FgAgSbfcYlt/wzh3APAhwh0AfIhwBwAfItwBwIcIdwDwIcIdAHyI\ncAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfIgpf4F+ZudO2ySpokK65prUng8G\nJsId6AfWr5dWrpT27JEyM225NslW9YnFpClTpEcflRYuTOlpYgAh3IEUchzphhukX/1KGjbMKvWM\nDAt3dwHmSESqrZWWLJGefFKqrk71WWMgINyBFGluloYOlQoLpenT7TU319bidBxblzMSkZqapPx8\nqb5e+uQT219fL+XkpPoK0J8R7kAKOI4Fe1mZbVdcIQ0fboEdCtniym6wNzTYQsuDBklpaVJNjVX5\n8XiqrwL9GeEOpMANN1ilXlYmlZdL48dLI0dKgwdbO6alRWpslLKzLdSDQWvVOI7tj0almTOl2bNT\nfSXorwh3oI+tX2899unTrWIfP16aNMm+zsiQwmFru6SnW6AnEhbobW3WqmltlQoKpB07pOLiKZo8\neXeqLwn9EOPcgT62cqW1VQoLrRUzcqQF+5Qp0pgxFtxDhkh5edamGTzYKnh3y8ry2jTvvTcr1ZeD\nfopwB/rYnj1SUZHdPHXDOyPDeuiJhPXcQyHrr6enJ28ZGd6WnS3V1BSl+nLQT9GWAfrQzp02jj0j\nwyrvUMhaLuGw3SgNhexGalubBb3jWGsmGLR97qsb/mlpHTp2rDjVl4V+iHDHadXV1aqqqkr1afja\nzp0zFAjcevrmaDxuN0/r6+37UMiCvaXFxrjH41J7uwV9IuEdJxBwt6C++GJU6i4I/RbhDknSI488\nIsdx5DhOqk/F1+x/Y6vW3XHsjY3Wcmlrs2o8kbBgD4fttbXVC/mODq+i5/8qnAvhjtNWrlyZ6lPw\nvV27pHff9Z48bWqy3nkgYD9LT/cqejfgI5HkkG9r84LecRIaNeqLVF8W+iHCHehD11yTHOwNDdZ7\nTyRs7Lo7/LG93avsw2Fr00SjXqvGq+RDGjnyaKovC/0Q4Q70sSlTbK6Y/Hwb0hgMWhUeiVi4B4MW\n9m713tKSHPDRqBf8RUUnUn056KcId6CPPfqoTQJWX29VeyBgrZbsbAv3UMjC3a3eo9HkgI9E7GeJ\nhDR79nupvhz0U4Q70McWLrTZHT/5xG6gupOEZWfbEMngqadPOjq86t0N+OZmb2qC666TJk/+jBur\n6BLhDqRAdbVV7TU13siZrCwL91DI3pNIWEUfj1u4RyIW7E1Ntm/rVomRq+gO4Q6kSH29TUMQjdqU\nA53D3R0H395u4e62ZxobLdjr6lJ99ujvCHcgRXJyLLhnzrRJwAYNstZMWlryDJDt7Va1JxLWitm6\nNdVnjoGAuWWAFKuultaulUaPlo4dk06csNfjx6WTJ611U1IirVtHsKPnqNyBfuC222yT7Ebrjh32\ndUWFTQ0MXCjCHehnpk8n0HHpaMsAgA8R7gDgQ4Q7APgQ4Q4APkS4A4APEe4A4EOEOwD4EOEOAD5E\nuAOADxHuAOBDhDsA+BDhDgA+RLgDgA8R7gDgQ0z5CwCX4I03pM2bpQ8/lA4etJ+VlkrXXivNnSvN\nm5ea8yLcAeAi7N4t3XSTLYHoOLZs4ogRthzi0aPSSy9Jr7wi5ebaClrl5X17foQ7AFygpUulp5+W\nioutSs/NtcXNJamtzRYzD4elpiaptlaaNk168EFp1aq+O0fCHQAuwPLl0nPPSVdfLY0ZIxUVWbin\npdmC5m6wNzZKDQ1SVpZUVyc99ZQtgF5V5R3r5MmTOnz4sHJyclR+CaV9V8fhhioA9NCLL0pr1khX\nXSVNmiRNnSpdeaX0la/Y96Wl0tixttj5yJHS8OFSYaFUUGB/BJ54whY6d23ZskUzZszQ97///fN+\n9vvvv68f/OAH2rhx41n7ujoOlTsA9NA991jFXlIilZXZNmqUVefxuLVh6uul9HQpeKp0TiSk9nZr\n17S2SosXS3feafsKCgqSXruze/du3XrrrWpsbNR3v/tdLViwIGl/V8ch3AGgB+6+26rw4mLbxo6V\nJkywLRj0Ql1KDvNYzFo10ajddG1utj8SL7zQs3A/ceKEbrnlFkUikW7f09VxaMsAQA9s3GjhPnSo\nNGyYvbotl9xc66dnZ1sVn5UlZWZ2vWVnS6+9Zsc8X7hHo1HNmzdPDQ0NeuyxxyRJgUDgrPdRuQPA\nRTh61Crx/HwL7owMu4Ha1mY3TSWr1h1HCgSskg+F7D3uNmiQbRkZdrP1+PFzh7vjOFq8eLF27dql\nt99+W/F4vNvzo3IHgIuwc6cX1oGA9dFjMQvpmhobDdPcbK2XeNxGzSQSFvaShX3n0A8GpR07pJyc\nHK1atUqVlZVnfeby5cu1YcMGPfPMM5o9e/Y5z6+r41C5AwNYZWWlKisrVdV5fB0ug+9p8ODVchyr\n0KNRC/O6OvveHQYZidi+1lar6tvbLeQTieSjOY60f799vWzZsrM+7dlnn9U//uM/avny5fr2t7/d\nozM88ziEOzCAzZw5U45bHuKyeestaeHChOLx4Olx7A0NVoW3tlq7JZGwqj0SkVpavJCPx72gdyv6\nQMBG2nTl3Xff1X333advfvOb+ru/+7uLPmfCHQDOY8YMyXGCSQ8oZWbavljMwl2yAI/FLNzdgI/F\nLODjcdvvVvMVFWd/zpEjR7Rw4UJ99atf1b/92791efO0pwh3ADiP4mIb5tjU5I1l79yKSU+3aryj\nwxv+GA7b1rlV425ZWTYPzZkOHjyoxsZGDR48WPfcc0/SvqNHj0qyB5YWLVqkqVOn6oc//GG350y4\nA0APLFggvf66NGSIjXgJBq0Kd8M9GLSK3B3f7rZnWlq8gHfHvC9a1PVnpKena8SIEfr000/P2tfa\n2ipJ+t3vfqfm5uZzjp6RpIBDww4AeiQUsidShw+3se55eTZuvfMTqe3t1oLp3J5pbraKv6HBvm5r\nu/DPfvvttzV37lzdd999+pd/+Zfzvp/KHQB6aO1a6d57LeQTCavQ3XAPhew9HR1e9R6NWtXutmhi\nMenll/vmXAl3AOihu+6S9u61qXs7Oiys3XBPO5WmjmPh7lbv0ahXwT/8sHT77X1zroQ7AFyAqioL\n7NWrrcWSl+eFeyBg4d65eg+HLegfekhasaLvzpOeOwBchH37pBtvtGGR7lj39PTkyj0ry27Abtsm\njR/ft+dHuAPAJdi0ydZQ/eAD6cAB+9mECdL119saqjffnJrzItwBwIeYOAwAfIhwBwAfItwBwIcI\ndwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfIhwBwAfItwBwIcI\ndwDwIcIdAHyIcAcAHyLcAcCHCHcA8CHCHQB8iHAHAB8i3AHAhwh3APAhwh0AfChYXV2t7du3p/o8\nLplfrgMAekNw+/btvghFv1wHAPQG2jIA4EOEOwD4EOEOAD6UdujQIR06dEhVVVWpPpdLsn37dpWU\nlKT6NACgX0iTJMdx5DhOqs/lkvnhGgCgN6SVlJRo3LhxWrlyZarP5ZJUVVUR7gBwCj13APAhwh0A\nfIhwBwAfSqusrEzqVbe2tiojI+OSDxwOh5WWlqbMzMxLPlZPnHkdAPBlFpw5c6YqKyu1a9cu/eEf\n/qHy8vJUXl6ud955p9tfGj16tIYMGXJ6u/baa0/vO378uObOnatx48YpJydHixcv7ovrkHsdAIBT\nQyFjsZi+9a1vacmSJdq+fbs2bdqkhQsXavfu3briiiuSfuHYsWNqa2vThx9+ePpnnavz++67T2Vl\nZXrrrbf061//WnPnztX777+v66+/vo8uCQCQJkn/+Z//qdraWi1btkzBYFDz58/XqFGj9Oqrr2rp\n0qVJv7B//35NnjxZU6dOPetg+/bt08aNG3XkyBGFQiFdffXV+o//+A/l5eX1zdUAACSduqH62Wef\n6fd+7/cUDHr3V8vLy7V3796zfmH//v2KRCKaM2eOvvKVr2jJkiX63e9+J0nau3evSktL9c477+jG\nG2/Ud7/7XUWjUZWWlvbR5QAApFPh/tvf/lZDhgxJ2jFu3DiFw+GzfuHQoUNqbm7Wgw8+qCeeeEKH\nDx/WjTfeqFgspn379uno0aN69dVXdf/99ystLU2VlZV6++23++ZqAACSTrVlBg8erBMnTiTtCIfD\nmjx58lm/sHLlyqSnWSdOnKiJEyfqo48+UkdHh9ra2rR27VqNGDFCt912mz7//HOtX79ec+bMubxX\nAgA4LShJEyZM0IEDB5J27Nu3T9OmTTvrF959910dOXLk9PcjR45URkaGRo4cqbFjxyo/P1+FhYWn\n9xcXF6u1tfVynT8AoAtBSfrGN76h2tpabdiwQY7jaN26dfr000/1ta99TZK0fv16/c///I8ku/m6\nYMECtbe3q7W1VT/60Y80ZcoUTZ48WTfddJMSiYSeffZZJRIJ7d+/X2+88YZuvvnm1F0hAHwZOaf8\n9Kc/dbKyspzi4mInJyfH2bBhg7vLGTNmjPPjH//YcRzHqaurcyorK50xY8Y4RUVFzsiRI51t27ad\nfu+bb77pDB061CkuLnaCwaDzwAMPOIlEwgEA9J2A43iPdba2tuqzzz7TlClTNGjQoHP+UaitrVV9\nfb1KS0sVCoWS9rnHGT16tAoKCi7PXyUAQLeSwh0A4A9MHAYAPkS4A4APEe4A4EOEOwD4EOEOAD5E\nuAOADxHuAOBDhDsA+BDhDgA+RLgDgA8R7gDgQ4Q7APgQ4Q4APkS4A4APEe4A4EOEOwD4EOEOAD5E\nuAOADxHuAOBDhDsA+BDhDgA+RLgDgA8R7gDgQ4Q7APgQ4Q4APkS4A4APEe4A4EOEOwD40P8HNlIz\nrMj5u9IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ete3 import Tree, TreeStyle , NodeStyle, AttrFace, faces\n",
    "\n",
    "t = Tree(str(T)+\";\")\n",
    "\n",
    "ts = TreeStyle()\n",
    "#ts.rotation = 90\n",
    "ts.scale = 90\n",
    "ts.branch_vertical_margin = 50 # 10 pixels between adjacent branches\n",
    "ts.show_leaf_name = False #\n",
    "\n",
    "nstyle = NodeStyle()\n",
    "nstyle[\"shape\"] = \"sphere\"\n",
    "nstyle[\"size\"] = 20\n",
    "nstyle[\"fgcolor\"] = \"blue\"\n",
    "\n",
    "for n in t.traverse():\n",
    "    n.set_style(nstyle)\n",
    "\n",
    "for node in t.iter_leaves():\n",
    "    node.add_face(AttrFace(\"name\", fsize=20), column=0)  \n",
    "    \n",
    "t.render(\"%%inline\", tree_style = ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the tree $T$ is computed, one can compute the code $C_{i}$\n",
    "associated to each symbol $i$. This requires to perform a deep first\n",
    "search in the tree and stop at each node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "codes = {}\n",
    "\n",
    "def huffman_gencode(T,codes,c):\n",
    "    if type(T) == str: #test if T is a leaf\n",
    "        codes[T] = c\n",
    "    else:\n",
    "        huffman_gencode(T[0],codes, c + \"0\")\n",
    "        huffman_gencode(T[1],codes, c + \"1\")\n",
    "\n",
    "huffman_gencode(T,codes,\"\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Code of token 1: 101\n",
      "Code of token 3: 110\n",
      "Code of token 0: 100\n",
      "Code of token 2: 0\n",
      "Code of token 4: 111\n"
     ]
    }
   ],
   "source": [
    "for e in codes:\n",
    "    print(\"Code of token \" + e + \": \" + codes[e])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We draw a vector $x$ according to the distribution $h$.\n",
    "\n",
    "Size of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 1024"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import random\n",
    "\n",
    "def rand_discr(p, m = 1):\n",
    "    \"\"\"\n",
    "        rand_discr - discrete random generator \n",
    "        \n",
    "          y = rand_discr(p, n);\n",
    "        \n",
    "          y is a random vector of length n drawn from \n",
    "          a variable X such that\n",
    "              p(i) = Prob( X=i )\n",
    "        \n",
    "          Copyright (c) 2004 Gabriel Peyré\n",
    "    \"\"\"\n",
    "\n",
    "    # makes sure it sums to 1\n",
    "    p = p/np.sum(p)\n",
    "    \n",
    "    n = len(p)\n",
    "    coin = random.rand(m)\n",
    "    cumprob = np.append(0,+ np.cumsum(p))\n",
    "    sample = np.zeros(m)\n",
    "    \n",
    "    for j in range(n):\n",
    "        ind = [(coin > cumprob[j]) & (coin <= cumprob[j+1])]\n",
    "        sample[ind] = j\n",
    "        \n",
    "    return sample\n",
    "\n",
    "x = rand_discr(h, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 1__\n",
    "\n",
    "Implement the coding of the vector $x$ to obtain a binary vector $y$, which corresponds to replacing each sybmol $x(i)$ by the code $C_{x(i)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "run -i nt_solutions/coding_2_entropic/exo1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the length of the code with the entropy bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy bound = 2197.95\n",
      "Huffman code  = 2228.00\n"
     ]
    }
   ],
   "source": [
    "e = - np.sum(h*np.log2([max(e,1e-20) for e in h]))\n",
    "print(\"Entropy bound = %.2f\" %(n*e))\n",
    "print(\"Huffman code  = %.2f\" %len(y)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Decoding is more complicated, since it requires to iteratively parse the tree $T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial empty decoded stream."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1 = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform decoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T0 = T\n",
    "for e in y:\n",
    "    if e == '0':\n",
    "            T0 = T0[0]\n",
    "    else:\n",
    "            T0 = T0[1]\n",
    "    if type(T0) == str:\n",
    "        i = i+1\n",
    "        x1 += T0\n",
    "        T0 = T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We test if the decoding is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error (should be zero) : 0.000000 \n"
     ]
    }
   ],
   "source": [
    "from numpy import linalg\n",
    "\n",
    "err = linalg.norm(np.subtract(x,[float(e) for e in x1]))\n",
    "print(\"Error (should be zero) : %f \" %err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huffman Block Coding\n",
    "--------------------\n",
    "A Huffman coder is inefficient because it can distribute only an integer\n",
    "number of bit per symbol. In particular, distribution where one of the\n",
    "symbol has a large probability are not well coded using a Huffman code.\n",
    "This can be aleviated by replacing the set of $m$ symbols by $m^q$\n",
    "symbols obtained by packing the symbols by blocks of $q$ (here we use $m=2$ for a binary alphabet). This breaks\n",
    "symbols with large probability into many symbols with smaller proablity,\n",
    "thus approaching the Shannon entropy bound.\n",
    "\n",
    "\n",
    "Generate a binary vector with a high probability of having 1, so that the\n",
    "Huffman code is not very efficient (far from Shanon bound).\n",
    "\n",
    "Proability of having 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t = .12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Probability distriution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h = [t, 1-t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import random\n",
    "\n",
    "n = 4096*2\n",
    "x = (random.rand(n) > t) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For block of length $q=3$, create a new vector by coding each block\n",
    "with an integer in $\\{1,...,m^q=2^3\\}$. The new length of the vector is\n",
    "$n_1/q$ where $n_1=\\lceil n/q\\rceil q$.\n",
    "\n",
    "Block size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "q = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximum token value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n1 = (n//q+1)*q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1 = np.zeros(n1)\n",
    "x1[:len(x)] = x\n",
    "x1[len(x):] = 1\n",
    "x1 = x1 - 1\n",
    "x2 = []\n",
    "\n",
    "for i in range(0,n1,q):\n",
    "    mult = [m**j for j in range(q)]\n",
    "    x2.append(sum(x1[i:i+q]*mult))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate the probability table $H$ of $x_1$ that represents the probability\n",
    "of each new block symbols in $\\{1,...,m^q\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H = h\n",
    "for i in range(q-1):\n",
    "    Hold = H\n",
    "    H = []\n",
    "    for j in range(len(h)):\n",
    "        H = H + [e*h[j] for e in Hold]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simpler way to compute this block-histogram is to use the Kronecker product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H = h\n",
    "for i in range(1,q):\n",
    "    H = np.kron(H, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise 2__\n",
    "\n",
    "For various values of block size $k$, Perform the Huffman coding and compute the length of the code.\n",
    "Compare with the entropy lower bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy bound = 0.529361\n",
      "---\n",
      "Huffman(block size = 1) = 1.000122\n",
      "Huffman(block size = 2) = 0.675171\n",
      "Huffman(block size = 3) = 0.572876\n",
      "Huffman(block size = 4) = 0.544434\n",
      "Huffman(block size = 5) = 0.537476\n",
      "Huffman(block size = 6) = 0.532593\n",
      "Huffman(block size = 7) = 0.543579\n",
      "Huffman(block size = 8) = 0.539795\n",
      "Huffman(block size = 9) = 0.537231\n",
      "Huffman(block size = 10) = 0.537109\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFRCAYAAADgqHO9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8HFWd///XOzcJWSEhEJaYEFYFMQNhSQQCV3EwoLIp\ng8giiA7iKPgddx0lOs6o89NxGcVRWWSRRREBFZRxuYQlCRB2SIAAgQDKmoRA9uTz++PUTTrNXbrv\n7e7q5f18POpxq6uqqz7dnfSnzzl1zlFEYGZm1goG5B2AmZlZrTjpmZlZy3DSMzOzluGkZ2ZmLcNJ\nz8zMWoaTnpmZtQwnPasqSWdKek7SK5JGSzpQ0qOSlkk6MufY1kvaqZt9p0q6uZ/nb5e0qD/n6MM1\nN4k7e58n1jKGvpJ0oaSXJc3OOxZrXk561i1JCyUdWrSt5GQgaRDwHeDQiNg8IhYDXwN+EBEjI+K6\nykdthbL3eWHecfRG0jTgHcD2ETE173iseTnpWU8iW/pqW2AIMK9g2wTgof4EZc1FUhuwA7AwIlb2\n4fkDKx+VNSsnPSvXJkmwuIpQ0s8l/bukXYH52eYlkv4saQGwE/DbrLpzsKSO7Phbs6q46yRtJekX\nkpZKul3SDgXn/76kp7J9d0o6qGDfDEm/lHRRdv4HJO3Ty+t5l6THJL0g6b8kqauDJB0g6Q5JS7KY\n3lqwb8usau6ZrHruN92c4yxJD0ravpv9H5H0UBb7g5L2zrbvnr1Pi7PX9J6C54zJ3rOlkuYAOxed\nc8Pnk302P5L0u+was4s+u8MkPZy9xh9JuknS6d3EOkPSVZKuyM41V9Kkgv3bS/q1pOclPS7pE108\n9xJJS4HTgZ8Bb83+DZxT8H48KuklSddK2q7odX1M0qPAw5IOkfS0pM9k13xW0tGSjpD0SHaOzxc8\nf39Js7L39FlJ/5PVTBSe/4zsuYsl/bDEz6rb1211IiK8eOlyAZ4gVU0WbjsVuLng8Xpgp4LHFwJf\ny9Z3yPYPKDrn2wsedwCPADsCmwMPAo8CbwfagIuACwqOPxEYTfrB9q/A34DB2b4ZwApgOiDgP4FZ\nPby+9cCfgVHAeOBh4PTi1wlsCSzOrj0AeD/wMjA62/974HJgC2AgMC3b3g4syta/AtwJjOkmluOA\np4F9ssc7k0rFg4AFwOezc78NeAXYLTvuimwZCrw5O8fMrj4f4OfAi8C+2Xt7KXB5tm8rYClwdPYa\nzwJWAx/qJt4Z2f5js3N9Cng8Wx8AzAX+LYt5R+Ax4LCi5x6ZPR4CfJBN/129HXgB2AsYDPwAuKno\ndf0x++w2y97rNdk124APZ6/1F8BwYA9gObBD9vzJwP5ZrDuQah/OLjr/daR/k+OB54F39vJZ9fi6\nvdTHknsAXup3ARYCy0hf+J3La919qWaPLwT+PVufSO9J76/AFwoefxv4fcHjdwN39xDjy8BbsvUZ\nwI0F+/YAlvfw3PWFX0jAmcCfsvVT2Zj0TgZmFz33tuyLejtgHbBFF+dvz74c/xuYCYzsIZY/Ap/o\nYvs04G9F2y4Dzsm+3FeTJcBs33/QzY+S7LP5acG+w4F52fopwK1F13mKnpPebQWPBTwLHARMAZ4s\nOv4LZD9esud2FO3f8H5nj88HvlnweHj2WicUvK72ovd6OaDs8cjsmP0KjrkTOKqb1/NJ4Oqi9+2A\ngsdXAp/t5bPq8XV7qY/FdeHWkyB9Sfylc4OkD5J+RVfScwXrK0m/qgsfjyi4/qeBDwHbZ/FtTiql\ndHWu5cAQSQMiYn031y68u/Kp7LzFts/2FXoy2/4G4OWIWNrN+UeR3q/3R8Sybo4hO89j3Vy7+A7Q\nzmtvRSpRFL+GnhS+PyvY+N5uT0rQhYofF9uwPyJC0tNs/Fy2l7S44Ng2UuIv9dzbkZJU5/lfk/QS\nMI6Nr7H4fXkpskxDem3w+tc7HEDSbqQfI/sAw0jv451s6u8F68vZ+F5191ntQO+v23LmNj0rV3Gb\n13LSl0an7ejfzS/dPlfpDr/PAMdFxKiIGE2qkuuyHa5EE4rWn+nimGdIX2iFdsi2LwK2lLRFN+df\nTCqtXijpgB7iWATs0sX2Z4HxRW2Nndd+AVjbxWvoi2dJX+YAZNd7Q/eHA6nar/P4Adnxne/JExEx\numDZPCLenR1eyg1Sz5JqCjrPPxwYw6afT3/+nf2YVKW5S0RsAXyJ0r8Pu/usnqLn1211wEnP+use\n4ERJbZKmAwf34RzqZr3YSNKX/ItKN8F8hVTS649PSxolaTypHevKLo65AdhN0gmSBko6HngT8LuI\n+Hu2/9zsPIMkbfIeRMRMUnvg1ZL26yaO87JYJivZRdIEYDbph8Vns3O3k5LoFVnp9WpghqShkvYg\nVbl2p6f39nrgLZKOUrob8l9Id9/2ZB9Jx2THf5JUKp8N3AEsk/TZLK42SXtK2reEODpdDpwm6R8k\nbUZqn50dEb2VZEs1glR1v1zSm0hV2z0RG+Pu7rO6nZ5ft9UBJz0rV/Gv9LOB95BKNB8Aiu9cLOXX\neBStFz+n8/EfsuURUnvjCjatzuvpud25lnTzwd3A70htSZucKyJeIiWaT5Fujvg08O6IeDk79mTS\nTRTzSdVpZxVfPyL+RKqW/a2kvYqDiIirSO1xl5FuVLmadKPMGtL7ezipZPdD4OSIeCR76sdJX+B/\nBy7IluL3s3C9y/cnIl4k3aDxX9lr3J1U3bfqde/YxuddCxxPalc9ETg2ItZFxLrs/dqLdHPLC8BP\n2fgDpbs4NmyLiD8DXwZ+TSr17Ui6gair19Xdtp4++0+T/r2+ksV2Bd2/b5vE18NntZ6eX7fVAW2s\nAq/BxVJJ4Hukeu7zIuJbRfvbSf+RHs82/Toivl6zAM0M2FBduQj4QETc1MX+c0hVgyfXPDizfqjZ\njSxKHVB/SBp14RngDknXRcS8okNviohch6cya0WSDiNV0a0gtZ1Cqq7s8vCaBGVWYbWs3twfWBAR\nC7MqmyuAo7o4zv+ZzPLxVlKfwBeAdwFHR0RP1Zu1qyYyq5BadlkYx6a3GD9N6tdSKIADJN1LKg1+\nOiI8ZJVZDUTEV4GvlnGsWcOpZdIr5VfhXcD4iFgu6XDgGmC36oZlZmatopZJ7xkK+vVk65t0UC3s\nvBsRN0g6V9KWBXfJASDJ1SpmZvY6EdFjE1kt2/TuBHaVNFHSYNKtzptMLSNpm85OuJL2J91d+vLr\nT9U4w6edc845ucfQrPE2UqyNFq9jdbyNFmtEaWWhmpX0ImKtpI+Txq1rA86PiHmSzsj2/wR4H3Cm\npLWkDrnv7/aEZmZmZarp2JsRcQNp9IrCbT8pWP8R8KNaxmRmZq3DI7JUWXt7e94hlKWR4m2kWKGx\n4nWs1dNI8TZSrKWq6YgslSIpGjFuMzOrHklEHd3IYmZmlisnPTMzaxlOemZm1jKc9MzMrGU46ZmZ\nWctw0jMzs5bhpGdmZi3DSc/MzFqGk56ZmbUMJz0zM2sZDZv0Fi/OOwIzM2s0DZv0br017wjMzKzR\nNGzSu/nmvCMwM7NG46RnZmYto2GnFho0KFiyBIYNyzsaMzOrB009tdCaNTBnTt5RmJlZI2nYpAeu\n4jQzs/I46ZmZWcto2DY9CIYPT/31Bg3KOyIzM8tbU7fp7borvPYa3H133pGYmVmjaNikd/DB6a+r\nOM3MrFQNm/SmTUt/nfTMzKxUDdum99hjwc47w5ZbwgsvwICGTd9mZlYJTd2mt+OOsP328PLLMG9e\n3tGYmVkjaNikJ7ldz8zMytOwSQ/crmdmZuVpiqQ3cyY0YNOkmZnVWMPeyBIRrF8PW22VOqg/8QRM\nnJh3ZGZmlpemvpEF0h2bBx2U1l3FaWZmvWnopAdu1zMzs9I1TdKbOTPfOMzMrP41dJsewOrVMHo0\nLF8Ozz0HY8fmHJyZmeWi6dv0AAYPhqlT0/ott+Qbi5mZ1beGT3rgKk4zMytNUyU938xiZmY9afg2\nPUjz6o0aBevXpz57m2+eY3BmZpaLlmjTAxg+HPbZJyW9WbPyjsbMzOpVUyQ9cLuemZn1rumSntv1\nzMysO03RpgdpXr0xY2CzzWDp0vTXzMxaR8u06UGaQX3PPWHVKrjjjryjMTOzetQ0SQ/crmdmZj1r\nqqTnmdTNzKwnTdOmB/DMM/CGN8DIkam/XltbDsGZmVkuWqpND2DcONhxR1i2DO69N+9ozMys3tQ0\n6UmaLmm+pEclfa6H4/aTtFbSseVew10XzMysOzVLepLagB8C04E9gBMk7d7Ncd8C/gD0WEztitv1\nzMysO7Us6e0PLIiIhRGxBrgCOKqL4z4BXAW80JeLFJb0GrC50szMqqiWSW8csKjg8dPZtg0kjSMl\nwh9nm8pOW7vumiaSff55eOSRvoZqZmbNqJZJr5QE9j3g89mtmaIP1ZuSqzjNzKxrA2t4rWeA8QWP\nx5NKe4X2Aa6QBLAVcLikNRFxXfHJZsyYsWG9vb2d9vb2DY+nTYOrrkpJ78MfrlT4ZmZWTzo6Oujo\n6CjrOTXrpydpIPAwcCjwLHA7cEJEzOvm+AuB30bE1V3s67KfXqe774bJk1P3hccfr0j4ZmZW5+qq\nn15ErAU+DvwReAi4MiLmSTpD0hmVvNakSWki2SeegKeLy5JmZtaymmpElkJHHAE33ACXXQYnnFCj\nwMzMLDd1VdKrNXdSNzOzYk2f9DzjgpmZdeq1elPSpwoednYl6FwnIv67OqH1GFOv1ZurVsEWW6S/\nL76YJpg1M7PmVanqzZHACFJ3gjOB7Umdyj8KTO5vkNWy2WYwZUpav/XWfGMxM7P60GvSi4gZEfFV\nUr+6yRHxqYj4V1IS3KHaAfaH2/XMzKxQOW16Y4E1BY/XZNvqltv1zMysUDkjslwM3C7palK73tHA\nRVWJqkIOOAAGDIC77oLXXoPhw/OOyMzM8lRySS8i/gM4DVgCvAycGhH/Wa3AKmHkSNh7b1i7FmbP\nzjsaMzPLW1ldFiJibkR8LyK+HxF3VyuoSnK7npmZdSq5elPSEOC9wMSC50VEfK0KcVXMtGnwve+5\nXc/MzMpr07uWVLU5F1hZnXAqr7OkN3s2rF4NgwfnG4+ZmeWn5LE3JT0QEXtWOZ6SlNI5vdDuu8P8\n+TBrFkydWsXAzMwsN5Uee/M2SZP6GVMu3HXBzMygvKQ3DZgr6RFJ92fLfdUKrJJ8M4uZmUF51Ztd\njr4SEU9WNKLSYimrevPJJ2HiRBg9Oo3DOaBph9k2M2tdpVRvlpP0zqHrAadrfvdmuUkPYMIEWLQI\n7rsP3vKWKgVmZma5qXSb3mvZ8iqwDjiC1H2hIbhdz8zMyhmR5dsR8Z1s+TpwCLBz9UKrrIMPTn/d\nrmdm1rr607o1nDTFUEMovJmlzJpRMzNrEuWMyHJ/wcMBpBkW6no0lkK7754mkn32WXjiCdhpp7wj\nMjOzWitnRJb3ZH8DWAs8HxFreji+rkhw0EFw7bWpXc9Jz8ys9ZTTprcQGAUcCRwD7FGlmKrG7Xpm\nZq2t5KQn6WzgUmBrYBvgUklnVSuwanAndTOz1lZOP737gakR8Vr2eDgwOyJq3uutL/30IM2rN2pU\nmlD2b3+DbbetQnBmZpaLSvfTA1jfzXpDGDgwzaYOLu2ZmbWicpLehcAcSTMkfRWYDVxQnbCqx1Wc\nZmatq+S7NyPivyXdBByYbTq1UWZPL+SkZ2bWusrppzcUaAcOInVbaJM0LyIaZkJZgClTYNAguPde\nWLIktfGZmVlrKKd682JSN4UfAD8E3gxcUo2gqmnoUNhvvzQqy2235R2NmZnVUjlJ780RcXpE/DUi\n/hIRHyYlvobjKk4zs9ZUTtK7S9JbOx9ImgrMrXxI1eekZ2bWmnrtp1cw5uZA4I3AIlKb3gTg4YjY\nvaoRdh1Tn/rpdVqyBLbcMnVhWLo0VXmamVljK6WfXik3srynh30NOV/BqFEwaVK6meX22+GQQ/KO\nyMzMaqHXpJeNudl0pk1LSe/mm530zMxaRX/m02tobtczM2s9JY+9WU/626YHaezN7beHESNg8eLU\nvmdmZo2romNvSjpL0uj+h1UfttsOdtkFXn0V7rkn72jMzKwWyqne3Aa4Q9IvJU2X1GM2bQSu4jQz\nay3lTCL7JWA30iDTpwKPSvpPSTtXKbaq60x6M2fmG4eZmdVGWTeyRMR64O/Ac8A6YDRwlaT/rwqx\nVV3nTOq33JKGJTMzs+ZWziSyZwOnAC8B5wG/iYg1kgYAj0ZEzUp8lbiRBVKiGzcu3dTy0EOwe827\n2ZuZWaVUehLZLYFjI+KwiPhlRKyBDaW/njqw1y3J7XpmZq2knDa9cyLiyW72PVS5kGrL7XpmZq2j\n3Pn0PsbG+fRuBn7caPPpFets13NJz8ys+ZXTpvcr4BXgUkDAB4AtIuK46oXXbSwVadMDWL8exoxJ\ng1A/+SRMmFCR05qZWY1Vuk2vaebTKzRgABx4YFp3ac/MrLm15Hx6xTqrON2uZ2bW3MpJevsCt0p6\nUtJC4DZgX0n3S7qvlBNkI7nMl/SopM91sf8oSfdKulvSXElvLyO+PvMdnGZmraGcNr2J2WrnEzap\nN+1tCiJJbcDDwDuAZ4A7gBMiYl7BMcMj4rVs/S2kvoC7dHGuirXpAaxenebYW7ECnn8ett66Yqc2\nM7MaqWibXpbURgFHkvrlbRERCzuXEk6xP7AgO34NcAVwVNE1Xit4OAJ4sdT4+mPwYJgyJa3fckst\nrmhmZnkoZ5aFs0l3bm5NGnz6UklnlXGtccCigsdPZ9uKr3O0pHnADUA55+8Xd10wM2t+5cwi92Fg\nSkH14zeB2cAPSnx+SfWREXENcI2kacAlwBu7Om7GjBkb1tvb22lvby8xjK65Xc/MrLF0dHTQ0dFR\n1nPKadO7H9g/IlZkj4cCt0fEW0p8/lRgRkRMzx5/AVgfEd/q4TmPZdd8qWh7Rdv0IM2rN2pUWl+8\nGEaOrOjpzcysyirdT+9CYI6kGZK+SirlXVDG8+8EdpU0UdJg4HjguqKAd+6cp0/SZIDihFctI0bA\n5Mmwbh3MmlWLK5qZWa2VlPSyRHQVcBqwmDTTwqkR8d1SLxQRa4GPA38EHgKujIh5ks6QdEZ22HuB\n+yXdDXwfeH/Jr6QC3K5nZtbcSqrezJLe/RGxZ/VD6l01qjcBrr0Wjj4aDjkEyqwmNjOznFWsejPL\nMHMl7V+RyOrUQQelv3PmwKpV+cZiZmaVV06b3lRglqTHs1FYSh6JpVGMGQNvfjOsXAl33pl3NGZm\nVmnldFk4jKJRWCixG0IjmTYNHnwwtet1DkRtZmbNoZyS3scKR2DJRmH5WJXiyo3765mZNa9ykt5h\nXWw7olKB1IvOpHfLLan7gpmZNY9ek56kM7OO6W8saMu7P5tpoana9ADGj4eJE+GVV+D++/OOxszM\nKqmUNr3LSONgfhMonA5oWUS8XJWocjZtGixcmKo499or72jMzKxSei3pRcTSrP3uVOAg4KRs/ROS\nvlLN4PLidj0zs+ZUzt2b1wJLSLOlr6xOOPWhcCb1CFCPXR3NzKxRlDPg9APNPiJLpwjYdts0oewj\nj8Cuu1btUmZmViGVHnD6NkmT+hlTQ5A2js7iKk4zs+ZRTtKbRhqK7JFmHZGlkNv1zMyaTzlteod3\nsa3pRmTpVNiuZ2ZmzaGUfnqfBcju4NyvaESWM3p6biP7h39IE8k+/jg8+2ze0ZiZWSWUUr15QsH6\nF4v2dVX6awptbXDAAWndVZxmZs2hnDa9ltPZrucqTjOz5uCk1wPPpG5m1lx67acnaR2wPHs4FFhR\nsHtoRJRzM0xFVLufXqeVK2GLLWDNGnjpJRg9uuqXNDOzPqpIP72IaIuIkdkysGB9ZB4Jr5aGDIH9\n90+d1W+9Ne9ozMysv1y92Qt3XTAzax5Oer1wJ3Uzs+ZR8tib9aRWbXqQ5tUbPRoGDIClS2HYsJpc\n1szMylTpsTdb0uabp47qa9fCnDl5R2NmZv3Rp6Qn6X2VDqSeuV3PzKw59LWkN7yiUdQ5t+uZmTUH\nV2+WoHOaoVmzUp89MzNrTE56JdhmG9htN1i+HO66K+9ozMysr5z0SuQhyczMGl9fk17j9XPoJ7fr\nmZk1vj7105M0LCKW935kddSyn16nJ56AnXaCLbeEF15I/fbMzKx+VK2fXp4JLy8TJ8Ib3gAvvwwP\nPZR3NGZm1hcur5RIchWnmVmjc9Irg5OemVljKynpKRlf7WDqXeFM6g04ZKmZWcsrp6R3Q9WiaBB7\n7JFuZHnmGVi4MO9ozMysXCUlvexWybmS9q9yPHVtwICNo7O4itPMrPGUU9KbCsyS9Lik+7PlvmoF\nVq/crmdm1rgGlnHsO7O/na1ZPfaFaFaF7XpmZtZYyuqcLmkvYBop8d0cEfdWK7Be4qh55/ROa9bA\nqFFpHM6//z2Ny2lmZvmraOd0SWcDlwJbA9sAl0o6q38hNp5Bg+Ctb03rt9ySbyxmZlaectr0PgxM\niYivRMSXSW18H6lOWPXN7XpmZo2p3M7p67tZbymeSd3MrDGVcyPLhcAcSVeTbmI5GrigKlHVuSlT\nUjXnvffCK6/A5pvnHZGZmZWi5BFZgKuA04DFwEvAqRHx3SrGVreGDYN99oH16+G22/KOxszMSlVO\n9eb1ETE3Ir4fET+IiLurFlUDcLuemVnj8YgsfeR2PTOzxlNyPz1JDwO7AE8Cr2WbIyImVSm2nmLJ\nrZ9ep8WLYcyY1La3dCkMGZJrOGZmLa9i/fSyNr2PADsDbwfeky1H9iGo6ZLmS3pU0ue62H+ipHsl\n3SfpVkk1T6qlGD0a9twTVq+GO+7IOxozMytFOW1650bEwuKlnItJagN+CEwH9gBOkLR70WGPAwdn\nJch/B35azjVqyVWcZmaNpdZtevsDC7KEuQa4Ajiq6FqzImJp9nAO8IZ+XrNqfDOLmVljKaef3lTg\nJEn9adMbBywqePw0MKWH408Hri/j/DXVmfRuuw3WrYO2tnzjMTOznvV1loW+zrBQ8t0nkt4GfAg4\nsKv9M2bM2LDe3t5Oe3t7H0Pqu+23h512gscfTx3VJ0+ueQhmZi2ro6ODjo6Osp5Tzt2bA4ATgR0j\n4muSJgDbRsTtJV9MmgrMiIjp2eMvAOsj4ltFx00CrgamR8SCLs6T+92bnU47DX7+c/jud+GTn8w7\nGjOz1lXRWRaAc4G3Ah/IHr+abSvHncCukiZKGgwcD1xXeECWTK8GTuoq4dUbt+uZmTWOcqo3p0TE\n3pLuBoiIlyUNKudiEbFW0seBPwJtwPkRMU/SGdn+nwBfAUYDP049JVgTEXXbKb4w6UWAWnJqXTOz\nxlBO9eYc4ADgziz5bQ3cGBF7VzPAbmKpm+rNCNhuO3juOZg/H974xrwjMjNrTZWu3vwf4DfAWEn/\nCdwKfKMf8TUFyf31zMwaRclJLyIuBT5HSnTPAkdFxC+rFVgjcbuemVljKKdNj4iYB8yrUiwNy0nP\nzKwxlNymV0/qqU0PUsf0MWPSwNNPPQXjx+cdkZlZ66l0m551o60NDsy60Lu0Z2ZWv5z0KsRVnGZm\n9a/kNj1JQ4D3AhMLnhcR8bUqxNVwnPTMzOpfOf30/ggsAeYC6zq3R8R3qhNaj7HUVZsewKpVMGoU\nrFwJL76Y2vjMzKx2SmnTK+fuzXER8c7eD2tNm20GU6bATTfBLbfAUUf1/hwzM6utctr0bqvXWczr\nhas4zczqWzklvWnAaZKeAFZl28qdT6+pOemZmdW3ctr0JnaxOSLiyUoGVIp6bNMDWLYMRo9O60uW\nwIgR+cZjZtZKKtpPLyIWAkuBscCEbNmhPwE2m5EjYe+9U2f12bPzjsbMzIqVnPQkfQSYCdwIfJU0\nPdCM6oTVuFzFaWZWv8q5keVsYH9gYUS8DdibVPKzAp5xwcysfpWT9FZGxApIHdUjYj7g2eOKHHRQ\n+jt7NqxenW8sZma2qXKS3iJJo4FrgP+TdB2wsCpRNbCttoLdd0+d1OfOzTsaMzMrVM6NLMdExOKI\nmAF8GTgPOLpagTUyt+uZmdWnPg04HREdEXFdRLgCrwtu1zMzq0/l9NPbD/girx9wuuad0+u1n16n\np56CHXZIY3G+9BIM8FwWZmZVV0o/vXKS3iPAp4EHgPWd27P+ezVV70kPUtJ76im4916Y5DFrzMyq\nrtKTyL6QVWk+HhELO5f+hdi8Oqs43a5nZlY/ykl6X5V0vqQTJL03W46tWmQNrvNmFrfrmZnVj3IG\nnP4gqV/eQAqqN4GrKxpRkyi8gzMC1GOB28zMaqGcNr2HgTfVQ2NaI7TpRcDYsWlC2QULYOed847I\nzKy5VbpN7zZgj/6F1DokaG9P6+eck5KgmZnlq5yS3nxgZyD3+fQaoaQHMG8e7LsvLF8OP/kJ/PM/\n5x2RmVnzqkiXBUmXRMTJkj5JGoJsE+6y0LNLL4WTT4bNNoNZs9LUQ2ZmVnmVqt7cR9L2wIeAV7pY\nrAcnnZRKeKtWwXHHwVLPS2FmlptSSnpnAWcCOwHPFu2OiNipSrH1FFPDlPQgDT59wAFw991wzDHw\n61/7bk4zs0qr9IgsP46IMysSWT81WtIDeOwxmDwZXnkFvvtd+OQn847IzKy5VKR6U9Il2erDFYmq\nRe28M1x4YVr/zGdS+56ZmdVWWW16krYsXqodYDM59lj4f/8P1q6Ff/qn1IfPzMxqp79tekTEjtUJ\nrceYGq56s9OaNXDIIamkN306/P73noXBzKwS3KZXpxYtSl0XXnoJvv51+NKX8o7IzKzxVTrpnQO8\n7uCI+Frfwuu7Rk96AH/4Axx+eCrl/elP8La35R2RmVljq/QwZK8VLOuBI0gTylofTJ+eSnjr18MJ\nJ8Df/pZ3RGZmza/kkt7rnihtBtwYEYdUNqSSrt3wJT2AdevgH/8R/vrX1M73pz/BwHLmvTAzsw0q\nXdIrNhwY14/nt7y2NrjsMth2W7jppjQwtZmZVU85bXr3FzwcAIwFvhYR/1ONwHqJpSlKep06OuDQ\nQ1NV5+9/D0cckXdEZmaNp9I3suxQ8HAt8FxErO1HfH3WbEkP4BvfgC9+EbbcMg1XNmFC3hGZmTWW\nSs2yUFiqTu3wAAARkElEQVSSC6DwhBERZ/U9xL5pxqS3fj285z1w/fUwZQrMnAmDB+cdlZlZ46hU\nm95c4M7s71EF652LVcCAAXDxxamEN2cOfPazeUdkZtZ8yrp7U9LdEZH7jHDNWNLrNGcOTJuWRm75\n1a/gfe/LOyIzs8ZQ7bs3rQqmTIFvfzutf+hDsGBBvvGYmTUTJ7069IlPwHvfC8uWpZLeihV5R2Rm\n1hxKmVroVUnLJC0D3tK5ni1lzZwuabqk+ZIelfS5Lva/SdIsSSslfaqcczcTCc4/H3bZBe69F86q\n+a1CZmbNqc8jspR9IamNNCffO4BngDuAEyJiXsExWwM7AEcDiyPiO92cq2nb9Ardcw9MnQqrVsFF\nF8Epp+QdkZlZ/aq3Nr39gQURsTAi1gBXkO4G3SAiXoiIO4E1NYyrbu21F/zwh2n9zDPhwQfzjcfM\nrNHVMumNAxYVPH4aD2PWq9NPTyW85ctT+96rr+YdkZlZ46pl0mv++sgqkODcc+HNb4b58+GMM6AF\nanbNzKqilmP6PwOML3g8nlTa65MZM2ZsWG9vb6e9vb2vp6p7w4enPnv77ZcGqJ42DT760byjMjPL\nV0dHBx0dHWU9p5Y3sgwk3chyKPAscDtFN7IUHDsDWNbqN7IUu+wyOPHENDzZrFkweXLeEZmZ1Y+K\nDjhdCZIOB74HtAHnR8Q3JJ0BEBE/kbQt6a7OzUkT1S4D9oiIV4vO05JJD9INLf/7v7DjjnDXXTBq\nVN4RmZnVh7pLepXSyklv5Uo48MCU8I4+Gq6+OrX7mZm1unrrsmAVMGRIat/bYgu45hr47nfzjsjM\nrHG4pNegrrkGjjkGBg5Ms64fcEDeEZmZ5cslvSZ29NHwqU/B2rXwT/8EL7yQd0RmZvXPJb0GtmYN\ntLfDbbfBYYfBDTekefnMzFqRS3pNbtAguPJK2GoruPFG+I//yDsiM7P65pJeE7jxRpg+Pa3/3//B\noYfmG4+ZWR5c0msRhx0G//ZvaXiyD3wAnn0274jMzOqTS3pNYt06eOc74c9/hoMPTn8H1nKQOTOz\nnLmk10La2uAXv4DttoOZM+HLX847IjOz+uOSXpOZORPe/vZU8vvtb+Hd7847IjOz2nBJrwUdfPDG\nuzhPOQUWLsw1HDOzuuKSXhNavx6OOgp+97s0HdEtt6SZGczMmplLei1qwAC46CLYYQe44w749Kfz\njsjMrD64pNfEbr8dDjoojdzyy1/CccflHZGZWfW4pNfi9t8fvpNNw3v66fDII/nGY2aWN5f0mlwE\nHH98mo5o0iSYPRuGDs07KjOzynNJz5DgvPNg113hvvvgE5/IOyIzs/w46bWAzTeHq65KE9Cef366\nycXMrBU56bWISZPgRz9K62eeCfffn288ZmZ5cNJrIaedBh/8IKxYke7kXLYs74jMzGrLSa+FSHDu\nubDnnvDww/DP/5xudDEzaxVOei1m2LDUvjdiBFxxBfz4x3lHZGZWO+6y0KKuuAJOOCENT3brrbDv\nvnlHZGbWP+6yYN16//vhYx+D1atT+97ixXlHZGZWfS7ptbBVq+DAA2HuXDjySLjmmtTuZ2bWiFzS\nsx5ttlkaqWXUKLjuuo1DlpmZNSsnvRa3444bO6t//vPpJpd16/KNycysWpz0jCOPhM98JiW7446D\nCRPgc5+DBx7IOzIzs8pym54BsHYtfPvb8LOfweOPb9y+995w8snwgQ/ANtvkF5+ZWW9KadNz0rNN\nRMCsWXDxxXDllbBkSdre1gbvfGdKgEcd5ZkazKz+OOlZv6xcCb//fUqA11+fSoOQBrA+7riUAKdN\nSzO1m5nlzUnPKuaFF1LJ7+KL4Y47Nm7fYQc46aSUAN/4xvziMzNz0rOqmD8fLrkELr0Unnpq4/Yp\nU1Lye//7YcyY/OIzs9bkpGdVtX49zJyZSn9XXbVx1oZBg+Bd70oJ8F3vSv0BzcyqzUnPamb5crj2\n2pQAb7wxJUSA0aPh+OPhlFNg6lSP+GJm1eOkZ7n429/g8stTFeg992zcvssuqfR30kmw0075xWdm\nzclJz3J3330p+f3iFykZdjrooFT6O+64NAyamVl/OelZ3Vi3Dv7851T9+ZvfpOpQSO19Rx6ZEuA7\n35naA83M+sJJz+rSsmVw9dWpBPiXv2ycvX3rrdMcf6ecApMnu/3PzMrjpGd1b9GiVPV58cUwb97G\n7bvvnpLfiSfC+PH5xWdmjcNJzxpGBNx1V0p+l1+eOsNDKu297W0pAR57LIwcmW+cZla/nPSsIa1Z\nA3/8Y6r+vPbaNNktwLBhcMwxKQEeemgaD9TMrJOTnjW8JUvSRLeXXAI337xx+3bbwbvfDePGpdkf\nipcRI/KL2WpnxQp4/vn0g2jkyHRjlNuCW5eTnjWVJ55IQ59dfDEsWNDzscOGwdixXSfEwmXs2NRl\nwl+U9eu11+Cxx+DRR9PnXrg8/fSmx7a1peQ3YsTGv4Xr5f4dMcI1Co3ESc+aUgTMng133gnPPdf1\nsnJl6ecbPLj3BNm5f8wYzypRDa+8khLbggWvT26F/TuLDRyYPpcVK+DVV2H16srHNnRo/xNn4bYh\nQ/wjq1qc9KwlRaRuEc89l6q+ukuMncurr5Z+7ra21LWiu6RYuGy9dfpStmTJkk2TWWFye/757p83\neHAawWeXXTZddt0VJkzY9D1evTp9nq++mv4N9PS31GOq8VUzZEhahg7d9G9X28o5prd9zV5qddIz\nK8Hy5aUlx+ee2zipbimkVDLsTIpbbpmqUntaRo9Of4cNa7zSQAS8/PLrqyA7k9tLL3X/3CFDYOed\nX5/YdtkldVnJ68s6Iv37qETyXLYsLZ03ZuVh0KDykuagQRv/HRb+e+xpW2/7q3meb36zzpKepOnA\n94A24LyI+FYXx/wAOBxYDpwaEXd3cYyTnuVi9ereE2Tn/hdf7HspYeDA3hNkcaIsXIYOrU7SjEjd\nSYoTW2dy6+lHwdChG0toxYlt3LjWqTZevz4lvpUrU7VsV3+rsW/FiuqUWutLHSU9SW3Aw8A7gGeA\nO4ATImJewTFHAB+PiCMkTQG+HxFTuzhXwyS9jo4O2tvb8w6jZI0Ub73HunZtSnydifDWWzsYN66d\nJUtg8eKUILpbVqzo37UHDSo/URYu11/fwdix7V0mt1de6f66I0ZsWv1YmNi22646ibje/x0Uyyve\niNQdqJxk+eCDHey2W/uGZFn4tVvJbZU6x5e+1HvSq2WLw/7AgohYCCDpCuAooGAcDo4ELgKIiDmS\nRknaJiKeq2GcFeX/kNVT77EOHAjbbpsWSEnvIx9pL+m5q1bB0qUbk2BvSbJ4Wbkylcg6O/mXrwPo\nOtbNN980oRWujx1b+2rZev93UCyveKXUPjp4cPoMSzFjRgcf/Wh7VeOqpC99qfdjapn0xgGLCh4/\nDUwp4Zg3AA2b9Mz6YrPNUgIZO7Zvz1+5srwkWbgsXpyqGidNev2NI7vsktopG6290axTLZNeqfWR\nxf+dGqMe06yODBmyaSmzXDNmpMWs2dSyTW8qMCMipmePvwCsL7yZRdL/Ah0RcUX2eD5wSHH1piQn\nQjMze516atO7E9hV0kTgWeB44ISiY64DPg5ckSXJJV215/X2oszMzLpSs6QXEWslfRz4I6nLwvkR\nMU/SGdn+n0TE9ZKOkLQAeA04rVbxmZlZ82vIzulmZmZ90VDdQSVNlzRf0qOSPpd3PD2RdIGk5yTd\nn3csvZE0XtJfJT0o6QFJZ+UdU08kDZE0R9I9kh6S9I28Y+qNpDZJd0v6bd6x9EbSQkn3ZfHennc8\nPcm6NV0laV72b+F1/XrrgaQ3Zu9n57K0Af6ffSH7Trhf0mWSNss7pu5IOjuL8wFJZ/d4bKOU9Erp\n3F5PJE0DXgUujoi35B1PTyRtC2wbEfdIGgHMBY6u1/cWQNKwiFguaSBwC/DpiLgl77i6I+lfgX2A\nkRFxZN7x9ETSE8A+EfFy3rH0RtJFwE0RcUH2b2F4RCzNO66eSBpA+g7bPyIW9XZ8HrJ7L/4C7B4R\nqyRdCVwfERflGlgXJO0JXA7sB6wB/gB8NCIe6+r4RirpbejcHhFrgM7O7XUpIm4GFucdRyki4u8R\ncU+2/ippwIDt842qZxGxPFsdTGojrtsvaElvAI4AzuP1XXLqVd3HKWkLYFpEXADpvoF6T3iZdwCP\n1WvCy7xCSiDDsh8Tw0iJuh69CZgTESsjYh1wE3Bsdwc3UtLrquP6uJxiaVrZL7y9gTn5RtIzSQMk\n3UMauOCvEfFQ3jH14LvAZ4D1eQdSogD+JOlOSR/JO5ge7Ai8IOlCSXdJ+pmkYXkHVYL3A5flHURP\nslL+d4CnSHfbL4mIP+UbVbceAKZJ2jL7/N9FGtSkS42U9BqjHraBZVWbVwFnZyW+uhUR6yNiL9I/\n7oMlteccUpckvRt4Phs4ve5LT5kDI2Jv0sDv/5JV1dejgcBk4NyImEy64/vz+YbUM0mDgfcAv8o7\nlp5I2hn4JDCRVOszQtKJuQbVjYiYD3wLuBG4AbibHn5gNlLSewYYX/B4PKm0ZxUgaRDwa+DSiLgm\n73hKlVVn/R7YN+9YunEAcGTWTnY58HZJF+ccU48i4m/Z3xeA35CaFurR08DTEXFH9vgqUhKsZ4cD\nc7P3tp7tC9wWES9FxFrgatK/5boUERdExL4RcQiwhHT/R5caKelt6Nye/Vo6ntSZ3fpJkoDzgYci\n4nt5x9MbSVtJGpWtDwX+kfTrru5ExBcjYnxE7Eiq1vpLRJySd1zdkTRM0shsfThwGFCXdyBHxN+B\nRZJ2yza9A3gwx5BKcQLpx0+9mw9MlTQ0+354B1C3TQiSxmZ/JwDH0EP1ccPM69xd5/acw+qWpMuB\nQ4AxkhYBX4mIC3MOqzsHAicB90nqTB5fiIg/5BhTT7YDLsrughsAXBIRf845plLVezX9NsBv0vcc\nA4FfRMSN+YbUo08Av8h+CD9GHQ9okf2IeAdQz+2kAETEvVmNxJ2kqsK7gJ/mG1WPrpI0hnTzzcci\notsJsBqmy4KZmVl/NVL1ppmZWb846ZmZWctw0jMzs5bhpGdmZi3DSc/MzFqGk56ZmbUMJz2zCpK0\nLps65h5JcyW9Nds+sa/TTGVT/WzZh+d9VdKhfbmmWbNqmM7pZg1ieTZuJZIOA74BtPfznH3qTBsR\n5/TzumZNxyU9s+rZgi6mPMomwb0wm6j1rs7BsrOJZr+dTYZ5r6R/KXreUEk3SDq9aHubpJ9nz7uv\ncxLNbNt7Je1TMHnp/ZLWZ/t3zs53p6SZkt7YRaxjJN2YTc75s76WOs3qhUt6ZpU1NBvKbQhpuLS3\nd3HMvwDrImJSlmhuzMaP/BAwAfiHiFgvaXTBc0YCVwIXRcSlRefbC9i+c7JiSZtn2wOIiJhLmi4K\nSf8FXJ/t/ylwRkQskDQFOBcorg49B5gZEV+XdARwOmYNzEnPrLJWFFRvTgUuBvYsOuZA4AcAEfGw\npCeB3UgJ58cRsT7b1zkJsYBrgW9FRFeDFT8G7CTpB6QZJwrHytwwnZGk40mzEPxjNo3UW4FfZeNs\nQpqQt9g00gC+RMT1khpiYmSz7rh606xKImI2sJWkrbrY3d3cel1tD+AW0rQ0XV1nCTAJ6AA+Spqh\nfdOTSnuSSm3HRxpwdwBpYtC9C5Y3lxGTWUNy0jOrEklvIs0I8lLRrpuBE7NjdiNVac4H/g84Q1Jb\ntq+wevMrwGJJP+riOmOAgRFxNfBlsqrMTGTTMF0OnBwRLwFko9A/Iel92TkkaVIXL2Mm8IHsmMOB\n0V0cY9YwnPTMKmto500jwBXAKbFxKpPOv+cCAyTdlx3zwYhYQyqhPUWa4uke0txrG0TE2dn5v1V0\nzXHAX7NrXgJ8oWj/kaTEel4W213Z9hOB07NrPZAdV+yrpJnpHyBVcz5V8jthVoc8tZCZlSybAX6f\niHjdXalmjcAlPTMrh38lW0NzSc/MzFqGS3pmZtYynPTMzKxlOOmZmVnLcNIzM7OW4aRnZmYtw0nP\nzMxaxv8PiFQTwmPwzewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11027af98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "run -i nt_solutions/coding_2_entropic/exo2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Insert your code here."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}