residual_entropy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
      "[1 1 1 0 0 1 0 1 0 0 0 0 1 0 0 0 0 1 1 1 0 0 0 1 0 0 1 0 0 1 1 0 1 1 1 1 1\n",
      " 0 1 1 1 0 0 1 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1\n",
      " 1 1 0 1 1 1 0 0 0 1 1 0 1 0 1 1 0 0 0 0 0 1 0 0 0 0 1 1 0 1 0 1 1 1 1 1 1\n",
      " 0 1 1 1 1 0 1 1 1 0 0 1 0 1 1 0 1]\n",
      "01001001100001101001101100100000100110110111110110010010111000110011011110111111010001001010101100010001011111001100010010111111\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "from gen_bits_with_probability import gen_bits_with_probability\n",
    "import numpy as np\n",
    "from numpy import random\n",
    "\n",
    "length  = 128\n",
    "probabilities = [0.5, .5] # probability of  0 and 1 respectively\n",
    "np.random.seed(1527)\n",
    "random.seed(1527)\n",
    "bit_numpy = np.random.randint(0, 2, size=128)\n",
    "bit_pyrandom = [random.randint(0, 1) for _ in range(128)]\n",
    "bit_custm = gen_bits_with_probability(128, probabilities)\n",
    "print(bit_pyrandom)\n",
    "print(bit_numpy)\n",
    "print(bit_custm)\n",
    "print(np.array_equal(bit_numpy, bit_pyrandom))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Test with *n* item's prob**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribution name: zipf\n",
      "Probabilities:  [0.365 0.159 0.098 0.069 0.053 0.043 0.035 0.03  0.026 0.023 0.021 0.019\n",
      " 0.017 0.015 0.014 0.013]\n",
      "Sum of probabilities:  1.0\n",
      "Max probability:  0.365\n",
      "Shannon entropy :  3.134\n",
      "Min entropy:  1.453\n",
      "Remaining items: 15\n",
      "Ousted prob: 0.365\n",
      "New Probabilities:  [0.251 0.154 0.109 0.083 0.067 0.056 0.047 0.041 0.036 0.032 0.029 0.027\n",
      " 0.024 0.022 0.021]\n",
      "Sum of new probabilities:  1.0\n",
      "New Max probability:  0.251\n",
      "New Shannon entropy:  3.445\n",
      "New Min entropy:  1.997\n",
      "Calculated remaining entropy: 3.445\n",
      "True\n",
      "Remaining entropy estimate: 3.614\n"
     ]
    }
   ],
   "source": [
    "# %%time\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "from math import log2\n",
    "from numpy import zeros\n",
    "\n",
    "''' set number of bits so that it we can have 2^b items, here items represents each bit sequence.\n",
    "For example 4 bits would have 2^4 = 16 bit sequences'''\n",
    "\n",
    "number_bits = 4\n",
    "\n",
    "'''different ways to generate the items'''\n",
    "\n",
    "# items = zeros(2 ** number_bits, dtype=np.int8) # 2 represents 0 and 1 of a bit string. initialize the array with smallest possible data type to save RAM!\n",
    "\n",
    "\n",
    "items = []\n",
    "for i in range(2 ** number_bits):\n",
    "    items.append(i)\n",
    "\n",
    "# items = np.linspace(start=0.0, stop=2 ** number_bits, num=2 ** number_bits, endpoint=True)\n",
    "# items =  np.arange(0,(2**number_bits), step=1)\n",
    "\n",
    "num_items = len(items) # values of the random variable X\n",
    "probs = []\n",
    "\n",
    "\n",
    "'''choose one of the following distribution names'''\n",
    "    ## dist_name: custom, zipf, zipfian, logser, boltzmann, dlaplace, geom, betabinom,uniform, yulesimon, skellam\n",
    "    ## poisson, planck, nchypergeom_fisher, nchypergeom_wallenius, uniform\n",
    "\n",
    "dist_name = 'zipf' #  insert as string\n",
    "\n",
    "\n",
    "match dist_name:\n",
    "    case 'custom':\n",
    "        '''chosse one of the following probability distributions'''\n",
    "        probs = [.4, .4, .1, .1]\n",
    "        # probs = [1/4, 3/4]\n",
    "        # probs = [7/8, 1/16, 1/16]  # probabilities of the values manually defined, 3 items\n",
    "        # probs = [1/2, 1/4, 1/8, 1/8]  # probabilities of the values manually defined, 4 items\n",
    "        # probs = [3/4, 1/8, 1/16, 1/16]  # probabilities of the values manually defined, 4 items\n",
    "        # probs = [1/2, 1/4, 1/8, 1/16, 1/16]  # probabilities of the values manually defined, 5 items\n",
    "        # probs = [1/4, 1/8, 1/8, 1/8, 1/8, 1/8, 1/16, 1/16]  # probabilities of the values manually defined, 8 items\n",
    "\n",
    "        ##  paper and book example probabilities:\n",
    "        # probs = [2 ** (-3/4), (1 - 2 ** (-3/4))/7,(1 - 2 ** (-3/4))/7, (1 - 2 ** (-3/4))/7, (1 - 2 ** (-3/4))/7, (1 - 2 ** (-3/4))/7, (1 - 2 ** (-3/4))/7, (1 - 2 ** (-3/4))/7 ]\n",
    "\n",
    "        probs = np.array(probs) # convert manually set probs to numpy array\n",
    "        print('Nr. Items:', len(probs))\n",
    "    case 'zipf':\n",
    "        from scipy.stats import zipf # Generate the probabilities using the Zipf distribution\n",
    "        \n",
    "        _param = 1.2 # _param > 1.0\n",
    "        probs = zipf.pmf(np.arange(1, (num_items+1)), _param)\n",
    "        probs /= probs.sum() # Normalize the probabilities so their sum is 1.0\n",
    "    case 'zipfian':\n",
    "        from scipy.stats import zipfian\n",
    "        a, n = 1.25, 10\n",
    "        probs = zipfian.pmf(np.arange(1, (num_items+1)), a,n)\n",
    "        probs /= probs.sum()\n",
    "    case 'logser':\n",
    "        from scipy.stats import logser\n",
    "        p = 0.99 #for smaller p values the first item has a higher probability, 0<p<1\n",
    "        probs = logser.pmf(np.arange(1, (num_items+1)), p)\n",
    "        probs /= probs.sum()\n",
    "    case 'boltzmann':\n",
    "        from scipy.stats import boltzmann\n",
    "        lambda_, N = 1.4, 19\n",
    "        probs = boltzmann.pmf(np.arange(1, (num_items+1)),  lambda_, N)\n",
    "        probs /= probs.sum()\n",
    "    case 'dlaplace':\n",
    "        from scipy.stats import dlaplace\n",
    "        a = 0.1\n",
    "        probs = dlaplace.pmf(np.arange(1, (num_items+1)),  a)\n",
    "        probs /= probs.sum()\n",
    "    case 'geom':\n",
    "        from scipy.stats import geom\n",
    "        p = .2\n",
    "        probs = geom.pmf(np.arange(1, (num_items+1)),  p)\n",
    "        probs /= probs.sum()\n",
    "    case 'betabinom':\n",
    "        from scipy.stats import betabinom\n",
    "        n, a, b = 5, 2.3, 0.63\n",
    "        probs = betabinom.pmf(np.arange(1, (num_items+1)),  n,a,b)\n",
    "        probs /= probs.sum()\n",
    "    case 'uniform':\n",
    "        from numpy import random\n",
    "        probs = np.random.random(size=num_items)\n",
    "        probs /= probs.sum()\n",
    "    case 'yulesimon':\n",
    "        from scipy.stats import yulesimon\n",
    "        a = .01 # 11 gives first item 0.9 probability and 0.01 gives 0.299 for the first item\n",
    "        probs = yulesimon.pmf(np.arange(1, (num_items+1)),  a)\n",
    "        probs /= probs.sum()\n",
    "    case 'skellam':\n",
    "        from scipy.stats import skellam\n",
    "        mu1, mu2 = 15, 8\n",
    "        probs = skellam.pmf(np.arange(1, (num_items+1)),  mu1, mu2)\n",
    "        probs /= probs.sum()\n",
    "    case 'poisson':\n",
    "        from scipy.stats import poisson\n",
    "        mu = 0.6\n",
    "        probs = poisson.pmf(np.arange(1, (num_items+1)),  mu)\n",
    "        probs /= probs.sum()\n",
    "    case 'planck':\n",
    "        from scipy.stats import planck\n",
    "        alpha = 0.5\n",
    "        probs = planck.pmf(np.arange(1, (num_items+1)),  alpha)\n",
    "        probs /= probs.sum()\n",
    "    case 'nchypergeom_wallenius':\n",
    "        from scipy.stats import nchypergeom_wallenius\n",
    "        M, n, N, odds = 140, 80, 60, 0.5\n",
    "        probs = nchypergeom_wallenius.pmf(np.arange(1, (num_items+1)),  M, n, N, odds)\n",
    "        probs /= probs.sum()\n",
    "    case 'nchypergeom_fisher':\n",
    "        from scipy.stats import nchypergeom_fisher\n",
    "        M, n, N, odds = 140, 80, 60, 0.5\n",
    "        probs = nchypergeom_fisher.pmf(np.arange(1, (num_items+1)),  M, n, N, odds)\n",
    "        probs /= probs.sum()\n",
    "    case _:\n",
    "        print('No distribution name given or wrong distribution name.')\n",
    "\n",
    "\n",
    "'''Calculate the minimum entropy'''\n",
    "min_entropy = - math.log2(max(probs))\n",
    "\n",
    "''' Calculate the Shannon entropy '''\n",
    "entropy_shannon = stats.entropy(probs, base=2) # by python stats module\n",
    "\n",
    "\n",
    "print('Distribution name:', str(dist_name))\n",
    "print('Probabilities: ', np.round(probs, decimals=3))\n",
    "print('Sum of probabilities: ', np.round(sum(probs), decimals=3))\n",
    "print('Max probability: ', np.round(max(probs), decimals=3))\n",
    "print('Shannon entropy : ', np.round(entropy_shannon, decimals=3))\n",
    "# print('Calculated Shannon entropy: ', np.round(entropy_formula(probs), decimals=3)\n",
    "print('Min entropy: ', np.round(min_entropy, decimals=3))\n",
    "\n",
    "'''\n",
    "Problem : what would be the entropy of the remaining items if we take out THAT one item with max probability,\n",
    "given that the rest items DOES NOT have UNIFORM prob distribution ?'''\n",
    "item_remove = 'one_max' # choose one of the following: 'all_max', 'one_max'\n",
    "# item_remove = 'all_max' # choose one of the following: 'all_max', 'one_max'\n",
    "match item_remove:\n",
    "    case 'all_max':\n",
    "        '''remove ALL the item with max probability and renormalize probability for the rest of the items'''\n",
    "        prob_max = probs.max()\n",
    "        new_probs = probs[probs != prob_max]\n",
    "        prob_adjust = 1 - prob_max* (len(probs)- len(new_probs))\n",
    "        # print(prob_adjust)\n",
    "        # new_probs = [item / prob_adjust for item in new_probs]\n",
    "        # print(new_probs)\n",
    "    case 'one_max':\n",
    "        '''remove only ONE item with max probability'''''\n",
    "        prob_max_val = probs.max()\n",
    "        prob_max = np.argmax(probs)\n",
    "        probs =list(probs)\n",
    "        del probs[prob_max]\n",
    "        new_probs = np.array(probs)\n",
    "        prob_max = prob_max_val\n",
    "        prob_adjust = 1 - prob_max\n",
    "        # new_probs = [item / prob_adjust for item in new_probs]\n",
    "        # print(new_probs)\n",
    "\n",
    "new_probs = [item / prob_adjust for item in new_probs]\n",
    "\n",
    "'''# Calculate the new minimum entropy'''\n",
    "min_entropy_new = - log2(max(new_probs))\n",
    "\n",
    "''' Calculate the new Shannon entropy using the stats module'''\n",
    "entropy_shannon_new = stats.entropy(new_probs, base=2)\n",
    "entropy_shannon_new = np.round(entropy_shannon_new, decimals=3)\n",
    "\n",
    "'''with shannon entropy from original probs distribution, calculate the remaining entropy'''\n",
    "bin_entropy_h = - (prob_max * log2(prob_max) + (1-prob_max) * log2(1-prob_max))# binary entropy of the prob max\n",
    "\n",
    "'''remaining entropy calculation'''\n",
    "entropy_remaining = (entropy_shannon - bin_entropy_h)/ (1 - prob_max)\n",
    "entropy_remaining = np.round(entropy_remaining, decimals=3)\n",
    "\n",
    "entropy_remaining_test = entropy_shannon + (prob_max * log2(len(new_probs))) - bin_entropy_h\n",
    "\n",
    "# print('New Probabilities: ', np.round(new_probs), decimals=3)\n",
    "print('Remaining items:', len(new_probs))\n",
    "print('Ousted prob:' , np.round(prob_max, decimals=3))\n",
    "print('New Probabilities: ', np.round(new_probs, decimals=3))\n",
    "print('Sum of new probabilities: ',np.round(sum(new_probs), decimals=3))\n",
    "print('New Max probability: ', np.round(max(new_probs), decimals=3))\n",
    "print('New Shannon entropy: ', entropy_shannon_new)\n",
    "# print('Updated Calculated Shannon entropy: ', entropy_formula(new_probs))\n",
    "print('New Min entropy: ', np.round(min_entropy_new, decimals=3))\n",
    "print('Calculated remaining entropy:', np.round(entropy_remaining, decimals=3))\n",
    "print(np.array_equal(entropy_shannon_new, entropy_remaining))\n",
    "print('Remaining entropy estimate:', np.round(entropy_remaining_test, decimals=3))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "latex"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "#Distribution name: logser\n",
    "Probabilities:  [0.215 0.106 0.07  ... 0.    0.    0.   ]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.215\n",
    "Shannon entropy :  5.282\n",
    "Min entropy:  2.218\n",
    "Remaining items: 1048575\n",
    "Ousted prob: 0.215\n",
    "New Probabilities:  [0.136 0.089 0.066 ... 0.    0.    0.   ]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.136\n",
    "New Shannon entropy:  5.772\n",
    "New Min entropy:  2.883\n",
    "Calculated remaining entropy: 5.772\n",
    "True\n",
    "Remaining entropy estimate: 8.831\n",
    "\n",
    "# zipf\n",
    "Nr. Items: 16\n",
    "Probabilities:  [0.365 0.159 0.098 0.069 0.053 0.043 0.035 0.03  0.026 0.023 0.021 0.019\n",
    " 0.017 0.015 0.014 0.013]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.365\n",
    "Shannon entropy :  3.134\n",
    "Min entropy:  1.453\n",
    "Remaining items: 15\n",
    "Ousted prob: 0.365\n",
    "New Probabilities:  [0.251 0.154 0.109 0.083 0.067 0.056 0.047 0.041 0.036 0.032 0.029 0.027\n",
    " 0.024 0.022 0.021]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.251\n",
    "New Shannon entropy:  3.445\n",
    "New Min entropy:  1.997\n",
    "Calculated remaining entropy: 3.445\n",
    "True\n",
    "\n",
    "#Book example\n",
    "Nr. Items: 8\n",
    "Probabilities:  [0.595 0.058 0.058 0.058 0.058 0.058 0.058 0.058]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.595\n",
    "Shannon entropy :  2.112\n",
    "Min entropy:  0.75\n",
    "Remaining items: 7\n",
    "Ousted prob: 0.595\n",
    "New Probabilities:  [0.143 0.143 0.143 0.143 0.143 0.143 0.143]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.143\n",
    "New Shannon entropy:  2.807\n",
    "New Min entropy:  2.807\n",
    "Calculated remaining entropy: 2.807\n",
    "True\n",
    "\n",
    "# Nr. Items: 5\n",
    "Probabilities:  [0.5   0.25  0.125 0.062 0.062]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.5\n",
    "Shannon entropy :  1.875\n",
    "Min entropy:  1.0\n",
    "Remaining items: 4\n",
    "Ousted prob: 0.5\n",
    "New Probabilities:  [0.5   0.25  0.125 0.125]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.5\n",
    "New Shannon entropy:  1.75\n",
    "New Min entropy:  1.0\n",
    "Calculated remaining entropy: 1.75\n",
    "True\n",
    "\n",
    "#Nr. Items: 4, probs = [3/4, 1/8, 1/16, 1/16]\n",
    "Probabilities:  [0.75  0.125 0.062 0.062]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.75\n",
    "Shannon entropy :  1.186\n",
    "Min entropy:  0.415\n",
    "Remaining items: 3\n",
    "Ousted prob: 0.75\n",
    "New Probabilities:  [0.5  0.25 0.25]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.5\n",
    "New Shannon entropy:  1.5\n",
    "New Min entropy:  1.0\n",
    "Calculated remaining entropy: 1.5\n",
    "\n",
    "#Nr. Items: 4, probs = [1/2, 1/4, 1/8, 1/8]\n",
    "Probabilities:  [0.5   0.25  0.125 0.125]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.5\n",
    "Shannon entropy :  1.75\n",
    "Min entropy:  1.0\n",
    "Remaining items: 3\n",
    "Ousted prob: 0.5\n",
    "New Probabilities:  [0.5  0.25 0.25]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.5\n",
    "New Shannon entropy:  1.5\n",
    "New Min entropy:  1.0\n",
    "Calculated remaining entropy: 1.5\n",
    "True\n",
    "\n",
    "#Nr. Items: 3 , probs = [7/8, 1/16, 1/16]\n",
    "Probabilities:  [0.875 0.062 0.062]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.875\n",
    "Shannon entropy :  0.669\n",
    "Min entropy:  0.193\n",
    "Remaining items: 2\n",
    "Ousted prob: 0.875\n",
    "New Probabilities:  [0.5 0.5]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  0.5\n",
    "New Shannon entropy:  1.0\n",
    "New Min entropy:  1.0\n",
    "Calculated remaining entropy: 1.0\n",
    "True\n",
    "\n",
    "#Nr. Items: 2\n",
    "Distribution name: custom\n",
    "Probabilities:  [0.25 0.75]\n",
    "Sum of probabilities:  1.0\n",
    "Max probability:  0.75\n",
    "Shannon entropy :  0.811\n",
    "Min entropy:  0.415\n",
    "Remaining items: 1\n",
    "Ousted prob: 0.75\n",
    "New Probabilities:  [1.]\n",
    "Sum of new probabilities:  1.0\n",
    "New Max probability:  1.0\n",
    "New Shannon entropy:  0.0\n",
    "New Min entropy:  -0.0\n",
    "Calculated remaining entropy: 0.0\n",
    "True\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK= zipf, logser, poisson, zipfian, yulesimon, geom\n",
    "NOT OK = planck, , betabinom,[hypergeom, ficher noncentralhyper, Wallenius (new sum is >1)]\n",
    "\n",
    "for boltzman, the new distribution keeps the values of the original distribution\n",
    "\n",
    ">> Min entropy does not capture the  guessing complexity which is needed for the brute force attack\n",
    " to break a key. For example, if we have a very high random bit sequence then we should have  Min-entropy ~1.0 \n",
    " but that does not change with the bit string's length. That is uniformly generated 2^16 bit sequence and 2^128 bit\n",
    " sequence have Min-entropy ~ 1.0.  This does not give us the full picture of the key's strength.\n",
    " Min-entropy can only hints about the underlying distribution. SO it tells about the goodness of the key's \n",
    " randomness but not the hardness of the key.\n",
    "\n",
    ">> Duality of the distribution.\n",
    "Random number generator produces bitstrings from some entropy source. Which gurantees that the distribution of the bitstring is near to uniform. But then those raw bitstring is to clip of to a certain (desired length) that can be used a key for cryptopgraphic application. What if one of this bit sequences is having non-uniform probability?\n",
    "If we can pop it out and ommit ths from the avaiable pool of keys, can we asses the remaining probability of the key pool, without reassesing all the possible keys?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "vscode": {
     "languageId": "latex"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8112781244591328 1.0\n",
      "0.4150374992788438\n"
     ]
    }
   ],
   "source": [
    "## Proof of Lemma 5 when |X|= 2 and 0 < p < .5\n",
    "\n",
    "\n",
    "x1 = 1/4 # <0.5\n",
    "x2 = 1 - 1/4\n",
    "X = [x1, x2]\n",
    "\n",
    "leamma_5 = - 2 * x1 * log2(x1) # -2plog(p) when 0<p<0.5\n",
    "\n",
    "H_X = stats.entropy(X, base=2)\n",
    "H_min = - log2(max(X))\n",
    "\n",
    "print(H_X, leamma_5)\n",
    "print(H_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the  original and new distributions after ousting the item that has the highest probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwkAAAIgCAYAAADOTyXGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwRklEQVR4nO3deVxU9f7H8fcwbOKCioZairibmPu+gN00l/JqWZrlUmpZmduv0rKb1q1su2WLy20xu2Vp+2omlZpe82qmlommSWomkqig7Azf3x8nBg4MiiwO4uv5eMyDmXO+55zPHEHmzfd7ztdhjDECAAAAgL/4eLsAAAAAAOULIQEAAACADSEBAAAAgA0hAQAAAIANIQEAAACADSEBAAAAgA0hAQAAAIANIQEAAACADSEBAAAAgA0hAQAKsWnTJkVFRcnhcKhFixaKiopS165d1bp1a/373/8u8f53797t3v+aNWuKvZ8BAwaoevXqmjNnTqFt0tPTVb9+fW3evFmS9Nxzz6lFixZq2LBhoW0kac2aNVqyZEmxa/Nk/fr16tq1qxwOh3777Tf38qFDh2revHlF2seJEyc0Z84cnThxosjHbd++vT744ANJ0vLly9W2bVs5HI6zqNwzT+fI07kEgPMJIQEACtG5c2f3h/eZM2dqzZo12rhxoxYtWqQ777xTy5YtK9H+mzdvXqJwkOOLL75Q27ZtT9vG6XSqefPmqlq1qiRpypQpmjlz5mnbSGUTEnr27Onx3DVs2FChoaFF2seJEyf00EMPnVVIaNasmWrWrClJGj58eJEDyZl4OkeeziUAnE98vV0AAJxvevTooYiICL3//vsaMWKEt8spEl9fX3311VclblOWnn322TLdf0lD3dnw9rkEgJKiJwEAiiEzM1P+/v7au3eve8jQK6+8omHDhql169aqXr26JCkrK0szZ85URESEOnXqpD59+mj79u0F9rd3714NHTpU7du3V7t27bRlyxb3um3btmngwIHq1auXevbsqaFDh+r3338vsI+UlBTddttt6tGjhxo3bqw33njDva5fv35nHJKUv80zzzyjJUuWaNu2bYqKilJUVJTeeOMN1alTRzVr1tTNN98sSfrmm28UERGh8PBwRUdHe9z37t271aNHD7Vu3VpXXXWVNm3aZFt/7733qmHDhoqKinIv++GHHxQZGamoqCh169ZNt9xyi+Li4rRz5053OBsxYoSioqL04Ycf6sUXX3QPoVqyZIkGDhyomjVraurUqRo9erTq1KmjsWPHFqjt888/18CBA9W0aVMNHjxYR48elSQ9/vjjtpoSExMLDA/zdI5SU1M9nu/TfS+cOnVKUVFRCgwM1JNPPqlRo0apU6dO6tatm2JjYwv9NwOAMmMAAKclybz22mvu18uXLzcOh8OsWrXK1ubKK680aWlpxuVyme7duxtjjLnvvvtM27ZtzcmTJ40xxvz73/82tWvXNidOnLBtGxkZaVJSUowxxsyZM8dcfPHFJi0tzRhjzAsvvGCmTZvmbv/www+bPn362GqMjIw0derUMb/++qsxxpg1a9YYHx8fs337dlub2bNnu1+/9tprJiwsrMB+8raZPXu2iYyMtLV55ZVXTM2aNd31GWPM9OnTzddff+3x/LlcLtOyZUszadIkY4wxWVlZZsSIEUaSiY2NLfRYLVu2NK+++qoxxpjMzEwTGRlpVq9ebYwxJjY2tsD2Oe+pUqVKZsGCBcYYY7755hszc+ZMY4wxY8aMMWPGjHG3Xb16tZFk7rnnHvcx+vbta6699trTvn9J7joKa2NMwXNZlO+FsLAw06lTJ3eboUOHmtGjRxfYNwCUNXoSAKAIHn/8cfeFy4sXL9aKFSvUt29fW5sbbrhBAQEB8vHx0X//+1+lpqbq2Wef1Z133qkqVapIksaNG6fs7Gy99NJLtm1HjRqlSpUqSbKuF/jjjz/0/vvvS7L+Wp73L9LXX3+91qxZo9TUVNs+rrjiCjVq1EiSFBkZqTZt2ujFF18s1fOQc/z09HR9/PHHkqxelfXr16tPnz4e20dHRysmJkbTp0+XZI3XHz9+/BmPc+jQIR08eFCSNXxn0aJFat269Rm3c7lc7v336dNHc+fOPW372267zX2MSZMm6YMPPvDYU1MSZ/O9MHjwYHebqKgobdu2rVRrAYCi4JoEACiCmTNnehyqktcll1xie713716lpaWpadOm7mVOp1MNGzbUjh07bG3DwsLcz6tXr66QkBDFxMRIkhwOh2bPnq1NmzbJ19dX6enpMsYoPj7etl3e55LUuHFj9z5KU9WqVTVs2DAtXrxY119/vT777DMNGDCg0DsF7dq1S06n01ZfgwYNznicuXPnatq0aXr77bc1cuRIjR8/XiEhIWfc7qKLLpKfn1+R30/euho3bixjjHbv3l3g37MkzuZ7oW7duu7nVatWVVJSUqnVAQBFRU8CAJQSp9Npe22MKbRt/g/U+dsaY9xtRo8erY0bN2rlypVau3at+wJcT9sUto/Sdssttyg6Olq///67lixZojFjxhTa9nTn4XTuuOMOHThwQBMmTNDbb7+tFi1aFOmWovn/Hc4kb305z3POW/7z53K5zmrfno6RX/5j5K3f4XAU+/wBQEkQEgCgjDRt2lSBgYHas2ePe5nL5dJvv/2miIgIW9sDBw64n584cULHjh1TixYtJElr167VwIED3bfTzMjI8Hi8vPuQpH379rn3UVw+Prm/JtLS0pSZmSlJ6t27txo1aqQnnnhCqampCg8PL3Qfl156qVwul/bv319orZ689957Cg0N1f/93//pp59+UqtWrfSf//ynQF3Z2dlKTk4+6/fmqZZ9+/bJ4XCoWbNmkqy/5J86dcq9/tChQwW2L+wc5XU23wsAUB4QEgCgjFSqVEnTpk3TggUL3B9iX331Vfn4+GjChAm2tosXL3ZfY/Dcc8+pXr16uvbaayVJrVq10po1a5SVlSVJ7msV8luxYoX7Tjhr167V9u3bNWnSpBK9h9q1a+v48eOSpOnTp2vVqlXudWPHjtWLL76oUaNGnXYfV1xxhVq2bKlnnnlGkvXheOHChWc89oQJE3TkyBFJ1l/UXS6XmjdvLkkKCQmRj4+Pjh8/ru+///6MQ8FOZ/78+ZKsuw+9+OKLuuaaa9xDjdq2bauYmBj3OXj77bcLbH+6c5TjbL4XAKBc8Mrl0gBwHvjf//5nIiMjjSTTvHlz079//wJtDh8+7G7Tpk0bM2vWLNv6zMxMM2PGDNOqVSvTsWNHExkZabZu3WqMMWbXrl3ubZ999lkzYMAA065dO9O2bVuzefNm9z527NhhevbsaZo3b24GDx5s7r33XiPJdOnSxfz444+mf//+Jjg42Nx1113mxhtvNN27dzfh4eHm9ddfd++jb9++Jjg42ISFhZlZs2aZefPmmebNm5uAgAATGRlpTp48WaCNMcYcOXLEdOrUyfTo0cMMHDjQdkejgwcPmmrVqpnk5OQznstdu3aZ7t27m1atWpm+ffual19+2f0e1q1bZ+655x4TFhZmgoODzaBBg4wx1t2AOnToYKKiokynTp3M9OnTTVZWlnuf9957r2nVqpXp0qWLWb9+vXnttdds72ndunXutqNGjTKhoaEmNDTUjBs3zixbtsy0adPGSDJLly41ffv2NY0bNzZXXXWV+fPPP22133HHHaZZs2Zm0KBB5uOPP3b/W7/77ruFniNP5/J03wvGWHdDCggIMM2bNzdLly41y5Ytc7+fyy+//IznGABKk8MYBjsCAM7exo0b9eqrr+rll1/2dikAgFLGcCMAwFl57LHHJEmLFi0q0q1MAQDnH3oSAABnpX79+qpdu7Z69eql5557ztvlAADKACEBAAAAgI1Xhxt9++23uvrqq1WvXj05HA599NFHZ9xm7dq16tChgwIDA9WoUSMtWrSo7AsFAAAALiBeDQnJyclq06aNXnzxxSK1j42N1cCBA9WrVy9t3bpV999/vyZPnlzo7QABAAAAnL1yM9zI4XDoww8/1JAhQwptM2PGDH3yySeKiYlxL5s4caK2b9+u77777hxUCQAAAFR8vt4u4Gx899136tevn23ZlVdeqVdffVWZmZny8/MrsE16errS09Pdr7Ozs3Xs2DGFhITI4XCUec0AAADAuWCM0cmTJ1WvXj3bbPDFcV6FhLi4OIWGhtqWhYaGKisrS0ePHlXdunULbDN37lw99NBD56pEAAAAwKsOHjzonjm+uM6rkCCpwF//c0ZLFdYrcN9992n69Onu14mJiWrQoIEOHjyoatWqlV2hAAAAwDmUlJSk+vXrq2rVqiXe13kVEurUqaO4uDjbsvj4ePn6+iokJMTjNgEBAQoICCiwvFq1aoQEAAAAVDilMaT+vJpxuVu3boqOjrYtW7VqlTp27OjxegQAAAAAZ8+rIeHUqVPatm2btm3bJsm6xem2bdt04MABSdZQodGjR7vbT5w4Ufv379f06dMVExOjxYsX69VXX9Xdd9/tjfIBAACACsmrw42+//579enTx/0659qBMWPGaMmSJTp8+LA7MEhSeHi4VqxYoWnTpmn+/PmqV6+enn/+eV177bXnvHYAAACgoio38yScK0lJSQoODlZiYiLXJAAAUI4YY5SVlSWXy+XtUoByyel0ytfXt9BrDkrzc+55deEyAAComDIyMnT48GGlpKR4uxSgXAsKClLdunXl7+9fpschJAAAAK/Kzs5WbGysnE6n6tWrJ39/fyY8BfIxxigjI0N//vmnYmNj1bRp0xJPmHY6hAQAAOBVGRkZys7OVv369RUUFOTtcoByq1KlSvLz89P+/fuVkZGhwMDAMjvWeXULVAAAUHGV5V9FgYriXP2c8NMIAAAAwIaQAAAAAMCGkAAAAFBMqampmjNnjrp166aoqCh17dpVY8aM0W+//VbkfcybN09Dhw49q+Omp6erfv362rx581lWXLjdu3crKipKDodDa9asKbX9no21a9eqa9eucjgcZ3UO80pMTFRUVJQCAwO1ZMmSQtsdOnRIoaGhOnTokCRrEt+GDRsqKiqq0DaS9NFHH+mjjz467b4qAkICAABAMaSlpelvf/ubTp06pW+//VZr1qzRxo0b1b9/f3Xr1k0///xzkfZz0UUXqWHDhmd1bKfTqebNm6tq1arFqNyz5s2bey0c5IiMjNSyZctKtI/g4GCtWbNGderUOW27wMBANW/e3H3x79y5czV27NjTtpE8hwRP7c533N0IAACgGGbPnq309HQ99dRTtlu23nDDDdq4caPGjBmj77///oz7GTlypEaOHHlWx/b19dVXX3111jUjV0hIiL799tsStzmbducTehIAAEC5YoyRK9nllYcxpkg1ZmVladGiRRo+fLjHOR1uuOEGbdmyRZs3b9bevXvdw3heeeUVDRs2TK1bt1b16tX11ltvqW3btgX28cknn6h58+bq0qWLhg4dqjlz5igwMFBRUVE6deqU+vXrp+rVq2vOnDmSpPfee8+9n88++0yDBw9W06ZNddddd9n2++6776p79+7q06ePOnfurOnTpys9Pb3I/zZ5hwO98sor6tevn1q1aqXevXtr3759kqT169e727zzzjv6+9//rsaNG6tt27aSpFOnTunWW29V69at1b59e1199dUehxZt2rRJV111lSIiImz7l6RvvvlGffr0UVRUlLp166axY8fqxIkTBfYRHx+vkSNHqkuXLmrZsqW+/PJLSdKxY8fOOCTJU5t7771XK1eu1MqVKxUVFaW///3vhe7r5MmTGjdunNq1a6fIyEgNGTJEBw4ccK+Pjo5Wt27d1KdPH3Xp0kWTJ09WcnJykf8tyho9CQAAoFzJTsnWuirrvHLsXqd6yVnZecZ2u3fvVlJSklq0aOFxfcuWLSVJmzdv1h133KE1a9bI4XDovffe08cffyw/Pz/16tVLI0eOVL169dSnTx/3tvv379d1112nt956S9dee62OHj2qLl26qE6dOu7hQKtWrbKNnR82bJhq1aqlPn36aOfOnfrkk08UHx+vBg0a6JprrnHvf/ny5Zo5c6YGDx6szMxMDR48WE888YQefPDBIp2fnOFA4eHhWrFihb744gs5nU7dcsstuv766/X999+rZ8+e7jYrV67Uhx9+qJSUFA0ZMkSSdOutt+rEiRPaunWrfH19df/992vQoEHavn27fH1zP5p++OGH+vjjjwvsX5JWrFihoUOHavLkyTLG6NZbb9X06dO1ePFiW72vv/661q5dq1q1aun111/XNddco3379ik0NFRr1qw57TCvmjVrFmjz5JNPKj4+XpJsgcDTvsaPHy+n06ktW7bIx8dHjz76qAYMGKAff/xRxhgNGzZMH374oS6//HIlJyerffv2mj59uipXrlykf4uyRk8CAADAWcr5q3WVKlU8rs9Zfvz4cdvyG264QQEBAfLx8dF///tfj9suWrRIderU0bXXXitJqlWr1lkNR8ppe9FFF+nSSy/Vtm3b3OuefvppXXXVVZIkPz8/DRkyRF988UWR953XnXfeKafTClRTp07Vli1btGHDBlubMWPGyMfHR1WqVNFXX32lffv26e2339b06dPdgeDuu+9WTEyMPvzwQ9u2OR+yPe1/+vTpuvXWWyVJDodDw4YN8/g+rrnmGtWqVUuSNGrUKAUFBenVV18t1vs9G/v27dM777yj6dOnu+c1uO2227Rz506tWbNGJ0+eVFJSkrtnoXLlylq2bJlCQ0PLvLaioicBAACUKz5BPup1qpfXjl0U1atXl6RCh4ecOnVKklSjRg3b8ksuueSM+961a5caNWpkW9agQYMi1SVJdevWdT+vWrWqkpKS3K+Tk5N14403av/+/fL391dcXNxZDTfKKywszP28cePGkqSYmBh1797dvTz/+825mLtp06buZTVr1lTNmjW1Y8cOXXfddUXaf1ZWliZNmqSdO3fK399fJ06cUFxc3Glr9PHxUcOGDRUTE1Os93s2duzYIUmaMmWK/Pz8bPX8+eefqlGjhu677z6NGzdOL774okaOHKmbb75ZlSpVKvPaioqQAAAAyhWHw1GkIT/elHNnoZiYGF199dUF1u/atUuS1LVrV9vynL+Mn44xxuN1DkWV9xgOh8N9ncWpU6d0+eWXa/jw4Vq6dKl8fHy0ZMkS93UNZ8vT9Rv5687/fk93zUf+bU+3/wEDBqhFixZavXq1AgICtGbNGtuQrcL2UdJze7befPNNhYeHe1z32GOP6dZbb9Xrr7+uefPm6cknn9TGjRvP+k5XZYXhRgAAAGfJ19dXEydO1Lvvvutx/dtvv62uXbuqffv2Z73vSy+9VL/++qttWd4LXotr165dio+P13XXXeceApORkVHs/eWtKafewq7RyBERESGHw6E9e/a4lx07dkzHjh1TREREkfZ/9OhR7dy5U0OGDFFAQMBp30fefWRnZ2v//v1nrPFMcs6dJKWkpMjlchVok/M+d+/ebVv+4IMPateuXTp58qS+/PJLNWzYULNnz9auXbsUGBio999/v0S1lSZCAgAAQDE8/PDD8vPz04wZM5SVleVe/s477+j999/Xf/7zn2Lt97bbblNcXJz7A2NCQkKB8frF0ahRI1WqVMl961SXy6WPP/642PtbvHixsrOzJVkTwnXo0ME21KiwGkaMGKFnnnnGfc6efvpptWzZ0n1hc45FixZ53H9ISIjq1Kmjr7/+2t32gw8+8Hi8ZcuWKSEhQZL0xhtvKCUlRePGjSvW+81Ru3Zt97Umw4YNc/caeXqfTz75pNLS0iRJGzZs0Pvvv68mTZooISFBd955p224msvlUvPmzUtUW6kyF5jExEQjySQmJnq7FAAAYIxJTU01O3fuNKmpqd4u5aylpKSYOXPmmG7dupnIyEjTtWtXM27cOPP777+72xw+fNhERkYaSaZNmzZm1qxZ7nVLly41bdq0MZJMZGSk2bNnjzHGmI8//tg0a9bMdOnSxQwfPtw89thjpkmTJu7t+vbta4KDg01YWJiZNWuW+eKLL2z7SUhIMGPHjnW3efLJJ40xxnz44YemWbNmpnPnzmbIkCHm5ptvNgEBAebyyy83u3btstX57rvvenzPsbGxRpJ5++23zcCBA82ll15qevbsaX799VdjjDFbt241Xbp0MZJMly5dzAsvvGDb/uTJk2bChAkmIiLCtGvXzgwaNMjExsYaY4xZs2aNe9uXX37Z9OvXr8D+jTFm3bp1pm3btuayyy4zgwcPNnfddZf7vf/2228mMjLSBAQEmDlz5pghQ4aYTp06mebNm5uVK1caY4xJSEhwt2nevLl54YUXzMyZM01YWJgJDg42gwYN8tjGGGNiYmJMRESE6dmzp7npppsKbXfy5Elz6623mubNm5uoqChz1VVXuf99T506Ze666y7ToUMHExUVZTp27Gjmzp1bpO+50/28lObnXIcxRbwhcAWRlJSk4OBgJSYmqlq1at4uBwCAC15aWppiY2MVHh5eoWasLS5jjI4dO6aQkBD3sscee0yrV69WdHS0Fyuz/PbbbwoPD1dsbGy5GT9/ITndz0tpfs5luBEAAEA5kpycrF69eik1NVWSdbvVpUuX6qabbvJyZbiQEBIAAADKkYCAAHXr1k09e/ZUVFSU+vXrp/Hjx2v06NHeLk1r167ViBEjJEkjRozQd9995+WKUFa4BSoAAEA54ufnd04m/CqOyMhIbdy40dtl4BygJwEAAACADSEBAAAAgA0hAQAAAIANIQEAAACADSEBAAAAgA0hAQAAAIANIQEAAOAsfffdd4qKipLD4VCLFi00e/bsAusCAwMVFRWl2NjYYh3jvffeU9u2beVwOM5quzVr1mjJkiXFOmaO9evXq2vXrnI4HPrtt99KtK/iWr58ebHef167d+92/zutWbOm0HabN29W/fr1lZ6eLkm6+eabVadOHY0dO7bQNpK0ZMmSAvv11O58REgAAAA4S926dXN/OJw5c6YeeuihAuvq1KmjNWvWKDw8vFjHGDZsmObNm3fW25VGSOjZs6eWLVtWon2U1PDhw4v1/vNq3rz5acNBjqpVq6p58+by9bWmEHvttdfUv3//07aRPIcET+3OR+d39QAAAEAJtWjRQl999VWJ25xNu/KOngQAAFAuJScX/khLK3rb1NSitS1LOcNXbrjhBk2cOFG9e/dWw4YN9dRTT9nafffdd2rTpo06dOigoUOHau/evbb1WVlZuuuuu9SpUydFRUWpa9eu+uSTT9zrn3nmGS1ZskTbtm1TVFSUoqKilPrXCdi4caN69eql7t27q1u3bvrnP/8pl8vl3nb37t3q0aOHWrdurauuukqbNm067XvKOxxo6dKl6tevn5o2barBgwfr6NGjkuxDpj7//HNdffXVqlevnoYMGSJJOnLkiIYPH662bduqTZs2uvHGG3Xs2LECx/r88881cODAAvuXpHfffVfdu3dXnz591LlzZ02fPt3jUJ+9e/dq6NChat++vdq1a6ctW7ZIknbu3HnGIUme2owePVrbtm3TkiVLFBUVpdtuu63QfR0+fFjDhg1Tx44d1bNnT40ZM8b2Pt966y116tRJffr0UdeuXXX//fef9tyfE+YCk5iYaCSZxMREb5cCAACMMampqWbnzp0mNTXVtlwq/DFwoH0fQUGFt42MtLetVctzu+KQZF577TWP68LCwmyvx4wZY/z8/MzGjRuNMcbExMSYwMBA8+mnnxpjjElKSjIhISHm6aefNsYYk5ycbCIjI03ej2snT540YWFhJikpyRhjzO7du01wcLDZs2ePu83s2bNNZL43feTIEVOtWjWzYsUK937atm1r5s6da4wxxuVymZYtW5pJkyYZY4zJysoyI0aMMJJMbGxsoe9/9erVRpKZPHmyMcaYzMxM07dvX3PttdcWaDN79mxjjDF79+41I0eONMYY0717dzNhwgRjjDHZ2dlm5MiRpm/fvgW2veeeewrd/7XXXms+/vhjY4wxGRkZpn///uahhx6y1SnJREZGmpSUFGOMMXPmzDEXX3yxSUtLs7VZvXq1+/WYMWPMmDFjCuwnb5vIyEj3+zpdu65du5oZM2a43+eECRPMlVdeaYwx5tChQ8bpdJpff/3VGGNMXFycqVGjRoF95ijs58WY0v2cS08CAADAOdKlSxd16dJFkjUsZcCAAXr++eclWX9NPnXqlO644w5JUlBQkEaNGmXbPigoSGvWrFHVqlUlSc2aNVPLli319ddfn/a4L774ourXr68BAwZIkqpUqaIbb7xR8+fPlyRFR0crJiZG06dPlyQ5nU6NHz++yO9r8uTJkiRfX19NmjRJH3zwgX7//Xdbm5tvvlmS1LhxYy1dulSrV6/Whg0bdO+990qSHA6H7rnnHkVHR2vz5s22bW+77bZC9//000/rqquukiT5+flpyJAh+uKLLwrUOGrUKFWqVEmSNGXKFP3xxx96//33i/wei+ubb77Rxo0bdffdd0uy3uett96qL7/8Ur/++quOHDkil8ulAwcOSJJCQ0P16aeflnldZ8I1CQAAoFw6darwdU6n/XV8fOFtffL9SbQ0b9bjcDhk/eHYzhgjn/wHlhQWFmZ73bhxY/cFwrt27VLdunXdH2QlqUGDBrb2Pj4+Wrt2rW655RZlZmbK6XQqJiZGcXFxp61zx44dOnz4sKKiotzLTp06JT8/P2VmZmrXrl1yOp22+vIf+3Tybte4cWMZY7R7925dcskl7uV5n+fU5HQ61ahRI/eyJk2auNd16tSpSPtPTk7WjTfeqP3798vf319xcXEehxvl3Uf16tUVEhKimJiYIr/H4tqxY4d8fHw0bNgw97KsrCyFhYXp8OHD6tGjh0aNGqXLL79cvXr10o033qibbrqpzOs6E0ICAAAolypX9n7bM6latapOnjxZYHlSUpKqVatWYHn+QGGMcd/iM+/zwrz77rsaP368vv32W3Xr1k2SFBUV5TGo5BcREVHomPuibH86ebfPeZ7/vTjzJTtPx8zZJv+2he3/1KlTuvzyyzV8+HAtXbpUPj4+WrJkiebMmXPaGnNel+T2qmfr66+/LnAOcvznP//RjBkztGTJEs2aNUv/+te/tGnTJgUHB5+z+vJjuBEAAEAxRURE6Keffiqw/Mcff9Sll15aYHnOkJIc+/btU4sWLSRJl156qf744w/3hcae2q9du1YNGjRwBwRJysjIsLXJ24ORlpamzMxMtW7dWnv27FF2drZ7XXx8vCZNmuQ+tsvl0v79+ws99unkbbtv3z45HA41a9bstNu0bt1aLpdL+/btcy/bs2ePJOu8FmX/u3btUnx8vK677jr3+85/Pjzt48SJEzp27Jj73BdX3nN96tQpj8GndevWys7Odr+3HLfffrsSEhJ06NAhfffdd2rVqpWeeuop/fzzz/r999+9fockQgIAAEAxjRs3TsuWLdN///tf97Ljx4/rH//4h20irhzbt293j7ffvXu3vvjiC/d4/pEjR6pKlSpasGCBJCk1NVWLFy+2bd+qVSsdPHhQu3fvliT9+uuv2r59u61N7dq1dfz4cUnS9OnTtWrVKk2aNEkpKSl65ZVXJFl/Rf/nP/+p2rVrS5KuuOIKtWzZUs8884wkyeVyaeHChUU+Dy+//LIkaxjNiy++qGuuuabA8KL8+vTpo+7du+vJJ5901/TUU0+pX79+6tixo61tzrUT+fffqFEjVapUyf2B2uVy6eOPP/Z4vMWLF7sD2HPPPad69erp2muvLfJ79CTvue7SpYtOeRgjl/M+H3nkEXdIe/fdd7Vr1y6FhIRoz549mjFjhrKystznQZKaNm1aotpKrMSXPp9nuLsRAADly+nu1nI+mD9/vmndurWJiIgw3bt3N7169TJvvvlmgXZjxowxN910k5k+fbqJjIw0DRo0ME888YStzYYNG8xll11m2rVrZwYMGGCeffZZ95159uzZYzIzM83tt99u6tevb6688kozfvx4c9lll5mwsDDz5JNPGmOsOxl16tTJ9OjRwwwcONB9B59NmzaZnj17mnbt2pmePXua++67z2RlZbmPvWvXLtO9e3fTqlUr07dvX/Pyyy8bSaZLly5m3bp1Ht97zt2HPvnkE9O/f3/TpEkTc9VVV5k///zTGGPMF198Ydq0aeN+D++++65t+7i4OHPdddeZyy67zFx22WVm5MiR5ujRo8YYY5YtW+bedunSpaZv376mcePGtv0bY8yHH35omjVrZjp37myGDBlibr75ZhMQEGAuv/xys2vXLvcdop599lkzYMAA065dO9O2bVuzefNmY4wxP//8s7tNmzZtzLvvvmvGjh1rQkNDTWhoqBk3bpzHNsYYs27dOtO8eXPTvXt3M3PmzELbxcXFmeHDh5uWLVuaqKgoM3z4cHPkyBFjjDGHDx82Y8eONR07djRRUVGmU6dOZvHixYV+v52ruxs5jCnhILTzTFJSkoKDg5WYmOhxrCAAADi30tLSFBsbq/DwcAUGBnq7nDKT07NQ0tmQy5M1a9aoT58+Jb6mAUV3up+X0vycy3AjAAAAADaEBAAAgDJ28803a+XKlVq5cqX7nv/nu+XLl2vq1KmSrDssxcbGercglCpugQoAAFDGXnvtNW+XUOqGDx+u4cOHe7sMlBF6EgAAAADYEBIAAEC5wMWvwJmdq58TQgIAAPAqPz8/SVJKSoqXKwHKv5yfk5yfm7LCNQkAAMCrnE6nqlevrvj4eElSUFCQHA6Hl6sCyhdjjFJSUhQfH6/q1avL6XSW6fEICQAAwOvq1KkjSe6gAMCz6tWru39eyhIhAQAAeJ3D4VDdunV10UUXKTMz09vlAOWSn59fmfcg5CAkAACAcsPpdJ6zD0EACseFywAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbAgJAAAAAGwICQAAAABsCAkAAAAAbLweEhYsWKDw8HAFBgaqQ4cOWrdu3WnbL126VG3atFFQUJDq1q2rm2++WQkJCeeoWgAAAKDi82pIWL58uaZOnapZs2Zp69at6tWrlwYMGKADBw54bL9+/XqNHj1a48aN088//6x3331Xmzdv1vjx489x5QAAAEDF5dWQ8Mwzz2jcuHEaP368WrZsqXnz5ql+/fpauHChx/YbN25Uw4YNNXnyZIWHh6tnz5667bbb9P3335/jygEAAICKy2shISMjQ1u2bFG/fv1sy/v166cNGzZ43KZ79+76/ffftWLFChljdOTIEb333nsaNGhQocdJT09XUlKS7QEAAACgcF4LCUePHpXL5VJoaKhteWhoqOLi4jxu0717dy1dulTDhw+Xv7+/6tSpo+rVq+uFF14o9Dhz585VcHCw+1G/fv1SfR8AAABAReP1C5cdDofttTGmwLIcO3fu1OTJk/Xggw9qy5YtWrlypWJjYzVx4sRC93/fffcpMTHR/Th48GCp1g8AAABUNL7eOnCtWrXkdDoL9BrEx8cX6F3IMXfuXPXo0UP33HOPJOmyyy5T5cqV1atXLz3yyCOqW7dugW0CAgIUEBBQ+m8AAAAAqKC81pPg7++vDh06KDo62rY8Ojpa3bt397hNSkqKfHzsJTudTklWDwQAAACAkvPqcKPp06frlVde0eLFixUTE6Np06bpwIED7uFD9913n0aPHu1uf/XVV+uDDz7QwoULtW/fPv33v//V5MmT1blzZ9WrV89bbwMAAACoULw23EiShg8froSEBD388MM6fPiwIiIitGLFCoWFhUmSDh8+bJszYezYsTp58qRefPFF/d///Z+qV6+uyy+/XE888YS33gIAAABQ4TjMBTZOJykpScHBwUpMTFS1atW8XQ4AAABQKkrzc67X724EAAAAoHwhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALDxekhYsGCBwsPDFRgYqA4dOmjdunWnbZ+enq5Zs2YpLCxMAQEBaty4sRYvXnyOqgUAAAAqPl9vHnz58uWaOnWqFixYoB49eujf//63BgwYoJ07d6pBgwYet7n++ut15MgRvfrqq2rSpIni4+OVlZV1jisHAAAAKi6HMcZ46+BdunRR+/bttXDhQveyli1basiQIZo7d26B9itXrtSIESO0b98+1axZs1jHTEpKUnBwsBITE1WtWrVi1w4AAACUJ6X5Oddrw40yMjK0ZcsW9evXz7a8X79+2rBhg8dtPvnkE3Xs2FFPPvmkLr74YjVr1kx33323UlNTCz1Oenq6kpKSbA8AAAAAhfPacKOjR4/K5XIpNDTUtjw0NFRxcXEet9m3b5/Wr1+vwMBAffjhhzp69KjuuOMOHTt2rNDrEubOnauHHnqo1OsHAAAAKiqvX7jscDhsr40xBZblyM7OlsPh0NKlS9W5c2cNHDhQzzzzjJYsWVJob8J9992nxMRE9+PgwYOl/h4AAACAisRrPQm1atWS0+ks0GsQHx9foHchR926dXXxxRcrODjYvaxly5Yyxuj3339X06ZNC2wTEBCggICA0i0eAAAAqMC81pPg7++vDh06KDo62rY8Ojpa3bt397hNjx499Mcff+jUqVPuZb/88ot8fHx0ySWXlGm9AAAAwIXCq8ONpk+frldeeUWLFy9WTEyMpk2bpgMHDmjixImSrKFCo0ePdrcfOXKkQkJCdPPNN2vnzp369ttvdc899+iWW25RpUqVvPU2AAAAgArFq/MkDB8+XAkJCXr44Yd1+PBhRUREaMWKFQoLC5MkHT58WAcOHHC3r1KliqKjo3XXXXepY8eOCgkJ0fXXX69HHnnEW28BAAAAqHC8Ok+CNzBPAgAAACqiCjFPAgAAAIDyiZAAAAAAwKZYIeHrr78u7ToAAAAAlBPFCgm33nqr1q5dW9q1AAAAACgHihUSgoKC9Omnn6pr1656+OGHtX///tKuCwAAAICXFOsWqM8//7z69OmjrKwsffbZZ5oyZYrS0tI0evRoXXPNNQoMDCztOgEAAACcI8XqSejTp48kydfXV0OGDNH8+fPVrFkz3XTTTapbt65uu+02bd68uVQLBQAAAHBuFCskjBs3TpmZmXrvvfc0cOBANWzYUCtXrtQjjzyiHTt26JZbbtGCBQs0bdq00q4XAAAAQBkr1mRqISEh8vHxUVpamoYNG6ZbbrlFvXr1KtCuY8eO+v7770ul0NLCZGoAAACoiErzc26xrknw9/fXo48+quuvv15VqlTx2Gbu3LmKj48vUXEAAAAAzr1i9SS8+eabuummm2zLsrKyFB0drSuuuEJ+fn6lVmBpoycBAAAAFVFpfs4t1jUJixcvLrDM5XLps88+0zXXXFOiggAAAAB4V7FCgicBAQGaP3++EhMTS2uXAAAAALygyNckvP7663r99dclSdu2bdPll19eoM3x48cVEBBQetUBAAAAOOeKHBIaNmyoyMhISVJsbKz7eQ4fHx/Vrl1b1157belWCAAAAOCcKnJIiIyMdAeDatWqMQcCAAAAUEEV65qE0wWEBx98sNjFAAAAAPC+IvckbNmyRVWrVlWzZs308MMPF9ruzTffPO16AAAAAOVbkUPC0KFD1aJFC61atUrPPvus2rZt67HdiRMnSqk0AAAAAN5Q5JCwd+9eOZ1OSVLnzp315Zdfemx35ZVXlk5lAAAAALyiWDMun05KSoqCgoJKc5elihmXAQAAUBF5fcbl07nqqqtKe5cAAAAAzqEiDzdq1KhRkdrFxcUVuxgAAAAA3lfkkBAQEKCZM2eeto0xRk888USJiwIAAADgPUUOCbfffrvGjBlzxnZJSUklKggAAACAdxX5moTJkycXqd3Ro0eLXQwAAAAA72MyNQAAAAA2TKYGAAAAwIbJ1AAAAADYlPpkauUdk6kBAACgIirNz7lF7knIb8eOHXrxxRcVExMjh8OhFi1a6K677lKrVq1KVBAAAAAA7yrWjMtvvvmm2rZtq82bNys0NFS1a9fW5s2b1bZtWy1durS0awQAAABwDhWrJ+Ghhx7SN998o969e9uWr127VuPGjdONN95YKsUBAAAAOPeK1ZMQGhpaICBIUmRkpOrUqVPiogAAAAB4T7FCQkREhPbv319g+f79+xUeHl7iogAAAAB4T5GHG+WdIC04OFjdu3fXFVdcobCwMBljdODAAa1YsYKhRgAAAMB5rsi3QK1Ro0ahE6jltX37dh07dqykdZUZboEKAACAisgrt0Dt1q2bVqxYccZ2AwcOLFFBAAAAALyryNckFCUgSNJTTz1V7GIAAAAAeF+xJ1OTpPT0dMXHxyvviKUJEyZow4YNJS4MAAAAgHcUKyQcOnRIo0eP1tq1a1XESxoAAAAAnCeKFRImT56sPn36aMGCBRo7dqyWLVum9PR0vf/++0pKSirtGgEAAACcQ8UKCfHx8XrggQckSYGBgQoLC5Mk3XfffRo6dGjpVQcAAADgnCvWZGo+PrmbZWZmKiUlRZLkcrm0e/fu0qkMAAAAgFcUKyRUqVJFd999t06dOqXOnTurb9++evTRR9W/f3/Vrl27tGsEAAAAcA4VKyQ89thjCg8PV0ZGhh544AFVr15djz/+uE6ePKmFCxeWdo0AAAAAzqEiz7hcUTDjMgAAACoir8y4nF9mZqbef/99xcTEyOFwqGXLlrrmmmvk5+dXooIAAAAAeFexQsIPP/ygv//97/rjjz8UEhIiY4yOHTumiy++WB9//LHatWtX2nUCAAAAOEeKdU3ChAkTNGHCBB07dkzx8fH6888/dezYMY0bN07jx48v7RoBAAAAnEPFCgmS9OCDDyo4ONj9Ojg4WLNnzy6VogAAAAB4T7FCQlhYmDIyMgosT0tLU506dUpcFAAAAADvKfI1Cd9++637+aBBg9S/f3+NHj1aDRo0kDFGBw4c0OLFi3XdddeVSaEAAAAAzo0i3wI17yzLp92hwyGXy1WiosoSt0AFAABARVSan3OLPNwoMjJS2dnZZ3z07t27RAUBAAAA8K4ih4R///vfpdoOAAAAQPlU5GsSmjVrZnt98uRJvfTSS9qxY4ccDodat26t8ePHF2gHAAAA4PxS5GsS8vrxxx/Vt29fZWdnq2HDhpKk3377TU6nU9HR0WrdunVp11lquCYBAAAAFVFpfs4t1ozL06ZN0yOPPKJx48a5L2jOzs7WK6+8oilTpuibb74pUVEAAAAAvKdY8yScOnVKEyZMsN3xyMfHR7feequSk5NLrTgAAAAA516xQkJKSopSU1M9Lk9JSSlxUQAAAAC8p1jDjQYNGqSePXtq0qRJatKkiSRpz549WrBgga6++upSLRAAAADAuVWskPDoo4/Kx8dHd955p9LT02WMUWBgoKZNm6aHH364tGsEAAAAcA4V6+5GOdLS0rR3714ZY9S0aVMFBgaWZm1lgrsbAQAAoCLy+t2NatSooYiICK1bt04RERElKgAAAABA+VKsC5fr16+vb7/9trRrAQAAAFAOFCskNG/eXC6Xy+O6Bx98sEQFAQAAAPCuYg036t+/vwYPHqyRI0fq4osvltPpdK/77LPPuHgZAAAAOI8V68LlvJOoFdihw1FoL0N5wIXLAAAAqIhK83NusYYbRUZGKjs72+Ojd+/eJSoIAAAAgHedVUh45513NHz4cF1yySX65ptvPLZZtmxZqRQGAAAAwDuKHBJeeukl3Xjjjfrll1+0detW9evXT9HR0QXahYaGlmqBAAAAAM6tIoeE+fPna+3atdq6dat27Niht956S88++2xZ1gYAAADAC4ocEoKCgtS9e3f36+uvv17Hjx8vk6IAAAAAeE+RQ0KlSpWKtGzQoEElqwgAAACAVxV5noTDhw/rjTfeUN47psbFxRVYFhsbW7oVAgAAADinijxPwunmRrDtkHkSAAAAgHPOK/MknG5uBOZJAAAAACqOIoeEJ598slTbAQAAACifihwSOnXqVKrtAAAAAJRPZzXjMgAAAICKz+shYcGCBQoPD1dgYKA6dOigdevWFWm7//73v/L19VXbtm3LtkAAAADgAuPVkLB8+XJNnTpVs2bN0tatW9WrVy8NGDBABw4cOO12iYmJGj16tP72t7+do0oBAACAC0eRb4FaFrp06aL27dtr4cKF7mUtW7bUkCFDNHfu3EK3GzFihJo2bSqn06mPPvpI27ZtK/IxuQUqAAAAKiKv3AK1tGVkZGjLli3q16+fbXm/fv20YcOGQrd77bXX9Ouvv2r27NlFOk56erqSkpJsDwAAAACF81pIOHr0qFwul0JDQ23LQ0NDFRcX53GbPXv2aObMmVq6dKl8fYs2WfTcuXMVHBzsftSvX7/EtQMAAAAVmdcvXHY4HLbXxpgCyyTJ5XJp5MiReuihh9SsWbMi7/++++5TYmKi+3Hw4MES1wwAAABUZEX7c3wZqFWrlpxOZ4Feg/j4+AK9C5J08uRJff/999q6dasmTZokScrOzpYxRr6+vlq1apUuv/zyAtsFBAQoICCgbN4EAAAAUAF5rSfB399fHTp0UHR0tG15dHS0unfvXqB9tWrV9NNPP2nbtm3ux8SJE9W8eXNt27ZNXbp0OVelAwAAABWa13oSJGn69OkaNWqUOnbsqG7duumll17SgQMHNHHiREnWUKFDhw7pP//5j3x8fBQREWHb/qKLLlJgYGCB5QAAAACKz6shYfjw4UpISNDDDz+sw4cPKyIiQitWrFBYWJgk6fDhw2ecMwEAAABA6fLqPAnewDwJAAAAqIgqxDwJAAAAAMonQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAAAbQgIAAAAAG0ICAAAAABtCAgAAAACbCzYkvPWWtysAAAAAyqcLNiTcfrv0+uvergIAAAAofy7YkCBJN99MjwIAAACQ3wUbEsaOlYyRRo2S3n3X29UAAAAA5ccFGxKefdbqScjOlm64QfrwQ29XBAAAAJQPF2xI8PGRXn5ZuukmyeWShg+X1q71dlUAAACA9/l6uwBvcjql116TsrKkuDipY0dvVwQAAAB43wUdEiTJ11d64w0pI0MKCvJ2NQAAAID3XbDDjfLy9bUHhLlzpTVrvFYOAAAA4FWEhHzefFO6/37pqquk9eu9XQ0AAABw7hES8hk2TOrbV0pOlgYOlDZu9HZFAAAAwLlFSMgnMFD66COpTx/p5Empf3/p+++9XRUAAABw7hASPAgKkj79VOrVS0pMlPr1k7Zt83ZVAAAAwLlBSChE5crS559L3bpJx49bQ5COH/d2VQAAAEDZu+BvgXo6VatKX3whXXmlNHasVKOGtysCAAAAyh4h4QyCg627HPlypgAAAHCBYLhREeQNCAkJ0uDB0q+/eq8eAAAAoCwREs7SpEnWRc2XXy799pu3qwEAAABKHyHhLM2bJ7VoIR04YAWFgwe9XREAAABQuggJZyk0VPr6a6lJEyk21goKhw55uyoAAACg9Hg9JCxYsEDh4eEKDAxUhw4dtG7dukLbfvDBB+rbt69q166tatWqqVu3bvryyy/PYbWWevWkb76RwsOlvXulv/1Nios752UAAAAAZcKrIWH58uWaOnWqZs2apa1bt6pXr14aMGCADhw44LH9t99+q759+2rFihXasmWL+vTpo6uvvlpbt249x5VL9etbQaFBA2n3bmnMmHNeAgAAAFAmHMYY462Dd+nSRe3bt9fChQvdy1q2bKkhQ4Zo7ty5RdpHq1atNHz4cD344INFap+UlKTg4GAlJiaqWrVqxao7r19/tQLC669LjRuXeHcAAABAsZTm51yv3f0/IyNDW7Zs0cyZM23L+/Xrpw0bNhRpH9nZ2Tp58qRq1qxZaJv09HSlp6e7XyclJRWv4EI0biytWyc5HLnLjLG/BgAAAM4nXhtudPToUblcLoWGhtqWh4aGKq6IA/z/9a9/KTk5Wddff32hbebOnavg4GD3o379+iWq25O8geCzz6TevaXExFI/DAAAAHBOeP3CZUe+P7kbYwos8+Ttt9/WnDlztHz5cl100UWFtrvvvvuUmJjofhwsw3uWpqZKt91mzdDcv7908mSZHQoAAAAoM14LCbVq1ZLT6SzQaxAfH1+gdyG/5cuXa9y4cXrnnXd0xRVXnLZtQECAqlWrZnuUlUqVpM8/l2rUkDZulAYOlJKTy+xwAAAAQJnwWkjw9/dXhw4dFB0dbVseHR2t7t27F7rd22+/rbFjx+qtt97SoEGDyrrMs9a2rRQdLQUHWz0KV10lpaR4uyoAAACg6Lw63Gj69Ol65ZVXtHjxYsXExGjatGk6cOCAJk6cKMkaKjR69Gh3+7ffflujR4/Wv/71L3Xt2lVxcXGKi4tTYjm7AKBDB+nLL6WqVaU1a6S//11KS/N2VQAAAEDReDUkDB8+XPPmzdPDDz+stm3b6ttvv9WKFSsUFhYmSTp8+LBtzoR///vfysrK0p133qm6deu6H1OmTPHWWyhUly7SF19IlStLX30l5bnLKwAAAFCueXWeBG8o7XkSzmTtWumdd6Tnn5eczjI/HAAAAC5QFWKehAtFZKT1yOFySdnZkp+f92oCAAAATsfrt0C9kGRlSaNHS6NGWc8BAACA8oiehHNo61bp3XelzEyrJ2HJEoYgAQAAoPyhJ+Ec6tTJuj7B11d6801p/Hhr6BEAAABQnhASzrEhQ6S337Z6EJYskSZOJCgAAACgfCEkeMGwYdJ//iP5+EgvvyxNnixdWPeYAgAAQHlGSPCSkSOlxYslh0N65RVp505vVwQAAABYuHDZi8aMsW6JesklUqtW3q4GAAAAsBASvOyWW+yvExKkkBDv1AIAAABIDDcqV3bvli67TPrnP71dCQAAAC5khIRy5JtvpD/+kB58UHr8cW9XAwAAgAvVBRsSjn1zzNslFHD77dLcudbz++6T/vUv79YDAACAC9MFGxJ+Hvqzdt2yS5nHM71dis3MmdJDD1nP775beuEF79YDAACAC88FGxIkKe61OG2+dLP+/OBPb5di849/SLNmWc8nT5bmzfNqOQAAALjAXLAh4bJVlymoRZAy4jL087U/a8ewHUqPS/d2WZKsuRP++U/pnnus15s2ebceAAAAXFgu2JAQ3CVYHbZ2UINZDeTwdejo+0e1+dLNins9TqYcTH/scEhPPCE9/bR1fUKOPXuk556TkpO9VxsAAAAqNocpD5+Iz6GkpCQFBwcrMTFR1apVkySd3HZSu8ft1qkfTkmSavSroWb/bqZKDSt5s1SPxo2zZmquVUuaOlW6806penVvVwUAAABv8/Q5t7gu2J6EvKq2rar2/2uvRk80kk+gj46vOq7NEZv1+wu/y2SXrwzVo4fUqJF09Kj0wANSgwZWT8ORI96uDAAAABUFIeEvPr4+anBvA3Xc3lHBvYKVnZytvZP3amvvrUreVX7G9txyizXp2tKlUkSEdPKkNadCw4bWBc8AAABASRES8glqFqS2a9qq6YKmclZxKum/Sfq+zffa/9h+ZWdme7s8SZKvrzRypLR9u/Txx1KXLlJamuR0ersyAAAAVASEBA8cPg5dfPvF6vRzJ9UcUFMmwyh2Vqy2dNqikz+c9HZ5bj4+0uDB0nffWbM1T56cu27lSum666QffvBefQAAADg/ERJOI7BBoFp/3lot32wp3xBfJW9P1pbOW/TrzF/lSnV5uzw3h0Pq00eqWTN32WOPSe+9J3XoIA0YIK1b5736AAAAcH4hJJyBw+FQ6I2h6ryzs2oPry25pINPHNT3bb/XiXUnvF1eoRYulG66yRqCtHKl1Lu31KuX9MUX0oV1PysAAACcLUJCEflf5K9Wy1op4qMI+dfzV+ovqdrWe5t+ufMXZSVlebu8Alq1kt54Q/rlF2niRMnfX1q/Xho4ULrxRm9XBwAAgPKMkHCWav29ljr93El1J9SVJP2x4A9tjtishC8SvFyZZ40aWb0KsbHS3XdLlStLV1+duz49XcrI8F59AAAAKH8ICcXgV91PzV9qrjZft1Fgo0ClH0zXTwN/UszoGGUmZHq7PI/q1ZOeekrav9+6oDnHwoVSkybS889LKSneqw8AAADlByGhBGpcXkOdfuykS6ZfIvlIR944ok0tNyn+nXiV14msQ0KsW6jmeOst6eBBacoUa66Fxx6TTpzwVnUAAAAoDwgJJeSs7FSTfzVR+w3tFdQqSJl/Zmrn8J3aMXSH0v9I93Z5Z/Ttt9KiRVJ4uPTnn9KsWVJYmHT//VJ8vLerAwAAgDcQEkpJtS7V1PGHjgqbHSaHn0MJHydo06WbdPjVw+W2V0GSAgOl226zLnB+803rguekJGnuXKt3AQAAABceQkIp8vH3UficcHXY0kFVO1WVK9Gl3eN3a/sV25W6L9Xb5Z2Wr69116Mff5Q++kjq3Nm60DnHoUPS7t1eKw8AAADnECGhDFRpXUXtv2uvxv9qLJ9KPjrxzQltjtisg88elHGV314FyZrF+e9/lzZutCZiy/HII1LLltL110tbt3qvPgAAAJQ9QkIZcTgdqj+9vjr91EnV+1RXdmq2fp3+q37o8YNO7Tjl7fLOyOHIfW6MlJBgfX33Xal9e2u+hfXrvVcfAAAAyg4hoYxValxJbb5uo2YvNZOzmlMn/3dSW9pv0W8P/absjGxvl1ckDof0zjvWUKSRI63ehi++sGZw7t1b+uorb1cIAACA0kRIOAccDofqTainzjs7K+TqEJlMo9/m/KYtHbYoaXOSt8srstatpaVLrYucb73VmsV53Tpp9WpvVwYAAIDSREg4hwIuDlDExxG6dNml8qvtp+Qdyfqh6w/ae/deuVJc3i6vyBo3lv79b2nfPun//s9+F6QNG6SHH7auacjK8l6NAAAAKD6HKc/35ywDSUlJCg4OVmJioqpVq+a1OjKOZmjv1L2KX2pNRhDYKFDNX2muGn1qeK2m0nD55bk9C9WrW6/79bMe4eFeLQ0AAKBCK83PuYQEL0v4PEG/TPxF6b9bE6/VvbWuGj/ZWL7BvmfYsvzJzpZefVX68kvp668LztzcsqV1XYPv+ffWAAAAyj1CQgmUt5AgSVlJWdo3c5/+WPiHJMm/nr+aLWqmWlfX8nJlxedySd9/L0VHS6tWSd99J3XvLq1dm9tm3DjpkkusXobOnSU/P+/VCwAAcL4jJJRAeQwJOU6sPaHd43crda818dpFIy5Sk+ebyL+2v5crK7mTJ6UjR6QmTazXR49KF11k3VZVkqpWtYYm9e1rhYYmTey3YQUAAMDpERJKoDyHBElypbr025zfdPDpg1K25Bviq/B/hqv2sNoVIizkOHnSuq1qdLR1C9WEBPv6O+6Q5s/3Tm0AAADnI0JCCZT3kJAj6fsk7R63W8k/JlsLHFK1btUUcnWIal1dS0GXBslRQf7Unp1tzeK8apUVGtavlxYtkm65xVq/c6c0dqzVw9C3r9Stm3X7VQAAAOQiJJTA+RISJCk7M1uHnj+kI0uP6NRW+yzNgY0CrcAwuJaCewXLx6/i3M02+a9cVLmy9fXZZ6Xp03PXV64sRUXl3jWpeXOGJgEAABASSuB8Cgl5pR1MU8JnCUr4JEHHvzkuk5H7z+YMdipkQIhCrg5RzQE15VejYl0BHBdn3TEpOtp6xMfb10dHS1dcYT03hsAAAAAuTISEEjhfQ0JeWaeydDz6uBI+SVDC5wnK/DMzd6VTqt6rukKuDlHI4BAFNQnyXqFlIDtb+umn3KFJmzZJf/whBf31NmfMsJbnDE3q0UMKDPRuzQAAAOcCIaEEKkJIyMu4jJL+l6SETxN09NOjSvk5xbY+qEWQQgZbvQzB3YLlcFasP7NnZtpvndq2rbR9e+7rSpWk3r1zQ0NEBD0NAACgYiIklEBFCwn5pe5LtQLDJ0eV+G2iTFbuP69viK9CBv01LOnKmvKtWvFmNYuLs+6WlDM/Q1xc7ro6daxeh5yQcOqUVKWKd+oEAAAobYSEEqjoISGvzBOZOrbymBI+TdCxFceUdSLLvc7h51D1PtXdd0sKDKt4Y3KMkX7+OXdoUoMG0r//ba3Lzpbq1ZNq1pTatLF6GFq3tr42bCj5VJzrwAEAwAWCkFACF1JIyCs7M1uJ/01UwqfWxc85E7blqHxZZdUaXEshV4eoaseqcvhU7DE5u3dLLVp4XlelinTbbdLTT1uvjZH+/NOa/A0AAKC8Ks3PuRVvvAk88vHzUY2oGqoRVUONn26slN0p7sCQuCFRyT8mK/nHZO1/ZL/86/gr5CprWFKNK2rIGeT0dvmlrnlz6y5JmzZJO3ZYj59+kmJirGFIeedhOHJEqltXql07t7ch52urVtZs0QAAABUJPQlQxtEMHVvx17CklcfkOuVyr/MJ9FGNK2pYFz9fFaKAugFerLTsZWVJe/ZYczE0aGAtW7dOioy0ehQ8mTNHmj3bep6aKu3da4UQJnwDAADnEsONSoCQcHrZ6dk6sfaE++Ln9APptvVVO1Z13y2pSpsqFWbW5zNJTrZ6GXJ6HHK+Hj4svfpq7uzQ69dLvXpJvr5WUMjpdeB6BwAAUNYICSVASCg6Y4ySf0rW0U+OKuHTBJ3cdNK2PqB+gDUfw9UhqtGnhnwCLrxPvwkJVo9BzpCjDz6wAkNiouf28+dLd9xhPT982LqwunVr63qHCyRvAQCAMkJIKAFCQvGlH05XwucJSvg0Qcejjys7Ndu9zlnFqRr9aijk6hBV61pNlZpUko/vhRcaJGtY0u+/F+x1iImx7rTUu7fVbskS6eabree1ahW83qFNm9xJ4gAAAM6EkFAChITS4Upx6fg3f836/FmCMg5n2NY7/B0KahmkyhGVVTmisqq0rqLKEZUV0CDgghmilF/WX3eg9f3rdgFLlkiPPWZdw+Dpp3DlSunKK63n27dbYaNpU+taidBQeh4AAIAdIaEECAmlz2QbnfzhpBI+SdCxVceU/FOyslOyPbZ1VnWqcqvK7vCQ8/APvXCv8k1J8Xy9w6ZN0sUXW20eeEB69NHcbQICrLDQoIEUFmZdOJ1zoXXO3Zm4cBoAgAsLIaEECAllz2Qbpf2WpuQdybZHyq4UmUzP325+tf0KBIfKrSrLN5i79ErWJHBvvinFxlqzRuf/qY2NtS6KlqRZs6S5c63J4sLCcoNEzvOoKOvuTQAAoGIhJJQAIcF7sjOzlbonNTc4/GR9Tf01VSrkuzCgfkBuaGhtfQ1qESRnpYo3d0NRZWRIhw5J+/fnPu67T/Lzs9aPHSu9/nrh2+cNFAsWWMOackJE3jDBkCYAAM4vhIQSICSUP64Ul1JiUgr0PKT/nu55Ax+pUpNKBXoeKjW9cC+WzitnhuicAHHggP35//6XGyhGjbJ6KDwJCLDmjKhf33q9erW1j5wgccklDGkCAKA8ISSUACHh/JF5IrNAcEj+KVlZx7I8ts9/sXTOI7BBoBw+/Enck40bpa1bCwaJQ4esXoS0tMIDhcNhzUSdExoWLpSqV7fW/fmnFTKqVqU3AgCAc4WQUAKEhPObMUYZRzIKhocdycpOLuRi6SpOBbUKst1lqXJEZfld5HfB3mnpTDIzrXkcci6GlqSnnrJu4XrggPVIS8td53BI6ekFA0VgoDUHRGio9TXnMXu2VKmS1fbwYWuCuZCQ3Ds/AQCAs0dIKAFCQsVkso3S9nu4WDrmNBdL1/JTUIsgBdQPsB6X5HnUD5D/Rf70QBTCGCk+PrcH4s8/pdtvz10/aJC0YoXnbR0O67qKnEBwww3SsmXW8lq1coNETrB4/PHcQLF/v5SdbS3n4msAAOwICSVASLiwZGdmK3Wvh4ul9xZ+sXQOh69D/hf7u4NDYP3AgkEi1F8OJ0HCk5QUK0gcOWL/evKkdfelHNdcI330kee5IhwOq1fD+dd16iNGSMuXW88rVy4YKJ57LncCuj17rDBy0UVWL4UPl6sAACo4QkIJEBIgSa5U62Lp1D2pSv893f1IO5im9N/TrcnhPI9esnNKAfXswaFAkKjjzwXVZ+BySUePWiEif6D4xz9y291wgxUo8g51yuHjY4UCT4HCx0eqXdseKF5+OTdQ/PyzdawaNXIfOUOnAAA4XxASSoCQgKLIzspWRlyG0g+m20JE+u/pucv+SJdcRdiZj+Rf199zb0ROkKjrLx8/gkRRGGNNGJc/UCQlSXffndtu7Fjps8+khISC+3A6rUCR07swfLj0zjv2NpUr5waG777LHd707rvSrl32QJH3Ubs2vRYAAO8ozc+5XCYIeODj66PASwIVeElgoW2My7qIOn+QyOmNSP89XRmHMmSyjDIOZSjjUIZO/u+k5505JP86/p57I3KW1QuQjz+fPh0O665JVatKjRsX3m7JEutrZmZuL0XeQJH3g3z16tbcEcePS4mJ1rLkZOvxxx+510RIVph4773Cj3vihBQcbD1/+GEpOtoeImrWzH1+7bW5+z51yrpOI7DwbzkAAM4ZehKAMmSyjTLiPfRI5Htd2MXV+fnW8JVfLT/5hfjJr5affEP+ep1nWd6vvjV96aE4Sy6XFRSOH7ceSUnS5Zfnrn/lFWnTptz1eR+nTll3ecoJINddV/RAceut1hCowEDPPRTPPWd9laQffrAuGs8JS9Wq5T6vUoWeDAC4UDHcqAQICShvTLZR5p+ZBYY15e2RSP89XSa9eD+qzmBnoSHCFjTyLKfHoniys+0f0H/4Qfr1V8+B4sQJa7br4gSKCROssFKYAwdyJ8F77jnp4489B4qqVaXRo3PDx++/W8fJ24bb0gLA+YPhRkAF4vBxyD/UX/6h/qraoarHNsYYZR7NVOafmcpMyLSe//U1KyHL9tq9/HiWZCRXokuuRJfSfvVwtW8hnFWdpw0RnkKGM9BZWqfkvJX/L/jt21uPoli+3Lp4OidEHDtmDxVV83xrhIVJ3bpZvRwnT1qPpCSrF0Syt/35Z2u27MIMGZIbEp57Tnr6afv6wMDc0PDFF1LTptbyDz+0XnsKHtWqST17Wl8lKTXVqi0oiF4OADhfEBKA84DD4ZB/bX/51/Yv8jbGZZR5vPAQ4TFkJGRK2ZLrpEuuky6l/Vb0YOET5GMPDzV95RtsPZzVnPKtlvu1wLJgX/lU8rmgJ7fz8bF6CoKDresjTueBB6xHXsZYd31KSsrtcZCkiROlPn3sYSLnec4dnXL4+1tzVZw8aQ2bkqx9pqVZ13LkvePTd99Zw6MKs327dNll1vOnn5YefNB6XqmSFRYqV7YeQUHSa69JrVtb61etsu5glbMub7vKla33Uru21TanRyZnHSEEAEoPIQGooBxOh/xr+cu/lr/UvGjbmGyjrBNZRe6tyEyw1pkso+yUbKUfSFf6gfTiFezU6YNENV85g/O9/itg5F3mE3Rhhg2Hw/oAnvcia+nsejMefdR6SNYF3/lDRb16uW2vvNLqNcgfOnJe5w0fKSm5z1NTrUfeu0658twlbPNmaeHCwmv89tvckPDGG9KUKfb1gYG5weI//5EiI63lX30lLVpUMHTkPB80SAoPt9rGx0t791r7yjmnOY/AQCssXYDfYgAuMIQEAG4OH4f8avrJr6af1LRo2xhj5EpyFQwRCZlyJbmUlZRl/5qYVWC5siW5pKzjWco6nqV0FTNoSAXCRlFDh7OKU87KuQ+fIB85g5wX7GR5fn7WnZhq1vS8/m9/sx5F8eijVs9HcrIVGHLuHJXzPO9dqnr1snodPLVLTrZ6OnIYY33AzxtCcno+EhLsE/Tt3i29/37hNTZsmBsSVq6UxowpvO1bb1lzduS0nTzZHiLyPp84Uerd22q7d691dyxP7SpVklq1ki6+OPd9HDuWuy4wkGAC4NwiJAAoEYfD4R5WVKlxpTNvkI8xRq5kV5EDhadlrsQyCBt/8Qn0kU9lKzA4Kzut55WdcgbleZ43WFTO0+4M21woQ6x8fHL/an8mvXvnfqg+kylTrEd2ttU7kT9YNM/Tg9a7t/Tii4UHlZwLvSXrQ3njxrm9Hmlp1te863MkJFizexdm4MDc5zt3SrNmFd524UIrVEjShg0FQ1hAQG6wmDPHuiNWzn4nTbKCREBA7iPn9dVXS337Wm2PHrVCTt71eduHh+eGpcxM6fBh+3p/f4Z0ARcKQgIAr3I4HPKt4ivfKr4KqBdQ7P0YYw15yhswihw6El1WUPnrkZ2SLf31V+jstGxlp2UrKyGrlN6xnS1Y5A8ZnoLFX218KvnIJ9BHzkpOK8j89dr2PE+bitwjkjeE5AxFyq9169zrHs7kuuusR17GWNdppKbmztQtSf36SevW5QaKvKEiNVXq2DG37SWXSOPGeW6XlmbNBp4jZ/bwvEOx0tOtx4kT1vocR4+e/uL0evVyQ8L+/QWHaOV1//25Q85iY+1BK4efnxUaJk/ObXvkiNS/vz145H1+xRW5vTOpqdLjj1uBw9OjaVOpa1erbXa2tH594W2rVMm9QB5A6SIkAKgQHA6H+4O06pZsX8YYZadlW4Eh2frqSsnzPE+YcD9PzpYrxWXfJl+bnNfZqdnuY2WnZCs7JVuZf2aW8AycnsPX4TFA5A8TRQkcha3z1M7h76gQvSUOh/WhN/9kd7VrFx5M8mvf/vS3rs2rf38pK8v6a37+MJGaar8+pEUL6e23rXU5QSI9Pfd19+65bYODrRnG867P2z7vfjMzrQ/46fk65DIzrUdWntycnCxt21b4+6lWLTcknDxpTTRYmDFjckNCamrudSWeDBtmzYIuWUEu55oRT4GiTx9p/vzcbYcMsUJIzvqAgNznLVpId96Z2/aFF6z37O9v7T/v46KL7HOpbNpkhbucOvK2rVTJap/D5bJCbgX4EUEFREgAgHwcDoeclZxyVnJKtc7c/myZbGMPHSlnCBYeAkhOD0d2qvXVlVpwmcnIHZRvsoxcp1xynXKdprIy4JB8AnzkCHDIJ8DHeu6f+zz/a0eAQz7+PmW/zXnSs5Lz4bKq57sjS7I+dI4YUbT9NWkiLVtWtLatWlnBwZjcsJI3VOStKTTUuiWup5CSni61aZPb1t/f+gCekWE90tNzn2dkSBERuW1dLusDe971ebfxz3PDt6ys3PXJyQXfT7Nm9tcrVljvy5O//c0eEh580OrB8aRLF2njxtzX11wjHTrkuW1EhPTTT7mvL71U+uWX3H/nvKGicWPrQv0cY8ZY17XkDR057WvVsl/wP2+e1Wvk65vbNud5UJDVC5Rj1SqrJyhnfd6vOeEqx2+/WefWU1s/PyuEEngqDkICAJxjDh9riJWqlO1xjMsoOz03NNjCRKo9UJwpcHhqX2BfeZ7nDNeSsYZsKU1y6RwHlNNxqmCw8M8NHA5/h3z8/uoJ8ftrmV/uugLL/Byn3ybf+rPexs8hh493Pn05HLl/YS9M5cpW70dRVK9uXR9SFNWqSTExha/Pe3G6r681IWD+QJHzqF7dvu0rr3gOHhkZuddl5BgxwuoBycjI7UnJebRsaW/bsKHVK5G3bc7z/D1ROSElp13ei/Cr5Pv/YetWe8DI6+KL7SFh+XJ7cMmrenV7SHjqKevuX574+dmHtk2ZIn3yiee2ktU251bJN91kXajvKUz06WPdfQzlGyEBACooh9NhXTwddG4nujPGyGQYW5gwGX8FlvRsmXSj7Iw8z9OzlZ2R53l6tq19/tfF3cbGJSvQ5Bn6Ve79FWzyhxWHb54g4fvXsrzL8y0rdJu/vpZ4mzPV4PyrTd5HMQNQ3r9aOxy5d4cqitGji972dLflzW/9+qK3/eEHK5zkDxOZmQVnO3/+eWtukLwBJad9/lsfjx0rRUXlDlnL+zV/UOnY0br+Je8wspyvznz/dVSrZg2vy98u59qZvO3zvpf88t4CGeWXwxhjztys4ijN6aoBAOcHY4w1n0f6GYJFulF2prXcZFphJu9Xk1H09cXZxr0844L61Sw5lBsYPIWIPA85ddr1Z9q+SG2c+Y7lzLON86/leV7b2udbftp9FGXfXupBOhvZ2VZYyNvbdPy4NTTJU1CpUqXg8C+UjtL8nEtPAgCgwnM4cv/iXdbDvEqDMUbGZYoWLLKsNibrr2WZeZblXZ5vWaHb/PW1xNt4qKHQEWdG7rYoqMgBJG/IcDokHxW6Lv/6s1rntIZNFmWd0+mQb/79+DiUsNdfIQNDvH1qcRqEBAAAyhmH46+/TFew39Lu8JNlJJd1QX2hD9dp1pXG+qLW4Mp97m5fyst1hlFvJstIWZJRxQlRwb2CCQnlXAX77wcAAJRXFTX8lFROeCpW2HDZA4172V/hw73NWazLv/6s1rmsO7idaV3lVkWYXRFexY8pAACAF9nCU/HnlARKldcnV1+wYIHCw8MVGBioDh06aN26dadtv3btWnXo0EGBgYFq1KiRFi1adI4qBQAAAC4MXg0Jy5cv19SpUzVr1ixt3bpVvXr10oABA3TgwAGP7WNjYzVw4ED16tVLW7du1f3336/Jkyfr/fffP8eVAwAAABWXV2+B2qVLF7Vv314L89yAuGXLlhoyZIjmzp1boP2MGTP0ySefKCbPzCoTJ07U9u3b9d133xXpmNwCFQAAABVRhbgFakZGhrZs2aKZM2falvfr108bNmzwuM13332nfv362ZZdeeWVevXVV5WZmSm/nGn+8khPT1d6err7dWJioiTrJAIAAAAVRc7n29LoA/BaSDh69KhcLpdCQ0Nty0NDQxUXF+dxm7i4OI/ts7KydPToUdWtW7fANnPnztVDDz1UYHn9+vVLUD0AAABQPiUkJCg4OLhE+/D63Y0cDvtMgsaYAsvO1N7T8hz33Xefpk+f7n594sQJhYWF6cCBAyU+eSgoKSlJ9evX18GDBxnOVco4t2WHc1t2OLdlh3Nbtji/ZYdzW3YSExPVoEED1axZs8T78lpIqFWrlpxOZ4Feg/j4+AK9BTnq1Knjsb2vr69CQjxPyBEQEKCAgIL3EwsODuYbswxVq1aN81tGOLdlh3Nbdji3ZYdzW7Y4v2WHc1t2fHxKfm8ir93dyN/fXx06dFB0dLRteXR0tLp37+5xm27duhVov2rVKnXs2NHj9QgAAAAAzp5Xb4E6ffp0vfLKK1q8eLFiYmI0bdo0HThwQBMnTpRkDRUaPXq0u/3EiRO1f/9+TZ8+XTExMVq8eLFeffVV3X333d56CwAAAECF49VrEoYPH66EhAQ9/PDDOnz4sCIiIrRixQqFhYVJkg4fPmybMyE8PFwrVqzQtGnTNH/+fNWrV0/PP/+8rr322iIfMyAgQLNnz/Y4BAklx/ktO5zbssO5LTuc27LDuS1bnN+yw7ktO6V5br06TwIAAACA8serw40AAAAAlD+EBAAAAAA2hAQAAAAANoQEAAAAADYXXEhYsGCBwsPDFRgYqA4dOmjdunXeLum8N3fuXHXq1ElVq1bVRRddpCFDhmj37t3eLqtCmjt3rhwOh6ZOnertUiqMQ4cO6aabblJISIiCgoLUtm1bbdmyxdtlnfeysrL0wAMPKDw8XJUqVVKjRo308MMPKzs729ulnXe+/fZbXX311apXr54cDoc++ugj23pjjObMmaN69eqpUqVKioqK0s8//+ydYs8zpzu3mZmZmjFjhlq3bq3KlSurXr16Gj16tP744w/vFXweOdP3bV633XabHA6H5s2bd87qO98V5fzGxMRo8ODBCg4OVtWqVdW1a1fbXUPP5IIKCcuXL9fUqVM1a9Ysbd26Vb169dKAAQPO6oShoLVr1+rOO+/Uxo0bFR0draysLPXr10/JycneLq1C2bx5s1566SVddtll3i6lwjh+/Lh69OghPz8/ffHFF9q5c6f+9a9/qXr16t4u7bz3xBNPaNGiRXrxxRcVExOjJ598Uk899ZReeOEFb5d23klOTlabNm304osvelz/5JNP6plnntGLL76ozZs3q06dOurbt69Onjx5jis9/5zu3KakpOiHH37QP/7xD/3www/64IMP9Msvv2jw4MFeqPT8c6bv2xwfffSR/ve//6levXrnqLKK4Uzn99dff1XPnj3VokULrVmzRtu3b9c//vEPBQYGFv0g5gLSuXNnM3HiRNuyFi1amJkzZ3qpooopPj7eSDJr1671dikVxsmTJ03Tpk1NdHS0iYyMNFOmTPF2SRXCjBkzTM+ePb1dRoU0aNAgc8stt9iWXXPNNeamm27yUkUVgyTz4Ycful9nZ2ebOnXqmMcff9y9LC0tzQQHB5tFixZ5ocLzV/5z68mmTZuMJLN///5zU1QFUdi5/f33383FF19sduzYYcLCwsyzzz57zmurCDyd3+HDh5f4/9sLpichIyNDW7ZsUb9+/WzL+/Xrpw0bNnipqoopMTFRklSzZk0vV1Jx3HnnnRo0aJCuuOIKb5dSoXzyySfq2LGjrrvuOl100UVq166dXn75ZW+XVSH07NlTX3/9tX755RdJ0vbt27V+/XoNHDjQy5VVLLGxsYqLi7P9bgsICFBkZCS/28pAYmKiHA4HvY2lIDs7W6NGjdI999yjVq1aebucCiU7O1uff/65mjVrpiuvvFIXXXSRunTpctohX55cMCHh6NGjcrlcCg0NtS0PDQ1VXFycl6qqeIwxmj59unr27KmIiAhvl1MhLFu2TD/88IPmzp3r7VIqnH379mnhwoVq2rSpvvzyS02cOFGTJ0/Wf/7zH2+Xdt6bMWOGbrjhBrVo0UJ+fn5q166dpk6dqhtuuMHbpVUoOb+/+N1W9tLS0jRz5kyNHDlS1apV83Y5570nnnhCvr6+mjx5srdLqXDi4+N16tQpPf744+rfv79WrVqloUOH6pprrtHatWuLvB/fMqyxXHI4HLbXxpgCy1B8kyZN0o8//qj169d7u5QK4eDBg5oyZYpWrVp1duMIUSTZ2dnq2LGjHnvsMUlSu3bt9PPPP2vhwoUaPXq0l6s7vy1fvlxvvvmm3nrrLbVq1Urbtm3T1KlTVa9ePY0ZM8bb5VU4/G4rW5mZmRoxYoSys7O1YMECb5dz3tuyZYuee+45/fDDD3yfloGcG0T8/e9/17Rp0yRJbdu21YYNG7Ro0SJFRkYWaT8XTE9CrVq15HQ6C/xlJT4+vsBfYFA8d911lz755BOtXr1al1xyibfLqRC2bNmi+Ph4dejQQb6+vvL19dXatWv1/PPPy9fXVy6Xy9slntfq1q2rSy+91LasZcuW3MygFNxzzz2aOXOmRowYodatW2vUqFGaNm0aPWKlrE6dOpLE77YylJmZqeuvv16xsbGKjo6mF6EUrFu3TvHx8WrQoIH7d9v+/fv1f//3f2rYsKG3yzvv1apVS76+viX+/XbBhAR/f3916NBB0dHRtuXR0dHq3r27l6qqGIwxmjRpkj744AN98803Cg8P93ZJFcbf/vY3/fTTT9q2bZv70bFjR914443atm2bnE6nt0s8r/Xo0aPA7Xp/+eUXhYWFeamiiiMlJUU+PvZfMU6nk1uglrLw8HDVqVPH9rstIyNDa9eu5XdbKcgJCHv27NFXX32lkJAQb5dUIYwaNUo//vij7XdbvXr1dM899+jLL7/0dnnnPX9/f3Xq1KnEv98uqOFG06dP16hRo9SxY0d169ZNL730kg4cOKCJEyd6u7Tz2p133qm33npLH3/8sapWrer+i1ZwcLAqVark5erOb1WrVi1wbUflypUVEhLCNR+lYNq0aerevbsee+wxXX/99dq0aZNeeuklvfTSS94u7bx39dVX69FHH1WDBg3UqlUrbd26Vc8884xuueUWb5d23jl16pT27t3rfh0bG6tt27apZs2aatCggaZOnarHHntMTZs2VdOmTfXYY48pKChII0eO9GLV54fTndt69epp2LBh+uGHH/TZZ5/J5XK5f7/VrFlT/v7+3ir7vHCm79v8gcvPz0916tRR8+bNz3Wp56Uznd977rlHw4cPV+/evdWnTx+tXLlSn376qdasWVP0g5To3kjnofnz55uwsDDj7+9v2rdvz206S4Ekj4/XXnvN26VVSNwCtXR9+umnJiIiwgQEBJgWLVqYl156ydslVQhJSUlmypQppkGDBiYwMNA0atTIzJo1y6Snp3u7tPPO6tWrPf4fO2bMGGOMdRvU2bNnmzp16piAgADTu3dv89NPP3m36PPE6c5tbGxsob/fVq9e7e3Sy70zfd/mxy1Qz05Rzu+rr75qmjRpYgIDA02bNm3MRx99dFbHcBhjzNnnFwAAAAAV1QVzTQIAAACAoiEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJAAAAACwISQAAAAAsCEkAAAAALAhJADAORYdHa2oqCg5HA5NmDDBti4xMVFRUVGqXr26unbtqjfeeKPM6ti9e7e7jjVr1pTZcYrjmWeeUdu2bdWpUyf16NGjwPo77rhDDRs2VGBgoKKiohQbGytJ+uijj/TRRx+d42oBoOJhxmUA8BKHwyFJ+vTTT3XVVVfZ1kVFRWnJkiVq2LDhOalj9erVioqKKvNjFcVvv/2m8PBwxcbGqmHDhnriiSc0Y8aMAu3mzJmjJUuW6LfffnMvGzt2rCRpyZIl56ZYAKig6EkAAC8JCwtT//79NWHCBB09etTb5ZQb+/fvlyR3QPIUEAAAZYuQAABe9NprrykrK0u33357oW2WLFmiFi1a2HoVBgwYoMDAQPdfzNevX6+uXbvK4XBo2bJlGjp0qJo0aaJJkyYpLS1NU6dOVdeuXdW1a1fbX95z7N27V0OHDlX79u3Vrl07bdmypUAN7dq1U69evdS9e3d9+OGH7nVXXXWVqlevrnvvvVe33367evXqJYfDoW3btnl8P3v27NGAAQPUoUMHtW7d2l2jJL333nuaMmWKJKs35Wx6N+69916tXLlSK1euVFRUlP7+97+7133++efq3LmzevbsqW7dumnRokXudePHj1edOnU0atQozZgxQ7169VKrVq20efNmffPNNxoyZIiaNGmixx9/3Ha8efPmqWPHjurTp4+6d++uefPmFblWACj3DADAK8LCwowxxnz44YdGknnjjTfc6yIjI01sbKz79WuvveZun3f71157zf06NjbWSDJ33XWXMcaYY8eOmcqVK5vhw4eb+Ph4Y4wxI0aMMGPHjrXtR5KJjIw0KSkpxhhj5syZYy6++GKTlpZmjDFmxYoVJiQkxBw8eNAYY8wvv/xigoKCzIYNG2z11q9f3xw4cMAYY8z48ePNjz/+WOA9p6WlmfDwcPPII48YY4xJT083kZGRZsKECe42q1evNkX59TR79uwC52TMmDFmzJgxtmU//fSTqVSpktm2bZsxxpgjR46YevXqmbffftu2Xc2aNc2ePXuMMcbcf//9JiIiwsyfP98YY0xMTIzx8fEx+/btM8YY87///c9UqVLFnDhxwr2+cePGZ6wZAM4X9CQAgJcNGTJE48aN01133aXff/+9xPu7/vrrJUk1atTQpZdeqipVqqh27dqSpJ49e2rr1q0Fthk1apQqVaokSZoyZYr++OMPvf/++5KkuXPnasSIEbrkkkskSU2bNlWfPn20YMEC2z7+9re/qX79+pKkl19+Wa1bty5wnLfeekt//PGHpk6dKkny9/fX1KlT9eqrr+rIkSMlfu+ePPnkk7r88svVpk0bSdJFF12koUOHav78+bZ27dq1U5MmTSRJPXr00I4dOzR48GBJUosWLVSzZk1t375dknTo0CFlZmYqLi7Ovf71118vk/oBwBt8vV0AAEB67rnntHbtWt1yyy368ssvS7SvunXrup8HBQXZXleuXFmJiYkFtgkLC3M/r169ukJCQhQTEyNJ2rFjhw4dOmQb+nP06FEFBgba9pETIk5nx44dqlu3ripXruxe1qRJE2VnZ2vnzp0KDQ098xs8Szt27NCRI0ds9Z84caJA/fnPW/5lec/dgAED1KtXL0VERKhfv3666aabNGzYsFKvHQC8hZAAAOVA5cqV9eabb6pnz54F/sIt5d4JKS+Xy+VxX06n87SvjYeb2uVfZoxxH9PhcOimm27SQw89dNr3kP84nng6dg5P77G0XHHFFWf8S7+n+gs7d4GBgYqOjtb//vc/LVmyRLfeeqvmz5+vNWvWyNeXX60Azn8MNwKAcqJLly564IEHNGPGDB04cMC2rmrVqjp16pT7dWZmpuLj40vt2HmPd+LECR07dkwtWrSQJEVERGj37t229qtXr9bChQvP+jitW7fW4cOHlZyc7F62d+9eOZ1OtWzZspjV5/Lxyf21lpKSIpfLpdatWxeof8eOHXr44YeLfZxdu3Zpx44d6tKlixYuXKiNGzfqv//9r3s4EgCc7wgJAFCOPPDAA7rsssvck4PlaNOmjY4dO+b+sLt06VLbB+KSWrx4sVJTUyVZQ5/q1auna6+9VpI0a9YsffLJJ+4PwMnJybr//vvdIeJsjBw5UvXq1XPfCSgjI0Pz5s3TuHHjSmWoUe3atXX8+HFJ0rBhw7Rr1y7NmDFDP/zwg1atWiXJClj/+Mc/bEOsztbGjRv12GOPuXsWXC6XAgICSrRPAChXvHnVNABciFatWmUiIyNNQECAiYyMND///LNt/d69e02VKlVsdzcyxphHHnnENGnSxPTr18+88sorJiwszDRv3ty88MILZuvWraZLly5GkunSpYv5+eefzahRo0xwcLAJCwsz//rXv8zSpUtN8+bN3cfdtWuXiYyMNJLMs88+awYMGGDatWtn2rZtazZv3mw79htvvGFat25tunXrZnr06GHefPNN97rhw4e7jzNo0KAzvv9ffvnFXHnllaZ9+/amVatW5o477nDfWendd981bdq0cd9xad68eR73cfvtt5uwsDD3e8m561BMTIyJiIgwPXv2NDfddJO7/cqVK02HDh1Mp06dTI8ePcwzzzzjXjdlyhQTGhpqQkNDzT333GO+/vprWw0JCQmmb9++JiAgwDRv3ty89tprZvfu3WbYsGGmS5cuJioqynTt2tV8+umnZ3zvAHC+YMZlAAAAADYMNwIAAABgQ0gAAAAAYENIAAAAAGBDSAAAAABgQ0gAAAAAYENIAAAAAGBDSAAAAABgQ0gAAAAAYENIAAAAAGBDSAAAAABgQ0gAAAAAYPP/16cJZmBrDr0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.font_manager import FontProperties\n",
    "\n",
    "# from math import log2\n",
    "\n",
    "font = FontProperties()\n",
    "font.set_family('serif')\n",
    "# font.set_name('Times')\n",
    "font.set_size(10)\n",
    "\n",
    "# plt.rcParams.update({\n",
    "#     # 'font.size': 12,\n",
    "#     'text.usetex': True,\n",
    "#     'font.serif': 'Times New Roman',\n",
    "#     })\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(9,6 ))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(probs, 'm', label='Original probabilities')\n",
    "ax.plot(new_probs, 'b--', label='Updated probabilities')\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_xlim([0, num_items])\n",
    "# ax.set_xscale('log', base=2)\n",
    "# ax.set_xticks([0, 2, 4, 8, 16])\n",
    "# ax.set_xticks([0, 2, 4, 8, 16, 32, 64, 128, 256])\n",
    "\n",
    "\n",
    "ax.legend(loc='upper right')\n",
    "ax.set_xlabel('Number of Items', fontproperties=font)\n",
    "ax.set_ylabel('Probability', fontproperties=font)\n",
    "# ax.set_ylabel(r'$\\int\\ Y^2\\ dt\\ \\ [V^2 s]$')\n",
    "ax.set_title('Probability distribution', fontproperties=font)\n",
    "plt.legend(prop=font, loc='upper right')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# fig.savefig('plot.pdf', dpi=300)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.75, 1.75)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''### Example 2.1.2[Cover, Thomas, p:14]'''\n",
    "\n",
    "from entropy_formula import entropy_formula\n",
    "num_items = [0, 1, 2, 3] # values of the random variable X\n",
    "probs = [1/2, 1/4, 1/8, 1/8]  # probabilities of the values\n",
    "\n",
    "entropy_shannon = stats.entropy(probs, base=2) # by python stats module\n",
    "entropy_shannon, entropy_formula(probs)\n",
    "\n",
    "## this means that to encode the information of the random variable X = number of items in (xk), we need at least 2 bits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Updated probabilities after removing index 1 : [0.6666666666666666, 0.16666666666666666, 0.16666666666666666]\n"
     ]
    }
   ],
   "source": [
    "class RenormalizeDistribution:\n",
    "    def __init__(self, probabilities):\n",
    "        self.probabilities = probabilities\n",
    "\n",
    "    def remove_and_renormalize(self, index_to_remove):\n",
    "        # Check if the index is valid\n",
    "        if index_to_remove < 0 or index_to_remove >= len(self.probabilities):\n",
    "            print(\"Invalid index to remove.\")\n",
    "            return\n",
    "\n",
    "        # Remove the probability at the specified index\n",
    "        removed_probability = self.probabilities.pop(index_to_remove)\n",
    "\n",
    "        # Calculate the normalization factor\n",
    "        normalization_factor = sum(self.probabilities)\n",
    "\n",
    "        # Renormalize the probabilities\n",
    "        if normalization_factor > 0:\n",
    "            self.probabilities = [p / normalization_factor for p in self.probabilities]\n",
    "\n",
    "        # Print the updated probabilities\n",
    "        print(\"Updated probabilities after removing index\", index_to_remove, \":\", self.probabilities)\n",
    "\n",
    "# Example usage:\n",
    "initial_probabilities = [1/2, 1/4, 1/8, 1/8]\n",
    "distribution = RenormalizeDistribution(initial_probabilities)\n",
    "\n",
    "# Remove the first probability (index 0) and renormalize\n",
    "distribution.remove_and_renormalize(1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n",
      "findfont: Font family 'Times New Roman' not found.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEbCAYAAABUV7o5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIRUlEQVR4nO3deVzU1f748dfMsG/DjiDIoriksohKmKZpWlZmm1mW2XItb/kzo+WmfcuyunbL27VuZZnlzbJr3WvdNjOt3ModAXfEREHZ932AYX5/DDNJoDI4CwPv5+Mxj+DDZz6f98d03pxz3ucchU6n0yGEEELYIaWtAxBCCCE6S5KYEEIIuyVJTAghhN2SJCaEEMJuSRITQghhtySJCSGEsFuSxIQQQtgtSWJCCCHsloOtA+jqmpubyc3NxdPTE4VCYetwhBCi29PpdFRVVRESEoJSeeG2liSxi8jNzSUsLMzWYQghRI+Tk5NDaGjoBc+RJHYRnp6egP4P08vLy8bRCCFE91dZWUlYWJjx8/dCJIldhKEL0cvLS5KYEEJYUUeGcKSwQwghhN2SJCaEEMJu2VUS27ZtG1OmTCEkJASFQsH//ve/i75n69atJCQk4OLiQlRUFO+++67lAxVCCGEVdpXEampqiI2N5a233urQ+VlZWVx33XWMGTOG1NRUFi5cyLx581i3bp2FIxVCCGENdlXYMXnyZCZPntzh899991369OnDsmXLABg0aBD79u1j6dKl3HrrrRaKUgghhLXYVUvMVDt37mTSpEmtjl1zzTXs27ePxsZGi9+/rkHL5oxC1u7Jtvi9hBCiJ7Krlpip8vPzCQoKanUsKCiIpqYmiouLCQ4ObvMejUaDRqMxfl9ZWdnp+/9WVM19q/bi7qTi1oRQHFXd+ncGIYSwum7/qfrHeQY6na7d4wZLlixBrVYbX5eyWsdlwV74ujtR06AlLae809cRQgjRvm6dxHr16kV+fn6rY4WFhTg4OODn59fuexYsWEBFRYXxlZOT0+n7K5UKRvXV32d7ZnGnryOEEKJ93TqJJSUlsWnTplbHNm7cyPDhw3F0dGz3Pc7OzsbVOcyxSseYaH8AfsksuqTrCCGEaMuuklh1dTVpaWmkpaUB+hL6tLQ0srP1hRMLFizgnnvuMZ4/Z84cTp8+TXJyMkePHuXDDz/kgw8+4IknnrBazKOjAwBIP1NBZb3li0mEEKInsasktm/fPuLj44mPjwcgOTmZ+Ph4nnvuOQDy8vKMCQ0gMjKS9evXs2XLFuLi4njxxRd58803rVpe39vblUh/d7TNOnb9VmK1+wohRE+g0BkqHUS7KisrUavVVFRUdLpr8dn/HeLjXae5JymcxVOHmDlCIYToXkz53LWrlpi9Gm0cF5PiDiGEMCdJYlaQ1NcPpQJOFtdwtrzO1uEIIUS3IUnMCrxcHIkN8wbgV2mNCSGE2UgSs5Ix/fRdittPSBITQghzkSRmJYZS+19PFNPcLLU0QghhDpLErCS+jzfuTipKaxrILKy2dThCCNEtSBKzEkeVkv69PAH9wsBCCCEunSQxK4r0dwcgq7jGxpEIIUT3IEnMiqJaktjJIkliQghhDpLErCjS3wOArGLpThRCCHOQJGZF0p0ohBDmJUnMiiL83QAoq22kvLbBxtEIIYT9kyRmRW5ODvTycgGkNSaEEOYgSczKpEtRCCHMR5KYlUUGSBITQghzkSRmZcYye0liQghxyTqUxFasWMHu3bupra21dDzdnrE7UeaKCSHEJXPoyElz5sxBqdTnu759+xIbG9vqFRYWZtEgu5Nzx8R0Oh0KhcLGEQkhhP3qUBIDyMzMpKioiAMHDpCWlsaGDRt47bXXqKqqQq1WExMT0yqxDRs2TD6g2xHm64ZKqaCuUUtBpYZeahdbhySEEHarQ0ns1VdfJSgoiMjISEaOHNnqZ1lZWaSnpxtf3377LadOnSIzM5OoqCiLBG3PHFVK+vi6kVVcw8niakliQghxCTqUxJ544gl+/vlnxo8f3+ZnkZGRREZGctNNNxmPeXp6mi3A7ijS352s4hqyimsY1dff1uEIIYTd6nB14tVXX820adPIzs6+6Lk33XQTrq6ulxRYdybFHUIIYR4dTmK//PILp0+fZtCgQbzwwgvU19ef99yPP/6Y4OBgswTYHcmEZyGEMI8OJ7FRo0axZ88eli9fzsqVKxk4cCDr1q2zZGzdVpQkMSGEMAuTJzvfc889ZGRkcM8993DPPfcwYcIEDh8+bInYuq2IliSWXVpLo7bZxtEIIYT96tSKHW5ubixevJhjx44RFBREfHw8jz76KOXl5WYOr3vq5eWCi6OSpmYdZ8rqbB2OEELYrUtadiosLIxPP/2ULVu28O233zJgwABzxdWtKZUKIvwMXYqyQaYQQnRWhyc7G+Tn53Po0CEOHTrEwYMHOXToEEeOHKGmpgYXF5nz1FGR/u4cy6/idIks5SWEEJ3V4SQ2btw4Dh8+TGlpKaCfHzZ06FCuvfZannzySYYOHUp0dLTFAu1uQrz1UxDyKs5f5SmEEOLCOpzEYmNjmTlzJkOHDmXIkCG4ublZMq5uL7hlpY7cchkTE0KIzupQEps9ezZvvPGGJC4zkpaYEEJcug4Vdnz44YfU1MicJnMytMTypCUmhBCd1qEkptPpmDhxIvfffz9vvPEGW7ZsuWA5/apVq6irkw/nCwlW61tiBVUatM06G0cjhBD2qUPdiZs2beLAgQOkp6fz0UcfceTIERobGwkNDSUuLo7Y2Fji4uKIi4sjKiqKRx99lLFjx8oq9hcQ4OmMg1JBU7OOwqp6Y1ITQgjRcR1KYhMmTGDChAnG77VaLceOHTNuv7J3715WrlxJQUEBHh4esgN0B6iUCoK8XDhbXkduuSQxIYToDJPniQGoVCoGDx7M4MGDmTFjhvF4cXExaWlpHDx4EB8fH7MF2V0Fq/VJLK+iDpA/LyGEMFWnktj5+Pv7c/XVV3P11Veb87LdVrC3K5wuI69cKhSFEKIzLmnZKXFpQgxzxSqkCEYIITpDkpgN/V5mLy0xIYToDEliNhRsnPAsLTEhhOgMSWI2FNJSkZgrq3YIIUSnmFzYUVpayttvv01RURH9+vUzzg/z8vKyRHzdWrC3vjuxuFpDQ1MzTg7yO4UQQpjC5E/N2267jRUrVpCfn8/q1au55ppr8PHxISoqiltuucUSMXZbfu5OODko0emgoFJaY0IIYSqTW2K7d+9m+/btDBs2DICmpiaOHDlCamoqaWlp5o6vW1MoFASrXThdUktueR1hvrLAshBCmMLkJNa3b18cHH5/m4ODAzExMcTExDBr1iyzBtcTGJKYrGYvhBCmM7k78bnnnmPp0qWWiKVH+r24QyoUhRDCVCYnsQULFvDNN98wduxYPvroIzIyMiwRV49hKO6QuWJCCGE6k5PY008/zcyZMwGYP38+gwYNwtPTk9GjRzNv3jyzB9jdGRb+lbliQghhOpPHxB544IFW3586dcq4TUt6errZAuspQlpaYrnSEhNCCJN1egHguro6GhsbiYiIICIightvvNGccfUY0hITQojOM7k7MT8/n4kTJ+Lp6YmPjw/h4eEsXLiQqqoqS8TX7RkKO8pqG6lr0No4GiGEsC8mJ7GHH36YxsZGfvzxR7Zv387//d//8c033zBq1CjKy8stEGJr77zzDpGRkbi4uJCQkMD27dvPe+6WLVtQKBRtXseOHbN4nB3l5eqAm5MKkNaYEEKYyuQk9vPPP7N8+XLGjRvHqFGjmD17Nunp6URHR7NgwQJLxGj02WefMX/+fJ555hlSU1MZM2YMkydPJjs7+4Lvy8jIIC8vz/iKjo62aJymMEx4BmSumBBCmMjkJObi4kJjY2PriyiVLFq0iC+//NJsgbXn9ddf54EHHuBPf/oTgwYNYtmyZYSFhbF8+fILvi8wMJBevXoZXyqVyqJxmiqkZTX73HJpiQkhhClMTmJTp07l+eefp7m5udVxR0dH6ust15JoaGggJSWFSZMmtTo+adIkduzYccH3xsfHExwczIQJE9i8ebPFYuwsaYkJIUTnmFyd+MorrzBmzBgSEhJ4+OGHGTx4MHV1dfz1r39l3LhxFghRr7i4GK1WS1BQUKvjQUFB5Ofnt/ue4OBgVqxYQUJCAhqNho8//pgJEyawZcsWrrzyynbfo9Fo0Gg0xu8rKyvN9xDnIRWKQgjROSYnMR8fH3bs2MHixYtZtGgR+fn5KJVKxo8ff9FuPXNQKBStvtfpdG2OGQwYMIABAwYYv09KSiInJ4elS5eeN4ktWbKEF154wXwBd4ChJVZYqbnImUIIIc7Voe7EmpqaVt97eXmxdOlScnNzKS4upqamho0bNxIcHGyRIAH8/f1RqVRtWl2FhYVtWmcXcvnll5OZmXneny9YsICKigrjKycnp9Mxd1SQlz6JFVRJd6IQQpiiQ0nMz8+P0tJSAF588UW+/vprY0Wgr68vzs7OlouwhZOTEwkJCWzatKnV8U2bNjFq1KgOXyc1NfWCydbZ2RkvL69WL0sL9NL/+RVIS0wIIUzSoe7ELVu24O7uDsBPP/3EsmXLKCsrw9vbm9jYWGJjY4mLiyM2Npb4+HiLBZucnMzMmTMZPnw4SUlJrFixguzsbObMmQPoW1Fnz55l9erVACxbtoyIiAgGDx5MQ0MDn3zyCevWrWPdunUWi7EzDC2x4moNTdpmHFSyw7MQQnREh5LY5Zdfbvx6y5YtAGRnZ5OWlkZ6ejppaWl8++23nDp1iqamJosECjB9+nRKSkpYvHgxeXl5DBkyhPXr1xMeHg5AXl5eqzljDQ0NPPHEE5w9exZXV1cGDx7Md999x3XXXWexGDvD180JB6WCpmYdxdUN9GoZIxNCCHFhCp1OpzPlDXfeeScrV640tszOVV1djYeHh9mC6woqKytRq9VUVFRYtGtx1JKfyK2o56tHriA2zNti9xFCiK7OlM9dk/utPv/88zaFHqAv/vjoo49MvZxoEWgo7qiU4g4hhOioDiex5ORk1q5dC+jL2v+ourpa9hO7BEGG4o4qKe4QQoiO6vA8saysLD7//HMAIiMjGTJkiLGoIyYmhv379+Pn52exQLs7Q3FHobTEhBCiwzqcxAzrIiqVSj755BNycnJITU1l5cqVHD16FJ1OxyuvvGKxQLu7IOlOFEIIk5m8YkdWVha9e/fGweH3t2q1WnQ6XatjwjSBnjJXTAghTGVy1jGUs5+rq60Kb4+kJSaEEKYzOYkdOnSI2bNnc+zYMUJDQ4mLiyMuLo74+Hji4uLw9fW1RJzdnnFMTAo7hBCiwzq1s3NQUBCfffYZv/32GzU1Nbz44otMnDiRgIAAS8TYIxiqE0trGtA0aW0cjRBC2AeTW2L79+9n9+7dDB48GJVKxdKlS3FxceG2227j0UcftUSMPYLa1REnByUNTc0UVWkI9XGzdUid9tPRAgoqNYyJ9ifM136fQwjR9ZmcxNzc3HB11e9/5eXlRWVlJVFRUSxatIjXX3+d6dOnmz3InkChUBDk5UxOaR0FlfabxL4/mMef1+w3fh/l784NMcH8vwnROMqakEIIMzP5UyUuLo5Dhw4BEBERweHDhwH93l2//vqreaPrYYI87Xuu2NG8SpI/Twcg0t8dlVLByeIa3vz5BCu2nbRxdEKI7sjkltjChQuNe3rdddddLF68GLVazTfffENgYKDZA+xJ7LlCsbSmgdmr91HXqGV0P3/+dd8Iahu1fLo7m1e+P8abP2UyJSaEPn722cIUQnRNJiexcePGGb9+6KGHSElJ4ZZbbsHb25tVq1aZM7YeJ9BOl55q0jbz8JoUzpTVEe7nxlsz4nFQKfFSKXnoyii2Zxbx64kSnv3qEP+6b8R5d+IWQghTdag7MTU1lcbGxjbHVSoVH3zwAfX19RQWFnL99debPcCexF5bYj8cLmDXyVLcnVS8f89wvN2cjD9TKBS8OHUITiolW48Xsf5g/gWuJIQQpulQEhs+fDgeHh7ExsYya9Ys/vGPf7B582bKysr0F1HKgL05GMrsC+1s1Y4vU88AMGtUBP2DPNv8PCrAgz+P6wvAC98cpqq+7S9EQgjRGR3KPhUVFfz444/Mnj0bR0dH1qxZw3XXXYe/vz/h4eFMnTqVRYsWGddXFJ1jKOywp5ZYSbWGLRlFANwc3/u85/15XF8i/d0prNKw6tdTVopOCNHddWhMzMPDgzFjxjBmzBjjMa1Wy5EjR4w7O//6668sX76cm2++2WLBdnf2uKfYtwfyaGrWMaS3F9HttMIMXBxVzJvQj8c+S+c/KTnMvaofSqWMjQkhLk2HktiXX35JXFwckZGRxmMqlYqhQ4cydOhQ7r77bosF2JMYuhMr65uoa9Di6tT116T8MvUsADfHh1703GsG98LD+TA5pXXsPVVKYpRs3SOEuDQd6k588MEH6devH97e3owdO5ZHH32UVatWnbfgQ3SOh7MDbi2Jq7Cq67fGThZVk5ZTjkqp4MbYkIue7+bkwPVDgwH4b8oZS4cnhOgBOpTEioqKOH36NJ988gkTJ07k7NmzvPzyy+0WfIjO06/aYehS7PrFHf9raYWN7udPQMtWMhdz23B9i+27g3nUaJosFpsQomfo8Dyx0NBQQkNDueGGG4zHqqurSU1NNY6L/fvf/+axxx6zSKA9RaCnM1nFNV1+XEyn0/Flmj6J3TLs/AUdfzQ83IdwPzdOl9Sy4VA+tyZcvBtSCCHOp0MtsZUrV6LVtl1Z3VDwMXfuXFauXMmePXvMHmBPYy9zxVJOl5FTWoe7k4pJl/Xq8PsUCgW3DdMnLulSFEJcqg4lsaVLl1JaWgpAVFQUt9xyC4sXL+arr77i9OnTFg2wpzHOFeviq3ZsOloAwKTBvUwuQLklIRSFAnaeLCGntNYS4QkheogOJbFjx44Z9wp77rnnCA8PZ8uWLdx///1ERkbi4+PDuHHjZCsWM7CXltieLP0vNaP7+Zv83t7erozqq69M/GL/WbPGJYToWUxeO/Hee+9t9X1OTg5paWmkp6eTnp5urrh6LHuYK1bXoOXgmQoARkZ2bifvW+JD+fVECd8fyuPRq6PNGZ4QogcxOYmtW7cOHx8fxo8fD0BYWBhhYWFMmTLF7MH1REEtVX5duToxNbuMpmYdwWoXQn1cO3WN8QMDUSrgWH4VZ8pq7Xb/NCGEbZm86OH//d//UVxc3Ob4rl27yMrKMktQPZmhO7Er7ym255S+K3FEhG+nV6T3cXciIdwHgM3HCs0WmxCiZzE5iWVlZTFixIg2x3/77Tfuv/9+swTVkxm2Y6lp0FLdRedRGcbDOtuVaDB+YBAAP0kSE0J0kslJLCAggLNn2w7GJyYmkpqaapagejI3Jwc8XfS9vF1xXKyhqZn92frdCy41iU0YpN9EdcdvJdQ2dM2ELYTo2kxOYnfeeSfPP/98m+WmampqZEsWMzFWKFZ0vSR2KLeC+sZmfNwc6RfgcUnXig70IMzXlYamZn7JbNtFLYQQF2Ny1lm0aBG5ubnExsby9ddfU1BQwPHjx1m4cCEjR460RIw9TpBxh+eul8T2tnQlDo/wveRV6BUKBRNauhR/li5FIUQnmJzE3N3d2blzJ+PGjWPGjBmEhIQwcOBAMjIyZO1EM/l9X7GuV6FoGA9LvMSuRIPxA/Vdij8fK6S5WWeWawoheo5O9f+p1WreeecdSkpKSE1N5fDhw2RkZDBo0CBzx9cjddW5Ys3NOvaeU5loDolRvrg7qSis0nA4t9Is1xRC9BwmJ7GFCxfy3//+FwBnZ2fjHmMyHmY+xqWnulhLLKOgisr6JtycVAwO8TLLNZ0dVIyJ1q8G82PLUlZCCNFRJmeeTz75BH9//VJDtbW1xMTEMGTIEGJiYsjNzTV7gD1RV116ytCVmBDug4PKfL+0jB/0e5eiEEKYwuRPosLCQqKiogD44osvAP3csfj4eBYuXGje6HqorlrYYe6uRINxA/QtsYNnKyip7lqtTyFE12ZyEgsODiY/Px/QJ7H77ruPPn368NRTT7Fp0yazB9gTBZ5T2KHTdZ1iB8OYVVyYt1mvG+jpwqBgfffkLyek1F4I0XEmJ7HbbruNv/zlL7z33nt899133HTTTYB+fKysrMzc8fVIhlU7GpqaqahrvMjZ1lGtaSKruAaAy8w0HnauK6P1XdTb7Xy+WEVtI//Zl8PfNhzjq7SznCqu6VK/iAjR3Zi8APDixYuZN28ef/3rX3nooYeIiYkBID09neDgYLMH2BM5O6jwdXeitKaBgkoN3m5Otg6JY3n6VliQlzP+Hs5mv/6Y6ADe23aS7ZlF6HS6Tq/JaCtbMgpZvfM02zOLaNS2Tlo+bo48clU/7rsiEtUlzq0TQrRmchJzdXXl/fffb3M8JSWFa6+91ixBCQj0dG5JYvUM6OVp63CMXYmDQ9QWuf7wCB+cHZQUVGrILKymf5Dtn7mjVm4/yUvfHTV+PyDIk/g+3mQUVHE4t5Ky2kZe+u4oGw8XsHRaLH38ZMV+IczF5CQG0NzcTGZmJhqNhgEDBuDs7MySJUvMHVuPFuTlwrH8KvK7SIXikZYkdlmw+bsSAVwcVSRG+bHteBHbjhfZRRLT6XS88v0x3tt2EoA7RoTxwOhIos+JvaGpmf+k5PDyd0fZc6qUa9/Yxiu3xnBjbIitwhaiWzF5TCwjI4PBgwczaNAg4uLi8Pb2ZsaMGe0uCiw67/e5Yl0jiR3O02+Caa75Ye2xp3Gx5mYdT/zngDGBPT15IEtuGdoqgQE4OSi5KzGcDY9eychIX2obtMxfm8q3B2Q6ihDmYHISe+SRR4iJieHkyZOcOXOGr7/+mrKyMhITEzlz5owlYuyRfp8rZvuS80ZtM8fzqwHLdScCxknPu7NKqG/UWuw+5vDBL1ms238GlVLB0mmxzBnb94LjeH383Fg7+3LuGBFGsw7mr03jxyMyuVuIS2VyEtu9ezcvv/wyERERhISEMHHiRL7//nuuv/56nnzySUvE2CN1paWnThRW06BtxtPZodM7OXdE/yAPgrycqW9sZt+prlvpeiy/ktd+yADgxalDuC0htEPvUyoVvHzzUG6KC6GpWcfDa/az7XiRJUMVotszOYn5+PhQWlra5vj8+fPZuHGjWYISEORpmPBs+5aYoahjUIjXJa9cfyEKhcLYGtue2TU/3DVNWuavTaNB28zVgwK5c2SYSe83tNyuHdyLBm0zcz5JIbOgykLRCtH9dWo/sSeeeIKamppWx+vq6uyuLLorM3QndoUxsSPGykTLjYcZjGkZF9vWRcfFXt94nGP5Vfi5O7HklphO/Z13UCl58854kqL8qG3Q8tDHKVTWd435gELYG5OT2AsvvIBKpaJ///4sXLiQTz/9lA8++ID77ruPG264wRIx9kjGJFalsfkWJYdz9UUdlqpMPNfofvokdjSvksIutuxWyulSVmzXF3K8cmsMAZ6dny/n5KDkrRnxhKhdOFlcw+Ofp9v8/7MQ9sjkJObi4sKmTZtYsGABmzdvZs6cOTz33HMkJibyxhtvWCLGHsnfwwmlArTNOkpqGmwWh06n40ieZeeIncvPw5mhvfX32X6867TGdDodL313FJ0OpiWEMvGyoEu+pp+HM8vvTsBJpWTTkQKWb/3NDJEK0bN0ap6Yg4MDc+fOZe7cueaOR7RwUCnx93CmsEpDQWX9Jf3WfynOlNVRVd+Eo0pBv0APq9xzbP8ADp6tYOvxIm7tYNGEpW08UkBqdjmujiqevHaA2a4bG+bN4qmDefqLgyzdmEFcmDdXtLRG7UVueR1bMopIzS6jpKaBkpoGquoa8XF3IljtQm9vV4b0VjMm2r9LrD4jupdOJTFhHUFeLhRWaVq61SzfCmqPoSuxf5AnTg7W2TNu7IAA3tp8gu2ZRWibdTZfqqlJ22ysRnxgdKRxgWZzuWNkH1Kzy/lsXw7zP0tj/bwxNvulpaOq6hv5eNdpvk7L5Vj+eQpTiluPmysV+qQ98bIgpiWEdflnFPZBklgXFuTlzMGzkF9huwpFaxZ1GMSFeePp7EBZbSOHzlYQa+ZV8021bv8ZThRW4+PmyINjoyxyj+dvHExqThnHC6pJ/jyNj+4badFK0M6qqGvkX7+e4oNfTlJZ3wTok1N8Hx+u6OdPiNoFX3cnvFwdKaluIK+ijpzSWnadLCWjoIrU7HJSs8t5feNxrhnci7su70NSlJ8UhYlOs7sk9s477/Daa6+Rl5fH4MGDWbZsGWPGjDnv+Vu3biU5OZnDhw8TEhLCU089xZw5c6wYced1hblihy283FR7HFVKrujnz4bD+Ww9XmTTJFbfqOUfmzIBeOSqfni5OFrkPq5OKt6eMYwpb/3C9sxilm/9jUeu6meRe3WGTqfjf2lneeGbI5TX6isp+wa489CVfZl4WRA+7hfvJsyrqGPzsSL+k5JDanY53x3M47uDecSGqvnzuH5MuiyoSyZu0bXZVRL77LPPmD9/Pu+88w5XXHEF7733HpMnT+bIkSP06dOnzflZWVlcd911zJ49m08++YRff/2Vhx9+mICAAG699VYbPIFpgjwNFYq2S2KGoo7LrFDUca6xAwKMSWzehGir3vtcH+88TX5lPSFqF+6+PNyi94oO8mTxjUN4at0BXt90nBERvoyMNO8GpJ2RX1HPM18e5KeWnbejAz2YNyGa64YGm9TVG6x2ZUZiH2Yk9uFIbiVrdp9m3f4zpJ+pYM4nKUQHejB3fD9uiAmxeRfy+eh0OgqrNBwvqOJ4QTWnimsorKqnuLqBkmoNDU3NNDXr0DbrcHZQ4uHigIezA77uTgSrXQn21o8R9g3wINLfHXdnu/oI7pIUug5sdhQSEkJcXBxxcXHExsYSFxdH//79rd4FkJiYyLBhw1i+fLnx2KBBg7jpppvaXYD4L3/5C19//TVHj/6+wvicOXNIT09n586dHbpnZWUlarWaiooKvLys1xoBWLsnm6e/OMj4gYF8eO8Iq94b9F1HsS/oJ7CnL5qE2tUyrZD2nC2v44pXfkapgNRnJ6F2s969DeobtYx5dTNFVRpevTWG20eYNrG5M3Q6HY99lsb/0nIJ9HTm23mjzT4GZ4oNh/J58r/pVNU34aRSMm9CPx4a2xdHlXnGR0uqNXz4axard5ymSqPvnozyd2fu+H7cGBuCg5nu01nNzfrq3F0nS0g5XcbeU2UUm3H38RC1C/17eTKwlxeDgvX/jQpwN9ufrznVNWgpqdFQVtNIlaaR6vomqjVN1Dc2o2nS0tDUzB9niUwfEYZvB1rpf2TK526Hfg14/fXXOXDgAOnp6axevZq8vDxcXV0ZMmRIq+QWExODu7u7yQF3RENDAykpKTz99NOtjk+aNIkdO3a0+56dO3cyadKkVseuueYaPvjgAxobG3F0bPvBqNFo0Gh+/0taWVlphug7J8jG3YknCvUD9r28XKyawAB6e7sSHehBZmE1v5wo5voY6+9V92XqWYqqNASrXbgpvrdV7qlQKPjrLUM5klfJ8YJq5n6aypo/JVr9Q81QzGJY4Dg2VM1r02LNvruAn4czT14zkIfG9mX1jlO8vz2Lk8U1JH+ezrIfM5k9JpJpw8NwcVSZ9b4XUlrTwJaMQrYeL+KXzOI2U1xUSgXhfm70D/Skb6A7vbxc8Pdwxs/DGVdHFUql/pyGpmaq6puoqm+kuGV8MK+8npyyWk4W1VBS00BuRT25FfVsyfh9hRonlZJ+gR4M7OVJvyAP+gd60i/Qg94+rhb5e9DQ1ExBZT0FlfpY8ivqyKvQf19YqaGgqp6iKg31jc0mX3viZYGdSmKm6FASu+OOO7jjjjsA2LZtGw899BA333wzjo6O7N+/n5UrVxo3MmxqarJIoMXFxWi1WoKCWs/PCQoKIj8/v9335Ofnt3t+U1MTxcXF7W7iuWTJEl544QXzBX4JDDs822oR4OMF+kV/o4OsU1r/R2P7B5BZWM3W44VWT2LaZh3vtczbemB0pNUqMwHcnBxYfncCU9/6lT1Zpbz2QwYLrxtktfsXVWl45NP97MnSLy83e0wkT1070KKJ1MvFkbnjo5k1KoLVO0+zcvtJsktrefarwyz7MZO7Evtw+4gwQn3MvxebTqfjeEE1Px8r5KejBezPLmvVonB30m8TNDzChxERvgztrTZLUi2raeBEUTXH8qs4llfJsfwqMvKrqNY0cSSv0tiVb6BSKgjxdiHMx40gLxcCPZ0J8HTG08UBd2f9y1GpRKEAhQKatDrqGrXUN2qp0Wgpr2ugvLaRspqGlqpnDUUtXaEd5aRS4uPuiJeLo7Gr1MVRhZODEmeVsk03sKeFxpDPZXKH7IMPPsi7777LuHHjjMcOHz7MtGnTmDVrljlja9cfuzAvtgtwe+e3d9xgwYIFJCcnG7+vrKwkLMzy3UjtMbTESmo0NGqbrf7b+PGWNf0G2Ghvr7EDAlj5SxZbj1t/t+cfDudzqqQWtasjd45sO95qaX0DPHjtthj+vGY/K7adJCZUzQ0xlt+D7OCZCh78eB95FfW4O6l4bVos1w213i8Qni76XbDvvyKSz/fl8P72k5wpq+PNn0/wz80nGN3Pn9sSQhnbP+CS5pyV1jSwJ6uE7ZnFbMko4mx5XaufDwr24qoBAVzZP4BhfXws8kuMj7sTI9x9GRHx+7inTqfjTFkdR/IqySyoIrOwmuMF1WQVV1Pf2ExOaR05pXUXuGrnOKmUBKmd9eN2ahd6qV0I8tT/15As/TyccXdSdblKUpOTWE5ODr17t+5aGTx4MG+//TbPPvssf/nLX8wW3Ln8/f1RqVRtWl2FhYVtWlsGvXr1avd8BwcH/Pz82n2Ps7Mzzs5dY/6Kr5sTTiolDVp9c98Sv4VeiCGJ2WqDyhERvrg46nd7PpZfxSArVUjqdDrebWmFzUoKt9ng++ShwTx4ZRQrtp0k+bN0/D2cuTyq/b+35vB1ei5P/icdTVMzUQHurJg53GoT3P/I1UnFrFER3JXYh+8P5bN2bza/ntAnne2Zxcay/jHR/gzs5UV0kAfhvm5txtC0zTrKahs4WVTDsfxKjuZVkZpd1mZum7ODklF9/Rg/KIgJAwMJ8bbcbg0XolAoCPN1I8zXjWsG9zIeNxSUnC6p5UxZrb4lVamhuFpDtUY/NlWjaULbrEOng2adDgeVEldHJS6OKtycHPB2c8THzRFvNycCPJ2NySlY7YqPm2OXS04dZfK/zquvvpqXX36Zf/3rX62Oh4aGcujQIXPF1YaTkxMJCQls2rSJm2++2Xh806ZNTJ06td33JCUl8c0337Q6tnHjRoYPH97ueFhXo1QqCFI7k1Oq76O2fhKzbXeii6OKUX39jd081kpiO38r4cCZClwclcwaFWGVe57PX64dyOmSGn44XMDsj/bx+Zwks/85NGmbWbrxuDFxXzUggDfujLfYdAJTOKiUTIkNYUpsCNkltfwnJYeNhwvIKKgi5XQZKad/37JHpVTg5qjCxUmFs4OSak0TFXWNnK90rX+QB5dH+TFuQABJUf64Ollv3M1UCoWCIC8XgrxcukTFaldichJbvnw5iYmJjBs3juTkZOLj46mrq2PhwoX06tXr4he4BMnJycycOZPhw4eTlJTEihUryM7ONs77WrBgAWfPnmX16tWAvhLxrbfeIjk5mdmzZ7Nz504++OAD/v3vf1s0TnMKVruSU1pHbrn5uxAupLy2gaKWbWD+uFuxNU26LIifjxXyw+EC5o63Tqm9YQ3D24eH4edh21a5SqngjTviueeDPew5Vcq9q/aw7s+jzPYLTXG1hnn/TmXHbyUAzBnblyevGdAlS9z7+Lnx+KQBPD5pAGfL69iSUUjKqTJOFFWTWVBNXaOWKk2TscrxXCFqFwYFezEw2JPBIWpGRvrib+P/t8I8TE5iISEh7N+/n6effpo777yT+vp6dDodAQEBrFmzxhIxGk2fPp2SkhIWL15MXl4eQ4YMYf369YSH6+fv5OXlkZ2dbTw/MjKS9evX89hjj/H2228TEhLCm2++aRdzxAxC1PpxsbwK61YoGlphvb1d8bDhXJarLwtC+eVBDp6t4ExZrcVbo0fzKo3dVbPHWGZ1DlO5OKp4/57h3P7eTjIKqpj+3i4+uHc4A3tdWoss5XQZj6zZT35lPW5OKv52awxTYi0/7mYOvb1duSsxnLsS9f/2m5t1FFVrqG3QFzLUNWrxcHbAx80JHzdHm5fqC8vp0Dyx89FoNGRkZKBUKomOju4yY0nmZMt5YgB/23CM5Vt+Y1ZSOC9MHWK1+3686zTP/u8QVw0IYNV9I6123/bc/t5O9mSV8twNl3H/6EiL3uvxz9NZt/8M18cE8/aMYRa9l6nyK+q58/1dZBXX4O6k4p8z4hk/0PTV9OsatPzjx+Os3H6SZp1+5Y13706waYtbiHOZ8rnbqV9Pfv31V1566SWWLVtGVFQUQ4YM6ZYJrCswtMRyrdwSy7RxUce5rm0Z4N5wuP2pFOZSUFnP1+lnga7TCjtXL7ULXz48iqQoP2oatDzw0T7e3nwCTZO2w9fY8Vsxk9/Yxopt+gR2U1wIX80dLQlM2C2Tk9h7773HpEmT2LNnD88++yyFhfqlaD766CP27dtn9gB7umC1vkoqr8K6Y2KGysSu8OE2abC+tbH3VKlZV0v4o3/tOEWjVsfICF/ibLzo8Pl4uzmx+oGR3DkyDJ0OXvshg/FLt/LZ3mwate1PRm3SNvPdgTxuXb6DGe/v5lRJLb28XFh5z3CW3RFv0+5iIS6VyX97X331Vf71r38xbdo0PD1//4DTarW88MILbaoBxaUJ9m4ZEyu3dktMPyZmqzli5wr1cWNobzUHz1bw45EC7rDAvK0aTRNrdp0G4E9jLNtleakcVUr+evNQ4sK8eX3Tcc6W1/GXdQd5fdNxhvb2ZkAvD0K89QVBvxVVc+BMuXHCvKNKwfQRYTx17cAuUX0oxKUyOYnl5uYycmTbMZLhw4fzzDPPmCUo8buQlpZYSU0D9Y1aqyy/U1ytoaSmAYUCm80T+qNrBgdx8GwFGw7nWySJfb4vh8r6JiL93bl60KXv2mxpCoWC6SP6MDWuN5/sOs3yLb9RUKmhoLKAH48WtDnf192JuxP7cPfl4cbdEYToDkxOYkOGDCE1NdVYEWjg6elJRUWF2QITet5ujrg4KqlvbCa/op4If8usTXkuQ1dimI9bl5k7c+2QXizdeJwdJ0qorG80ayuiSdvMh79mAfolpuxpOxAXRxV/GhPFXYnhpGaX6VdXL6wmr7yOUB83+gV60DfAg+ERPlZdf1AIazE5iS1atIh58+YxePBg4PflmzZs2NDudiji0igUCkLUrpwsriG3os4qSczQldjfRpOc29Mv0JO+Ae78VlTD5mOFTI0z34K83x7II6e0Dl93J24dFmq261qTq5OKUf38GdXP39ahCGFVJiexG264gczMTIYNG0ZdXR1/+9vfqKmp4fPPP+dvf/ubJWLs8YK9XThZXGO1cTFbLzd1PpOHBPPW5hOs23/WbEmsuVnHW5tPAPpWWFdpeQohOqZTJfaPPfYYBw4c4Nlnn6WoqIjGxkbef/995s+fb+bwBFi/QrGrJrHbh+sXYt52vIjTJTVmueYPh/M5UViNl4sD9yRZdtNLIYT5dbq2NjIykkWLFpkzFnEe1pwrZtiWAmy3ZuL59PFz48r+AWw7XsSnu7NZcInbk+h0Ov75s74Vdu8VkVbZNkIIYV4mt8Ty8/OZNWsW0dHRREdHM23aNL7//ntLxCZaBLesqJ1nhfUTi6o0VNQ1olTotwPpau5O1I+7fr4vx6RJvu3ZnFHIkbxK3J1U3GfjhX6FEJ1jchKbOXMmv/76K3/+85+ZP38+Xl5e3Hbbbdx+++1otZf2oSLaF2zF9RMNrbBwP/cuWc02fmAgwWoXymob+f5g51fw0Ol0vPmTvhV2d1I4PhbefVYIYRkmJ7EdO3awZs0akpOTeeSRR/jggw/47bffOHr0aJfZEbm7MextZI2V7DMLW1bq6CLzw/7IQaXkjhH61tgnLZOTO+PnY4Wk5ZTj7KDkT6O73hJTQoiOMTmJBQUFtVqpA/SbT65YsYL33nvPbIGJ3xlaYpX1+o3vLOm4sby+axV1nOuOkWGolAr2nS7jWH7lxd/wB/WNWl745ggA914RQYCnrPsphL0yOYnNnTuXp556iubm1uu0+fr6Ul9v3aWRegpPF0c8W9a3s3SFYqZxzcSu2RIDCPJyYdJl+lU1Ptphemvsva0nyS7Vrx84z0p7lAkhLMPkJLZ06VJ++OEHYmJieP/999mzZw87d+7kscce48EHH7REjILf11DMteBcMZ1OR2ZhS2ViYNdtiQHGHZc/25vNobMdXykmp7SWd7box8KeuX4Q7rL4rRB2zeQk9o9//IMnnniCPn36sGjRIi6//HJGjx7Njh070Gg0rFy5kj179lBbW2uJeHssa8wVO7cyMSrA8iuDXIrLo/y4ISaYZh0s+OIgTedZwf2PXvjmCJqmZkb11b9fCGHfTP41dPr06UyfPt34fVFREWlpacbXsmXLyMjIQKfT0dRk2fGbniTECi2xrl6Z+EfPTbmMbceLOHi2gn/tOMWfLrIH2Nfpufx4tAAHpYIXbhxsXDJNCGG/LrkvJSAggIkTJzJx4kTjMY1Gw6FDhy710uIc1miJdfXKxD8K9HRhwXWDWPDFQf6+8TjXDulFqI9bu+duzijk8c/TAPjTmKgusU+aEOLSmdyduH//fqZMmcKIESO4/fbb+etf/8r69evJzc01nuPs7ExCQoJZA+3prDFXzB4qE/9o+vAwRkb4Uteo5fHP0ylpZ9PMXSdLmPNxCo1aHdfHBPPkNQNsEKkQwhJMTmIzZsygvLycadOmERgYyPfff8+MGTMIDQ0lKKjr78Nkr6wxV8weKhP/SKlU8NdbhuDsoGR3VilXv76Vr9LOotPpOFNWy9o92Tzwr71ompoZPzCQf9weh8qOtloRQlyYyd2JZ8+e5dtvv6Vfv36tjp88eZIDBw6YLTDR2rktMZ1OZ/bxHP2aiYbuRPtpiYF+m5b/zhnFk/9N51h+FY+uTePFb49SfE6rLCnKj3fuGoaTQ6fWvBZCdFEm/4ueMGECx48fb3M8KiqKm266yRwxiXYYxsRqG7RU1pm/YKaoSkNlfZNdVCa2Z2iomq/njubxif1xUikprtagUioY1seb+VdHs3LWcLsoVhFCmMbkltibb77JvffeS3h4uHFjTGF5rk4qfNwcKattJLeiDrWbeVdct7fKxPY4OSj5fxOiuXlYb06X1BITqpaV6YXo5kxOYseOHSM7O5vY2FhGjhxJYmIi8fHxxMXFcdlll+HgIJNHLSXE25Wy2kbOlNUxKNjLrNf+vSvRfsbDzifUx+28VYpCiO7F5Iwze/Zshg4dynPPPUd2djapqal89dVXnDp1CicnJ1l6yoIi/Nw5nFtptg0hz2VcqcOOijqEEMLkJFZWVsabb75JVFTriaUVFRWkpaWZKy7RjnA/fevidIn5V0PJ7KK7OQshxIWYXNgxdepUUlNT2xxXq9WMHTvWLEGJ9kX46QsuTpm5JWbPlYlCiJ7N5CQWExPDM888w9atWy0Rj7iACH99EjN3S8zeKxOFED2Xyd2J77zzDqdPn2b8+PGEh4dz+eWXEx8fz7Bhw4iPj8fX19cScQogoqU78UxZLQ1NzWab89QdKhOFED2TyZ+CWVlZlJWVsXnzZubPn4+rqytr167lhhtuICAgwBIxihYBns64Oqpo1ukTmbkYNpbsL0UdQgg706l6eLVazZVXXsmVV15pPKbVajl27JjZAhNtKRQKwv3cOJZfxemSWqICzJN0juTpk5i5y/aFEMLSOtUfdezYMb766itSUlLQaPRL+6hUKpn8bAWWKO44mqcv6pAkJoSwNya3xN59913mzp1LUFAQ+fn5qFQqBgwYYJzwnJycbIk4RYtwf/OW2Tc0NXOiZQuWyySJCSHsjMktsVdffZV//vOfnD17Fjc3Nz788EMGDhzIV199xaeffmqJGMU5IltaYlnF5mmJ/VZUTaNWh6ezA6E+rma5phBCWIvJLbG8vDyuv/56AJRKJYmJidx99908+OCDJCYmmj1A0Vq4n6HM3jxJ7GjLeNjAYE/Z6VgIYXdMbokFBwdTUVEBQFBQEMXFxQDMmzeP1157zbzRiTYi/A1l9nU0apsv+XpHpahDCGHHTE5i48aNY8uWLQDExsayYcMGABobG8nJyTFrcKKtIE8XnB2UNDXrzLJBphR1CCHsmcndiW+88Qb5+fkAJCcnM378eA4fPkxqaipXXHGF2QMUrSmV+jL74wXVnCqpNXYvdoZOp5OWmBDCrpncEvP09CQ6OhqApKQk1q9fj1qt5pZbbmHNmjVmD1C0ZUhcpy6xuKOoSkNJTQNKBQyQhX+FEHbokjf/uuqqq7jqqqvMEYvooEh/88wVM0xyjvB3x9VJlpsSQtgfk1piGzZsICEhAScnJ1xdXYmPj+fll182FnoI6zDXliwyHiaEsHcdTmI7d+5kypQpeHt7s3jxYl566SVGjhzJP//5T4YMGSJ7iVmRuVbtMIyHySRnIYS96nB34muvvcatt97K2rVrWx1/4403mDNnDjfccAOHDh3C29vb3DGKPzC0xHJKa9E261ApOze/6/eiDhkPE0LYpw63xHbt2sXDDz/c5riLiwurVq2iT58+LF++3KzBifYFq11xUilp1Ha+zL6+UcvJlsIQ6U4UQtirDiexoqIiIiIi2v2ZQqFg3rx5rF+/3lxxiQtQKRX0aWmNdbZLMbOgGm2zDh83R3p5uZgzPCGEsJoOJzGtVouLy/k/7BISEmQrFisybJDZ2TL7c+eHyXJTQgh7ZVJ14urVq9m9ezf19fVtfubl5SVVilbUN1C/l9jR/KpOvf9Qrv7/lXQlCiHsWYcLO0aPHs2LL75IVVUVDg4ODBw4kGHDhjFs2DASEhIICgpCq9VaMlZxjthQbwDSc8o79f69p8oASAj3MVNEQghhfR1OYtu2bQMgMzOTlJQU9u/fT0pKCl9//TXl5eXSJWVlMaFqADLyq6hv1OLi2PHJypX1jRzL13cnDpckJoSwYyav2BEdHU10dDR33HGH8VhWVhb79u0jNTXVrMGJ8+vt7Yq/hxPF1Q0cyatkWJ+OJ6P9p8vQ6fSl+oFS1CGEsGOXvOwUQGRkJJGRkUybNs0clxMdoFAoiAn15udjhaTnlJuUxPaeKgVgeLivpcITQgirMHkBYFspKytj5syZqNVq1Go1M2fOpLy8/ILvuffee1EoFK1el19+uXUCtgLDuNiBM6YV1BjGw0ZGSleiEMK+maUlZg0zZszgzJkzxv3LHnzwQWbOnMk333xzwfdde+21rFq1yvi9k5OTReO0ppgw/bhY+pnyDr9H06Q1FoMMj5CWmBDCvtlFEjt69CgbNmxg165dJCYmAvD++++TlJRERkYGAwYMOO97nZ2d6dWrl7VCtSpDS+xkUQ0VdY2oXR0v+p5DZyvQNDXj6+5ElH/n9yITQoiuwC66E3fu3IlarTYmMIDLL78ctVrNjh07LvjeLVu2EBgYSP/+/Zk9ezaFhYUXPF+j0VBZWdnq1VX5ujsR5usK6JNTRxi6EoeH+0hFqRDC7tlFEsvPzycwMLDN8cDAQOMu0+2ZPHkya9as4eeff+bvf/87e/fuZfz48Wg0mvO+Z8mSJcZxN7VaTVhYmFmewVJiDPPFOtiluK+lqGNkpHQlCiHsn02T2PPPP9+m8OKPr3379gG022rQ6XQXbE1Mnz6d66+/niFDhjBlyhS+//57jh8/znfffXfe9yxYsICKigrjKycn59If1IJiW+aLdWTSc3Ozjn2nW1piMh4mhOgGbDomNnfu3FbzzdoTERHBgQMHKCgoaPOzoqIigoKCOny/4OBgwsPDyczMPO85zs7OODs7d/iatmZKheKJomrKaxtxcVQyOESWmxJC2D+bJjF/f3/8/f0vel5SUhIVFRXs2bOHkSNHArB7924qKioYNWpUh+9XUlJCTk4OwcHBnY65qxnSW41SAXkV9RRW1l9w8rJhflh8mA+OKrvoSRZCiAuyi0+yQYMGce211zJ79mx27drFrl27mD17NjfccEOrysSBAwfy5ZdfAlBdXc0TTzzBzp07OXXqFFu2bGHKlCn4+/tz88032+pRzM7d2YHoQP2mlukXaY3tzdInsREyHiaE6CbsIokBrFmzhqFDhzJp0iQmTZpETEwMH3/8catzMjIyjCvpq1QqDh48yNSpU+nfvz+zZs2if//+7Ny5E0/P7rWTsWEdxQMXKO6o1jSx6Yi+S/aKvn7WCEsIISzOLuaJAfj6+vLJJ59c8BydTmf82tXVlR9++MHSYXUJcX28+U/KGX46WkjyxP7tFrt8nZZLTYOWKH93qUwUQnQbdtMSE+d33ZBgXByVHMmrZHdLl+G5dDoda3afBmBGYh+ZHyaE6DYkiXUDPu5O3DIsFIAPf8lq8/MDZyo4nFuJk4OSW1vOE0KI7kCSWDdx/xURAGw6WsDpkppWPzO0wq4fGoyPe/dZO1IIISSJdRP9Aj0Z2z8AnQ5W/XrKeLyirpGv03MBuCuxj42iE0IIy5Ak1o08MDoSgP/sy6GyvhGA/6Wepb6xmQFBniTILs5CiG7GbqoTxcWNifYnOtCDzMJqnl53AJ0OfsksBuCuy6WgQwjR/UhLrBtRKBTc39IaW38wn+8P5VOlaaKXlws3xfe2cXRCCGF+0hLrZm6O783ukyVUa7QkhPuQEO5DTKgaF0eVrUMTQgizkyTWzbg4qlh2R7ytwxBCCKuQ7kQhhBB2S5KYEEIIuyVJTAghhN2SJCaEEMJuSRITQghhtySJCSGEsFuSxIQQQtgtmSd2EYaNNisrK20ciRBC9AyGz9tzNzo+H0liF1FVVQVAWFiYjSMRQoiepaqqCrVafcFzFLqOpLoerLm5mdzcXDw9PTu1gG5lZSVhYWHk5OTg5eVlgQi7tp78/D352aFnP39Pfna49OfX6XRUVVUREhKCUnnhUS9piV2EUqkkNPTSd0P28vLqkX+ZDXry8/fkZ4ee/fw9+dnh0p7/Yi0wAynsEEIIYbckiQkhhLBbksQszNnZmUWLFuHs7GzrUGyiJz9/T3526NnP35OfHaz7/FLYIYQQwm5JS0wIIYTdkiQmhBDCbkkSE0IIYbckiQkhhLBbksQs6J133iEyMhIXFxcSEhLYvn27rUOymm3btjFlyhRCQkJQKBT873//s3VIVrNkyRJGjBiBp6cngYGB3HTTTWRkZNg6LKtYvnw5MTExxkmuSUlJfP/997YOy2aWLFmCQqFg/vz5tg7FKp5//nkUCkWrV69evSx6T0liFvLZZ58xf/58nnnmGVJTUxkzZgyTJ08mOzvb1qFZRU1NDbGxsbz11lu2DsXqtm7dyiOPPMKuXbvYtGkTTU1NTJo0iZqaGluHZnGhoaG88sor7Nu3j3379jF+/HimTp3K4cOHbR2a1e3du5cVK1YQExNj61CsavDgweTl5RlfBw8etOwNdcIiRo4cqZszZ06rYwMHDtQ9/fTTNorIdgDdl19+aeswbKawsFAH6LZu3WrrUGzCx8dHt3LlSluHYVVVVVW66Oho3aZNm3Rjx47VPfroo7YOySoWLVqki42Nteo9pSVmAQ0NDaSkpDBp0qRWxydNmsSOHTtsFJWwlYqKCgB8fX1tHIl1abVa1q5dS01NDUlJSbYOx6oeeeQRrr/+eq6++mpbh2J1mZmZhISEEBkZyR133MHJkyctej9ZANgCiouL0Wq1BAUFtToeFBREfn6+jaIStqDT6UhOTmb06NEMGTLE1uFYxcGDB0lKSqK+vh4PDw++/PJLLrvsMluHZTVr165l//797N2719ahWF1iYiKrV6+mf//+FBQU8NJLLzFq1CgOHz6Mn5+fRe4pScyC/rh1i06n69R2LsJ+zZ07lwMHDvDLL7/YOhSrGTBgAGlpaZSXl7Nu3TpmzZrF1q1be0Qiy8nJ4dFHH2Xjxo24uLjYOhyrmzx5svHroUOHkpSURN++ffnoo49ITk62yD0liVmAv78/KpWqTaursLCwTetMdF//7//9P77++mu2bdtmlu187IWTkxP9+vUDYPjw4ezdu5c33niD9957z8aRWV5KSgqFhYUkJCQYj2m1WrZt28Zbb72FRqNBpVLZMELrcnd3Z+jQoWRmZlrsHjImZgFOTk4kJCSwadOmVsc3bdrEqFGjbBSVsBadTsfcuXP54osv+Pnnn4mMjLR1SDal0+nQaDS2DsMqJkyYwMGDB0lLSzO+hg8fzl133UVaWlqPSmAAGo2Go0ePEhwcbLF7SEvMQpKTk5k5cybDhw8nKSmJFStWkJ2dzZw5c2wdmlVUV1dz4sQJ4/dZWVmkpaXh6+tLnz59bBiZ5T3yyCN8+umnfPXVV3h6ehpb5Gq1GldXVxtHZ1kLFy5k8uTJhIWFUVVVxdq1a9myZQsbNmywdWhW4enp2Wbs093dHT8/vx4xJvrEE08wZcoU+vTpQ2FhIS+99BKVlZXMmjXLcje1ai1kD/P222/rwsPDdU5OTrphw4b1qBLrzZs364A2r1mzZtk6NItr77kB3apVq2wdmsXdf//9xr/zAQEBugkTJug2btxo67BsqieV2E+fPl0XHBysc3R01IWEhOhuueUW3eHDhy16T9mKRQghhN2SMTEhhBB2S5KYEEIIuyVJTAghhN2SJCaEEMJuSRITQghhtySJCSGEsFuSxIQQQtgtSWJCmEl0dHSbXW3P9/r3v/8NQFJSkvHY0aNHbfwEbYWHh7eJ3dQNLm+77bY211i9erWFIhY9jUx2FsIMTF3cOTs7m4CAANRqNQ0NDfj5+VFUVNSldjkoLS1ts32Gi4sL1dXVJq0BOGDAAI4fP97q2IEDBxg6dKhZ4hQ9m7TEhDCDwMBAdDqd8bVx40bjzyIiIlr9TKfTERYWhouLCxqNBp1OR3FxcZdKYACpqanGr5988kl0Oh11dXUmL2KbkZGBTqfjqquuAvSJcNCgQWaNVfRcksSEsIC0tDTj1/Hx8e2e8/rrrxu71xYvXmw8fuONNxqPf/XVVzz66KMEBQXh5eXFPffcQ21tLenp6Vx//fV4enrSq1cvnn322TbXr6mp4dVXX2XkyJF4e3vj5eVFfHw8y5cvp6mp6aLPcG4SO3drEYOUlBTuvvtuoqKicHV1xdPTk8jISG6++WZ279593usNHToUBwdZe1yYiUVXZhSih5oxY4Zx4d/Fixe3e85dd91lPOfbb781Hu/du7fxeGhoaJuFhK+77jqdi4tLm+NffPGF8RoZGRm6yMjI8y5G/MADD5j0DJmZma1+9ssvv+gcHBzOe/0ff/yx1fknTpww/uzBBx805Y9SiAuSlpgQFnBuK+Z8LbH9+/cbvx42bBgARUVFnD171nh8zJgxpKSkMH36dOOx9evXM2/ePPbv38/EiRONx9PT0wGoqqri2muvJSsri/79+7N582bq6+vJzc2lf//+AKxatYrS0tIOPYNaraZv376tfvbRRx/R1NSEg4MDe/bsoaGhgby8PH7++WfmzZvX5plTUlIu+uchRGdIm14IM6urq2tVyNDeh3ZtbS0ZGRkABAcHGzcNPPfDfsCAAXz00Uc4OjoyduxYPvvsMwAmTpzI3/72NwBiYmKMm6/27t0bgH/+859kZWUBcPz4ceNY1Lmam5spKSnB19e33Weora01PsOwYcPajNd5e3sD0NTUxOOPP87o0aMZNWoUEyZMaPd+5z6XIWELYQ7SEhPCzA4cOIBWqwXA39/fmFzOlZaWRnNzM9B6vOnc1tmdd96Jo6Mj8HsrC+Duu+82ft1ei++bb765aIxOTk6Eh4d36BnaGw97/PHHGTlyJADbt29nyZIlTJkyhdDQUH744Yc25xuSmEqlkqpEYVaSxIQws44UdZzbMjk3SZx7fNSoUcav9+7d2+5xQxJzcnIiJiYGgIKCAuPP8/Pz21RG6nQ66uvrcXJyOu8zXKyoIygoiN27d5OamsqyZcu46aabAH1ZfnJycpvzDcl54MCB3X53a2FdksSEMLOOJLH2xsP+eNyQPDQaDYcOHQLAx8eHfv36AfDbb79RVlYG6Cv+DEkpKirKeI3nn3+ewsJCGhoaOHHiBMuXLycxMZHMzMwLPsOFktgPP/zAU089xZ49e+jduzePPPIIzz33HC4uLgBtuh5PnjxpjFO6EoW5yZiYEGZmalGHIUmUlZVx6tQpQJ+IDONV6enpNDQ0ADB8+HDj+/bt29fmGgALFy5k69atNDU18e677/Luu++2urenpyfR0dEdegYvLy9j0jT46aefeO2113jttdfavE+hUPDkk0+2OiZFHcKSpCUmhBk1Nzdz8OBB4/dxcXFtztFoNBw5cgTQT5I2jJmdm9jOl6zOPX5ucjj3+Pjx49mzZw/Tpk0jIiICJycn/P39GTp0KH/+85/573//e8GJ1U1NTcaWX3tFHYmJidx2223069cPT09PVCoVvXr14sYbb2TTpk3MmjWr1fmSxIQlybJTQohWDh48aBxfS05O5u9///slXW/ixIn8+OOPAJSXl6NWqy85RiEMpCUmhGjl3O5Qw6oi51ZEdsS+ffuMq44YElhUVJQkMGF2ksSEEK2cm8QMDC2zjjq3S9VAuhKFJUh3ohBCCLslLTEhhBB2S5KYEEIIuyVJTAghhN2SJCaEEMJuSRITQghhtySJCSGEsFuSxIQQQtgtSWJCCCHsliQxIYQQduv/A1hAM24C59RPAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.font_manager import FontProperties\n",
    "\n",
    "x1 = np.linspace(0.0, 5.0, 100)\n",
    "y1 = np.cos(2 * np.pi * x1) * np.exp(-x1)\n",
    "\n",
    "font = FontProperties()\n",
    "font.set_family('serif')\n",
    "font.set_name('Times New Roman')\n",
    "font.set_style('italic')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 3))\n",
    "fig.subplots_adjust(bottom=0.15, left=0.2)\n",
    "ax.plot(x1, y1)\n",
    "ax.set_xlabel('Time [s]', fontsize='large', fontweight='bold', fontproperties=font)\n",
    "ax.set_ylabel('Damped oscillation [V]', fontproperties=font)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import ga\n",
    "\n",
    "'''Set number of items that needs to be generated'''\n",
    "num_items = 27\n",
    "num_seq = 2 ** num_items # values of the random variable X\n",
    "probs = np.random.random(size=num_seq)\n",
    "probs /= probs.sum()\n",
    "\n",
    "\n",
    "print(probs)\n",
    "print(sum(probs))\n",
    "print(max(probs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construct probabilities from a chosen distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from entropy_formula import entropy_formula\n",
    "from scipy import stats\n",
    "\n",
    "'''Set number of items that needs to be generated'''\n",
    "num_items = 8\n",
    "num_seq = 2 ** num_items # values of the random variable X\n",
    "\n",
    "\n",
    "\n",
    "# '''Generate the probabilities using the Normal/Gaussian distribution (randint func)'''\n",
    "\n",
    "probs = np.random.random(size=num_seq)\n",
    "\n",
    "'''Generate the probabilities using the discreet Laplace distribution'''\n",
    "\n",
    "# from scipy.stats import dlaplace # Discrete Laplace distribution\n",
    "# '''Define  the scale parameter for the Laplace distribution'''\n",
    "# scale_param = 1.0\n",
    "# probs = dlaplace.pmf(np.arange(1, num_seq+1), scale_param)\n",
    "# probs /= probs.sum()\n",
    "\n",
    "'''Print some metrics'''\n",
    "prob_max= sorted(probs, reverse=True)\n",
    "entropy_min = - (np.log2(prob_max[0]))\n",
    "entropy_shannon = stats.entropy(probs, base=2) # by python stats module\n",
    "\n",
    "\n",
    "from typing import Counter\n",
    "\n",
    "print(probs)\n",
    "# print(Counter(probs))\n",
    "print('Sum of probabilities: ',sum(probs))\n",
    "print('Shannon entropy: ', entropy_shannon)\n",
    "print('Min entropy: ', entropy_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bitstring generator Class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import zipf, planck\n",
    "\n",
    "class BitStringGenerator:\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def generate_uniform(self, size):\n",
    "        bit_string = np.random.randint(0, 2, size)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "    def generate_normal(self, size, mean=0, std_dev=1):\n",
    "        random_values = np.random.normal(mean, std_dev, size)\n",
    "        bit_string = (random_values >= 0.5).astype(int)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "    def generate_gamma(self, size, shape=2, scale=1):\n",
    "        random_values = np.random.gamma(shape, scale, size)\n",
    "        bit_string = (random_values >= 0.5).astype(int)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "\n",
    "    def generate_laplace(self, size, loc=0.5, scale=0.1):\n",
    "        random_values = np.random.laplace(loc, scale, size)\n",
    "        bit_string = (random_values >= 0.5).astype(int)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "    def generate_zipf(self, size, a=2):\n",
    "        zipfian_values = zipf.rvs(a, size=size)\n",
    "        max_value = max(zipfian_values)\n",
    "        bit_string = (zipfian_values == max_value).astype(int)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "    def generate_planck(self, size, lam=1.0):\n",
    "        planck_values = planck.rvs(lam, size=size)\n",
    "        bit_string = (planck_values >= max(planck_values) / 2).astype(int)\n",
    "        return ''.join(map(str, bit_string))\n",
    "\n",
    "# Example usage:\n",
    "size = 16\n",
    "\n",
    "bit_string_generator = BitStringGenerator()\n",
    "\n",
    "# Generate bit strings with different distributions\n",
    "uniform_bit_string = bit_string_generator.generate_uniform(size)\n",
    "normal_bit_string = bit_string_generator.generate_normal(size)\n",
    "gamma_bit_string = bit_string_generator.generate_gamma(size)\n",
    "laplace_bit_string = bit_string_generator.generate_laplace(size)\n",
    "zipf_bit_string = bit_string_generator.generate_zipf(size)\n",
    "planck_bit_string = bit_string_generator.generate_planck(size)\n",
    "\n",
    "print(\"Uniform Bit String:\", uniform_bit_string)\n",
    "print(\"Normal Bit String:\", normal_bit_string)\n",
    "print(\"Gamma Bit String:\", gamma_bit_string)\n",
    "print(\"Laplace Bit String:\", laplace_bit_string)\n",
    "print(\"Zipfian Bit String:\", zipf_bit_string)\n",
    "print(\"Planck Bit String:\", planck_bit_string)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example plot distribiution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm, gamma, laplace, zipf, planck\n",
    "\n",
    "# Generate random values for different distributions\n",
    "size = 1000  # Sample size\n",
    "uniform_values = np.random.rand(size)\n",
    "normal_values = np.random.normal(0.5, 0.1, size)\n",
    "gamma_values = np.random.gamma(2, 1, size)\n",
    "laplace_values = np.random.laplace(0.5, 0.1, size)\n",
    "zipf_values = zipf.rvs(2, size=size)\n",
    "planck_values = planck.rvs(1, size=size)\n",
    "\n",
    "# Create subplots for each distribution\n",
    "distributions = {\n",
    "    'Uniform': uniform_values,\n",
    "    'Normal': normal_values,\n",
    "    'Gamma': gamma_values,\n",
    "    'Laplace': laplace_values,\n",
    "    'Zipf': zipf_values,\n",
    "    'Planck': planck_values\n",
    "}\n",
    "\n",
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "for i, (name, values) in enumerate(distributions.items(), start=1):\n",
    "    plt.subplot(2, 3, i)\n",
    "    plt.hist(values, bins=30, density=True, alpha=0.6)\n",
    "    plt.title(name + ' Distribution')\n",
    "    x = np.linspace(min(values), max(values), 100)\n",
    "    \n",
    "    if name == 'Normal':\n",
    "        # Plot the probability density function for the normal distribution\n",
    "        mean, std_dev = np.mean(values), np.std(values)\n",
    "        pdf = norm.pdf(x, mean, std_dev)\n",
    "        plt.plot(x, pdf, 'r-', lw=2)\n",
    "    elif name == 'Gamma':\n",
    "        # Plot the probability density function for the gamma distribution\n",
    "        shape, scale = 2, 1\n",
    "        pdf = gamma.pdf(x, shape, scale)\n",
    "        plt.plot(x, pdf, 'r-', lw=2)\n",
    "    elif name == 'Laplace':\n",
    "        # Plot the probability density function for the Laplace distribution\n",
    "        loc, scale = 0.5, 0.1\n",
    "        pdf = laplace.pdf(x, loc, scale)\n",
    "        plt.plot(x, pdf, 'r-', lw=2)\n",
    "    elif name == 'Zipf':\n",
    "        # Plot the probability density function for the Zipf distribution\n",
    "        a = 2\n",
    "        pdf = zipf.pmf(x, a)\n",
    "        plt.plot(x, pdf, 'r-', lw=2)\n",
    "    elif name == 'Planck':\n",
    "        # Plot the probability density function for the Planck distribution\n",
    "        lam = 1\n",
    "        pdf = planck.pmf(x, lam)\n",
    "        plt.plot(x, pdf, 'r-', lw=2)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  }
 ],
 "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}