density-ratio-trick.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Density Ratio Trick\n",
    "### (Trik Ilorazu Funkcji Gęstości Prawdopodobieństwa)\n",
    "\n",
    "work based on Louis Tiao article: https://tiao.io/post/density-ratio-estimation-for-kl-divergence-minimization-between-implicit-distributions/\n",
    "\n",
    "---\n",
    "\n",
    "Density Ratio Trick to technika określająca podobieństwo pomiędzy dwiema funkcjami gęstości prawdopodobieństwa\n",
    "(estymuje ich iloraz). Używana jest przy wyliczaniu wartości Dywergencji Kullbacka-Leiblera będącego miarą podobieństwa dwóch rozkładów.\n",
    "\n",
    "W prostym i akademickim przypadku mamy do dyspozycji dwie funkcje gęstości prawdopodobieństwa opisane\n",
    "równaniami. Policzenie ich ilorazu nie jest problemem. Możemy naiwnie wyliczyć ich stosunek, przechodząc po każdym punkcie X lub skorzystać\n",
    "ze wzorów, jeśli wiemy, że rozkłady są Gaussowskie i znamy ich wariancję i średnią.\n",
    "Zadanie jest utrudnione w przypadku gdy nie dysponujemy postacią analityczną rozkładu a jedynie próbkami (samplami) z tego rozkładu.\n",
    "W takiej sytuacji często niemożliwe jest wyznaczenie analitycznej formy pierwotnego rozkładu.\n",
    "\n",
    "Zacznijmy jednak od początku, czyli sytuacji, w której dwa rozkłady są określone funkcjami. Ich iloraz może zostać wyznaczony\n",
    "bez większych problemów."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9pElEQVR4nO3deVxU5f7A8c/DFiCI5oLmvivu4HpdE0PFLUszTW0xTdOym/3qmuVNq9tiVrfbVa920RbTlDJLRU3FXK77nhuSC+4LomwiAs/vjwMIyKbMcGaG7/v1Oi9m5pw58z0wfOeZ5zzP9yitNUIIIeyfk9kBCCGEsAxJ6EII4SAkoQshhIOQhC6EEA5CEroQQjgISehCCOEgCkzoSqlqSqlwpdQRpdQhpdSEXLbpqpS6oZTal75MsU64Qggh8uJSiG1SgIla6z1KKW9gt1LqN6314RzbbdJa97F8iEIIIQqjwBa61vqC1npP+u044AhQxdqBCSGEuDeFaaFnUkrVBFoC23NZ3V4ptR84D7ymtT6Uy/NHA6MBPDw8AqpVq3bPAZstLS0NJ6eSdepBjtnxlbTjBfs95oiIiKta6wq5rtRaF2oBvIDdwGO5rCsNeKXfDgaOF7S/gIAAbY/Cw8PNDqHYyTE7vpJ2vFrb7zEDu3QeebVQH09KKVfgR2CB1vqnXD4UYrXW8em3VwKuSqny9/jBI4QQoggKM8pFAf8FjmitP81jm0rp26GUapO+32hLBiqEECJ/helD7wAMBw4qpfalP/YmUB1Aaz0bGAiMVUqlADeBJ9O/GgghhCgmBSZ0rfVmQBWwzZfAl0UN5vbt25w9e5akpKSi7spqfHx8OHLkiNlhFBt3d3fSv3wJIWzcPY1ysbazZ8/i7e1NzZo1bTaJxMXF4e3tbXYYxUJrTXR0NKVKlTI7FCFEIdjUmJ2kpCTKlStns8m8pFFKUa5cOZydnc0ORQhRCDaV0AFJ5jZG/h5C2A+bS+hCCCHujyR0C6tZsyZXr169p+e88sorbNy4Md9tunfvTkxMTFFCE0I4OEnoedBak5aWZvXXuXbtGtu2baNz5875bjd8+HBmzpxp9XiE7dNaM3PmTH799VeuX79udjjChkhCz+LUqVM0atSIF198EX9/f86cOcP06dNp3bo1zZo14+9//3vmto8++igBAQE0btyYOXPmFLhvLy8vJk6ciL+/P4GBgVy5cgWA0NBQevbsCcCNGzdo0KABx44dA2DIkCHMnTsXgH79+rFw4UJLH7KwA7dv32bWrFlMnz4dMAYPjBs3jn79+uHr68vzzz/PhQsXTI5S2AKbGraYzSuvwL59lt1nixbw+ef5bnLs2DHmzZvHzJkzWbNmDcePH2fHjh1orenXrx9btmyhZ8+ehISE8OCDD3Lz5k1at27N448/Trly5fLcb0JCAv7+/syYMYNp06YxdepUvvzyS7Zs2cLAgQMBY4z7l19+yTPPPMOECROIiYlh1KhRAJQtW5Zbt24RHR2d7+sIxxIVFcWTTz7J1q1bGTJkCGDMDbh48SLHjh1j0aJFhISEEBoayvLly+nYsaPJEQsz2W5CN0mNGjVo164dAGvWrGHNmjW0bNkSgPj4eP78808AvvjiC5YuXQrAmTNnOH78eL6J1snJicGDBwMwbNgwHnvsMQAuXLhAhQp3Cqc98sgjLFmyhHHjxrF///5s+6hYsSLnz5+XhF5CHD58mMDAQBISEvj+++958sknAWPkka+vL76+vnTu3JlXXnmFjz/+mBYtWpgbsDCd7Sb0AlrS1pJ1Eo3WmkmTJvHCCy9kPhYXF8eGDRtYu3YtW7duxdPTk65du97z7NaM4YAeHh7ZnpuWlsaRI0fw8PDg2rVrVK1aNXNdUlISHh4e93towo5ER0fTtWtXnJ2d2bp1K40bN85z2/r16/PVV18BxjfBdevW0a9fv+IKVdgQ6UPPR48ePQgJCSE+Ph6Ac+fOceXKFW7cuEHZsmXx9PTk6NGjbNu2rcB9paWlERoaCsD333+f+dW4UaNGREZGZm732Wef0ahRIxYuXMhzzz3H7du3AePD5eLFi9SsWdPCRylsUbly5fjkk0/YsGFDvsk8p/fff59HH32UlStXWjE6YaskoecjKCiIoUOH0r59e5o2bcrAgQOJi4ujZ8+epKSk0KxZM95+++3MLpr8lCpVikOHDhEQEMD69euZMsW47Grv3r3ZsGEDABEREXz11VfMmDGDTp060blzZ9577z0Adu/eTbt27XBxsd0vVaLotNacOnUKgBEjRtCgQYN7ev7bb79Ns2bNGDZsWOZ+RAmSV6F0ay+5XeDi8OHDRav8XgxiY2Pv63mlSpXKc12HDh10TExMvs9/+eWX9dq1a+/rtYtqz549pryumcy6+MHs2bO1u7u73rdv333vIzIyUvv4+Oj27dvrlJSUQj3HXi/2UBT2eswU9QIXwrpmzJhBVFRUvts0adKEwMDAYopImOHEiRNMnDiRjh070rRp0/veT506dfjXv/7F1q1b+fLLIhdBFXZEvr8Xk4x++Ny0bdu2wOdnDF8Ujuull17CycmJkJCQIl/rctiwYURERNC7d28LRSfsgSR0IWzAihUrWLlyJTNmzMASF09XSvHuu+9aIDJhT6TLRQgbcPDgQRo3bsz48eMtut9Lly7x2GOPsWXLFovuV9gmSehC2IC//e1v7N69Gzc3N4vu18vLi+3bt/Pqq6+i5aqQDk8SuhAmunnzZuY8hgceeMDi+y9VqhTvvfceO3bsYPny5Rbfv7AtktAt7H7K537++ed88803+W7z5JNPcvz48aKEJmzQnDlzaN++PfssXbcoi2HDhlG7dm2mTZsmrXQHJwk9D7qYyuempKQQEhLC0KFD891u7NixfPzxx1aPRxSfxMREPvjgAx5++GGr1mFxdXXlzTffZNeuXYSFhVntdYT5JKFnYc3yufPmzaN+/fp06dKFUaNGZZ78Wr9+Pf7+/ri4uJCSkkLr1q0zZ45OmjSJyZMnA9CpUyfWrl1LSkqK5Q9cmOLrr7/m0qVLvPPOO1Z/rREjRjBt2jRatWpl9dcS5rHZYYsmVc+1SvncCxcu8Pe//53du3fj4+PDww8/nFnBccuWLQQEBADg4uLC/PnzGThwIF988QWrVq1i+/btgFGtsW7duuzfvz9ze2G/0tLS+Oyzz2jTpg2dOnWy+uu5urry9ttvW/11hLmkhZ5DXuVz/f39OXr0aLbyuc2bN6ddu3aZ5XPzsn37drp27UqFChVwc3PLLKMLd5fPbdy4McOHD6dv376EhIRkG/WQUT5X2L+IiAiuXLnCq6++WqwX4l69enVmfSDheGy2hW5S9Vyrlc/N6582Z/lcMMYklylThkuXLmV7XMrnOo6GDRty5swZ3N3di/V1w8PDmT59OsOGDZPKnQ5IWuj5sFT53LZt27Jhwwaio6O5ffs2S5YsyVyXs3zuTz/9RHR0NBs3buTll1/Ods3IiIiIeyqlKmxTfHw8Wmu8vLyKvXrmuHHjUEpJjRcHJQk9H5Yqn1u5cmXeeecd2rdvT/fu3fH3989c16tXLzZu3AjA1atX+dvf/sZ///tf6tevz/jx45kwYQJgzPjz8PCgcuXK1jtgUSz++te/EhAQUCyjqHKqVq0agwYN4quvviIxMbHYX19YWV5lGK29lLTyuVnNmzdPjxs3LvP+o48+qiMiIvJ9zqeffqq/+uqrIr/2/ZDyuZZz48YNXapUKf3cc89ZZf+F8fvvv2tAz58/P/Mxey0lWxT2esxI+Vzb9uGHHxZ41fYyZcrw9NNPF1NEwloWLlxIQkICo0ePNi2GTp060bt3b6vMTBXmstmToo7smWee4Zlnnsm836BBgwKvTPPss89aOSpRHObMmUOzZs1o06aNaTEopaQMgIOSFroQxWT37t3s2bOH0aNHF+tQxbwkJSWxe/dus8MQFiQJXYhi0rRpU5YsWcJTTz1ldigAjBkzhqCgoAKH3Ar7IQldiGLi5ubGwIEDKVOmjNmhADB06FCuXbsm3S8ORBK6EMXgt99+491337WpoYKBgYFUrlyZ7777zuxQhIVIQrew+ymf+8orr2SORc/L8uXLsxUHA6NsQFBQUIH7HzhwICdOnMhzfXJyMp07d5bCX1Y0c+ZMZs6caVMjS5ydnRk6dCgrV67kxo0bZocjLKDAhK6UqqaUCldKHVFKHVJKTchlG6WU+kIpFamUOqCU8s9tX/ZEF1P53GvXrrFt2zY6d+6c73a9e/fml19+ydbCW7VqFT169Mj3eYcOHSI1NZXatWvnuY2bmxuBgYH88MMP9xa8KJRr166xYsUKhg4dirOzs9nhZDNs2DBu375d4GxnYR8K00JPASZqrRsB7YBxSim/HNv0AuqlL6OBWRaNsphYs3yul5cXEydOxN/fn8DAQK5cuQJAaGgoPXv2zNxu1apVNGzYkI4dO/Lyyy/Tp08fwBhq1rVr12z9natWraJXr15orRk/fjx+fn707t2b4OBgQkNDAViwYAH9+/cH4PTp09SrV4+rV6+SlpZGp06dWLNmTebxLFiwoIi/QZGbJUuWcPv2bYYNG2Z2KHdp3rw5+/btK9Q3PWH7CkzoWusLWus96bfjgCNAlRyb9Qe+SZ/ItA0oo5Qq8hz1rl273rXMnDkTMC4OkNv6+fPnA8Y0+pzrCuPYsWOMGDGCvXv3cuzYsczyufv27WP37t2ZF9sNCQlh9+7d7Nq1iy+++ILo6Oh895uQkIC/vz979uyhS5cuTJ06FchePjcpKYlRo0bx66+/smnTJi5evJhtH61atWLTpk0ApKamcuzYMfz8/Fi6dCnHjh3j4MGDzJ07l//973+Zz8m6/xo1avDGG28wZswYZsyYgZ+fX+Y/cpMmTdi5c2ehfkfi3nz33Xf4+flZ9SIW90spRfPmzW1iGKUounuaWKSUqgm0BLbnWFUFOJPl/tn0x7JNf1RKjcZowePr65t5IYcMPj4+xMXFZd5PTU29K4akpCTi4uJITEzMd318fPxd67PuOzfx8fFUr16dxo0bExcXx/Lly1m9ejXNmzfPXH/8+HE6dOjA9OnTM1vLZ86cYd++fbRp0watNfHx8Xf1lTo5OREcHExcXByPPvoow4YNIy4ujjNnzuDp6UlcXBwHDhygevXqVKpUifj4eB5//HHmzZuXGbeXlxdRUVHExcWxfft2/P39iYuLY+3atQwYMIDExES8vb3p3LkzN2/eJC4ujnPnzuHh4ZG5j8GDB7Nw4UJmzpzJli1bsv1OXF1dOX/+PN7e3tli11rf9bdydPHx8RY55pSUFLTWdO7cmd9//73ogVlBamoq06dPJzQ0lIEDB5odTrGx1N/YlhQ6oSulvIAfgVe01rE5V+fylLsuXqi1ngPMAWjVqpXO2Wo+cuRItmSS0RrNjbe3d5HW58bLywsvL6/MGDIu3ZWzfO7u3bvZtGkT27dvzyyf6+zsjLe3N0qpbPvIGZOLiwteXl6Z23t7e2feLlWqFC4uLpnP9fDwyHbfyckp8zm///47ffv2xdvbGzc3Nzw8PDK3c3Fxybyfc5+JiYlcvHgRJycnlFLZ4kxOTqZChQq4urpmizuju6ck2bBhg8WOuXv37hbZjzW9+uqrnDlzpkRVYbTk39hWFGqUi1LKFSOZL9Ba/5TLJmeBalnuVwXs/koMliqfC8YVajL6tb///ns6duwIZC+f27BhQ06ePJl5EY2FCxdm20dERARNmjQBYN26dQQGBgLQuXNnFi1aRGpqKhcuXCA8PDzzOTnL877xxhs89dRTTJs2jVGjRmU+Hh0dnWsyF0WTs9vMVnXt2pXDhw9z+PBhs0MRRVCYUS4K+C9wRGv9aR6b/QKMSB/t0g64obXOv9qUHbBU+VwwLpxx6NAhAgICWL9+PVOmTAGM0SsZX/vc3d2ZM2cOvXv3pmPHjtSoUSPbPsLDw+nduzdXrlzB3d2d0qVLAzBgwADq1atH06ZNGTt2LF26dMl8Ttb9//777+zcuTMzqbu5uTFv3rzMfQcHBxf1VyayOHHiBJUrV+abb74xO5QCde7cGaVUZqND2Km8yjBmLEBHjO6TA8C+9CUYGAOMSd9GAf8G/gQOAq0K2m9JK59bqlSpPNd16NBBx8TE3PV4eHi47t27t9Za64sXL+pu3bpprbX+9ttv9QcffJDn/p5++mm9ZMkSrbXWiYmJum3btjolJSXf+AYMGKCPHj2a6zopn3t/PvroIw3okydPFnlf1hYeHq47deqkmzRpYnYoxcYRy+cW2Ieutd5M7n3kWbfRwLj7+kQRzJgxg6ioqHynhEdFRTFjxgyAexr+5uHhwdSpUzl37hzVq1fPdZvk5GQeffTRAis+insTGhpK69at7eZSb6NGjWLHjh0kJydnu5atsB9SPreYZPTD56Zt27a5Pp51uGXr1q0L/VoZQzczFDT5yM3NjREjRhR6/6Jgp0+fZufOnXz00Udmh1Jow4cPZ/jw4WaHIYrA5qb+G419YSvk73F/Mvqi7W0YYFpampTUtWM2ldDd3d2Jjo6WJGIjtNZER0fnOt5f5G/48OEsWrQo35ILtmjOnDm0atWKiIgIs0MR98GmulyqVq3K2bNnM6fF26KkpCTc3d3NDqPYuLu7k5CQYHYYdqdixYoMHjzY7DDuWd++fRk7diyhoaG8+eabZocj7pFNJXRXV1dq1apldhj52rBhAy1btjQ7jGJ1+vRps0OwK8uWLePSpUs8//zzODnZ1JfgAlWpUoU2bdqwbNkySeh2yL7ebULYgX/+85988cUXdpfMM/Tv358dO3Zw/rzdzw0scezzHSeEjYqJiWHjxo2ZFS7tUUbsv/76q8mRiHtlU10uQti7lStXkpqaSr9+/cwO5b75+fmxbt06OnToYHYo4h5JQhfCgn755Rd8fX3vad6ArVFK0a1bN7PDEPdBulyEsJCMYZ79+vWz2/7zDPHx8UyePJnffvvN7FDEPZAWuhAWopRi7dq1DjFu393dnf/85z+cPn2aRx55xOxwRCHZdzNCCBuSkcht7bqh98PFxYU+ffqwYsUKbt++bXY4opAkoQthAVprGjduzPvvv292KBbTv39/rl+/fs8XihHmkYQuhAUcOHCAY8eOUblykS+lazOCgoJwd3dn2bJlZociCkkSuhAW8Msvv6CUonfv3maHYjGlSpViwIABUlvJjshJUSEs4JdffqFdu3b4+vqaHYpFff/992aHIO6BtNCFKKJz586xa9cuu55MVJCkpCSzQxCFIAldiCJycXFh6tSpPP7442aHYhWDBg2S683aCelyEaKIfH19My/67Yjq1KnDzz//zI0bN/Dx8TE7HJEPaaELUQTJycn8+uuvDl0zvnfv3qSkpLB27VqzQxEFkIQuRBFs3ryZfv36OXSya9++PWXKlGHFihVmhyIKIAldiCIICwvD1dWVwMBAs0OxGhcXF3r06MHKlStJS0szOxyRD+lDF6IIwsLC6Ny5M15eXmaHYlVjx46lZ8+epKam2n3hMUcmCV2I+xQVFcWhQ4d49tlnzQ7F6rp06UKXLl3MDkMUQD5qhbhPGf3mJWVI35kzZ/jhhx/MDkPkQxK6EPfp2WefZf/+/TRs2NDsUIrF/PnzGTJkCJcvXzY7FJEHSehC3CelFM2aNUMpZXYoxaJ3795orVm1apXZoYg8SEIX4j7873//Y/To0Vy8eNHsUIpNixYtqFSpkgxftGGS0IW4Dz/99BNff/21w49uycrJyYng4GBWr14tF72wUZLQhbgPJWW4Yk7BwcHExsZy4MABs0MRuZBhi0Lco6ioKA4fPszIkSPNDqXYBQcHc/HiRSpWrGh2KCIXktCFuEdhYWFAyRmumJWHhwceHh5mhyHyIF0uQtyjtLQ0OnXqRIMGDcwOxRS7d++mZ8+enD171uxQRA6S0IW4R2PHjmXjxo0lZrhiTg888ACrV6+W4Ys2qMCErpQKUUpdVkr9kcf6rkqpG0qpfemL4xaGFiVeQkJCiS9Q1bhxY6pVq8bKlSvNDkXkUJgW+nygZwHbbNJat0hfphU9LCFs05QpU6hduzYpKSlmh2IapRTBwcH89ttvJCcnmx2OyKLAhK613ghcK4ZYhLB5YWFh1K9fHxeXkj2eIDg4mPj4eDZv3mx2KCILS70r2yul9gPngde01ody20gpNRoYDcZluzZs2GChly8+8fHxdhl3UcgxGy5evMiRI0d4+OGHHe73ca9/Y1dXV5o1a8bevXvttpyuQ76vtdYFLkBN4I881pUGvNJvBwPHC7PPgIAAbY/Cw8PNDqHYyTEbZs2apQF95MiR4g/IyuRvbD+AXTqPvFrkj1atdazWOj799krAVSlVvqj7FcLWrFy5klq1apXY4Yq5iY+PJz4+3uwwRLoiJ3SlVCWVPn5LKdUmfZ/RRd2vELZm/PjxfPDBByV2uGJOZ8+epVy5cnz77bdmhyLSFdiHrpRaCHQFyiulzgJ/B1wBtNazgYHAWKVUCnATeDL9a4EQDiUoKMjsEGxKlSpVqFq1KmFhYYwdO9bscASFSOha6yEFrP8S+NJiEQlhg9asWUO5cuUICAgwOxSboZSiV69ezJs3j6SkJNzd3c0OqcSzz9PTQhSzCRMm8Oabb5odhs0JDg4mMTGRjRs3mh2KQBK6EAU6deoUR48epVevXmaHYnO6du2Ku7t7ZsEyYa6SPTtCiELISFaS0O/m6enJ/PnzadasmdmhCCShC1GgsLAwatWqRf369c0OxSYNHjzY7BBEOulyESIfqampbNu2jV69eslwxTxorfnxxx9Zt26d2aGUeNJCFyIfzs7OnD59moSEBLNDsVlKKSZNmkTdunUJDAw0O5wSTVroQhTAw8OD8uVl8nN+goODCQ8PJzEx0exQSjRJ6ELkY/DgwXzzzTdmh2HzgoODSUpKcrxiV3ZGEroQeTh9+jSLFy/m6tWrZodi8zp37oynp6cMXzSZJHQh8lCSLwZ9r9zd3enWrRsHDx40O5QSTU6KCpGHsLAwatasKdUVC2nBggV4e3ubHUaJJi10IXKRnJzMunXrZLjiPShdurT8rkwmCV2IXMTFxREYGEj//v3NDsWuTJkyhaeeesrsMEosSehC5KJcuXIsW7aMHj16mB2KXUlMTCQ0NFQuemESSehC5OLaNbku+v0IDg4mOTmZ9evXmx1KiSQJXYgczpw5w+OPP868efPMDsXudOzYES8vLxm+aBJJ6ELkkJGM2rZta3Ik9sfNzY3u3buzcuVK5MJlxU+GLQqRQ1hYGL6+vjRq1MjsUOzSU089RdWqVUlKSsLDw8PscEoUSehCZJGcnMzatWvp1q2bDMG7TwMHDmTgwIFmh1EiSZeLEFls3ryZ+Ph42rRpY3Yodi01NZXDhw+bHUaJIwldiCyaNWvGV199hb+/v9mh2LVp06bRvHlzYmNjzQ6lRJGELkQW5cuXZ+TIkdL3W0Tdu3cnJSWFtWvXmh1KiSIJXYh0586dY9asWTIG3QLat2+Pj48PK1euNDuUEkUSuhDpli1bxosvvsiVK1fMDsXuubi4EBQURFhYmAxfLEaS0IVIt3z5curWrSsXg7aQ4OBgzp8/z4EDB8wOpcSQhC4EkJCQwPr16+nTp48MV7SQvn37smbNGho2bGh2KCWGjEMXAli3bh23bt2iT58+ZofiMMqVK8cjjzxidhglirTQhQAOHDiAj48PnTp1MjsUh3Ly5EkmT55MTEyM2aGUCJLQhQDeeustoqKicHNzMzsUh3LhwgX+8Y9/8Ntvv5kdSokgCV2IdKVLlzY7BIfTtm1bHnzwQRm+WEwkoYsS76OPPqJv376kpKSYHYrDcXZ2pkePHoSFhZGWlmZ2OA5PEroo8UJDQ4mOjsbFRcYIWEOvXr24fPkye/fuNTsUhycJXZRoFy5cYNeuXTK6xYp69OhB6dKliYyMNDsUhydNElGirVixAjDGTAvrqFixonwDKibSQhcl2vLly6levTpNmjQxOxSHlpHMpQyAdRWY0JVSIUqpy0qpP/JYr5RSXyilIpVSB5RSUndU2I2OHTsyYcIEmR1qZefPn6dp06YsXLjQ7FAcWmFa6POBnvms7wXUS19GA7OKHpYQxeO1117j1VdfNTsMh1epUiUuX77M8uXLzQ7FoRWY0LXWG4H86on2B77Rhm1AGaVUZUsFKIS17N27l4SEBLPDKBGcnJzo27cvK1asIDk52exwHJYlzlJUAc5kuX82/bELOTdUSo3GaMXj6+vLhg0bLPDyxSs+Pt4u4y4KRzzm1NRUBg0aREBAAJMnT75rvSMec36K43jr1KlDbGwsn3/+uU1c4s8h/8Za6wIXoCbwRx7rVgAds9xfBwQUtM+AgABtj8LDw80Oodg54jFv2rRJA3rRokW5rnfEY85PcRzvzZs3tZeXl37hhRes/lqFYa9/Y2CXziOvWqKFfhaoluV+VeC8BfYrhNUsXboUNzc3evXqZXYoJYa7uzvvvPMOtWrVMjsUh2WJhP4LMF4ptQhoC9zQWt/V3SKErdBa8/PPPxMYGCj1W4rZxIkTzQ7BoRVm2OJCYCvQQCl1Vik1Uik1Rik1Jn2TlcAJIBKYC7xotWiFsIBDhw5x4sQJBgwYYHYoJdLp06fZtGmT2WE4pAJb6FrrIQWs18A4i0UkhJU1btyYvXv3UqNGDbNDKZHGjBlDREQEkZGRMv7fwmSmqChxlFK0aNGCsmXLmh1KifToo49y4sQJDh48aHYoDkcSuihRTp06xciRI6VQlIn69++PUoqlS5eaHYrDkYQuSpSff/6ZkJAQ+apvokqVKtG+fXtJ6FYgCV2UKD/99BNNmjShTp06ZodSog0YMICDBw9y4YIMiLMkSeiixDh//jybN29m0KBBZodS4o0cOZILFy5QubJUCbEkKVAsSozQ0FC01jzxxBNmh1LiyQlp65AWuigxlFL06NGDhg0bmh2KAPbv30/37t05efKk2aE4DEnoosR46aWXWLVqldlhiHQ+Pj6sW7eOJUuWmB2Kw5CELkqEixcvylXnbUzNmjVp06YNP/zwg9mhOAxJ6KJEeOyxxwgODjY7DJHD4MGD2bNnj8wLsBBJ6MLhRUVFsXXrVrp06WJ2KCKHjBFHixcvNjkSxyAJXTi80NBQABmuaIOqVavGCy+8IHV1LESGLQqH98MPP+Dv70/dunXNDkXkYvbs2WaH4DCkhS4cWmRkJDt27JCx5zYuNjZWinVZgLTQhUOrWbMmK1eupEWLFmaHIvIxcOBATp8+zdGjR6XOThFIC104NBcXF3r16iVTzG3c4MGDiYiIYOfOnWaHYtckoQuHtXv3biZNmkR0dLTZoYgCDBw4EHd3d7799luzQ7FrktCFwwoJCeHzzz/H1dXV7FBEAXx8fOjfvz8LFy4kOTnZ7HDsliR04ZCSk5NZtGgR/fv3lwtB24kRI0YQHR3Nxo0bzQ7FbslJUeGQVq1axbVr1xg+fLjZoYhCCgoK4uDBgzRp0sTsUOyWtNCFQ/r222+pUKECQUFBZociCsnFxUWSeRFJQhcOR2uNt7c3zz77rPSf25nExESGDRvG/PnzzQ7FLkmXi3A4SilCQkLMDkPcBw8PD/bt28eff/7JM888Y3Y4dkda6MKhaK2JiIgwOwxxn5RSPPvss2zbto1Dhw6ZHY7dkYQuHMquXbto0KCBXDTBjj399NO4ubkxd+5cs0OxO5LQhUOZO3cunp6e9OjRw+xQxH0qX748jz32GN988w1JSUlmh2NXpA9dOIz4+HgWLlzIE088IWPP7dz48eOpVq0aSUlJuLu7mx2O3ZCELhzG4sWLiY+P5/nnnzc7FFFEHTp0oEOHDmaHYXeky0U4jAULFtCwYUP+8pe/mB2KsIC0tDRWr17NyZMnzQ7FbkhCFw5j6dKlLFy4UMqvOogrV67Qp08fvvzyS7NDsRuS0IXDKF26tNQ9dyC+vr4MGDCAefPmkZCQYHY4dkESurB7V69epW3btmzZssXsUISFvfTSS8TExPDdd9+ZHYpdkIQu7N7cuXPZsWMHZcuWNTsUYWEdO3akZcuWfPHFF2itzQ7H5klCF3YtJSWFWbNmERgYiJ+fn9nhCAtTSjFhwgRiY2M5c+aM2eHYvEIldKVUT6XUMaVUpFLqb7ms76qUuqGU2pe+TLF8qELcbdmyZZw5c4aXXnrJ7FCElQwdOpQTJ05QvXp1s0OxeQWOQ1dKOQP/Bh4BzgI7lVK/aK0P59h0k9a6jxViFCJXWmtmzJhBrVq16NNH3nqOKqNiZnJyMjdv3sTHx8fkiGxXYSYWtQEitdYnAJRSi4D+QM6ELkSxGzt2LO7u7jg7O5sdirCi5ORkGjZsSHBwsAxjzEdhEnoVIGvn1VmgbS7btVdK7QfOA69pre8qlaaUGg2MBmNI0oYNG+45YLPFx8fbZdxFYcvHXK1aNQCLx2fLx2wN9nC8jRo1Yu7cuXTr1o0HH3ywyPuzh2O+Z1rrfBdgEPBVlvvDgX/l2KY04JV+Oxg4XtB+AwICtD0KDw83O4RiZ4vH/Mcff+iPPvpIx8XFWWX/tnjM1mQPx3vs2DGtlNKTJk2yyP7s4ZhzA+zSeeTVwpwUPQtUy3K/KkYrPOuHQqzWOj799krAVSlVvigfNELk5+OPP2bq1KncunXL7FBEMalfvz6DBg3i3//+N9evXzc7HJtUmIS+E6inlKqllHIDngR+ybqBUqqSSp9vrZRqk77faEsHKwTA6dOn+f777xk1ahTlypUzOxxRjCZNmkRsbCw//PCD2aHYpAL70LXWKUqp8cBqwBkI0VofUkqNSV8/GxgIjFVKpQA3gSfTvxoIYXHvv/8+Tk5OTJw40exQRDFr0aIFO3bsoFWrVmaHYpMKVT43vRtlZY7HZme5/SUgp56F1Z04cYJ58+YxZsyYzBOiomRp3bo1ALdv35aLgOcgM0WFXUlKSiIwMJBJkyaZHYow0XfffUe9evWIjY01OxSbIgld2BU/Pz9WrVrFQw89ZHYowkSNGjXi9OnTfPrpp2aHYlPkikXCbsybN4+goCCqVKlidihFd+0a7NkDJ07AmTNw+TLcvg0pKVCqFJQrB5UqQaNG4OcHlSubHbFNCQgI4PHHH2fGjBmMHz+e8uVlUB1IQhd2YufOnTz33HO8+eabvP/++2aHc+9iY2HtWlixAsLDIetVeJycoHx5cHMDFxeIjzcSflranW1q1ICuXaF7d+jbF2T6O++++y5Lly7lww8/5JNPPjE7HJsgCV3YPK01r776KhUrVuSNN94wO5zCS001kvi8ebB0KSQnQ+nSEBgIL7wAAQHQsKHREnfJ8a+YlgaXLsGRI3DwIGzaBMuXw9dfG4m/Vy94+mkjued8bgnRqFEjRowYwcyZM5k0aZIMYUUSurADS5cuZfPmzcyePZvSpUubHU7BEhLgv/+FTz+F06fhwQdh9GgYOBD+8hcozMgMJyejm6VyZejWDSZMMJL8jh3www+weDEsWwbVqsHYscYHhAWmw9ub999/n5deekmSeTo5KSpsWnJyMq+//jqNGzdm5MiRZoeTv7g4mDbN6B6ZMMFItosXw/nz8K9/QZcuhUvmeXFygnbt4LPPICoKfv4Z6teHN9+EWrVgyhSIibHY4diDhx56CH9/fwBu3rxpcjTmk4QubNrNmzdp3749M2bMwMVWuxaSk+HLL6FOHfj736FDB9iyxegmGTQIHnjA8q/p7Az9+xtdOgcOQFAQvPsu1KwJ778PJSy5vf3227Rt25aUlBSzQzGVJHRh03x8fPj222/p0aOH2aHkLiwMGjeGl14yfu7YYXSF/OUvxRdD06awZAns3290z7z1ljEyJjQUSsiE7YCAAA4ePFjiS+tKQhc2SWvNG2+8wYEDB8wOJXcxMfDMMxAcbHSjrFgB69dD+ixGUzRrZpx8Xb/eOPk6aJCR4I8eNS+mYtK/f3969erF5MmTOXXqlNnhmEYSurBJixcv5uOPP7bNetUrVkCTJvDdd0ZreO9eI7Eb9enM9/DDxhj3WbOMVnvz5kZ3THKy2ZFZjVKKWbNm4eTkxAsvvFBiLygtCV3YnEuXLjF+/HhatWrFuHHjzA7njpgYY6hgnz7GxJ/t241EaY0+8qJydoYxY4xhjwMGGCdM/f1h2zazI7OaGjVq8MEHH/C///2PiIgIs8MxhSR0YVO01owaNYq4uDi+/vpr27m03K+/Gn3kCxbA22/Drl3GOHJb5+sLixYZ8d+4YfTtT5wISUlmR2YVL774IkePHqVBgwZmh2IKSejCpixevJhff/2VDz/8ED8/P7PDMWZsDh8O/fpBhQrGSc9p04zJPfakTx84fNhotX/6qdFa37XL7KgszsnJiSpVqqC1ZunSpaSmppodUrGShC5sSv/+/Zk5cyYvv/yy2aHAL78YrfJFi4zhiDt3GonQXnl7w8yZsHq1UYqgXTt45x2jhoyDWbduHY899hgfffSR2aEUK0nowiYkJCQQExODu7s7Y8eOxcnJxLdmdDSN3n/fGOft62sk8nfesb9WeV6CgoxyAkOGwNSp0L49ng42MiQwMJDBgwczZcoUtm7danY4xUYSujCd1pqxY8fSunVrEhMTzQ3m55+hcWMqhIcbSXzHDmjRwtyYrKFsWfj2W/jxRzh9mlajRxtdMQ7SRaGUYvbs2VSvXp3HH3+c8+fPF/wkByAJXZhu5syZfPvtt4wYMQJPT09zgrh6FYYONUaEVK7M7tmzjW4WR2mV5+Wxx+CPP7jWpo1xsrRrV4iMNDsqiyhTpgw///wzN27cYNCgQSViKKMkdGGqFStW8PLLL9O3b1/eeustc4L46Sejrzw01DjhuWMHCXXrmhOLGXx9+ePdd2H+fKMrpnlzo5RB1vK9dqpZs2YsWLCAt956C2Ur8wSsyEaLY4iSYN++fQwePJgWLVrw/fffW7XfPDXVKG+SmJjl5/kYEv/xOcnrN6HrDkC/+zppNWuj18O+fQ+SmGgM537ggdwXLy9jQqatjKwsEqWMMfbdu8OoUUYpgx9/hJAQo/CXHXv00Uczbx8+fNg2Rk9ZiSR0YZqqVavSp08fPvvsM7y8vAr9PK2N0YQXLhiFDC9cMJboaOPxa9eMOUBZb+feNV8WmGrcjAReyLquWaHj8fIyrjdRurTx08fHmHfk63tnqVgx+22b7cmpUsWYCRsSAn/9q1FO4JNPjPK/dt7C3bBhA4GBgXz22We2MYrKCiShi2IXFRVFxYoVKV++PIsWLbpr/c2bRhnxkyeN5dQpYzl79k4Cz20Wu4eHca7vwQeNpXZtaNXKeKx0afD0BM/k63gs+QbP/f/Ds15VPF4bh2v9Wjg5GdVplTJ+7t27h4AAf1JT4datu5ekJOPCQjduGCMAb9y4s1y7BhERxvUp8jrHW6kSVK+e91K+vIn5UykYORIeecT4OWaM0Vr/6isjODvVsWNH+vXrx4QJE/D29qaWnX/zyI0kdFGsIiMj6datG+3adWDKlIUcOmTMdzl+3EjaJ0/CxYvZn+PmZpQYr1YNOnUyrvnw0EPZf1aubFyKM09aGwnptdeMcdcz3jNqlufRX3LrVixt2xb9eBMSjMR+6ZJx2dBLl4wPpDNnjJLmBw8aDeKc1W5Ll4a6daFevew/69Y1WvjFkuyrV4c1a2DOHOOEqZ+fMczx5ZeLVtfdJC4uLixcuJD+/fvz3HPPMWHCBLp27Wp2WBYlCV1YVVISHDtmJO21a3ezYEEwKSmp/PjjGyxZYmzj5GQk7Fq1jBpXtWoZZb1r1TKWSpWMbe7bvn0wfrxRo7xrV5g718iMxaBUKeObQu3aeW+jtdFdFBVlLKdOGQNNIiONyZyhodlHE3p7G+FnXD8642edOlbIs0oZV0Pq2dPoV3/tNeMyeLNnF2+JYAtxd3dn2bJlPPHEE/zzn/9k6NChtGnTxuywLEYSurCIxESjSuvhw8aS0fI+cSJjsMRq4HFcXMrTrdtq2rZtgJ+fkYjq1ze6SywuJsaouzJrltEH89//GiVvzZy0lAuljC6W8uVzn4ianGx0QR0/fifRR0TA5s3w/fd3tnN1NVryWZN8o0bQoAG4uxcxyBo1jJmzy5YZib1DB3j+efjgAyNwO+Lu7k5oaCjvvfeeQyVzkIQu7lF8vFHALyNxZywnT965loKLi5GkW7SAp56CunWTmDhxJL6+dVi1KoyHHnrIukHevm2c1HvrLaND+8UXjeGIZcta93WtxM3NSNT16t29Lj7+zgdpxt9l/35jJGbGqEMnJ+Mbgp+fMToz42fDhvfxQdq/v3GR62nTjIlIixfD5MlGN0yRPzWKj5ubG926dQNg7969TJ06lZCQEB608+uySkIXuYqNNRLEoUOwenUdPv7YSBanT9/Zxs3NaP21bm2MeMtIFnXrGq3F69ev4+npiZubO82bh1GrVq17Gs1yz7Q2+icmTzaasx07GtfydMSZnum8vIwTv61aZX88Kcn4FeT84F25EjKu0qaUkegbNwYvr1qcPWv8/Ro2NE4g5/uiH38Mzz4Lb7xhLP/+t9Faf/JJm/sGVJCIiAhWrlxJixYtWLBgAZ06dTI7pPuntTZlCQgI0PYoPDzc7BAsKiZG6y1btJ47V+tXXtE6KEjrqlW1NrKjsbi5pegWLbQeOlTr997T+qeftD52TOvbt/Peb1hYmK5SpYqeMmWK9Q8iLU3rsDCtW7UyAm7cWOtly4zH75Oj/Z0z3Lql9aFDWi9erPU772g9aJDx63J2Ts38eyulde3aWvftq/Xf/qb1N99ovXu31gkJeex0/Xqt/f2NJzdrZuw8NbVYj+t+ZP0b79y5U9etW1c7OTnpKVOm6Fu3bpkXWAGAXTqPvCot9BJAa2OExZEjd5aMFtuFC3e28/Q0+lwffpjM/m0/Pzh9ehOBgV0L9VqXL19m0qRJhISE4OfnR9++fa1zUGCcKfzxR/jwQ+OqQdWqwbx5Rrlbh5jtY3lubnf+rlmtXbuJhx7qku38x6FDsGrVnWKMShknqXN23TRq+zCeO3fCwoXGBT+eeMJ4I02eDIMHG31wNq5Vq1bs2bOHcePGMW3aNCpUqMD48ePNDuue2f5vWhRaWpoxSiJn4j5yxDg/mMHLy/haHRR05x/Tz88475Xbt+WzZwv3+qGhoYwcOZLExERef/11pk6dirs1+lVv3DAKS33xhdGvUL++ccJz2DAbnrFj21xcdOb7YODAO4/fvm2chM2Z6Fevzp7oa9Z0onHjp/DrO4TGgTvxW/UpjYaNptSbb8LYscYJVBs/eert7c0333zDiBEj6Ny5MwAbN26katWq1M5vmJINkYRuh65du3vEw5EjxvDArBNZKlQwGkoZDaaMpWpVy41jTkpKIj4+nvLly1OnTh06derEJ598QsOGDS3zAhm0Nlrhs2YZQzsSE43O+yVLjIJa0iK3ClfXO++bxx+/8/jt2/Dnn9mT/OHDsGaNE8nJbYEfUEpT8/IF/CbtofFb82nUviz1hrSizoBm+FZSNjvxtHv37oDRHT1mzBgiIyMZNWoUb7zxBtVtfGKVJHQbpLVR/O/PP7Mn7ozbWVvbShnzPxo2hC5d7gxVa9TImH5uLefOnSMkJIRZs2bRo0cP5s2bR8uWLVm+fLllXygy0vgqv3Ch8anl4WFURRwz5u4zgaLYuLoa77mGDbMn+pSUrIlecejQQxzeU5bfInuQvNkVNgPjoJTrLerUhjp+D1C3rjGGPuNntWq28fmslGLt2rW8++67zJkzh//85z8MHDiQN998k2bNCl8aojhJQjdBaqoxhf306dyXqKjsLW0nJyNp16tnDCLIOnOwdu3ivUZxWFgYs2fPZvny5aSlpREUFMTTTz9tuRdITTVm04SFwfLlsHu38XinTsbVdoYMgTJlLPd6wqJcXIyRTw0aGJV5DR6kpMDJg/H8uWAbkcuP8uexFCKP1eHoST9WplbnVuqdGVGursbEsho1jPd9xs+MpWrV4hsh+dBDDzFr1iwmTZrEv/71L+bOncuQIUNo1qwZly9fJikpyaZa7ZLQLUhruH79Tr2RrIWjMm6fPWssGUPHMpQvb7xx/fygVy/jjZsx1btWLXMuLK+15uTJk6xatYpOnTrh7OzM6tWr2b59O6+//jqjRo0qet9iair88Ycxi3PzZmOqeXS08dWjbVuYPt04sVatmmUOSpjCxQXqtfSiXsvu8El3OHfO6C779XPSft/EOXyJLB3An3WCiCzbmhNOdYmKLUNYmMp24j6Dr++dBF+lijGbuHJl42fG7fLlLdfSr169OtOnT2fKlCmZNftnzpzJ1KlTadWqFb169SI4OJjWrVubemFzpU0q+t6qVSu9y8YvUpuWZpx/u3oVrlwxfm7efJTy5Rty9SrZlosXjaR969bd+/H2vlNvpEoVI3FnXapXL6AOSTE6e/YsS5cuZdeuXWzYsIGoqCgANm/eTIcOHYiNjaVUqVL396a9dcuYBXPwoLHs2wfbthmD3sH4T3zkEeMTLSjIun1GBdiwYYPD1fnIj6nHe/26cZZ11Sr4/XdjlhoY/zgtW3KraSvO1fgLUQ+2IErVIOq8S2aZhNOnjYbSjRt379bZ2ah7k5HgfX2Nt1TGcuHCH3Tp0oTy5e88di+lE06ePMmiRYtYsWIFW7duJS0tjcqVK3P27FmcnJw4cuQIPj4+VK5c2aK12JVSu7XWufY3FiqhK6V6Av8EnIGvtNYf5liv0tcHA4nAM1rrPfnt01oJXWtjUkVi4t1LXNyding5K+RlXTLWxcTkfUWuBx64M127XLk7CTu3wlHWnEtzr5KTkzl+/DinTp3izz//zFxefvllgoKCCA8Pp1u3blSoUIHOnTvz8MMPU6pUKZ5++umC35RJSXcqUGWUS8youHXypNEfnvHVJGP8XLt2xjTyDh2M79k2cqZMErqJzpyBjRuNb2x798KBA3eql2UU/sn4+lq3LlSrRuKDVbnkUoULab5cjHmACxeMRlbGN+SLF423ZXS08TbNi7e38f/84INGgbSsS0aJ5Nzu3759jZ071xATc4GJE/+Kuzu0ahXAnj17KFu2LH5+fjRu3JhOnToxbNiwIv16ipTQlVLOQATwCHAW2AkM0VofzrJNMPASRkJvC/xTa51vrbr7TegbN8I//pF7wk5IMP7uhf3S4eSUvYZ1zprWZctqypXTmYm7YkUnIiO30atXazw9NaDR6fOrlVK4po+3vXXrFjotzRjsn77eyckJ9/R+k4SEBNLSPyl0+vwdFxcXPNPnYd+4cYPbyckkJyeTfPs2t5OT8SpVisqVK4PWbNm61ViXviQkJFC3dm1a+ftzMzGRt957j7i4OOLi44m5fp1Lly8z+umnGTtyJCdPnaJ2y5aZv4NSnp7UqVqVv48axWOdO5MUF0f01as85OODun0bbt7k6M6dNKxSxfhEzLlcuWIk8cuXjfs5lSlzp8pWw4bQtKmx1K9v0xX7bCrBFQObPt7UVGMo1969xlCujNEBx48brfucypQxmuZlyxq3M36mL4nO3kSnlmHn8Wv4VGnG1VveRCd5Ep3oSXSCO9FxrlyPc+FGgjOx8c7EJjgRG+fEjThFSkrhGxwPPLARF5f9KHWY1NRDJCcfpkyZrvj7h/Lkk/Dcc/f368gvoRemD70NEKm1PpG+s0VAf+Bwlm36A9+kz2LappQqo5SqrLXOpferaFI2bGL/6llcYSmQkbk1Lrgw3uV9PJ0SCUtbzD79R7b1pXFjuUsbfNQN/i8lktU6EZ0GsdfhxnWoBexP70bonprK9zletylwAGiF8Ym1I8f6vwBb0m/7k/2XA9ADWJV+uxFwJsf6gUB68UFqAtdzrH8O+G/67S5Azi8OL6XH5gT8B/BOX3yAKsCDr70Gr73GQ8BCoDpQB6iYmIiKiID/+z8A3NO3zyrbAER3d6MZk7GULw9t2hj/QFmX6tWNFrecwBRF5ex8Z+hWTteuGf3xWU9cnTtn9INev26sP3HCuB0TAykpeAKewL2eldFAEu7EUppYSnMDn8zbiXjmWEqReMuTm8keJNKORLqRiCfx0a7Erd3OTefL8JzlJ90VpoU+EOiptX4+/f5woK3WenyWbZYDH2qtN6ffXwe8obXelWNfo4HR6XcbAMcsdSDFqDxw1ewgipkcs+MraccL9nvMNbTWFXJbUZgWem7fMXJ+ChRmG7TWc4A5hXhNm6WU2pXX1x1HJcfs+Era8YJjHnNhyqKdJfu3k6rA+fvYRgghhBUVJqHvBOoppWoppdyAJ4FfcmzzCzBCGdoBN6zRfy6EECJvBXa5aK1TlFLjMS454wyEaK0PKaXGpK+fDazEGOESiTFs8VnrhWw6u+4yuk9yzI6vpB0vOOAxmzaxSAghhGXZ16VFhBBC5EkSuhBCOAhJ6EWglHpNKaWVUrZdub+IlFLTlVJHlVIHlFJLlVJlzI7JWpRSPZVSx5RSkUqpv5kdj7UppaoppcKVUkeUUoeUUhPMjqm4KKWclVJ70+fROARJ6PdJKVUNoxxClNmxFIPfgCZa62YYZSAmmRyPVaSXufg30AvwA4Yopfzyf5bdSwEmaq0bAe2AcSXgmDNMAI6YHYQlSUK/f58Br5PLBCpHo7Veo7XOKPi7DWOegSPKLHOhtU4GMspcOCyt9YWMQnpa6ziMBJez+oPDUUpVBXoDX5kdiyVJQr8PSql+wDmt9X6zYzHBc0CY2UFYSRWyl9k5SwlIbhmUUjWBlsB2k0MpDp9jNMjSTI7DouQCF3lQSq0FKuWyajLwJhBUvBFZV37Hq7Velr7NZIyv6AuKM7ZiVKgSFo5IKeUF/Ai8orWONTsea1JK9QEua613K6W6mhyORUlCz4PWuntujyulmpJenDG9PnhVYI9Sqo3W+mIxhmhReR1vBqXU00AfIFA77uSFElnCQinlipHMF2itfzI7nmLQAeiXXvbbHSitlPpOa120QuU2QCYWFZFS6hTQSmttj1XbCiX9AiefAl201lfMjsdalFIuGCd9A4FzGGUvhmqtD5kamBWlX5zma+Ca1voVk8Mpdukt9Ne01n1MDsUipA9dFMaXGOXVf1NK7VNKzTY7IGtIP/GbUebiCLDYkZN5ug7AcKBb+t92X3rLVdghaaELIYSDkBa6EEI4CEnoQgjhICShCyGEg5CELoQQDkISuhBCOAhJ6EII4SAkoQshhIP4fxO969RJkeijAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from torch.distributions import normal\n",
    "\n",
    "def log_density_ratio(p, q):\n",
    "    def log_ratio(point):\n",
    "        return p.log_prob(point) - q.log_prob(point)\n",
    "\n",
    "    return log_ratio\n",
    "\n",
    "def density_ratio(p, q):\n",
    "    log_ratio = log_density_ratio(p, q)\n",
    "\n",
    "    def ratio(point):\n",
    "        return torch.exp(log_ratio(point))\n",
    "\n",
    "    return ratio\n",
    "\n",
    "def get_discrete_distribution_values(p, xp):\n",
    "    values = torch.exp(p.log_prob(xp))\n",
    "    values = torch.div(values, torch.sum(values))\n",
    "\n",
    "    return values\n",
    "\n",
    "# define two probability density distributions\n",
    "mu1, std1 = torch.tensor(1.0), torch.tensor(1.0)\n",
    "mu2, std2 = torch.tensor(0.0), torch.tensor(2.0)\n",
    "p_real = normal.Normal(mu1, std1)\n",
    "q_real = normal.Normal(mu2, std2)\n",
    "\n",
    "# calculate their ratios for given space X\n",
    "x = torch.linspace(-5, 5, 10000)\n",
    "ratio = density_ratio(p_real, q_real)(x)\n",
    "\n",
    "# and plot results\n",
    "plt.figure()\n",
    "plt.plot(x, torch.exp(p_real.log_prob(x)), c='r', label='real p(x)')\n",
    "plt.plot(x, torch.exp(q_real.log_prob(x)), c='b', label='real q(x)')\n",
    "plt.plot(x, ratio, c='k', linestyle='--', label='real p(q)/q(x)')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.ylim(0, 2.5)\n",
    "plt.xlim(-5.5, 5.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "Docelowa wartość Dywergencji KL (miara podobieństwa rozkładów) również może zostać obliczona na wiele różnych sposobów:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pytorch kl-divergence (0.443147)\n",
      "gaussians' based analytical kl-divergence (0.443147)\n",
      "plot divergence based on integral calculation: (0.430858)\n",
      "monte carlo kl divergence estimation: (0.442540)\n"
     ]
    }
   ],
   "source": [
    "from torch.distributions import kl_divergence\n",
    "\n",
    "mc_no_samples = 50000\n",
    "xp_mc_samples = p_real.sample((mc_no_samples,))\n",
    "\n",
    "p_values = get_discrete_distribution_values(p_real, x)\n",
    "q_values = get_discrete_distribution_values(q_real, x)\n",
    "plot_div = torch.dot(p_values, torch.log(torch.div(p_values, q_values)))\n",
    "real_div = kl_divergence(p_real, q_real)\n",
    "gauss_div = torch.log(std2) - torch.log(std1) - (1 / 2 * (1 - (((std1**2)+((mu1-mu2)**2)) / std2 ** 2)))\n",
    "monte_carlo = torch.mean(torch.log(torch.exp(p_real.log_prob(xp_mc_samples))/torch.exp(q_real.log_prob(xp_mc_samples))))\n",
    "\n",
    "print(f\"pytorch kl-divergence ({real_div:.6f})\")\n",
    "print(f\"gaussians' based analytical kl-divergence ({gauss_div:.6f})\")\n",
    "print(f\"plot divergence based on integral calculation: ({plot_div:.6f})\")\n",
    "print(f\"monte carlo kl divergence estimation: ({monte_carlo:.6f})\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Wyliczanie KL i porównywanie dwóch dystrybucji komplikuje się w momencie, gdy nie są one jawnie podane.\n",
    "Współczynnik KL jest dany równaniem (przybliżenie metodą monte carlo):\n",
    "\n",
    "<img src=\n",
    "\"https://render.githubusercontent.com/render/math?math=%5Ctextstyle+%5Cbegin%7Balign%2A%7D%0ADL_%7Bkl%7D%5Cbigg%5Bp%28x%29%7C%7Cq%28x%29%5Cbigg%5D%3D%5Cmathbb%7BE%7D_%7Bp%28x%29%7D%5B%5Clog%7Br%5E%7B%2A%7D%28x%29%7D%5D%3D%5Cfrac%7B1%7D%7BM%7D%5Csum_%7Bi%3D1%7D%5E%7BM%7D%5Clog%7Br%5E%7B%2A%7D%28x_%7Bp%7D%5E%7Bi%7D%7D%29%0A%5Cend%7Balign%2A%7D%0A\"\n",
    "alt=\"\\begin{align*}\n",
    "DL_{kl}\\bigg[p(x)||q(x)\\bigg]=\\mathbb{E}_{p(x)}[\\log{r^{*}(x)}]=\\frac{1}{M}\\sum_{i=1}^{M}\\log{r^{*}(x_{p}^{i}})\n",
    "\\end{align*}\n",
    "\">\n",
    "\n",
    "Aby wyliczyć wartość KL metodą Monte Carlo, jednocześnie nie znając analitycznej formy rozkładów p(x) i q(x) estymuje się wartość <img src=\"https://render.githubusercontent.com/render/math?math=r^{*}=\\frac{p(x)}{q(x)}\">\n",
    "z użyciem optymalnego klasyfikatora probabilistycznego. W modelu tym wyjściem jest prawdopodobieństw przynależności do klasy p(x).\n",
    "\n",
    "Ale czy że rzeczywiście istnieje związek pomiędzy wyjściem klasyfikatora a r*(x)?\n",
    "\n",
    "Mając zbalansowany dataset <img src=\"https://render.githubusercontent.com/render/math?math=X={\\{x_{p_1}, ..., x_{p_i}, x_{q_1}, ..., x_{q_i}\\}}\"> ze\n",
    "zbiorem etykiet <img src=\"https://render.githubusercontent.com/render/math?math=Y={\\{1, ..., 1, 0, ..., 0\\}}\"> dowodzimy, że:\n",
    "\n",
    "<img src=\n",
    "\"https://render.githubusercontent.com/render/math?math=%5Ctextstyle+%5Cbegin%7Balign%2A%7D%0Ar%5E%7B%2A%7D%3D%5Cfrac%7Bp%28x%29%7D%7Bq%28x%29%7D%3D%5Cfrac%7BP%28x%7Cy%3D1%29%7D%7BP%28x%7Cy%3D0%29%7D%3D%0A%5Cbig%28%5Cfrac%7BP%28y%3D1%7Cx%29p%28x%29%7D%7Bp%28y%3D1%29%7D%5Cbig%29%5Cbig%28%5Cfrac%7BP%28y%3D0%7Cx%29p%28x%29%7D%7Bp%28y%3D0%29%7D%5Cbig%29%5E%7B-1%7D%3D%0A%5Cfrac%7BP%28y%3D0%29%7D%7BP%28y%3D1%29%7D%5Cbig%28%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7BP%28y%3D0%7Cx%29%7D%5Cbig%29%3D%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7BP%28y%3D0%7Cx%29%7D%0A%5Cend%7Balign%2A%7D%0A\"\n",
    "alt=\"\\begin{align*}\n",
    "r^{*}=\\frac{p(x)}{q(x)}=\\frac{P(x|y=1)}{P(x|y=0)}=\n",
    "\\big(\\frac{P(y=1|x)p(x)}{p(y=1)}\\big)\\big(\\frac{P(y=0|x)p(x)}{p(y=0)}\\big)^{-1}=\n",
    "\\frac{P(y=0)}{P(y=1)}\\big(\\frac{P(y=1|x)}{P(y=0|x)}\\big)=\\frac{P(y=1|x)}{P(y=0|x)}\n",
    "\\end{align*}\n",
    "\">\n",
    "\n",
    "Użyto twierdzenia Bayesa <img src=\"https://render.githubusercontent.com/render/math?math=P(y|x)=\\frac{P(x|y)P(y)}{P(x)}\"> w trzecim rozwinięciu i założono, że zbiór danych jest zbalansowany, tj. jest tyle samo sampli z\n",
    "p(x) co z q(x).\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe1UlEQVR4nO3dfXBd9X3n8fdXj7YlP2JbNrGxTdaYAGsnyEBYSiM1acYQGrczZJZk62nIUJc2TprtMoGWye4faafJZjdtQpJ6GMJmMk2q2VISCPFCu41VJu0A5inGjjFx7GCEH+UHGcnYsqTv/vG9h3stS/de2deS7m8/rxmN7rnnd875nnPP+dzf/enY19wdERGpfjUTXYCIiFSGAl1EJBEKdBGRRCjQRUQSoUAXEUlE3URteO7cub506VL6+vpoamqaqDLOW7XWDdVbe7XWDdVbu+oef6Vqf+GFF7rdfd6IM919Qn5aW1vd3X3z5s1ejaq1bvfqrb1a63av3tpV9/grVTvwvI+SqxpyERFJhAJdRCQRCnQRkUQo0EVEEqFAFxFJRMnbFs3sYeA24JC7XzPCfAO+BtwKnAQ+6e4vVrpQGV87d8LAAFx9dTz+xS/g9ddhyRK47TZ44on89PLlsGJFPNfVBTffDHv3Qn8/1NbCc8/B0aNw8CCcOgVTpsC8ebB1a0w3NMDgINTVweLF4A4f+hA8+mhs94orYOlSuPZa+NSn4MyZWP6KK2LdBw/GMgsWwCc+AU8+Cdu3x7YXLoTVq+Hll2H/fpg6NZ6rqYnlMjNnwoEDse53vztq7+mBkyej7bx5sY3jx+O3OzQ3w5w50NsLZvG7uTmOWXd3LN/TE/Puvx/+5E/g8OFY14IFsc3jx2P+zJkwa1Y8njMH3v/+OL7d3VH3e98Ld98dtW7cGMf5ppvgsstie9/6VtR5881xHFesiGOQPR7pdc0Ua1dopHYyuZRzH/p3gG8A3x1l/i3A8tzPDcDf5H5LFRsYgN2789PPPAOvvRYhOjgIW7bkp5cti4v99dcjfDo7YWgonjtyJELqtdfi8dBQLN/fH6H49tsRRPX1Me9nP4P58+GNN+JN4fhxOHYsalmxAn71q1ju4EF49dVYz+Bg1LhnDzzwQIRwFkZ798KhQxG2Q0Nw4kSEpFksNzQU2z9wIL/M1q0RXKdPx3RtbdQwNBTT2X9Q2t0ddQ4Oxvrq6qKuri6YPj3exE6fjn3r6YFt26LtoUNxrOrqYptTpkT7zNSp8UY1a1bs/8mTsc3t22P+Sy/F8kNDcOON8JOfxLzs2La1xfTu3XD55aO/rldfXbpdfX1Mj9ZOJpeSge7uT5vZ0iJN1gLfzd0f+YyZzTKzhe6+v1JFyvjLenDZxb9wYfzu7obvfx9mz44wX7gw36atLX53dkZoHT0Kb70FM2ZEj9QsAvXIkQi6mpoIS/d4nAXmiRPRU16yJHqq27dHz3ZoKNYBEVwDA7H89OmxvjNnYtlMFtq9vTFdXx/PZcGXtYF8sGfzCrdVU3Nue/dok9VeVwdNTfk3q5Mno7Zsm4OD8VNTE9Nvvw2NjdGjnzs3lj15Mp7v6Yn9WLgQrrkmXouaGujoiG2vXBnPHzoEP/5xvNksX55/jbLX4/LLz+6JD39dy2m3axf86Eejt5PJxbyM/w89F+hPjDLk8gTwJXf/aW76n4F73f35EdquB9YDtLS0tHZ0dNDb20tzc/OF7cUEqNa6YWy19/ScPd3Xlx82Gf6P2WbOzC/T1xePs7bZ72PH4vn+/gjBoaF8T7Xw99AQXHJJbOPIkVhfS0sv+/c3vxOkEKFZl+uWnD6dfz57Lmub9aCHhs5+LlM4na3D7Ow3muGGtyvcZvYGkU0vWNDLvn3N7zwHMb+hId8LhnhTgjgOU6fGMbv00rOPaVNTHOt9+/LH+JJLzq0vez1GUvi6Fm/XCzSXbDfZpHx9tre3v+Duq0eaV4l/+m8jPDfiu4S7Pwg8CLB69Wpva2ujs7OTtqxrV0WqtW4ov/bt26OHDTF0sX9/9NCPHYse+ty50Su87LJoM29efrkDB2Ic9sSJ6KGfOBHDEcN76BDB2NgYIeYeITZ1aqx34cJ8D/3++zu55562d8LXPQJ96tR8D72wf5IFbWH4Zz30wuGTwmDPni98rq4ueteF2x0pzAt76I2NUc+ZM9HmL/6ik89/vu2dTyU1NaP30E+fjmWG99APHIhtLliQH0o6fPjcHnrh6zFSj7rwdS3Vrqenk5qatqLtJqP/H67PkVQi0LuAxQXTi4B9FVivTKDhY6bPPpsfM7/llvwYOpw91LJ/fwTL/PkRNlnwHTiQH0OfMiUftNkYemEvOBuief31GFNvbIwwyZZxjwCtq4vwLAynGTPyY+hZ4E+ZEsMuWbjW1sZPYUjX1OTH0LNefzaGng2/ZNvO9ikL52z4p68v/+Y0bVoMOZ05Ez3wwm3W1MSbUF1d1Hro0Nlj6DNnxv7OmhWvwbZtEaSf/WzM//rX49hcey185CPxd4Xt22M9q1blX7PCsfKRXtfCMfTR2jU2wpo1o7eTyaUSgf44sMHMOog/hvZo/Lz61dXlL/qdO2Mse+HC/F0u2R0kS5bk73zYsyeez+5yueyy/F0u8+df+F0uTU3xuxrvcpk5M3rbF3KXSxak73tfrCO7y+W3fmvku1yy13G01xXyv0drd/hw8XYyuZRz2+LfAW3AXDPrAv4bUA/g7huBTcQti7uI2xbvvFjFyvgpvDVtxYpzb1W77bZzlyl8rrAXN1Lbcvz+75893dl59p03o/nzPz+/7V1MnZ3wYoVu5s1uXyz0R3907nMj9aRHuuWwWLvOzuLtZHIp5y6Xj5eY78CnK1aRiIicF/1LURGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUlEWYFuZmvMbKeZ7TKz+0aYP9PMfmRmPzOz7WZ2Z+VLFRGRYkoGupnVAt8EbgGuAj5uZlcNa/Zp4OfuvgpoA/6nmTVUuFYRESminB769cAud9/t7v1AB7B2WBsHppuZAc3AUWCgopWKiEhR5u7FG5jdDqxx97ty0+uAG9x9Q0Gb6cDjwJXAdOA/uvuPR1jXemA9QEtLS2tHRwe9vb00NzdXan/GTbXWDdVbe7XWDdVbu+oef6Vqb29vf8HdV484092L/gAfAx4qmF4HPDCsze3AXwEG/DtgDzCj2HpbW1vd3X3z5s1ejaq1bvfqrb1a63av3tpV9/grVTvwvI+Sq+UMuXQBiwumFwH7hrW5E3g0t71duUC/sox1i4hIhZQT6FuA5Wa2LPeHzjuI4ZVCe4EPAphZC7AC2F3JQkVEpLi6Ug3cfcDMNgBPAbXAw+6+3czuzs3fCHwR+I6ZvUIMu9zr7t0XsW4RERmmZKADuPsmYNOw5zYWPN4HfLiypYmIyFjoX4qKiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikggFuohIIhToIiKJUKCLiCRCgS4ikoiyAt3M1pjZTjPbZWb3jdKmzcxeNrPtZvYvlS1TRERKqSvVwMxqgW8Cvwl0AVvM7HF3/3lBm1nAt4A17r7XzOZfpHpFRGQU5fTQrwd2uftud+8HOoC1w9p8AnjU3fcCuPuhypYpIiKlmLsXb2B2O9Hzvis3vQ64wd03FLT5a6AeuBqYDnzN3b87wrrWA+sBWlpaWjs6Oujt7aW5ublCuzN+qrVuqN7aq7VuqN7aVff4K1V7e3v7C+6+esSZ7l70B/gY8FDB9DrggWFtvgE8AzQBc4FfAFcUW29ra6u7u2/evNmrUbXW7V69tVdr3e7VW7vqHn+lagee91FyteQYOjFuvrhgehGwb4Q23e7eB/SZ2dPAKuC1MtYvIiIVUM4Y+hZguZktM7MG4A7g8WFtHgNuNrM6M5sG3ADsqGypIiJSTMkeursPmNkG4CmgFnjY3beb2d25+RvdfYeZPQlsBYaIIZptF7NwERE5WzlDLrj7JmDTsOc2Dpv+CvCVypUmIiJjoX8pKiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIolQoIuIJEKBLiKSCAW6iEgiFOgiIokoK9DNbI2Z7TSzXWZ2X5F215nZoJndXrkSRUSkHCUD3cxqgW8CtwBXAR83s6tGafdl4KlKFykiIqWV00O/Htjl7rvdvR/oANaO0O4zwD8AhypYn4iIlMncvXiDGD5Z4+535abXATe4+4aCNu8Cvg/8BvBt4Al3f2SEda0H1gO0tLS0dnR00NvbS3Nzc6X2Z9xUa91QvbVXa91QvbWr7vFXqvb29vYX3H31iDPdvegP8DHgoYLpdcADw9r8PfD+3OPvALeXWm9ra6u7u2/evNmrUbXW7V69tVdr3e7VW7vqHn+lagee91Fyta6MN4wuYHHB9CJg37A2q4EOMwOYC9xqZgPu/sMy1i8iIhVQTqBvAZab2TLgTeAO4BOFDdx9WfbYzL5DDLn8sHJliohIKSUD3d0HzGwDcfdKLfCwu283s7tz8zde5BpFRKQM5fTQcfdNwKZhz40Y5O7+yQsvS0RExkr/UlREJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBEKdBGRRJQV6Ga2xsx2mtkuM7tvhPn/ycy25n7+zcxWVb5UEREppmSgm1kt8E3gFuAq4ONmdtWwZnuAD7j7SuCLwIOVLlRERIorp4d+PbDL3Xe7ez/QAawtbODu/+bux3KTzwCLKlumiIiUYu5evIHZ7cAad78rN70OuMHdN4zS/h7gyqz9sHnrgfUALS0trR0dHfT29tLc3HyBuzH+qrVuqN7aq7VuqN7aVff4K1V7e3v7C+6+esSZ7l70B/gY8FDB9DrggVHatgM7gEtKrbe1tdXd3Tdv3uzVqFrrdq/e2qu1bvfqrV11j79StQPP+yi5WlfGG0YXsLhgehGwb3gjM1sJPATc4u5HyliviIhUUDlj6FuA5Wa2zMwagDuAxwsbmNllwKPAOnd/rfJliohIKSV76O4+YGYbgKeAWuBhd99uZnfn5m8E/itwCfAtMwMY8NHGeERE5KIoZ8gFd98EbBr23MaCx3cB5/wRVERExo/+paiISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkQoEuIpIIBbqISCIU6CIiiVCgi4gkoq6cRma2BvgaUAs85O5fGjbfcvNvBU4Cn3T3FytcK3R1wZYtcPgwzJsH110HixaVv0xNDbjHz2jLl7ONri44fhwefLD8Ooqte6Qad++Gl16CN96AI0dgaAgaG+OnqQlmzYLLLovlh4Zgzhy45BJ4+WX46U/h5ElYuBCWLoUdO+C11+I5M/jyl+Fzn4OWlpju6YE9e+Ctt2L78+bBe98LV14Jx47BwYOxLERtJ05EXb29cPp0rCMzbRo0NMCpUzGvvh4uvxxWroT9+2Nd3d2xr2fOxPrq6+NncBAGBmI9tbUwf37sn3sse889cO+98J73wHPPwa5dsQ6IGhoa4tgsWADNzdDfH/WfPBnra2zM12WWP2ZmcYwPHIh52TkCUdfMmXEsIY7R4CBMnQpLlsSxco/tzZ4NO3fCK6/ENqdMiZoaG+Huu+HrX4/jM3161PPKK/Dmm3HMIY5rf39+/xsaYv2rVsU+zZ4NK1ace75l58+zz8b2zaL9rFlw9Cj09UV9S5bE8tk5d/hwtD12LP/6z5gR7a64IrZTjrFcl+dzDafkIu9/yUA3s1rgm8BvAl3AFjN73N1/XtDsFmB57ucG4G9yvyunqwseeyxO0paWCJPHHoO1a4ufPNkytbXw9NNxAt98c1xww5cvZxtZm4ULy6+j2Lqvvz7CqbDG7u4IzOPH44V3j9AeGoqLv6kp5u3eHaG0alWsf88eePXVCMKpU2HrVvjXf4W6ugiKLKQGB2PenDmxzaNH8+vOau3pgWeeiZOusTHq6O2NYOjqiuM4OBjLFXr77fzj2tqYv307bNsW6+rvj/0rdOZMPpgzg4Owb1/81NVFOLrH/j333LnH1z3CuL8/9mfatHx92Ztkf388Novajh2LEB8cjKA3yx+jwjpOnYrAr62NNnV1Ebb79sVrsWRJHPMf/jDCc+rUWG7//qh79uxYx9NPww03xLH41a+iRvd4LU+fjvVDvKmZxfy334bXX4cbb4TFi+O5ffvy51t2XnV3xxv5lCmxn7t3x/6tWhWv28BAvGGcPg0dHXDTTVH7k0/C3r2wbBkcOhR1nTgR+7BvXwR7MWO5Ls/nGk7JOOx/OUMu1wO73H23u/cDHcDaYW3WAt/18Awwy8wWVqTCzJYtcSBmzIiLcsaMmN6ypbxlfvlLmDs3emS//OXIy5ezjaxNbW35dRRb9yOPnFvjoUPRWzpzJi7smpp8ENXXRzj098fjt96Ki/6tt+Iibm6OdlmAm8VFXFMTQZT1ps1imb6+aJdtI2vX2xvb7+uLgKyvjwDYty+/78PDvKbgdMq2U1sbYeIeb6LHjxc/ToWy9We99uGPR5Lty+nT+eXr6vLLZYFdeCyHv5kMZxbLDw3l3wyyN9jsU8WOHfkgrq2N9dbXx/LHjsXzzc3xqauvL9pkn3qy9UKsL3sdstewpibOjebmeAMqPN+y82rr1uj5X3JJbLuvL57fuTPOrTlz4vXeujXaHDgQ6zx9Oua9/nq0mzs3Aj3bTl9f8WMzluvyfK7hlIzD/psP75EMb2B2O7DG3e/KTa8DbnD3DQVtngC+5O4/zU3/M3Cvuz8/bF3rgfUALS0trR0dHfT29tLc3Fy60u7u/AVS6MyZOAlLLdPTExc2xIU3c+a5y5ezjVybXqB5tDZjqf/o0bigCms8duzswMh6jtnF7R6BUFsbj6dOjeV7e2Mb2Wt65szZoZtbT++iRTR3deWfg3z4Z88NDUUtWUibxc+pU2cPXZWSBX+2/PA3gXLV1tJ76aVRdznbzWou3J+R5pW7rsJlsv3PHjc05IePhobyb2IFPf7eRYtozoaZCt8QCt/wRlJTk++5z52bP3ez8y07rw4ezJ/fp09Hu/r6eFPIrq9sm01NcV5B/lPLqVPRu8/aNTTAzJn0Dg3RPHv26MdkLNfl+VzD56nsXBlPZe5/qdrb29tfcPfVI80rZwzdRnhu+NlXThvc/UHgQYDVq1d7W1sbnZ2dtLW1la7iBz+IwJoxI//ciRPxEXS05QuXefXVONEhhhBuvPHc5cvZRq5NZ3MzbYVjn8XqKLbuV1+N8eXCGrNhl1OnIuQhH+4NDfF41qx8kLe2xnLPPZcfR4b4yN/XFxd3NlwwOEjnV75C2+c/HwGQ9VKzTwIQbQcHo6bsjS/rMe7fnw+27HhmCnvt2fqmTIlPEDU1ESzZx/9yZTU1N9P5xS/S9oUv5Hu2xZYp3J+Ghqghe2PJPoVkAQdnv3EOlz1fX58fbsmGnKZNi79l9PTEm3N/fxz/Eyfyxx3o/Mu/pO2rX83PP3o0PxZ/9Gi0ranJ15HV2dQUv+fOhdtui3P36qvz51t2Xj31VByXGTNi6O3YsWhz6lT8zQFiX+vr4+P+rFnx3LZtESjHj8fzWbtrroGrr6azqan49TmW6/J8ruHzVHaujKcy9/9Cai9nyKULWFwwvQjYdx5tLsx118VJd+JEnHAnTsR0sT/cFC7z7ndHSB45Eo9HWr6cbWRtsvHZcuootu7bbz+3xvnz4+Nz1tvOgsc9Lr7sj21nzkS7qVPj9+WXxwnjng8d9wiBoaGze4LusUxTU74Hm/0MDET4ZsMsc+bkh18uvfTsselChb3gwvH6bKhn2rR8kJQjW39dQb+jrkQfJNuXxsb88gMD+eUKe+bZsRyp11TIPf9pxf3sYZHa2pj3nvfE72z+lCn5oZzZs+P53l543/vyPeRp0/L7mR2vwt579hoODcW50dsbf+wsPN+y82rlyhhSOXIktp39nWXFiji3jh6N13vlymizYEGss7Ex5i1ZEu26uyNwsu1knYPRjOW6PJ9rOCXjsP/l9NC3AMvNbBnwJnAH8IlhbR4HNphZB/HH0B5331+xKiH+aLB2bYw3HTwYf2D7wAeK/zGhcJm+vvhjaHbBTZt27vLlbCNr8/zz5ddRat0LF55b4/nc5dLefvZdLitXjnyXS21tzBvLXS7ZsJB7XOwXcpdL9kfWsd7lYhY1Veoul8WLK3uXy513nn2Xy7Jl+btcpkyBX//1fI/+134tf5fL9Omx7nLuclm8+Ow7IwrPq7ffzt/lcs01o9/lcuut+btcPvzh/F0udXX5u1yy7ezadf7n9YW0TdE47H/JQHf3ATPbADxF3Lb4sLtvN7O7c/M3ApuIWxZ3Ebct3lmxCgstWjT2nR/rMuW0X7QoTvTf/u3K1HI++3UhOjsj+KtNZ2fcmleNOjvhD/7g4qw7O39+53fKX6bcXmGpQC/cfjnG+1yfbC7y/pd1H7q7byJCu/C5jQWPHfh0ZUsTEZGx0L8UFRFJhAJdRCQRCnQRkUQo0EVEElHyX4petA2bHQZeB+YC3SWaT0bVWjdUb+3VWjdUb+2qe/yVqn2Ju88bacaEBfo7BZg9P9o/Y53MqrVuqN7aq7VuqN7aVff4u5DaNeQiIpIIBbqISCImQ6A/ONEFnKdqrRuqt/ZqrRuqt3bVPf7Ou/YJH0MXEZHKmAw9dBERqQAFuohIIiZNoJvZZ8xsp5ltN7P/PtH1jIWZ3WNmbmaV/dqVi8jMvmJmr5rZVjP7gZnNmuiaijGzNbnzY5eZ3TfR9ZTDzBab2WYz25E7r/94omsaCzOrNbOXct9IVjXMbJaZPZI7v3eY2Y0TXVM5zOw/586TbWb2d2Y2ZazrmBSBbmbtxPeSrnT3q4H/McEllc3MFhNfoL13omsZo38CrnH3lcBrwJ9OcD2jKvii8luAq4CPm9lVE1tVWQaA/+Lu7wHeD3y6SurO/DGwY6KLOA9fA5509yuBVVTBPpjZu4DPAqvd/Rrivyq/Y6zrmRSBDvwh8Z2kpwHc/dAE1zMWfwV8nhG+cm8yc/d/dPfsu+CeIb5larIq54vKJx133+/uL+Yev0UEy7smtqrymNki4CPAQxNdy1iY2Qzg14FvA7h7v7sfn9CiylcHTDWzOmAa5/Gtb5Ml0K8AbjazZ83sX8ysKr6Tysw+Crzp7j+b6Fou0KeA/zPRRRTxLuCNgukuqiQYM2a2FHgfUC3f0PHXREflPL/Ve8JcDhwG/lduuOghMyvxPXoTz93fJEYm9gL7iW99+8exrqesL7ioBDP7v8CCEWbdn6tjNvGx9Drgf5vZ5T4J7qksUfefAR8e34rKV6x2d38s1+Z+Ymjge+NZ2xiV9SXkk5WZNQP/AHzO3U9MdD2lmNltwCF3f8HM2ia4nLGqA64FPuPuz5rZ14D7gC9MbFnFmdls4lPnMuA48Pdm9rvu/rdjWc+4Bbq7f2i0eWb2h8CjuQB/zsyGiP+g5vB41Tea0eo2s39PHPyfWXyn5iLgRTO73t0PjGOJoyp2zAHM7PeA24APToY3zyIu/peQXyRmVk+E+ffc/dGJrqdMNwEfNbNbgSnADDP7W3f/3QmuqxxdQJe7Z5+EHiECfbL7ELDH3Q8DmNmjwH8AxhTok2XI5YfAbwCY2RVAA5P8f0pz91fcfb67L3X3pcSJdO1kCfNSzGwNcC/wUXc/OdH1lPDOF5WbWQPxx6LHJ7imkize6b8N7HD3r050PeVy9z9190W58/oO4CdVEubkrr83zGxF7qkPAj+fwJLKtRd4v5lNy503H+Q8/pg7bj30Eh4GHjazbUA/8HuTvMeYgm8AjcA/5T5hPOPud09sSSMb7YvKJ7isctwErANeMbOXc8/9We47euXi+Qzwvdyb/24u1pfWV1BueOgR4EViCPQlzuO/ANA//RcRScRkGXIREZELpEAXEUmEAl1EJBEKdBGRRCjQRUQSoUAXEUmEAl1EJBH/D0iHnvrz/MelAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_samples = p_real.sample((200,))\n",
    "q_samples = q_real.sample((200,))\n",
    "\n",
    "plt.scatter(p_samples, [1]*len(p_samples), c='b', marker='x', alpha=0.3)\n",
    "plt.scatter(q_samples, [0]*len(q_samples), c='r', marker='o', alpha=0.3)\n",
    "plt.grid()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Wiemy, że:\n",
    "\n",
    "<img src=\n",
    "\"https://render.githubusercontent.com/render/math?math=%5Ctextstyle+%5Cbegin%7Balign%2A%7D%0Ar%5E%7B%2A%7D%3D%5Cfrac%7Bp%28x%29%7D%7Bq%28x%29%7D%3D%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7BP%28y%3D0%7Cx%29%7D%0A%3D%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7B1-P%28y%3D1%7Cx%29%7D%0A%5Cend%7Balign%2A%7D%0A\"\n",
    "alt=\"\\begin{align*}\n",
    "r^{*}=\\frac{p(x)}{q(x)}=\\frac{P(y=1|x)}{P(y=0|x)}\n",
    "=\\frac{P(y=1|x)}{1-P(y=1|x)}\n",
    "\\end{align*}\n",
    "\">\n",
    "\n",
    ", zatem iloraz dwóch dystrybucji jest proporcjonalny do prawdopodobieństwa, że wylosowany\n",
    " x pochodzi z rozkładu p(x).\n",
    "\n",
    ">Dysponując analityczną formą p(x) i q(x) możemy wyznaczyć rozkład P(y=1|x) stosując kilka przekształceń:\n",
    ">\n",
    "> <img src=\"https://render.githubusercontent.com/render/math?math=%5Ctextstyle+%5Cbegin%7Balign%2A%7D%0Ar%5E%7B%2A%7D%3D%5Cfrac%7Bp%28x%29%7D%7Bq%28x%29%7D%3D%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7BP%28y%3D0%7Cx%29%7D%0A%3D%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7B1-P%28y%3D1%7Cx%29%7D%3Dexp%5Cbig%28log%5Cbig%5B%5Cfrac%7BP%28y%3D1%7Cx%29%7D%7B1-P%28y%3D1%7Cx%29%7D%5Cbig%5D%5Cbig%29%3Dexp%5Cbig%28%5Csigma%5E%7B-1%7D%5Cbig%28P%28y%3D1%7Cx%29%5Cbig%29%5Cbig%29%0A%5Cend%7Balign%2A%7D%0A\"\n",
    "alt=\"\\begin{align*}r^{*}=\\frac{p(x)}{q(x)}=\\frac{P(y=1|x)}{P(y=0|x)}=\\frac{P(y=1|x)}{1-P(y=1|x)}=exp\\big(log\\big[\\frac{P(y=1|x)}{1-P(y=1|x)}\\big]\\big)=exp\\big(\\sigma^{-1}\\big(P(y=1|x)\\big)\\big)\n",
    "\\end{align*}\">\n",
    ">\n",
    "> dochodzimy do formy, która zawiera funkcję logitową pozwalającą nam dalej dowieść następującą własność:\n",
    ">\n",
    ">\n",
    "> <img src=\"https://render.githubusercontent.com/render/math?math=\\frac{p(x)}{p(x)%2bq(x)}=P(y=1|x)\">\n",
    ">\n",
    "> zatem prawdopodobieństwo, że wylosujemy samplę pochodzącą z p(x) (co jest podstawą estymacji ilorazu dwóch rozkładów r*(x))\n",
    "> może zostać przedstawione za pomocą p(x) i q(x) co jest w zasadzie intuicyjne.\n",
    "\n",
    "Idealny klasyfikator jest postaci: <img src=\"https://render.githubusercontent.com/render/math?math=\\frac{p(x)}{p(x)%2bq(x)}\">\n",
    "i można go zatem wyrysować na wykresie:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxJklEQVR4nO3deXxU1f3/8ddnJnvCKpCAYV8VQRGEulBBLEXAIkgrS3FpFVGQUqUtblWL9uu3ooi4UESLWBYXUBT3KtGv+hNBUSQV2WSJgIQ1JCFkmfP74+QyQ8wygUnuzOTzfDzyyCzn3vueO/d+5uTkzIwYY1BKKRX5PG4HUEopFRpa0JVSKkpoQVdKqSihBV0ppaKEFnSllIoSMW5tuEmTJqZNmzbk5eWRnJzsVoyTFqm5IXKzR2puiNzsmrv2VZX9iy++2GeMaVruncYYV3569uxpjDFm5cqVJhJFam5jIjd7pOY2JnKza+7aV1V2YI2poK7qkItSSkUJLehKKRUltKArpVSU0IKulFJRQgu6UkpFiSqnLYrIs8BQYK8x5qxy7hdgFjAYyAeuNcZ8GeqgqnZ99x0UF0PXrvbypk2wfTu0bg1Dh8KKFf7rHTtC5872tqws6NsXduyAwkLweuHzz+HAAfjxRygogIQEaNoU1q2z1+PioKQEYmKgZUswBi69FJYts9vt1AnatIFzz4Xf/Q6KiuzynTrZdf/4o10mLQ3GjIG334bMTLvt5s2hVy/46ivYvRsSE+1tHo9dztGgAezZY9fdvr3Nfvgw5Ofbtk2b2m0cOmR/GwMpKdC4MeTmgoj9nZJi99m+fXb5w4ftfXfeCbfeCtnZdl1paXabhw7Z+xs0gIYN7eXGjeFnP7P7d98+m/ucc2DCBJt1zhy7ny+8EFq1stt78kmbs29fux87d7b7wLlc3vPqqKxdoPLaqfASzDz0+cDjwIIK7r8M6Fj60wd4qvS3imDFxbB1q//6Z5/Bxo22iJaUwOrV/utt29qTfft2W3wyMsDns7ft32+L1MaN9rLPZ5cvLLRF8ehRW4hiY+19X38NzZrBzp32ReHQITh40Gbp3Bm2bbPL/fgjbNhg11NSYjN+/z3Mnm2LsFOMduyAvXttsfX5ICfHFkkRu5zPZ7e/Z49/mXXrbOE6dsxe93ptBp/PXnc+oHTfPpuzpMSuLybG5srKgnr17IvYsWP2sR0+DOvX27Z799p9FRNjt5mQYNs7EhPtC1XDhvbx5+fbbWZm2vvXrrXL+3xw/vnwwQf2Pmff9utnr2/dCu3aVfy8du1adbvYWHu9onYqvFRZ0I0xH4lIm0qaDAMWlM6P/ExEGopIc2PM7lCFVLXP6cE5J3/z5vb3vn2waBE0amSLefPm/jb9+tnfGRm2aB04AEeOQP36tkcqYgvq/v220Hk8tlgaYy87BTMnx/aUW7e2PdXMTNuz9fnsOsAWruJiu3y9enZ9RUV2WYdTtHNz7fXYWHubU/icNuAv7M59gdvyeH7a3hjbxskeEwPJyf4Xq/x8m83ZZkmJ/fF47PWjRyE+3vbomzSxy+bn29sPH7aPo3lzOOss+1x4PLBkid129+729r174Y037ItNx47+58h5Ptq1O7EnXvZ5Dabd5s3w+usVt1PhRUwQn4deWtBXVDDksgJ40Bjzcen194G/GGPWlNN2PDAeIDU1teeSJUvIzc0lJSXl1B6FCyI1N1Qv++HDJ17Py/MPm5R9M1uDBv5l8vLsZaet8/vgQXt7YaEtgj6fv6ca+Nvng9NOs9vYv9+uLzU1l927U44XUrBFM6a0W3LsmP925zanrdOD9vlOvM0ReN1Zh8iJLzRllW0XuE3nBcK5npaWy65dKcdvA3t/XJy/Fwz2RQnsfkhMtPusRYsT92lyst3Xu3b59/Fpp/00n/N8lCfwea28XS6QUmW7cBPN52f//v2/MMb0Ku++ULz1X8q5rdxXCWPMXGAuQK9evUy/fv3IyMign9O1iyCRmhuCz56ZaXvYYIcudu+2PfSDB20PvUkT2yts1cq2adrUv9yePXYcNifH9tBzcuxwRNkeOtjCGB9vi5gxtoglJtr1Nm/u76HfeWcGU6f2O158jbEFPTHR30MP7J84hTaw+Ds99MDhk8DC7tweeFtMjO1dB263vGIe2EOPj7d5iopsmwceyODPf+53/K8Sj6fiHvqxY3aZsj30PXvsNtPS/ENJ2dk/7aEHPh/l9agDn9eq2h0+nIHH06/SduGoLpyf5QlFQc8CWgZcTwd2hWC9ykVlx0xXrfKPmV92mX8MHU4catm92xaWZs1ssXEK3549/jH0hAR/oXXG0AN7wc4Qzfbtdkw9Pt4WE2cZY2wBjYmxxTOwONWv7x9Ddwp+QoIddnGKq9drfwKLtMfjH0N3ev3OGLoz/OJs23lMTnF2hn/y8vwvTklJdsipqMj2wAO36fHYF6GYGJt1794Tx9AbNLCPt2FD+xysX28L6eTJ9v7HHrP75txzYcgQ+3+FzEy7nrPP9j9ngWPl5T2vgWPoFbWLj4dBgypup8JLKAr6a8AkEVmC/WfoYR0/j3wxMf6T/rvv7Fh28+b+WS7ODJLWrf0zH77/3t7uzHJp1co/y6VZs1Of5ZKcbH9H4iyXBg1sb/tUZrk4hbRHD7sOZ5bL5ZeXP8vFeR4rel7B/7uidtnZlbdT4SWYaYuLgX5AExHJAu4BYgGMMXOAN7FTFjdjpy1eV1NhVe0JnJrWufNPp6oNHfrTZQJvC+zFldc2GDfccOL1jIwTZ95U5P77T257NSkjA74M0WReZ/pioJtv/ult5fWky5tyWFm7jIzK26nwEswsl9FV3G+AiSFLpJRS6qToO0WVUipKaEFXSqkooQVdKaWihBZ0pZSKElrQlVIqSmhBV0qpKKEFXSmlooQWdKWUihJa0JVSKkpoQVdKqSihBV0ppaKEFnSllIoSWtCVUipKaEFXSqkooQVdKaWihBZ0pZSKElrQlVIqSmhBV0qpKKEFXSmlooQWdKWUihJa0JVSKkpoQVdKqSihBV0ppaKEFnSllIoSWtCVUipKaEFXSqkooQVdKaWihBZ0pZSKElrQlVIqSmhBV0qpKKEFXSmlokRQBV1EBonIdyKyWUSmlXN/AxF5XUS+FpFMEbku9FGVUkpVpsqCLiJe4AngMuBMYLSInFmm2UTgv8aYs4F+wMMiEhfirEoppSoRTA+9N7DZGLPVGFMILAGGlWljgHoiIkAKcAAoDmlSpZRSlRJjTOUNREYCg4wx15deHwf0McZMCmhTD3gN6ALUA64yxrxRzrrGA+MBUlNTey5ZsoTc3FxSUlJC9XhqTaTmhsjNHqm5IXKza+7aV1X2/v37f2GM6VXuncaYSn+AXwPzAq6PA2aXaTMSmAkI0AH4Hqhf2Xp79uxpjDFm5cqVJhJFam5jIjd7pOY2JnKza+7aV1V2YI2poK4GM+SSBbQMuJ4O7CrT5jpgWen2NpcW9C5BrFsppVSIBFPQVwMdRaRt6T86R2GHVwLtAAYAiEgq0BnYGsqgSimlKhdTVQNjTLGITALeAbzAs8aYTBGZUHr/HGA6MF9EvsEOu/zFGLOvBnMrpZQqo8qCDmCMeRN4s8xtcwIu7wIGhjaaUkqp6tB3iiqlVJTQgq6UUlFCC7pSSkUJLehKKRUltKArpVSU0IKulFJRQgu6UkpFCS3oSikVJbSgK6VUlNCCrpRSUUILulJKRQkt6EopFSWC+nAupeqSbdu2sWrVKr766it27drFbbfdRvfu3dmwYQPPPvssnTp14sILL6RLly7Yb11UKjxoQVeqVGZmJqNGjWL9+vUAxMbG0rx5c8aNGwfAhg0bmDVrFoWFhQC0aNGC3/zmN0ybNo3U1FTXcivl0CEXVadt3ryZ999/H4BWrVqRlpbGzJkz+fLLL8nNzWX79u1ceumlAFxxxRXk5+ezceNGnn76afr06cMzzzxDXFwcAEeOHHHtcSgF2kNXdVRubi733HMPs2fP5uyzz2b16tXUq1eP9957r9LlvF4vHTt2pGPHjlx//fXk5+eTlJSEMYZLLrmEdu3aMXPmTFq0aFFLj0QpP+2hqzrn008/5ZxzzmHmzJlcffXVvPZa2W9UDF5SUhIAJSUl/OpXv2L58uWcccYZLFq0KFRxlQqaFnRVp6xatYq+fftSUlLChx9+yLx582jevPkprzcmJoa7776b9evX061bN8aOHcv48eM5evRoCFIrFRwt6KpO6d27NzNnzuTrr7+mb9++IV9/hw4dyMjI4Pbbb2fVqlUcO3Ys5NtQqiJa0FXUO3r0KNdeey2bN29GRJg8eTL169evse3FxMTw97//ndWrV9OwYUOOHj1KdnZ2jW1PKYcWdBXV8vLyGDp0KAsWLGDVqlW1um1n9svNN9/MlClT+OGHH2p1+6ru0YKuolZRUREjRowgIyOD+fPnM3bsWFdyTJgwgUOHDjFgwAD27NnjSgZVN2hBV1HJGMP111/Pu+++y9NPP83VV1/tWpY+ffrw4IMPsnPnToYMGUJeXp5rWVR004KuolJ+fj5btmzhvvvu43e/+53bcejWrRsvvPACa9euZcKECW7HUVFK31ikolJycjIffPABsbGxbkc5bujQocydO5eePXu6HUVFKe2hq6iyadMmrrzySrKzs4mLiwu7D8+6/vrr6dGjBwA//vijy2lUtNGCrqJGQUEBI0eOJCMjg/z8fLfjVOrBBx/krLPO0pkvKqS0oKuocffdd7Nu3Tqef/55Wrdu7XacSl1xxRUcPXqUsWPH4vP53I6jooQWdBUVPvroIx5++GEmTJjA4MGD3Y5TpS5dujB79mw+/PBDHn/8cbfjqCihBV1FhTvvvJO2bdvy0EMPuR0laNdeey2DBw9m2rRpbNmyxe04KgoEVdBFZJCIfCcim0VkWgVt+onIVyKSKSIfhjamUpV77bXXWL58OSkpKW5HCZqI8M9//pNmzZqxceNGt+OoKFDltEUR8QJPAL8AsoDVIvKaMea/AW0aAk8Cg4wxO0SkWQ3lVeoEe/bsoUmTJjRq1IhGjRq5Hafa0tPT2bRpU1hNr1SRK5geem9gszFmqzGmEFgCDCvTZgywzBizA8AYsze0MZX6KZ/Px/Dhw7nsssvcjnJKYmNj8fl8/Otf/+LQoUNux1ERTIwxlTcQGYnteV9fen0c0McYMymgzaNALNAVqAfMMsYsKGdd44HxAKmpqT2XLFlCbm5uRP2Z7IjU3BC52cvmfuONN5gxYwbTpk3jl7/8pYvJqlbVPv/+++/5/e9/z4gRI5g0aVKF7WpbtBwrkaSq7P379//CGNOr3DuNMZX+AL8G5gVcHwfMLtPmceAzIBloAmwCOlW23p49expjjFm5cqWJRJGa25jIzR6YOycnxzRr1sxcdNFFxufzuRcqSMHs85tuusl4vV6zbt26mg8UpGg4ViJNVdmBNaaCuhrMkEsW0DLgejqwq5w2bxtj8owx+4CPgLODWLdSJ2XGjBns3buXhx9+OOzeDXqy7r//fho0aMCf/vQnt6OoCBVMQV8NdBSRtiISB4wCyn4J43Kgr4jEiEgS0Af4NrRRlbJ8Ph9vv/02V111Fb1793Y7Tsg0btyY22+/nXfeeYeVK1e6HUdFoCpnuRhjikVkEvAO4AWeNcZkisiE0vvnGGO+FZG3gXWADztEs74mg6u6y+Px8Mknn5Cbm+t2lJCbNGkS7733XtT81aFqV1CftmiMeRN4s8xtc8pcfwiInHd1qIj0448/kpSURL169WjYsKHbcUIuISGBd955x+0YKkLpO0VVRPnTn/5E165dKSoqcjtKjcrJyWHmzJmUlJS4HUVFEC3oKmJkZWWxcOFCfvOb30T9G3Heeecdbr31VpYtW+Z2FBVBtKCriLFw4ULi4uKYOnWq21Fq3IgRI+jcuTP333+/MzVYqSppQVcRYcuWLbz77rtMmDCBtLQ0t+PUOK/Xyx133MG6detYsWKF23FUhNCCriLCK6+8gtfrrVNztEePHk3btm154IEHtJeugqIFXUWEqVOnsmDBAlq0aOF2lFoTGxvLtGnTaNSoEUeOHHE7jooAWtBV2Dt27BhAnRhqKeuGG27grbfeon79+m5HURFAC7oKa/n5+bRr144nn3zS7SiucN5gtG3bNnbu3OlyGhXutKCrsPbcc8+xa9cuunXr5nYU1+Tl5dGtWzemT5/udhQV5rSgq7Dl8/mYOXMm5513HhdddJHbcVyTnJzMmDFjeP7558nOznY7jgpjWtBV2FqxYgWbNm3itttuq/OfbTJlyhQKCgqYM2dO1Y1VnaUFXYWtxx9/nJYtW3LllVe6HcV1Z5xxBoMGDeKJJ544/k9ipcrSgq7C1hNPPMGzzz5LTExQnyEX9f7whz+wf/9+1qxZ43YUFab0TFFhq2PHjnTs2NHtGGFj4MCB7Ny5s05O31TB0R66CjtHjx7lt7/9LWvXrnU7SljxeDzHi7l+CqMqjxZ0FXZeeuklFi5cyMGDB92OEnZ8Ph+//OUvue2229yOosKQFnQVdp566ik6d+5M//793Y4SdjweD40bN+a5554jPz/f7TgqzGhBV2Hlq6++4rPPPmPChAl1fqpiRSZMmMChQ4d48cUX3Y6iwowWdBVW/vnPf5KQkMA111zjdpSw9fOf/5wuXbrw1FNPuR1FhRkt6CqstG/fnsmTJ9OoUSO3o4QtEeHGG2/k888/Z/16/S525afTFlVYqQvfRhQKY8eOxev1cvrpp7sdRYUR7aGrsPHJJ59E/Zc/h0rTpk255ZZb9C8ZdQIt6Cos7Nixg759+/LQQw+5HSViFBYW8swzz/Dxxx+7HUWFCS3oKiwsWLAAYwyjR492O0rEEBHuuOMOHnnkEbejqDChBV25zhjD/Pnz6d+/P23btnU7TsSIjY1l3LhxvP766/qxugrQgq7CwMcff8yWLVu47rrr3I4Sca699lqKi4tZuHCh21FUGNCCrly3fPly6tWrx4gRI9yOEnHOOussevXqxfz5892OosKAFnTlun/84x+sWbOG5ORkt6NEpOuuu47Y2Fj97BulBV25z+Px0KlTJ7djRKybbrqJ1atX6xRGpQVdueu6667j/vvvdztGRHM+8yYnJ0c/VreO04KuXJOdnc3zzz/PkSNH3I4S8T799FNSU1PJyMhwO4pyUVAFXUQGich3IrJZRKZV0u48ESkRkZGhi6ii1csvv0xJSQljx451O0rE69GjB7GxsTrbpY6rsqCLiBd4ArgMOBMYLSJnVtDuf4F3Qh1SRaeFCxfStWtXunXr5naUiJeYmMiIESNYunQpBQUFbsdRLgmmh94b2GyM2WqMKQSWAMPKaXcLsBTYG8J8Kkpt27aNTz75hDFjxujnnofI2LFjycnJYcWKFW5HUS4J5tMWTwd2BlzPAvoENhCR04HhwCXAeRWtSETGA+OB4+N9ubm5ETnuF6m5ITyy7969m8GDB9O2bdugs4RD7pNVG9mdbzOaNWsWTZo0Cck6I3WfR2puOMXsxphKf4BfA/MCro8DZpdp8xLws9LL84GRVa23Z8+exhhjVq5caSJRpOY2JnKzR2puY2ov+xtvvGEyMzNDtr5I3eeRmtuYqrMDa0wFdTWYHnoW0DLgejqwq0ybXsCS0j+dmwCDRaTYGPPqyb3MqGi2Y8cO9u7dS8+ePXW4JcQGDx7sdgTlomDG0FcDHUWkrYjEAaOA1wIbGGPaGmPaGGPaAC8DN2sxVxV56qmn+NnPfsb+/fvdjhKVPvroI6ZPn+52DOWCKgu6MaYYmISdvfIt8KIxJlNEJojIhJoOqKKLz+dj8eLFDBw4MGTjvOpEGRkZ/PWvf2Xnzp1VN1ZRJah56MaYN40xnYwx7Y0xD5TeNscYM6ecttcaY14OdVAVHT799FO2b9/OmDFj3I4StZx5/YsXL3Y5iapt+k5RVasWLVpEYmIiV1xxhdtRolb79u3p06cPixYtcjuKqmVa0FWtMcbw9ttvM2zYMFJSUtyOE9XGjBnD119/TWZmpttRVC3Sgq5qjYiQmZnJww8/7HaUqHfVVVfRpk0bsrKy3I6ialEw0xaVCpnExEQSExPdjhH1UlNT2bp1q04LrWO0h65qRV5eHn369OGtt95yO0qdISIUFxdz+PBht6OoWqIFXdWK1157jc8//5ykpCS3o9QZxcXFdOjQgbvuusvtKKqWaEFXtWLx4sWkp6fTt29ft6PUGTExMfTu3ZsXXniB4uJit+OoWqAFXdW4AwcO8PbbbzNq1Cg8Hj3katOYMWPIzs7mP//5j9tRVC3Qs0vVuKVLl1JUVMTo0aPdjlLnXHbZZTRs2FDnpNcRWtBVjWvXrh033HADPXr0cDtKnRMfH8/IkSN55ZVXyM/PdzuOqmE6bVHVuAEDBjBgwAC3Y9RZU6ZMYdSoUcTHx7sdRdUwLeiqRn3xxRc0a9aMli1bVt1Y1YiuXbvStWtXt2OoWqBDLqpGTZgwgeHDh7sdo87btm0bU6dO5cCBA25HUTVIC7qqMZs2bWLNmjX6z9AwcODAAR5++GGWLl3qdhRVg7SgqxqzZMkSRISrrrrK7Sh1Xo8ePejcubPOdolyWtBVjTDGsGjRIvr27Ut6errbceo8EWHs2LF8+OGH+oFdUUwLuqoRW7duZePGjfpFFmFk9OjRGGNYsmSJ21FUDdFZLqpGtG/fnh9++IHk5GS3o6hSHTp0YMCAAeTl5bkdRdUQLeiqxqSlpbkdQZXx3nvv6UfqRjEdclEh98knn3DppZeyZcsWt6OoMpxivn//fpeTqJqgBV2F3HPPPcdnn31Gamqq21FUOW688UZ69+6NMcbtKCrEtKCrkDp69CgvvvgiV155pX5vaJg6//zz2bp1K6tWrXI7igoxLegqpF5//XUOHz7M1Vdf7XYUVYHhw4cTHx+vc9KjkBZ0FVILFiygZcuW9O/f3+0oqgINGjRg6NCh+sUXUUgLugqpgQMHMm3aNP0iizA3ZswY9u7dy/vvv+92FBVCOm1RhdTkyZPdjqCCMHjwYP79739z4YUXuh1FhZB2o1TIvPrqq+Tm5rodQwUhISGBsWPH6j+uo4wWdBUS69atY/jw4cyfP9/tKCpIhYWFzJgxg9dff93tKCpEdMhFhcSCBQuIiYlh1KhRbkdRQYqNjWXu3LmkpaVx+eWXux1HhYD20NUpKywsZMGCBVx++eU0adLE7TgqSCLC73//e/7v//6P7777zu04KgSCKugiMkhEvhORzSIyrZz7x4rIutKfT0Xk7NBHVeFq+fLlZGdnM378eLejqGq65ppr8Hq9PPvss25HUSFQZUEXES/wBHAZcCYwWkTOLNPse+BiY0x3YDowN9RBVfj68MMPadWqFb/4xS/cjqKqKS0tjaFDhzJ//nyKiorcjqNOUTA99N7AZmPMVmNMIbAEGBbYwBjzqTHmYOnVzwD9RoM65PHHH2fNmjV4vV63o6iTcMMNN9C9e3eys7PdjqJOkVT1AT0iMhIYZIy5vvT6OKCPMWZSBe2nAl2c9mXuGw+MB0hNTe25ZMkScnNzI3LqVKTmhtBmLykpqbVCrvu89mnu2ldV9v79+39hjOlV7p3GmEp/gF8D8wKujwNmV9C2P/AtcFpV6+3Zs6cxxpiVK1eaSBSpuY0JXfaioiLToUMHM2vWrJCsryq6z2vWzp07TXZ29gm3RULu8kRqbmOqzg6sMRXU1WCGXLKAlgHX04FdZRuJSHdgHjDMGKMftlwHrFixgs2bN9O6dWu3o6hTtG/fPtq3b8+jjz7qdhR1CoIp6KuBjiLSVkTigFHAa4ENRKQVsAwYZ4zZGPqYKhzNnj2b9PR0hgwZ4nYUdYqaNGnCoEGDmDt3LgUFBW7HUSepyoJujCkGJgHvYIdTXjTGZIrIBBGZUNrsr8BpwJMi8pWIrKmxxCosrF+/ng8++ICJEycSE6PvT4sGt9xyC9nZ2bz00ktuR1EnKagz0RjzJvBmmdvmBFy+HvjJP0FV9Jo9ezYJCQnccMMNbkdRITJgwADOOOMMZs+ezbhx49yOo06CvlNUnZRbb72VZ555htNOO83tKCpERIRJkybx5ZdfsmnTJrfjqJOgfyurk9K5c2c6d+7sdgwVYtdccw2XX345LVu2rLqxCjvaQ1fVUlRUxPjx41m7dq3bUVQNSE5OPl7MS0pKXE6jqksLuqqWRYsW8fTTT7Nr109mrqoo4fP5uPzyy/njH//odhRVTVrQVdB8Ph8PPvgg3bt3Z/DgwW7HUTXE4/HQtGlT5s2bx6FDh9yOo6pBC7oK2quvvsqGDRu4/fbbERG346ga9Oc//5mCggKWLVvmdhRVDVrQVVCMMfz973+nQ4cO/PrXv3Y7jqphXbp0Yfjw4bzyyivk5OS4HUcFSQu6CkpxcTFDhgzhnnvu0U9VrCNuv/12cnNzeeKJJ9yOooKk0xZVUGJjY7nvvvvcjqFqUa9evbjrrruYOHGi21FUkLSHrqr01ltv8eqrrzqfqKnqkAEDBlC/fn23Y6ggaUFXlSoqKmLy5Mncc889WtDrqNWrV9OrVy+dqhoBtKCrSs2dO5fNmzfzwAMP4PHo4VIXNWnShHXr1nHXXXe5HUVVQc9QVaHs7GzuuusuBgwYoB+RW4e1bduWKVOm8K9//YtVq1a5HUdVQgu6qtCdd95Jbm4ujz32mM47r+PuvvtumjdvzsSJE/UjAcKYFnRVoYsvvpj77ruPM8880+0oymX16tVjxowZfPHFF7z88stux1EV0GmLqkJjx451O4IKI6NHjyYhIYFhw4a5HUVVQHvo6iemT5/Oo48+qrNa1AlEhBEjRuD1etm3b58eH2FIC7o6wapVq7j33nv55ptvdNxclWvjxo107tyZp59+2u0oqgwt6Oq4AwcOcNVVV5Gens4jjzzidhwVpjp06ECPHj2YMmUK33zzjdtxVAAt6AqwH4179dVXs2vXLl566SUaNGjgdiQVpjweD//+979p2LAhV155JYcPH3Y7kiqlBV0BdqjlzTffZObMmfTu3dvtOCrMpaWl8eKLL7J161auvvpqncoYJrSgKwDOP/98Vq1axc033+x2FBUhLrroImbNmkViYqIW9DChBb2OW758Oe+++y4A5513nv4jVFXLxIkTWbx4MXFxceTm5rodp87Tgl6HvfTSS4wcOZK//e1v+Hw+t+OoCCUiHDhwgF69eumHuLlMC3odNW/ePEaNGkWfPn1488039YO31Clp0KABF1xwAX/729+49dZbdQjGJfpO0TrG5/Px+OOPs3TpUgYOHMjSpUtJSUlxO5aKcF6vl3nz5lG/fn0effRRNmzYwKJFi2jUqJHb0eoU7ZbVMR6Ph6KiIqZMmcIbb7yhxVyFjMfj4dFHH2XOnDm8//773HbbbW5HqnO0h14HFBYW8uSTT3LBBRfQu3dv/vCHP3DJJZe4HUtFqRtvvJFu3bqRnp4OwA8//EBKSoq+t6EWaA89ipWUlLBs2TK6devGH//4RxYsWACg4+Wqxl1wwQW0atUKgJtuuokOHTowe/ZsCgoKXE4W3fTMjlJLly7ljDPO4MorrwRgxYoVzJ492+VUqi6699576d69O5MnT6Z169bce++97N692+1YUUkLepTIzs5m8eLFx+cCb9u2jQYNGvDCCy+QmZnJkCFDdI65csW5557Lf/7zH95//3169+7NfffdxzPPPON2rKgU1Bi6iAwCZgFeYJ4x5sEy90vp/YOBfOBaY8yXIc4KWVmwejVkZ0PTpnDeeVA6ThfUMh4PGGN/Klo+mG1kZcGhQzB3bvA5Klt3eRm3boW1a2HnTti/H3w+iI+3P8nJHEhOZpnPx9qiIj7NyuKrvXsBWPiLXzBm714m5+Vxa/PmyBtvwIwZsHEj5OeDCPzv/8KUKZCaaq8fPgzffw9HjtjtN20K55wDXbrAwYPw4492WbDZcnJsrtxcOHbMrsORlARxcVBQYO+LjYV27aB7d9i9265r3z77WIuK7PpiY+1PSQkUF9v1eL3QrJndP8bYZadOhb/8Bc44Az7/HDZvtusAmyEuDpKTIS0NUlKgsNDmz8+364uP9+cSgcaN4bTT7OX9+2HPHnufc4yAzdWgATRvbq8fOWJzJiZC69Z2Xxljt9eoEXz3HXzzjd1mQoLNFB8PEybAY4/Z/VOvns3zzTfwww92n4Pdr4WF/scfF2fXf/bZ9jE1agSdO//0eHOOn1Wr7PZFbPuGDeHAAcjLs/lat7bLO8dcdrZte/Cg//mvX9+269TJbicYQZwzIsIll1zCJZ06sWnQIBocPQqvvBL8uRMtTqaGVUOVBV1EvMATwC+ALGC1iLxmjPlvQLPLgI6lP32Ap0p/h05WFixfbg/S1FRbTJYvh2HDKt4hgct4vfDRR/YA7tvXnnBllw9mG06b5s2Dz1HBugtefpmC7t0pWLWKgsRECnJzqf/117TIz+fY9u0s2LWL/UeOsN8YfvD52GYM13m93FCvHgcTE7lh927qeTz0bNqU+7t35xc+Hz1Xr4bUVGKTkmzB+PRTiImxhcIpUiUlsG6dLWherz3pfT5/YcnKskX+s8/sQRcfbw/A3FxbGLKy7H4sKbHLBTp61H/Z67X3Z2bC+vV2XYWFtqAHKiryF2ZHSQns2mV/YmJscTQGNmywxbwsY2wxLiy0jycpyZ/PeZEsLLSXRWy2gwdtES8psYVexL+PAnMUFNiC7/XaNjExttju2mVfQFq3ti8+r75qi2diol1u926bu1Eju46PPoI+fey+2LbNZjTGdg6OHbPrB/uiJmLvP3oUtm+H88+Hli3tbbt2+Y8357jatw8+/thur7DQdggOHrQvBrm5dp05OXY7S5bAhRfa7G+/DTt2QNu2sHevzZWTYx/Drl22sFemOudladuODRva+4I9d6LFydSwagqmh94b2GyM2QogIkuAYUBgQR8GLDD2LWKfiUhDEWlujAndQNnq1XZH1K9vrzu/V6+ucGeM/93v+L/MTIwI5uhRjM9Hx0aNeKNFCzj/fH7z2GOsmT4dkpIwxmDy8ujRqhWv3HMPAANnzOC/27dj/vpXSEjAGEPftm154eabwevlvKlT2Z6djfH54I47MLGxDBkyhOeeew6wHzN6/IsAiosxPh+jL76Yf06cCPXrU/+aaygq8waMSeecw+yYGCQ3l/EHDwIQD7QA2oiQ5PVCQQHt4uPZmJ5O+5gYPL162RP1889t70/EX8BF/MXCKcJgLx85YnugxthC5xR0r9cebB6PLVB5ebZdcrI9yZ3CU7YIezz+Au/02r1efyHNz7frDZZTjJ1eO5x4uTzOYzl2zP94YmL8LzTOPhGxjymYf9KJ+ItsTIz/hcrn8/9V8e23/jZer3+fgS2sIvbFcO1aux+9Xrs/EhPtepz9VVLif9EpLLTr8Hhgyxb7l8mePdC1q/+4d86LDz6wz339+ra3nZdnb//uO7sc2Od73TpbTPbssbcdO2Zf2Ldvt7eDLejOdvLyKt831TkvT+Icjiq18PiDKeinAzsDrmfx0953eW1OB04o6CIyHhgPkJqaSkZGBrm5uWRkZFSdoqTEnhCBPcKUFFtUKljed9pppLZpg4ggxcWIx0O9xo3J6NIFfD6S27enQ3y87dUAUlhIwxYtyCjdRvPOnYk57TTEGCQxEYAmTZqQkZJCLtClRw9aHDlilzUGEhJo2rTp8cfTp08fCgoK7PYLCsDjoWmHDsfXf/24cXgKCohr2JC42FjiiotpmZ5ORoMG4PPx4oED1IuPJyE+3j8M4PGQYYwtCF4vPxhjiwLAJZf4CzTYfRO4v0p7oLnp6WT84x/+28Bf6JzbfD5bvJyi6BTB0sdxwrBEZZwi7yx/sh8x4PWSm5ZGxvTpwW3XyRz4eMq7rzpvU3eWcR6/czkuzj985PPZ58Yp7qXtctPSyPjTn2y7wBcEp21FOTwe/wtokya2beBx75wXv/ylfb7AFuniYnssFBba+8G/zeRk/wt7u3Z2GwUFx88DfD77mFJSyPX5Kj8/q3NensQ5fLKCriu1KcjHfyrZgyno5f0nrezRF0wbjDFzgbkAvXr1Mv369SMjI4N+/fpVneKVV2zvznlVA9uTSEqCCpbv54yf1q8P/+//2QMd7J/t559PvyFD7PLDh/u3kZ9/vIj1GzvWv43ANrm5ZKSk8Py4cSfmcNo42w/M5aw7IH+/Sy+1Pabu3f0Z9++3bfftsyeZ81nTTm8wLs5ebtjQX8h79vT30J1xZLB/8ufl2ZM7oIee8dBD9Pvzn20BcHqpIif20EtKbCZn7rDHY4vD7t3+wubsT0fZHrrHY4vE0aP2ckqK/8//YDmZUlLImD6dfnff7R/Pr2yZwMcTF2czOC8sHo997E6BA/t4yxtycR6LM9bvDLc4f+0kJUGrVvZ5OnDA7qPkZHtMOPsdyPif/6HfI4/47z9wwD8Wf+CAbevx+HM4OZOT7e8mTWDoUDv81bWr/7h3zot33vEfX99/b/8qSEqyz63TQ/f57GNITbXHD9jhn6IiO+zj9NB9PjjrLOjalYzk5MrPz+qclydxDp+soOtKbQry8Z9K9mBmuWQBLQOupwO7TqLNqTnvPHvQ5eTYAy4nx16v7B83gcu0b2+L5P799nJ5ywezDaeNMz4bTI7K1j1y5E8zNmtm/3x2ettO4THGnnzOP9uKimy7xET7u107e8AY4y86xtgi4POd2BM0xi6TnOzvwTo/Ti/QGWZp3NhuKy8PWrQ4cWw6UGDPI3C8PibGPybsFJJgBA6ZOGKq6IM4jyU+3r98cbF/ucCeubMvnaGRihjj/2vFGP+wiNN7jomxRTMmxn9/QoJ/SKpRI3t7bi706OHvIScl+R+ns78Ce+/Oc+jz2WMjN9f+szPweHOOq+7d7ZDK/v1228nJ9vbOne2xdeCAfb67d7dt0tLsOuPj7X2tW9t2+/bZguNsx+kcVKQ65+XJnMPRpBYefzA99NVARxFpC/wAjALGlGnzGjCpdHy9D3A4pOPnYMeYhg2z400//mj/wXbxxZWPPQUuk5dn/xnqnHBJST9dPphtOG3WrAk+R1Xrbt78pxmrmOVCw4a2Z5iebu9r3Bj694evvrL/HMvPtydvmzZ2fDdwlovXa++rziyXxo3t4zDGnuynMsvF+SdrdWe5iNhMoZrl0rJlaGe5XHfdibNc2rb1z3JJSICf/9zfo7/oIv8sl3r17LqDmeXSsuWJMyMCj6ujR/2zXM46q+JZLoMH+2e5DBzon+USE+Of5eJsZ/Pmkz+uT6VtNKqFx19lQTfGFIvIJOAd7LTFZ40xmSIyofT+OcCb2CmLm7HTFq8LWcJA6enVf/DVXSaY9unp9kC/4orQZDmZx3UqMjJs4Y80GRl2al4kysiAG2+smXU7x0+ZIb9KBdsrrKqgB24/GLV9rIebGn78Qc1DN8a8iS3agbfNCbhsgImhjaaUUqo69J2iSikVJbSgK6VUlNCCrpRSUUILulJKRQlx6wtdRSQb2A40AfZV0TwcRWpuiNzskZobIje75q59VWVvbYxpWt4drhX04wFE1hhjerka4iREam6I3OyRmhsiN7vmrn2nkl2HXJRSKkpoQVdKqSgRDgV9rtsBTlKk5obIzR6puSFys2vu2nfS2V0fQ1dKKRUa4dBDV0opFQJa0JVSKkqETUEXkVtE5DsRyRSRf7idpzpEZKqIGBFp4naWYInIQyKyQUTWicgrItLQ7UyVEZFBpcfHZhGZ5naeYIhISxFZKSLflh7Xf3A7U3WIiFdE1orICrezVEfpV2C+XHp8fysi57udKRgi8sfS42S9iCwWkYTqriMsCrqI9Md+L2l3Y0xXYIbLkYImIi2xX6C9w+0s1fQecJYxpjuwEbjd5TwVCvii8suAM4HRInKmu6mCUgzcZow5A/gZMDFCcjv+AHzrdoiTMAt42xjTBTibCHgMInI6MBnoZYw5C/tR5aOqu56wKOjATcCDxphjAMaYvS7nqY6ZwJ8p5yv3wpkx5l1jjPNdcJ9hv2UqXB3/onJjTCHgfFF5WDPG7DbGfFl6+Qi2sJzubqrgiEg6MASY53aW6hCR+sDPgWcAjDGFxphDroYKXgyQKCIxQBIn8a1v4VLQOwF9RWSViHwoIhHxnVQi8ivgB2PM125nOUW/A95yO0QlKvoS8oghIm2AHkCkfEPHo9iOykl+q7dr2gHZwL9Kh4vmiUgV36PnPmPMD9iRiR3Abuy3vr1b3fUE9QUXoSAi/wHSyrnrztIcjbB/lp4HvCgi7UwYzKmsIvcdwMDaTRS8yrIbY5aXtrkTOzSwsDazVVNQX0IerkQkBVgKTDHG5LidpyoiMhTYa4z5QkT6uRynumKAc4FbjDGrRGQWMA24291YlRORRti/OtsCh4CXROS3xph/V2c9tVbQjTGXVnSfiNwELCst4J+LiA/7ATXZtZWvIhXlFpFu2J3/tdjv1EwHvhSR3saYPbUYsUKV7XMAEbkGGAoMCIcXz0rU/JeQ1xARicUW84XGmGVu5wnShcCvRGQwkADUF5F/G2N+63KuYGQBWcYY5y+hl7EFPdxdCnxvjMkGEJFlwAVAtQp6uAy5vApcAiAinYA4wvyT0owx3xhjmhlj2hhj2mAPpHPDpZhXRUQGAX8BfmWMyXc7TxWOf1G5iMRh/1n0msuZqiT2lf4Z4FtjzCNu5wmWMeZ2Y0x66XE9CvggQoo5peffThHpXHrTAOC/LkYK1g7gZyKSVHrcDOAk/plbaz30KjwLPCsi64FC4Jow7zFGg8eBeOC90r8wPjPGTHA3Uvkq+qJyl2MF40JgHPCNiHxVetsdpd/Rq2rOLcDC0hf/rdTUl9aHUOnw0MvAl9gh0LWcxEcA6Fv/lVIqSoTLkItSSqlTpAVdKaWihBZ0pZSKElrQlVIqSmhBV0qpKKEFXSmlooQWdKWUihL/H8csZWzUW1CnAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "optim_classifier_vals = torch.exp(p_real.log_prob(x)) / (torch.exp(p_real.log_prob(x)) + torch.exp(q_real.log_prob(x)))\n",
    "\n",
    "plt.scatter(p_samples, [1]*len(p_samples), c='b', marker='x', alpha=0.3)\n",
    "plt.scatter(q_samples, [0]*len(q_samples), c='r', marker='o', alpha=0.3)\n",
    "plt.plot(x, optim_classifier_vals, 'k--', label='optim classifier')\n",
    "plt.grid()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Czyli chodzi nam o to aby zaproksymować tą funkcję, że użyciem dostępnych danych. Tutaj wchodzi probabilistyczny klasyfikator, którego\n",
    "wyjściem będą wartości z przedziału [0,1]."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch idx(0) with avg loss: 0.644580066204071\n",
      "epoch idx(10) with avg loss: 0.5642979741096497\n",
      "epoch idx(20) with avg loss: 0.5631875991821289\n",
      "epoch idx(30) with avg loss: 0.5634828805923462\n",
      "epoch idx(40) with avg loss: 0.5629462003707886\n",
      "epoch idx(50) with avg loss: 0.5627261996269226\n",
      "epoch idx(60) with avg loss: 0.5630586743354797\n",
      "epoch idx(70) with avg loss: 0.5628475546836853\n",
      "epoch idx(80) with avg loss: 0.562985897064209\n",
      "epoch idx(90) with avg loss: 0.5630849003791809\n"
     ]
    }
   ],
   "source": [
    "from torch import nn\n",
    "from datetime import datetime\n",
    "import torch.nn.functional as F\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "from torch.optim import RMSprop\n",
    "from pathlib import Path\n",
    "\n",
    "\n",
    "class ProbabilisticClassifier(nn.Module):\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.fc1 = nn.Linear(1, 32)  # input x of shape (1,)\n",
    "        self.fc2 = nn.Linear(32, 64)\n",
    "        self.fc3 = nn.Linear(64, 1)\n",
    "\n",
    "    def forward(self, input_x):\n",
    "        output_y = F.relu(self.fc1(input_x))\n",
    "        output_y = F.relu(self.fc2(output_y))\n",
    "        return torch.sigmoid(self.fc3(output_y))\n",
    "\n",
    "\n",
    "p_samples_classification = p_real.sample((10000,))\n",
    "q_samples_classification = q_real.sample((10000,))\n",
    "\n",
    "data = torch.hstack((p_samples_classification, q_samples_classification))\n",
    "y = torch.hstack((torch.ones_like(p_samples_classification),\n",
    "                                         torch.zeros_like(q_samples_classification)))\n",
    "\n",
    "output_folder = Path('.', 'output', f'{datetime.now().strftime(\"%d-%m-%Y-%H-%M-%S\")}')\n",
    "output_folder.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "# First attempt using classic NN-based classifier\n",
    "dataset = TensorDataset(data, y)\n",
    "dataloader = DataLoader(dataset, batch_size=64, shuffle=True, drop_last=True)\n",
    "model = ProbabilisticClassifier()\n",
    "optimizer = RMSprop(model.parameters(), lr=0.0001)\n",
    "criterion = nn.BCELoss()\n",
    "\n",
    "for epoch in range(100):\n",
    "    losses = 0\n",
    "    for batch, labels in dataloader:\n",
    "        optimizer.zero_grad()\n",
    "        output = model(batch.unsqueeze(1))\n",
    "        minibatch_loss = criterion(output.squeeze(), labels)\n",
    "        minibatch_loss.backward()\n",
    "        optimizer.step()\n",
    "        losses += minibatch_loss\n",
    "\n",
    "    with torch.no_grad():\n",
    "        predictions = (model(x.unsqueeze(1)).squeeze())\n",
    "        fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
    "        axs[1].plot(x, optim_classifier_vals, 'k--', label='P(y=1|x)')\n",
    "        axs[1].plot(x, predictions, 'r:', alpha=0.8, label='estimated')\n",
    "        axs[1].grid()\n",
    "        axs[0].plot(x, ratio, 'k--', label='real p(q)/q(x)')\n",
    "        axs[0].plot(x, torch.exp(torch.logit(predictions)), 'r:', label='estimated')\n",
    "        axs[0].legend(loc='upper left')\n",
    "        axs[0].grid()\n",
    "        plt.savefig(output_folder / f'{epoch}.png')\n",
    "        plt.close(fig)\n",
    "\n",
    "    if epoch % 10 == 0:\n",
    "        print(f'epoch idx({epoch}) with avg loss: {losses/len(dataloader)}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "W końcu możemy wyestymować wartość KL z założeniem, że nie znamy p(x) ani q(x) a dysponujemy tylko samplami\n",
    "z p(x)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated monte carlo kl divergence: 0.4381503760814667\n"
     ]
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    monte_carlo_c = torch.mean(torch.log(torch.exp(torch.logit(model(xp_mc_samples.unsqueeze(1))))))\n",
    "    print(f'estimated monte carlo kl divergence: {monte_carlo_c}')\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "pycharm-8e6839ef",
   "language": "python",
   "display_name": "PyCharm (buzz-based-anomaly)"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}