Kullback_Leibler.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing two 1D discrete distributions. Example\n",
    "\n",
    "We would like to know which distribution d1 or d2 is closer to distribution u.\n",
    "The correct answer here should be that distribution d1 is closer to u, because it covers the same area the u covers, whereas d2 only exists in the middle of the u (explanation in very simple English :))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFJpJREFUeJzt3X+wnmV95/H3p/lhRKxQyDpdQki0EUlhJ+BpcNeR7lTAuOyA7OgQf3Rwx07sjhG7jrODyw6y6TD+qNP1H1plShymraaIMHNgUynrj647XSRBEAg0a4gUTrGSBlbdlQQC3/3juXEfDyec5yRP7kO83q+ZM7l/XPf9vZ5zcj7nPtdz39dJVSFJasMvzXcHJEn9MfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDVk43x2Y7sQTT6wVK1bMdzck6ahy1113/WNVLZ2t3Usu9FesWMH27dvnuxuSdFRJ8nejtHN4R5IaYuhLUkMMfUlqyEtuTF+SjrRnnnmGqakp9u3bN99dmbMlS5awbNkyFi1adEjHG/qSmjM1NcUrX/lKVqxYQZL57s7Iqoq9e/cyNTXFypUrD+kcIw3vJFmXZGeSXUkun2H/7ya5L8k9Sf5HktVD+z7WHbczyVsPqZeSNEb79u3jhBNOOKoCHyAJJ5xwwmH9hjJr6CdZAFwDvA1YDbxrONQ7X6yqM6pqDfBp4A+7Y1cD64FfB9YBf9SdT5Lm1dEW+M873H6PcqW/FthVVbur6mlgC3DRcIOq+vHQ6iuA5/8G40XAlqraX1XfB3Z155MkzYNRxvRPAh4dWp8Czp7eKMkHgY8Ai4HfGjr2jmnHnnRIPZWkI2TF5f91rOd7+JMXzKn9VVddxbHHHsspp5zCVVddxYMPPsidd97JxMTEWPsFo4X+TL9LvOCvqVfVNcA1Sd4N/Cfg0lGPTbIB2ACwfPnyEbp0cAf74s31i6Cjy4t90z685N0z77jqR0eoN9PrvGr+6h+sdl/1e/CiX/uj7Pv+9NNP56abbuIDH/jAEasxyvDOFHDy0Poy4LEXab8FePtcjq2qa6tqoqomli6ddeoISTrqXX311Zx66qmce+657Ny5E4DTTjuNU0899YjWHSX0twGrkqxMspjBG7OTww2SrBpavQD4Xrc8CaxP8rIkK4FVwJ2H321JOnrdddddbNmyhbvvvpubbrqJbdu29VZ71uGdqjqQZCNwG7AA2FxVO5JsArZX1SSwMcm5wDPAkwyGduja3QA8ABwAPlhVzx6h1yJJR4VvfetbXHzxxRxzzDEAXHjhhb3VHunhrKraCmydtu3KoeUPv8ixVwNXH2oHJekX0XzdMurcO5LUs3POOYebb76Zp556ip/85CfccsstvdV2GgZJzev7Lp+zzjqLSy65hDVr1nDKKafw5je/GYCbb76ZD33oQ+zZs4cLLriANWvWcNttt421tqEvSfPgiiuu4IorrnjB9osvvviI1nV4R5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXEWzYl6cVmIz2k881tBtPnp1b+4Q9/yC233MLixYt57Wtfyxe+8AWOO+64sXbNK31Jeok477zzuP/++7n33nt53etexyc+8Ymx1zD0JWkezDS18vnnn8/ChYMBmDe+8Y1MTU2Nva7DO5LUs+GplQ8cOMBZZ53FG97whp9rs3nzZi655JKx1zb0Jalns02tfPXVV7Nw4ULe8573jL22oS9J8+BgUytff/313HrrrXzta187ItMvO6YvST072NTKX/3qV/nUpz7F5OTkz34LGDev9CWp5z8Sf7CplTdu3Mj+/fs577zzgMGbuZ/73OfGWtvQl6R5MNPUyh/96EePeF2HdySpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDvGVTUvPOuP6MsZ7vvkvvG+v5xskrfUlqyEihn2Rdkp1JdiW5fIb9H0nyQJJ7k3wtySlD+55Nck/3MTnOzkvS0ejhhx/m9NNP/9n6Zz7zGa666qpeas86vJNkAXANcB4wBWxLMllVDww1uxuYqKqfJvl3wKeB5+cEfaqq1oy535KkQzDKlf5aYFdV7a6qp4EtwEXDDarqG1X10271DmDZeLspSRqHUUL/JODRofWpbtvBvB/4y6H1JUm2J7kjydtnOiDJhq7N9j179ozQJUk6ei1cuJDnnnvuZ+v79u3rrfYooT/ThM41Y8PkvcAE8AdDm5dX1QTwbuCzSV77gpNVXVtVE1U1sXTp0hG6JElHr1e/+tU8/vjj7N27l/3793Prrbf2VnuUWzangJOH1pcBj01vlORc4ArgN6tq//Pbq+qx7t/dSb4JnAk8dBh9lqSx6vsWy0WLFnHllVdy9tlns3LlSl7/+tf3VnuU0N8GrEqyEvh7YD2Dq/afSXIm8HlgXVU9PrT9eOCnVbU/yYnAmxi8yStJTbvsssu47LLLeq87a+hX1YEkG4HbgAXA5qrakWQTsL2qJhkM5xwLfLn7816PVNWFwGnA55M8x2Ao6ZPT7vqRJPVopCdyq2orsHXatiuHls89yHF/A4z3UTdJ0iHziVxJTaqa8X6Ul7zD7behL6k5S5YsYe/evUdd8FcVe/fuZcmSJYd8Didck9ScZcuWMTU1xdH4XNCSJUtYtuzQn3819CU1Z9GiRaxcuXK+uzEvHN6RpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JashIoZ9kXZKdSXYluXyG/R9J8kCSe5N8LckpQ/suTfK97uPScXZekjQ3s4Z+kgXANcDbgNXAu5KsntbsbmCiqv4ZcCPw6e7YXwE+DpwNrAU+nuT48XVfkjQXo1zprwV2VdXuqnoa2AJcNNygqr5RVT/tVu8AlnXLbwVur6onqupJ4HZg3Xi6Lkmaq1FC/yTg0aH1qW7bwbwf+MtDPFaSdAQtHKFNZthWMzZM3gtMAL85l2OTbAA2ACxfvnyELkmSDsUoV/pTwMlD68uAx6Y3SnIucAVwYVXtn8uxVXVtVU1U1cTSpUtH7bskaY5GCf1twKokK5MsBtYDk8MNkpwJfJ5B4D8+tOs24Pwkx3dv4J7fbZMkzYNZh3eq6kCSjQzCegGwuap2JNkEbK+qSeAPgGOBLycBeKSqLqyqJ5L8PoMfHACbquqJI/JKJEmzGmVMn6raCmydtu3KoeVzX+TYzcDmQ+2gJGl8fCJXkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQ0YK/STrkuxMsivJ5TPsPyfJd5IcSPKOafueTXJP9zE5ro5LkuZu4WwNkiwArgHOA6aAbUkmq+qBoWaPAO8DPjrDKZ6qqjVj6Ksk6TDNGvrAWmBXVe0GSLIFuAj4WehX1cPdvueOQB8lSWMyyvDOScCjQ+tT3bZRLUmyPckdSd4+U4MkG7o22/fs2TOHU0uS5mKU0M8M22oONZZX1QTwbuCzSV77gpNVXVtVE1U1sXTp0jmcWpI0F6OE/hRw8tD6MuCxUQtU1WPdv7uBbwJnzqF/kqQxGiX0twGrkqxMshhYD4x0F06S45O8rFs+EXgTQ+8FSJL6NWvoV9UBYCNwG/AgcENV7UiyKcmFAEl+I8kU8E7g80l2dIefBmxP8l3gG8Anp931I0nq0Sh371BVW4Gt07ZdObS8jcGwz/Tj/gY44zD7KEkaE5/IlaSGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGjBT6SdYl2ZlkV5LLZ9h/TpLvJDmQ5B3T9l2a5Hvdx6Xj6rgkae5mDf0kC4BrgLcBq4F3JVk9rdkjwPuAL0479leAjwNnA2uBjyc5/vC7LUk6FKNc6a8FdlXV7qp6GtgCXDTcoKoerqp7geemHftW4PaqeqKqngRuB9aNod+SpEMwSuifBDw6tD7VbRvF4RwrSRqzhSO0yQzbasTzj3Rskg3ABoDly5ePeOrxOOP6Mw66775L75u3+n3Unu/6L9XP/XzX92vv5/5IGuVKfwo4eWh9GfDYiOcf6diquraqJqpqYunSpSOeWpI0V6OE/jZgVZKVSRYD64HJEc9/G3B+kuO7N3DP77ZJkubBrKFfVQeAjQzC+kHghqrakWRTkgsBkvxGkingncDnk+zojn0C+H0GPzi2AZu6bZKkeTDKmD5VtRXYOm3blUPL2xgM3cx07GZg82H0UZI0Jj6RK0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSEjhX6SdUl2JtmV5PIZ9r8syV90+7+dZEW3fUWSp5Lc0318brzdlyTNxcLZGiRZAFwDnAdMAduSTFbVA0PN3g88WVW/lmQ98Cngkm7fQ1W1Zsz9liQdglGu9NcCu6pqd1U9DWwBLprW5iLg+m75RuAtSTK+bkqSxmGU0D8JeHRofarbNmObqjoA/Ag4odu3MsndSf46yZtnKpBkQ5LtSbbv2bNnTi9AkjS6UUJ/piv2GrHND4DlVXUm8BHgi0l++QUNq66tqomqmli6dOkIXZIkHYpRQn8KOHlofRnw2MHaJFkIvAp4oqr2V9VegKq6C3gIeN3hdlqSdGhGCf1twKokK5MsBtYDk9PaTAKXdsvvAL5eVZVkafdGMEleA6wCdo+n65KkuZr17p2qOpBkI3AbsADYXFU7kmwCtlfVJHAd8KdJdgFPMPjBAHAOsCnJAeBZ4Her6okj8UIkSbObNfQBqmorsHXatiuHlvcB75zhuK8AXznMPkqSxsQnciWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpISOFfpJ1SXYm2ZXk8hn2vyzJX3T7v51kxdC+j3XbdyZ56/i6Lkmaq1lDP8kC4BrgbcBq4F1JVk9r9n7gyar6NeC/AJ/qjl0NrAd+HVgH/FF3PknSPBjlSn8tsKuqdlfV08AW4KJpbS4Cru+WbwTekiTd9i1Vtb+qvg/s6s4nSZoHo4T+ScCjQ+tT3bYZ21TVAeBHwAkjHitJ6kmq6sUbJO8E3lpVv9Ot/zawtqo+NNRmR9dmqlt/iMEV/Sbgf1bVn3XbrwO2VtVXptXYAGzoVk8Fdo7htb2YE4F/PMI1rP/Sq916/ZZfewv1T6mqpbM1WjjCiaaAk4fWlwGPHaTNVJKFwKuAJ0Y8lqq6Frh2hL6MRZLtVTXRVz3rvzRqt16/5ddu/f9vlOGdbcCqJCuTLGbwxuzktDaTwKXd8juAr9fgV4hJYH13d89KYBVw53i6Lkmaq1mv9KvqQJKNwG3AAmBzVe1IsgnYXlWTwHXAnybZxeAKf3137I4kNwAPAAeAD1bVs0fotUiSZjHK8A5VtRXYOm3blUPL+4B3HuTYq4GrD6OPR0JvQ0nWf0nVbr1+y6/d+p1Z38iVJP3icBoGSWpIc6E/25QSR7j25iSPJ7m/z7pd7ZOTfCPJg0l2JPlwz/WXJLkzyXe7+v+5z/pdHxYkuTvJrfNQ++Ek9yW5J8n2eah/XJIbk/xt93/gn/dY+9TudT//8eMkv9dj/X/f/Z+7P8mXkizpq3ZX/8Nd7R19vu6DqqpmPhi8Ef0Q8BpgMfBdYHWP9c8BzgLun4fX/qvAWd3yK4H/1fNrD3Bst7wI+Dbwxp4/Bx8BvgjcOg+f/4eBE/uuO1T/euB3uuXFwHHz1I8FwD8wuKe8j3onAd8HXt6t3wC8r8fXezpwP3AMg/dQ/xuwar7+H1RVc1f6o0wpccRU1X9ncHdT76rqB1X1nW75J8CD9Ph0dA38n251UffR2xtKSZYBFwB/0lfNl4okv8zgguM6gKp6uqr+9zx15y3AQ1X1dz3WXAi8vHuG6BhmeFboCDoNuKOqflqD2Qr+Gri4x/ov0FroOy0E0M2CeiaDq+0+6y5Icg/wOHB7VfVZ/7PAfwCe67HmsAL+Ksld3RPofXoNsAf4Qje89SdJXtFzH563HvhSX8Wq6u+BzwCPAD8AflRVf9VXfQZX+eckOSHJMcC/4ucfWO1da6GfGbY1dftSkmOBrwC/V1U/7rN2VT1bVWsYPJm9NsnpfdRN8q+Bx6vqrj7qHcSbquosBrPVfjDJOT3WXshgWPGPq+pM4P8Cvb6fBdA93Hkh8OUeax7P4Lf5lcA/BV6R5L191a+qBxnMOnw78FUGQ8oH+qo/k9ZCf6RpIX5RJVnEIPD/vKpumq9+dEML32Qw3XYf3gRcmORhBkN6v5Xkz3qqDUBVPdb9+zhwM/3ONjsFTA39ZnUjgx8CfXsb8J2q+mGPNc8Fvl9Ve6rqGeAm4F/0WJ+quq6qzqqqcxgM736vz/rTtRb6o0wp8Qupm+r6OuDBqvrDeai/NMlx3fLLGXwz/m0ftavqY1W1rKpWMPiaf72qervaS/KKJK98fhk4n8Gv/b2oqn8AHk1yarfpLQyeku/bu+hxaKfzCPDGJMd03wNvYfB+Vm+S/JPu3+XAv6H/z8HPGemJ3F8UdZApJfqqn+RLwL8ETkwyBXy8qq7rqfybgN8G7uvG1QH+Yw2etu7DrwLXd39E55eAG6qq91sn58mrgZsHmcNC4ItV9dWe+/Ah4M+7i53dwL/ts3g3nn0e8IE+61bVt5PcCHyHwbDK3fT/ZOxXkpwAPMNgKpone67/c3wiV5Ia0trwjiQ1zdCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakh/w/B+RU2oPCawwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110d603d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# small value\n",
    "sm = 1e-09\n",
    "d1 = np.repeat(sm, 10)\n",
    "d1[[0,4,8]]= 0.3\n",
    "d2 = np.repeat(sm, 10)\n",
    "d2[[4,5,6]] = 0.3\n",
    "\n",
    "u = np.repeat(0.1, 10)\n",
    "\n",
    "## plotting\n",
    "ind = np.arange(10)\n",
    "w = 0.2\n",
    "d1_plot = plt.bar(ind - w, d1, width = w, align='center')\n",
    "d2_plot = plt.bar(ind, d2, width= w, align='center')\n",
    "u_plot = plt.bar(ind + w, u, width= w, align='center')\n",
    "plt.legend( (d1_plot[0], d2_plot[0], u_plot[0]), ('d1', 'd2', 'u') )\n",
    "\n",
    "objects = list(map(str, range(10)));\n",
    "plt.xticks(ind, objects);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kullback-Leibler divergence\n",
    "\n",
    "Compares two distribution. Non commutative, thus not a metric.\n",
    "The smaller the value is the more similar two distribution are.\n",
    "\n",
    "$$ D_{KL}(p || q) = \\sum_{i=0}^N p(x_i) \\mathrm{log}\\frac{p(x_i)}{q(x_i)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def KLDivergence(p,q):\n",
    "    logTerm = np.log2(np.divide(p,q))\n",
    "    d = np.sum(np.multiply(p,logTerm))\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before the KL divergence could be applied to the discreet distributions, they should be convolved with gaussian kernel. They should be made more smooth.Formally the process is called **Kernel density estimation**. Otherwise in the given example the value of the divergence is the same for <d1,u> and <d2,u>, since they have the same number of bins of the same height. However, since the location of the bins are not the same, the smoothed version will populate the difference nearby bins. I am more than sure that there is a sane mathematical explanation behind this smoothness action, but I don't know it yet :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolving with gaussian\n",
    "you need to create a Gaussian that's discretized at the same spatial scale as your curve, then just convolve.\n",
    "Specifically, say your original curve has N points that are uniformly spaced along the x-axis\n",
    "Then the point spacing along the x-axis will be (physical range)/(digital range) = (10 -0)/N\n",
    "\n",
    "**Side note**: Convolving a probability distribution with a normalized kernel leads to a probability distribution. (Simple english) If your value sum up to 1 and convolve them with a kernel (sum of elements of which sum to 1), you will automatically get a probability distribution (all values sum to 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXuUXFd1p79d1d16v7sly3q1LCRjydgYhMHY4AcQ7IHYQCCxs0KAIfEMgwMJzEpMMsNkkZAHkxAgy8nCJCTOBGI8Dg4KMZiJ7YDfSDYWtmQkS63uVltv9UPvVlfVnj9O3e7q7nrcqrqvqtrfWlqlqrp17qnedfbZd5/f3UdUFcMwDKO5SMXdAcMwDCN4zLkbhmE0IebcDcMwmhBz7oZhGE2IOXfDMIwmxJy7YRhGE2LO3TAMowkx524YhtGEmHM3DMNoQtriOnFnZ6d2d3fHdXrDMIyG5Nlnnz2mql2VjovNuXd3d7Nt27a4Tm8YhtGQiEifn+MsLWMYhtGEmHM3DMNoQsy5G4ZhNCHm3A3DMJoQc+6GYRhNiC/nLiI3isguEdkjIncWef/DInJURJ7P//u14LtqGIZh+KWiFFJE0sBdwDuAAWCriGxR1Z1TDv2Wqt4RQh8NwzCMKvETuV8J7FHVHlU9D9wL3BJut5LJjgMjPNc/5J4cehF2/yDeDhmh8eIrI2zfP+yeHHge9vx7vB0yjCrx49xXAPsLng/kX5vKL4jIT0XkfhFZVawhEbldRLaJyLajR4/W0N14+bOHdvEH390JQ31wz7vhm79og75J+ZPv/YzPP/gSHN8L/3AzfOMD0PPDuLtlGL7x49ylyGtTd9X+V6BbVS8D/h24p1hDqnq3qm5W1c1dXRXvnk0co5kcjJ2F+z4IuRx0XQz3fxSGeuPumhEwo5ksMnYGvvVBkBQsXgf3fwRGBuLummH4wo9zHwAKI/GVwIHCA1T1uKqO5p9+DXh9MN1LFplsjo+d/is4uB3edzfc9k+Awrd+xTl9o2nIZHN87MRX4MhO+IW/dbbOnIf7fhUyo5UbMIyY8ePctwLrRWStiHQAtwJbCg8QkeUFT28GXgqui8nhutPf5+fOPwzX/g5cfCMsvgje9zU49AI8+N/j7p4RID935rtcd/4/4Ibfg1e9DTrXw3v/Gl55Fr7/mbi7ZxgVqejcVTUD3AE8hHPa96nqDhH5nIjcnD/sEyKyQ0S2A58APhxWh+Pk5898mx2pi+HaAjXohnfCGz8Gz38TzgzG1zkjUG45+y+8kN4I13x64sVLfh42/2d47h4YPRlf5wzDB7507qr6oKpuUNV1qvr5/GufVdUt+f9/RlU3qerlqnq9qv4szE7HwsgAq7L7eTR9FaSm/Nk2vRc0B/t+FE/fjGAZ7OHC3EEebbumuK1zGeh9PJ6+GYZP7A5Vv+x9FICn5bXT31vxepgxH/Y+EnGnjFDI2/rHqcunv7fqjdA+22xtJB5z7n7peZTjspiXdeX099JtsPatzinoVCGR0XD0PMoh6aJXl09/r20GdF8zPgEYRlIx5+6HXA72Psq29OVkS/nuddfDSD8M9kTaNSNgshno+RFbU+VsfQMcfxmG+yPtmmFUgzl3PxzaDmcH+XH6tWRyJUb8uhvco12uNzYHnoPREZ6WMra+6Hr3aNG7kWDMufsh77Cf4TKypcK5xRfBom5z7o3O3kcA4RkuJVvKuXddDPMuNFsbicacux/2PgrLXsMxXVA6mgMX0e17DLJj0fXNCJa9j8KFVzCkc8lkc8WPEXFXavt+CLlstP0zDJ+Yc6/E+dPQ/zSsu55MTktHc+AG/PmTMGAbfzck50ZgYCusu8GHra+Hs0Nw8Pno+mcYVWDOvRK9T0BuDNbdQDaXI5MrEc2BU8xIyi7XG5V9j4Fm87bWCldp17lHs7WRUMy5V2LvI9A2E1ZfRSan5BRypQb9rIVO824DvjHZ+wi0z4GVbyCTy5WP3Od0wvLLbVHVSCzm3CvR93j+xpWZ44M9W07LftF1TnFx/nQk3TMCpO8JWPNmaOsYj9y1kq33/9gKiRmJxJx7ObJjcORncKG7K9W7TC8b0S2/3JUiONJ8FRiamrFzcOzlabYuZ2qWX+5Sdsd2R9BBw6gOc+7lOL7HDd6lm4AJp142F7t0o3s8siPs3hlBcmyXy7cv3Ugup+M3GpddY8n/Ljg8dcdJw4gfc+7lOJx30Ms2ojqhniipdQdYtNbVHrEB31h49lq2adLkXfYqbck6SHfYRG4kEnPu5Ti8A1Jt0Llh0iAvG82lUrD0Ejj8YgQdNALj8IuQngGL102xdRnnnm53NzQdNuduJA9z7uU4shOWrIe2Gf6jOXCpmSM7rYhYI3Fkp3PU6bZJk3fZqzRwqRm7SjMSiDn3chzeCctcDt13NAewbBOcOQ6njoTZOyNIDu90dqNaW2+EkwdsoxYjcZhzL8W5EVflMT/gq4rc85+xXGyDcPo4nDpUp60tejeShTn3UhzJbwO7tIZoblxFYc69IfAm4aXFrtLKrK+AKWaMxGLOvRQFShmYPMizlQb8nCUwd5kN+EahQCkDVUbu8y6AWYtsAd1IHObcS3Fkp9s6b8EqoMrIHZyjsLRMY3D4RZidn5CZvIha0dYisOxSS8sYicOceykO73CX6SIAZAoHfCUFBeQVMz9zO/sYyebIzsm2nnSV5tfWL7kduwwjIZhzL4bqJKUMTB7kvgb8sk2QHbVt95JOLuccs7cwypSrND8T+bKNcP4UDPeF0UPDqAlz7sU48QqMjkwa8Jla0jJgqZmkM7QPxs6UtLW/ifxS92ipGSNBmHMvhrfAtrR4NOdrwHdeDJK2RdWkc6S8rSuqZQC6Xu0ezdZGgjDnXgxP+bD0kvGXCge5rwHfPtPVHjE5ZLI5vAMQWPrq8ZeqjtxnzHX755pixkgQ5tyLcWSnU8nMWjj+UtWRO5hiphE4vAMWr4WOOeMvZSdN5H5tbYoZI1mYcy/G4Z0TpXvzVJ1zB3epP9QLo6cC7JwRKEeK2Dpbw0S+dKMrET12LsjeGUbNmHOfSi4Hg3uhc/2klydF7n4UFDDRhilmkkl2zE2+nRsmv1zLRN65wW3SMtQbXP8Mow7MuU/l5AHInHP58gIy1dzY4uG1Mbg3qN4ZQTLcD7nMdFtPSsH51K4vucg9mq2NhGDOfSpelL34okkv15RzX7R2cptGshjc5x7L2NqXzr2wDbO1kRB8OXcRuVFEdonIHhG5s8xx7xcRFZHNwXUxYo7nI6/FU6O5KtUy4FQUcy+A4zbgE8lgKVvXMJHPWgSzFk/8fgwjZio6dxFJA3cBNwEbgdtEZGOR4+YBnwCeCbqTkTLY43bkmb9i0ss1Re7gIjqL5pLJYA90zIW5Sye9XJNaBszWRqLwE7lfCexR1R5VPQ/cC9xS5Lg/AL4ANLZcYLDHSeNSk/80NallwOViLQ+bTI7vdbbO15TxqClyB5e7N+duJAQ/zn0FsL/g+UD+tXFE5Apglap+N8C+xcNgz7QcLNQZuZ86bHLIJOLD1lVH7iMDJoc0EoEf5y5FXhv/xYtICvgL4NMVGxK5XUS2ici2o0eP+u9lVORyJQd8zZG7l8+1iC5ZZDOu0NeUfDtM1blXUelx8TpATQ5pJAI/zn0AWFXwfCVwoOD5POBS4D9EpBd4E7Cl2KKqqt6tqptVdXNXV1ftvQ6LkwedDLJoNFe4aXI1A95UFIlkJC+DDDpyB7O1kQj8OPetwHoRWSsiHcCtwBbvTVUdUdVOVe1W1W7gaeBmVd0WSo/DZFw9USRyr0XnDi6nW9i2kQw8B7ykSORecwrObG0kh4rOXVUzwB3AQ8BLwH2qukNEPiciN4fdwUgpM+BrzrnPmOd2+LFoLlkcL34/A0xRy/jVuQPMXuwkkWZrIwG0+TlIVR8EHpzy2mdLHHtd/d2KieN7Id0xTQYJdeTcwTkQ07oni8EeaJ8zvrVeITVH7pC3tUXuRvzYHaqFDPa4u0pT6Wlv1Ry5g1tos2guWQzudY5YpusFas65Q97W++rtnWHUjTn3QkooZaDeyH0tnDpkcsgkMdgzUQ9mCjXVlvFYfBGM7Dc5pBE75tw9cjkXcRXJt8MUtUy1A95rc8giukSQzTi5YomJvK7IfUleDmn7qRoxY87d4+RByJydUDxMoe6cO1guNimM7M/LIItP5DXVc/cwWxsJwdeCakswXg2yROSeH/Az2lL+67l7mP45WXhSxTJXaSLQnkrVPpGbrY2Yscjdo4zGHSai9Y62Ggb8jHkwZ6npn5NCiVK/Hpmc0pYS0impPnKfvRhmLjRbG7Fjzt3Dk0EuWFn07WxOSaeE9nSq+gEPLko0OWQyOL63pAwSJmzdlpLqdO4eS9ZZWsaIHXPuHoM9bgf7IjJIcNFcOh/NVR25g5WDTRKeKqqIDBK8yD1FOi3VL55D3ta2eG7Eizl3j8Gekvl2cHnYtnw0V/OANzlkMhjcW1IGCVMi95om8nUmhzRix5w7gKqLtEooZSAfuUudkTuYHDJuclkY6pvYArEImVyOdEpISQ05d8jb2uSQRryYcwc4dcTJIBd1lzwkm1PSaS9yr2HAL1rjHodswMfKiQOQG6ts63oid7O1kQDMucNEhFVmwBcqKGob8Gsnn8uIBz+2zuZtna51Iu+efC7DiAFz7jARYS1cU/KQXMGCatU6d3DVAjvmWTQXN97ff1FpW2fVi9xrVEbNXQZtM23TDiNWzLkDDPe6x4WrSx4yrqCo5cYWcMqMRWssmoub4T6QFCxYVfKQbD06d3C2XrjabG3Eijl3cBHWnKXQMbvkIYV52JrUMuCuDCyai5ehXlfSOd1e8pDMpJy72dpoTMy5Q1490V32kLpz7uDOMdzv1DlGPPiwdTbrXaXVGLmDO8dQf22fNYwAMOcO7vK5TA4WnM59InKvdcCvgbEzcDqBm4O3CsN9ZddWYGrkXoetR0fg7FBtnzeMOjHnnh2DkYHKAz5b5x2qMHEOu1yPh7Gzrvqnj4m8LV1Hzh3M1kbsmHMfGQDNVb5UzyltaaGtVnkcTJzDFDPxMLzfPfpIwXlqmZpqyxSew2xtxIQ59+HK0jjwBnwdahmYUON46hwjWoYrS14hALUMTPyeTDFjxIQ5dx8ad5gY8HWpZTpmO1WORXPx4KVIfE3k7iqtZrXMzAWu9K/Z2ogJc+5DvSBpJ48rg1dvJF1rGViPRSaRi42hXndzUYlSvx7ZXABqGTBbG7Fizn24DxaugnT5TakmR+71DPhuu1SPC08pU6LUr0cgahkwWxuxYs59qLI0DibXc6/LuS9cAyOvOJWOES1DlSWvMFHeORBbD/e7zdcNI2LMuQ/1+hzwE5F7fdHcGtCsU+kY0eJ3Is8WqGXqtXX2vJNfGkbEtLZzHz0FZ45VlMaBN+CdWqbutAzY5XrUnB1yNxX5sLUne60/554/l9naiIHWdu7D+dvDfURzkyP3Oi6zx29usQEfKT6qQXpk87LX+m3dPfnchhEhLe7cK9f29sjkcqTTddT49pi/wqlzLJqLFp8ad5hcR6im8s4eC1cBYrY2YqG1nbsnU6s6cq9jwKfbYMFKk8hFjU+NOxRUAE3Xaeu2GTBvudnaiIUWd+590D4H5nRWPHSSWqaeaA7yFQMtmouUoT63YcrMBRUPzQSllgGztREbre3cvWqQFXTPEGDkDrZpRxz4qAbpMVG7v061DJitjdjw5dxF5EYR2SUie0TkziLv/1cReUFEnheRx0VkY/BdDQGf0jiYXFum7mhu4RpX9vf86fraMfzjU+MOU3LuQdj6xAHIjNbXjmFUSUXnLiJp4C7gJmAjcFsR5/1NVX2Nqr4W+ALwxcB7GjSqeY17t6/DA1PLgFUMjJpczimj/No6G5BaBvLn1ImKlIYREX4i9yuBParao6rngXuBWwoPUNUTBU/nAMnfaujMcRg77T+ay07UlskpaD27KZn+OVpOHYLsaFVXaYHo3KGgOmRvfe0YRpWUL6jiWAEUhh0DwBunHiQiHwc+BXQANxRrSERuB24HWL269GbUkeCzGqRHYeQ+/jxdOVdfFNO6R0sVGneYvF9u3Tl3s7URE34i92IebNovXlXvUtV1wO8A/6NYQ6p6t6puVtXNXV1d1fU0aLxIqoo8rKdz957XzJxOaJ9tkXtUjGvcu30dPqGWSaEKuXpsPW85pDvM1kbk+HHuA8CqgucrgQNljr8XeE89nYqEKjTuUDxyrxmRvESut/Y2DP8M9QKSv6moPLmcklPGde5Q50SeSrlNWszWRsT4ce5bgfUislZEOoBbgS2FB4jI+oKn7wJeDq6LITHUB7M7Ycbcioeq6iS1DNQ54MFNKnapHg1DfTD/QndTUQWy+bUUTy0DdU7kYLY2YqFizl1VMyJyB/AQkAa+rqo7RORzwDZV3QLcISJvB8aAIeBDYXY6EIb9S+O8sR1Y5A7u3L2POdWOD529UQdVatyBcbUMkFfMpGs//6I1cOC52j9vGDXgZ0EVVX0QeHDKa58t+P8nA+5X+Az1wYVX+DrUk8OlC6K5uiVyC9fA+VNwZhDmLKmvLaM8Q32w9i2+DvWuyAKP3M8OwbkRX3fIGkYQtOYdqrksjOyvSj0BIUTuYLnYsMmMwolX/EfuWS9yl4LIPShbW2rGiI7WdO4nXoFcxvdNLZncxIAfj9yDqC8Dpn8Om5EBQKuwtbsiczp3NzwCqS8DppgxIqU1nXu1GvdsQeSeDvBSvbAvRjhUUQ0SCnPuAUbuZmsjBlrTuY/XcfevcQdIpwNUy8yYC7OXWDQXNlXUcYcSOfd6r9JmLYIZ883WRqS0pnMf6gVJwYLKumcIKecOeYlcb/3tGKUZ6nU3Ec1b7uvwSWqZdECL5yJmayNyWtS598H8lZBu93V4KGoZsFrfUTDU524iSvn7qYeilgF3lWi2NiKkNZ17FRp3CDFyX7TGLfjlsvW3ZRSnCo07QLZgIg8s5w5uIh/ud/c1GEYEtKZzr6KOO5RQywSVlsmNuXrfRjhUUccdpkbuAallwNk6cxZOHam/LcPwQes597GzrgRsTZF7irYgB7xp3cPl3Ak4O1jdRB6Gzh3M1kbktJ5zH+53jz51zzB5wAemcy/sg6kowmFcFdXt+yPjE3m6MOce0PpKYZ8MI2Raz7lXqXGHKTn3oHTu4NQ6krKFtrCoso47FKbgCmrLBDGRL1w9uU+GETKt59yr1LhDgVomHbBaJt0O81dYNBcWVWrcYfJEHqhapn0WzF1mdyQbkdF6zn2oF9pmuoHmk9DUMmD65zAZ6nU3D81a5PsjhbLXQOq5F2Klf40IaU3nvnB1VWV2Q1PLgOmfw8RTRVVh62xYahkwWxuR0nrOfbivqgU2CFEtA64vpw45FY8RLFXezwCTJ/JA1TLgbH1iALJjwbRnGGVoLeeuWrXGHUKO3L2+DO8vf5xRHZ6tq53IsxMTeaBqGXC21ly+UqVhhEtrOfezQzB6oupozhvck3PuAQ140z+Hw6kj7qahOiby4CN3s7URHa3l3If2ucdFa6v6WGg698K+eH0zgsH7ey6uztbFde5ma6PxaC3nPugN+Iuq+ljhgA9U5w4wdym0z4HBnmDaMxze37NKW09Sy3jlnYOayOevgPQMs7URCa3p3KvMwxarFBjYpbqIc0CDFs0FyuA+kLTvss4ek9QyQU/kqZT77ZmtjQhoMefeA/MuhI7ZVX1sUo3voNUy4FIHFs0Fy2APLFgJbR1VfSzUnDvYRG5ERus59ypzsBBy5A6uT0O9Vvo3SAZ7qk7JwGTZa+BqGcjbep+V/jVCx5y7D7zBnQpDLQPOCeXGTCIXJDU6d2/STqUgLSFF7mNn4NTh4No0jCK0jnMfPQWnj9Q14MOL3PN9MhVFMJwZhHPDtUXuWU/2mgo+5w4TwYWl4YyQaR3nPlSbUgYKc+4FkXtQCorCPtmAD4YaZZAQUc4dzNZG6LSOc/cGU5UadyheKTDQAT/vQpPIBUmNkleAnIZUFdJjwSqn4rFFVSNkWs+515Rzn4jmRNygD3TAm0QuWMYn8u6qPzo5cg9BGZVud4XrbCI3QqaFnPs+mN0JMxdU/dFMgYIC3MAPNHIHk8gFyeA+dzXUPqvqj07UlhHygXtItjbnboRLCzn32tQTMDlyBzfwA1XLwMSAN4lc/dRh68yUq7TwbG1ySCNcWsi576spJQMTt597C2zhRO5rXaGrk4eCbbcVqVHyCm4i9xw7hGjr0RFXyM4wQsKXcxeRG0Vkl4jsEZE7i7z/KRHZKSI/FZGHRaS6UnxhM3YOTrxSR+SeQ8Tp3MGL3EMY8GCX6/UyerJmySu4yN27QoO8rYNURoEpZoxIqOjcRSQN3AXcBGwEbhORjVMO+wmwWVUvA+4HvhB0R+tiuA/QugZ8W8GAT6dS4eRhwbTu9eKV0605cs9NsXVI6ytgzt0IFT+R+5XAHlXtUdXzwL3ALYUHqOqjqnom//RpYGWw3ayTGisEemSjiOYWrIZUmw34eqnT1tMi93Qq+Ku0hWsAsQV0I1T8OPcVQOE2QQP510rxUeB79XQqcOrQuIMXuU/8qUKJ5tJtTgNtzr0+6rR1dtpVWgi2bp/pyv+arY0QafNxTLHdhYv+2kXkV4DNwLUl3r8duB1g9erVPrsYAIP7YMYCmL24po9Pi9zTISgowOSQQTAueZ1f08dd5D4xkYeilgGrBGqEjp/IfQAoLIq9Ejgw9SAReTvwe8DNqjparCFVvVtVN6vq5q6urlr6WxueekKKzVOVyUSRhwWTyAVBHTJIcDr3yGxt6ytGiPhx7luB9SKyVkQ6gFuBLYUHiMgVwFdxjv1I8N2sk3oHfLGce1gDfnTEFb4yamNwX122LqqWCcvWp4/CuRPBt20Y+HDuqpoB7gAeAl4C7lPVHSLyORG5OX/Y/wbmAv9XRJ4XkS0lmoue7BgM99esngCncw9dLQMmh6yXOiWvkFfLpKOI3G0/VSNc/OTcUdUHgQenvPbZgv+/PeB+BcdwP2i25gU2yEfu6YiiOXDOfdUbgm+/2RnqxUle65jIp0XuqeCVUTDZ1ssvD759o+Vp/jtUj73sHjs31NxEJGoZcBOQpOHY7uDbbgW8v9uSV9XcRCRqGYDF6wCZ+H0aRsC0gHPf5R67anfuxXPuISgo2jpc1On12agO7+9W50Q+SS0TljKqYzYsXAVHzdZGODS/cz+6G+YshVmLam6iqFomjEt1gM6LXZ+N6jm6G+avhBlza24issgdnK1tIjdCovmd+7Fd0HVxXU0U17mHNOC7NsDgXrcQbFTHsV11XaFBhGoZcL/LY3sgjCsDo+Vpbueu6qK5Oi7TIaLaMh6dF0MuYzczVUsu5/LXnfVO5BHd0wDud5k5CyP94bRvtDTN7dxPHXa68aAj91CjufxEZJfr1XFiAMbO1B+5Z4uoZcKM3MHScEYoNLdzP1r/Aht4OvcI1DIw0VdbaKsOz0HWHblrNDp3mLC1TeRGCDS3c/ekcaFE7iHlSWfMc0WlTA5ZHeOqqPpsHVltGXC1jmZ32kRuhELzO/eOeTBveV3NZKK6a9Gjc4M592o5thtmLYY5nXU1U1QtE5YyCvKLqmZrI3ia27kfzasnaiwY5hFpzh3yA/5lKyBWDUd31x21Q7F67iHbunOD+52arY2AaW7nfqx+pQyUUMuEGc11rofzp1ydFMMfx3a5v1udTFfLhLigCm5COjfsiogZRoA0r3M/NwInDwbi3COP3L1FQcvF+uP0cThzvO7FVCiucw89BQdmayNwmte5ezU7ArpUn6SWSYc84L0+Wy7WHwEtpkLxnHvokTuYYsYInOZ17uMyyGAGfGRqGYA5XTBzoUVzfglI8gqezn2yWiYTpq3nr4COuaZ1NwKneZ37sV2Q7oBF3XU3FdlOTB4ipqKohmO7oX2224O2TiKP3EXcWoFF7kbANK9zP7rblVVN+ypZX5bslLsW0xLygIcJFYVRmaO7XJnfVP0/50yR2v2hTuSQt7VN5EawNK9zD6CIlEdm6l2LYefcwUXuZ47Zlnt+OBaMDBJKqGXCVEaBc+4nD9iWe0agNKdzHzvnduUJIN8OMahlwBQzfhk9BSP7A7N1MZ17JBM52MYdRqA0p3Mf3AuaCyyam74Tk9M+a5g3nlgBMX8c91RRwVylRZ5zh4mJyWxtBEhzOvfDO9zj0ksCaa5Y5O69HhoLVjsVhfddjOKM23pjIM0Vqy0TqloG3H6qbTPN1kagNKdzP7jdDZbALtWnq2Xc6yE691QKLrjMfRejNAe3u0lw8bpAmisWuecUcmHaOt0GyzaZrY1AaV7nvmxTIEoZiClyB1h+GRx6AXLZcM/TyBzcDhe8JhCljKqWtnXYtV8uuAwO/tRqzBiB0XzOXdUNkuWXB9bk9NoyEUTu4L7D2Bk4vjfc8zQquSwcejEwW3uT9VS1TOF7obH8crexzFBvuOcxWobmc+5DvW6QXHBZIM3lcooq0/KwEMGA976DXa4X5/heGDsdmHP3JuupOvfC90LD+w5mayMgms+5H/qpewx4wE/Wuafy74W80NZ1MaRnwCEb8EXxHGFAE3nxyD0/kYetdV+6ESQ98fs1jDppPud+cDuk2gJTT3gDPpace7rdFtrKcWi7m/wClLzClKu0tBe5hzyRt8906i6ztREQzencu17tBksAeIO6aDQXtnMHt6h6cLsttBVjfOG8PZDmykbukdj6crO1ERjN5dxV3eAIcDE11sgd3Hc5NwLD/eGfq5EIwdbeRF7M1qHn3MF9l9NH4eSh8M9lND3N5dxPHnKDI2ClDBSP5iIb8GCX61MZ7nOT3vJg8u0Qs1oGbAHdCJTmcu4BL7BBYeReqJaJcMAv3eQW2mzAT+ZgsAvnwPjWibFF7hdcCojZ2giE5nLuh34KSH6QBEPZyD1sBQW4tYOuV5uKYioHt7tJb+mmwJrMFlNGjafgQl5QBZgxz5UuNlsbAeDLuYvIjSKyS0T2iMidRd5/q4g8JyIZEXl/8N30ycHtbnDMmBdYk9ky0VwkkTtMLLQZEwS8cA4l1DJRRu4wsYBuGHVS0bmLSBq4C7gJ2AjcJiJTdYb9wIeBbwbdwap1Pb0QAAAPdUlEQVQ4uD3QHCwUqGWm1HMvfC90ll8Gpw7bQlshh4K9CxnKq2UiuUoD951G9lsdf6Nu/ETuVwJ7VLVHVc8D9wK3FB6gqr2q+lMgIm9XhDODblCENOBjj9xhIs/c6pw85Ca7gG1dVC2TjsvWFr0b9eHHua8A9hc8H8i/VjUicruIbBORbUePHq2lidIcfN49BriYCglQy4ArjAUT37HVOZD/OwR8lVZOLROdrT3FjNnaqA8/zl2KvFbTL11V71bVzaq6uaurq5YmStP3FEgKVrw+0Ga9AZ+Swsg9QrUMuDWEpZug78lozpd0+p+EVDssf22gzXoOPDVlv1yI0NazF8OS9e73bBh14Me5DwCF28qvBA6E05066HvCXdLOnB9os0Vry0QduQN0Xw37n4HsWHTnTCq9T7hJvGN2oM2WzblHtb4Cztb9T1mpZ6Mu/Dj3rcB6EVkrIh3ArcCWcLtVJWNnYWArdF8TeNPZ8TxssaqQUQ74a1z53wM/ie6cSWT0pPsbhGDrojr3qHPuAN1vgdETJok06qKic1fVDHAH8BDwEnCfqu4Qkc+JyM0AIvIGERkAPgB8VUSi3S9sYCtkz7tBETDZvP+OVUEBsOZq99j7WHTnTCL9z4BmQ3HuOfUi98L9cmNw7uO2fjy6cxpNhy+du6o+qKobVHWdqn4+/9pnVXVL/v9bVXWlqs5R1SWqGtydJX7ofcLl21e/KfCmE6GgAJjTCV2XuO/ayvQ97qp+rroy8KYzSVBGAcxf7rYNbHVbG3XRHHeo9j7uVAYzFwTedLE8bOQ3tnh0XwP9T7d23r338Xy+fU7gTWfLVACNxdZ9T1re3aiZxnfuY+dCy7dD8Wgu0mJShXRf7XYeOtCiMrnRU/DKcxNpi4ApXlsmLltf43YUO/RCtOc1mobGd+6vbIPsaGjO3Ss/0BbnLekea/Lfsa9Fc7H7w8u3Q/naMtHbOj+B9VlqxqiNxnfuvY8DAquvCqX54pF7DGoZgLldrp5Kqy609Xr59jeG0nyxG9ZiUUYBLFgBiy9qXVsbddMczv2C18CshaE0Xyyaiy1yBxfR9T8N2Uz0546b3sfhwitgxtxQmi9W3jkWZZTHmqtd5G55d6MGGtu5j+fbg5dAehRTy8Qij/PovgbOn2q92iPnT8OB50JLyUCJyD0OZZRH91vchiSHo1UWG81BYzv3V56FzLlQB3xxtUy+3kgc0Zz3XXt/FP2542T/M5DLhGzr0hN5LFdp3XZvg1E7je3cd3/f5WDXvDm0UxTNuccZzc1d6urM7P5B9OeOk90PQXoGrAr+XgaP4jn3mNQyAAtWujozux+K/txGw9O4zl0Vdv4LXHR9aPl2KIzcE6CW8dh4i6s9cuJgPOePmlwOdn4H1r8jtHw7FC/vHGvkDs7WvY/B6WPxnN9oWBrXuR94Dob7YdN7Qz1NotQyHpveAyi8lKwSP6Gx/xk4eTB8W5eRvcZn6/eC5lrH1kZgNK5z3/GAK/v66v8U6mmy2SJ3LUrM0VzXxS41s+OBeM4fNTsecCmZDe8M9TTjkXsSdO4eyza51Eyr2NoIjMZ07qqw419g3fUwa1Gop8oUGfCplJCSmPKwHpvem0/NJK/6cqBMSskEtzduMcrq3ONYPAcQcbbufRxOBbzBjdHUNKZzf+VZt6VeyJfpUFwt456n4ovmIJ+aAXY2+eX6/qfh1KGIbJ0wtYzHpvdYasaomsZ07l5K5uJwUzJQPOfuPY81cu9cD8subf7L9R0PQNtM2HBj6Kcat3XBrlsiEr+tl26Ezg3Nb2sjUBrPuXspmVe9LVSVjEcxtYx7LvHo3AvZ9B4X2Y68Em8/wiKXjUQl45HNKSmZvM0euIk81sjdS830PQEnD8fXD6OhaDznPrANTgxEcpkOBftqTtlJNp2W+BQUHhvzf4Od34m3H2HR/xScOhypradO4uAm8thtbaoZo0oaz7nvfRjSHXDxTZGcLpvL0ZYSRKbm3GOO5gA6X+X2jX3uHrfw2Gw8ew+0z4H14apkPLI5nZZ+gwRE7gBLL3HpmWfvcVevhlGBxnPu1/4OfPyZUDbmKEamzICPNQ/rcdUdcPRnsPt7cfckWAb3wYv3w+aPRJKSAadzn7pwDl7kngRbfxwOvwB7/j3unhgNQOM5dxFXCjUisiUHfMxqGY9N74OFa+CxP2+uiO7Jr7jSElfdEdkps7ncJMmrRzoptn7NL8L8lfDYF+PuidEANJ5zj5jER+7pNrjmN508dF+TFBM7eQh+8o/w2l92+4lGhMu5l4jc4148B2jrgDf/BvQ/CX1Pxd0bI+GYc69ANqe0pYsvsiUimgO4/Jdh7jJ4vEkiuqfuchUgr/5kpKdNdM7d43W/CrOXNI+tjdAw516B8pF7QhYx22e69EXPf7gIvpE5OwTbvu7STRGm36CMWiYJyiiPjtnwpo/Byz+Agz+NuzdGgjHnXgFPLTOVdBJ07oVs/gjMXAiP/lFj596f+IrbjOSa34r81A0RuQO84dehY17j29oIFXPuFSgVubtoLkEDa8Y8uO4zTknx7N/H3Zva2P9jeOLLcPltcMGlkZ++bM49SbaetRCu/W2nkHr+m3H3xkgo5twrkC0x4BOjoCjkyttdffuHfheO7Ym7N9UxehK+fbvbGPqmL8TShWwuVyJyT6Ctr/q424bve7/tZKOGMQVz7hUoGbknLZoDSKXgPX8NbTPg278O2bG4e+Sf798Jw33w3q/CzPmxdCGTbSRbp52tJQ0P/JfW3DDdKIs59wo4nfv0P5PLwyZkka2Q+cvh3V9ym5k88odx98YfL37bSR+v+a1Qt0yshFNGNUDO3WPhKnj3F91mJj/807h7YySMtrg7kHTKRe5j2QQ6d3AFxXo+DE98yT1/+++7m7+SyPZvwXf+G6x8A1x7Z6xdcbZOaG2ZUrzm/bDnYfhRPpV1/e8m19ZGpJhzr0A2lysZzZ0dS2A05/GuL7pL9ie+BGeOwbu/7G54ShJP3eXWB9a+FX7pG+4mnRgpvb6SMGXUVG7+S5em+dEX4PRReNefu+dGS5Ow0Z48GirnXkgq7Qb5nC744Z/AuRH4wD+4vHzcnBmER/7A6dk33gLv+5pbJ4iZTIkF1ba0MDqW0Mgd3KR981/CnE54/C9g9AS872+SYWsjNnxZX0RuFJFdIrJHRKZdO4vIDBH5Vv79Z0SkO+iOxkVZtUySozlwl+fXfwZ+7g/hpX+NPy+bOQ9P/RV85Qon17zqDnj/3yXCsUODKaOmIuLSbzf8T3jxn+0OVqNy5C4iaeAu4B3AALBVRLao6s6Cwz4KDKnqq0TkVuBPgV8Ko8NR07CReyFX3QGHd7oIfsXrQt9ouiiHXoB//jVXwfKi6+Gdn3ebPyeITE6Z2d7gtn7Lp+HITreYfuEVblMboyXxE7lfCexR1R5VPQ/cC9wy5ZhbgHvy/78feJtMLYDeoGRL3JKeTidULVMMEaequOAyJ5Ec7Inu3Lmcy61/7QZXWuC2b8EHH0icY4cKOfdGce4iLkWz9BI3mQ73x90jIyb85NxXAPsLng8Abyx1jKpmRGQEWAIcC6KThdy3dT9feyw659Q/eIY3XrRk2uttKaF/8Azv+OIPI+tLvVyQu4O7zn0K/cu3MiThb1EIMENHWa5HeLLtSv6cT3Di32bBvyWzemXf4Bneur5z2uttKWHvkVMNZesV2U9w19lPMfblqxmWaPY+MPxz/PW/yevf9WuhnsOPcy8WgU8NY/wcg4jcDtwOsHr1ah+nns7C2e2sXxbN5g0A65fN5ZbXrpj2+gdevyq5UsiSvIq7z/4R1w0/QIpsZGd9ePav8uT8m1gmwrLIzlo965fN5Rdet3La67/4hlUNqC7cwFfPfp5rh7+D0Gi/0+anY+7i0M8hWqHwkIhcBfy+qr4z//wzAKr6xwXHPJQ/5ikRaQMOAV1apvHNmzfrtm3bAvgKhmEYrYOIPKuqmysd5yfnvhVYLyJrRaQDuBWYukvvFuBD+f+/H3iknGM3DMMwwqViWiafQ78DeAhIA19X1R0i8jlgm6puAf4W+D8isgcYxE0AhmEYRkz4uolJVR8EHpzy2mcL/n8O+ECwXTMMwzBqxW5hMwzDaELMuRuGYTQh5twNwzCaEHPuhmEYTYg5d8MwjCak4k1MoZ1Y5CjQV+PHOwmhtEED0IrfuxW/M7Tm927F7wzVf+81qtpV6aDYnHs9iMg2P3doNRut+L1b8TtDa37vVvzOEN73trSMYRhGE2LO3TAMowlpVOd+d9wdiIlW/N6t+J2hNb93K35nCOl7N2TO3TAMwyhPo0buhmEYRhkazrlX2qy72RCRVSLyqIi8JCI7ROSTcfcpSkQkLSI/EZHvxt2XKBCRhSJyv4j8LG/zq+LuUxSIyG/lf98visg/icjMuPsUBiLydRE5IiIvFry2WET+n4i8nH9cFMS5Gsq5F2zWfROwEbhNRDbG26vQyQCfVtVLgDcBH2+B71zIJ4GX4u5EhHwZ+L6qvhq4nBb47iKyAvgEsFlVL8WVFm/WsuF/D9w45bU7gYdVdT3wcP553TSUc8ffZt1NhaoeVNXn8v8/iRvs0/f9a0JEZCXwLuBv4u5LFIjIfOCtuP0RUNXzqjocb68iow2Yld/JbTZwIOb+hIKq/gi350UhtwD35P9/D/CeIM7VaM692GbdLeHoAESkG7gCeCbenkTGl4DfhpbZBPQi4Cjwd/lU1N+IyJy4OxU2qvoK8GdAP3AQGFHVH8Tbq0hZpqoHwQVzwNIgGm005+5rI+5mRETmAv8M/Kaqnoi7P2EjIu8Gjqjqs3H3JULagNcBf62qVwCnCegSPcnkc8y3AGuBC4E5IvIr8faq8Wk05z4ArCp4vpImvXwrRETacY79G6r67bj7ExFXAzeLSC8u/XaDiPxjvF0KnQFgQFW9K7P7cc6+2Xk7sE9Vj6rqGPBt4M0x9ylKDovIcoD845EgGm005+5ns+6mQkQEl4N9SVW/GHd/okJVP6OqK1W1G2fnR1S1qaM5VT0E7BeRi/MvvQ3YGWOXoqIfeJOIzM7/3t9GCywkF7AF+FD+/x8CvhNEo772UE0KpTbrjrlbYXM18EHgBRF5Pv/a7+b3tTWaj98AvpEPXnqAj8Tcn9BR1WdE5H7gOZw67Cc06d2qIvJPwHVAp4gMAP8L+BPgPhH5KG6iC2Q/artD1TAMowlptLSMYRiG4QNz7oZhGE2IOXfDMIwmxJy7YRhGE2LO3TAMowkx524YhtGEmHM3DMNoQsy5G4ZhNCH/H65x95OS1iZ7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d3ed510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def convolve1D(y, sigma, discr):\n",
    "    gx = np.arange(-3*sigma, 3*sigma, discr)\n",
    "    gaussian = np.exp(-(gx/sigma)**2/2)\n",
    "    z = np.convolve(y, gaussian, mode=\"same\")\n",
    "    return z\n",
    "\n",
    "x = np.arange(0,10, 0.1)\n",
    "y = np.repeat(sm, x.shape[0])\n",
    "y[[10,50]] = 0.5\n",
    "\n",
    "z = convolve1D(y, 0.5,0.1)\n",
    "\n",
    "plt.plot(x,y)\n",
    "plt.plot(x, z)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Divergence d1 to u 1.3495510787832004\n",
      "Divergence d2 to u  1.6078272976845267\n"
     ]
    }
   ],
   "source": [
    "sigma = 0.2\n",
    "d1_blurred = convolve1D(d1, sigma, 1)\n",
    "d2_blurred = convolve1D(d2, sigma, 1)\n",
    "\n",
    "\n",
    "#normalize to make it a probability distribution after the convolving\n",
    "d1_blurred = d1_blurred / np.sum(d1_blurred)\n",
    "d2_blurred = d2_blurred / np.sum(d2_blurred)\n",
    "\n",
    "print \"Divergence d1 to u\", KLDivergence(d1_blurred,u)\n",
    "print \"Divergence d2 to u\", KLDivergence(d2_blurred,u)\n",
    "print \"Divergence u to u\", KLDivergence(u,u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entropy \n",
    "\n",
    "Found the biggest source of my confusion about entropy.\n",
    "When thinking about entropy the **less information the process \"hides\", the better** because there is less uncertainty about the outcome. If we have a small amount of information, there is nothing there to guess. That is why the fair coin toss has (better hides) the most information since we don't know anything more about what to expect next if we have already seen an outcome. Probably should consider the term information as \"information you don't know\" rather than the one you know.\n",
    "\n",
    "This definition of information is super confusing if you got used to the term \"information\" as the inverse of \"uncertainty\", where the less information you have the more is your uncertainty.\n",
    "\n",
    "### Random process perspective\n",
    "Let's assume there exist a random process (black magic) that produced the distribution we need to work with.\n",
    "By wikipedia definition, entropy is an average amount of information produced by a random process.\n",
    "This means that given a distribution, we can compute the entropy as  $$H = -\\sum_i^n p(x_i)\\mathrm{log(p(x_i))}$$\n",
    "\n",
    "### Data sending perspective\n",
    "Let's assume we need to send some information (string/text) though a channel, so that we need to minimize the amount of sent information (bits).\n",
    "Then we might want to reduce or decrypt the string. And entropy tells us the average amount of bits per letter that are needed to encrypt this string. (Side note: appearently *information theory* will tell how to exactly map those letters to the bits)\n",
    "If all the symbols occur with the same frequence/probability than there is no better way to reduce the string in less bits than there already are in a string (?) If on the other hand, someone (probably the same black magic from above) tells us the more precise probability of the letters, than we would assign less bits to the letters that occur often and more bits for the rare ones and the average bit information is captured as $$H = -\\sum_i^n p(x_i)\\mathrm{log_2(p(x_i))}$$\n",
    "\n",
    "### Other thoughts\n",
    "The bigger the entropy the more unpredictable the states / variables / values in distribution are?\n",
    "The biggest entropy is when every state is the same possible, e.g. uniform distribution. Since then every state needs the same amount of bits to be encoded.\n",
    "The smallest entropy is when the same event occur always/never.\n",
    "\n",
    "**Intuition behing log**: The more probable the event is, the less information it brings. Since if we observe the same event (variable) all the time, it doesn't tell us about other possible outcomes. Thus, the more rare the event is the more information it brings.\n",
    "That's why **-log(x)** function is a good choice, since it is monotonically decreasing for increasing values of probability. \n",
    "Information should also have other 3 properties discribed on wikipedia under Rationale.\n",
    "\n",
    "*Question*: if entropy is an average amount of information per state ($x_i$), why can't we just average the -log(x_i) like  $\\frac{1}{N}\\sum_i^n(-\\mathrm{log}p(x_i))$ ?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## visualization of log just to remind how the log looks like :)\n",
    "prob = np.arange(0.0001, 1.0, 0.01)\n",
    "f = -np.log2(prob)\n",
    "plt.subplot(121)\n",
    "plt.plot(prob, f, 'r-');\n",
    "plt.xlabel('prob');\n",
    "plt.ylabel('log');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy of a distribution 1.5632692440310552\n",
      "Average log information  21.449236676040247\n",
      "Entropy of a uniform distribution 3.321928094887362\n"
     ]
    }
   ],
   "source": [
    "def computeEntropy(dist):\n",
    "    ent = 0\n",
    "    for el in dist:\n",
    "        ent += -el * np.log2(el)\n",
    "    return ent\n",
    "\n",
    "\n",
    "print \"Entropy of a distribution\", computeEntropy(d1)\n",
    "le = 0\n",
    "for el in d1:\n",
    "    le += -np.log2(el)\n",
    "print \"Average log information \", le/d1.shape[0]\n",
    "print \"Entropy of a uniform distribution\", computeEntropy(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connection between KL and Entropy\n",
    "\n",
    "Based on the wikipedia article (https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) and Kraft-McMillan theorem, we can relate entropy and Kullback-Leibler divergence.\n",
    "$$ D_{KL}(P|| Q) = -\\sum_x{p(x)\\mathrm{log}(q_x)} + \\sum_x{p(x)\\mathrm{log}(p_x)}$$ where the first term is $H(P,Q)$ cross entropy(?) and second term is an entropy $H(P)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}