MINE_in_TF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MINE minimal working example\n",
    "\n",
    "This notebook gives an example of computing the mutual information using Mutual Information Neural Estimator (see 1801.04062) using Tensorflow.\n",
    "\n",
    "It is based on the Torch implementation at: https://github.com/MasanoriYamada/Mine_pytorch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import tensorflow.contrib.layers as layers\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.feature_selection import mutual_info_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data of two dependent variables with non zero MI. Compute it analytically and numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6565591963509047\n",
      "0.6564880831569115\n"
     ]
    }
   ],
   "source": [
    "# data\n",
    "var = 0.2\n",
    "def func(x):\n",
    "    return x\n",
    "\n",
    "def gen_x():\n",
    "    return np.sign(np.random.normal(0.,1.,[data_size,1]))\n",
    "\n",
    "def gen_y(x):\n",
    "    return func(x)+np.random.normal(0.,np.sqrt(var),[data_size,1])\n",
    "\n",
    "data_size = 20000\n",
    "x=gen_x()\n",
    "y=gen_y(x)\n",
    "p_y_x=np.exp(-(y-x)**2/(2*var))\n",
    "p_y_x_minus=np.exp(-(y+1)**2/(2*var))\n",
    "p_y_x_plus=np.exp(-(y-1)**2/(2*var))\n",
    "mi=np.average(np.log(p_y_x/(0.5*p_y_x_minus+0.5*p_y_x_plus)))\n",
    "print(mi)\n",
    "mi_numerical = mutual_info_regression(x.reshape(-1, 1), y.ravel())[0]\n",
    "print(mi_numerical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the MINE network and train it. Plot convergence of the loss towards true MI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/zgubic/anaconda3/envs/tf/lib/python3.6/site-packages/matplotlib/figure.py:445: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
      "  % get_backend())\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xl8VOXZ//HPlckeyB62BEjY9wQIi4ggImsVrCu0tfKAxaXUqtVWW4tK2+cnWteqxbX6WATcjYoiCuKC7DthCwRI2LLv62Tu3x8zGUNIyCRkIZPr/XrlxZwz95y55mT45p57zrmPGGNQSinlXjxaugCllFKNT8NdKaXckIa7Ukq5IQ13pZRyQxruSinlhjTclVLKDWm4K6WUG3Ip3EVkqogcEJEkEXmghvufFpEdjp+DIpLT+KUqpZRyldR1EpOIWICDwCQgFdgMzDbGJNbS/nfAUGPM3EauVSmllIs8XWgzEkgyxhwBEJHlwEygxnAHZgMP17XR8PBwEx0d7WKZSimlALZu3ZphjImoq50r4R4JpFRZTgVG1dRQRLoDMcCaujYaHR3Nli1bXHh6pZRSlUTkmCvtXBlzlxrW1TaWMwt4zxhTUUtR80Vki4hsSU9Pd6U+pZRSDeBKuKcCXassRwEna2k7C1hW24aMMS8bY+KNMfEREXV+qlBKKdVAroT7ZqC3iMSIiDf2AE+o3khE+gIhwI+NW6JSSqn6qjPcjTFWYAGwCtgHvGOM2Ssii0RkRpWms4HlRucQVkqpFufKF6oYY1YCK6utW1ht+ZHGK0sppdSF0DNUlVLKDWm4K6WUG9JwV0qpRrQrNYdNyVktXYaGu1KqeRSXVZCSVdTSZTS5+97dyR/e3dHSZWi4K9UWHM8sorzC1qI1/PPLA0x/9jvKrC1bR0Ok5Zfw9OqDddaelJbPwTMFpGQVk1VYdtZ9uUXlzfraNdyVaqAyq42CUmtLl1GnL/ac4vJ/ruX/rdzfYjUYY/hiz2nyS63sPZnb6Nt/c/1RXlibVK960vJKXP79LduYwrNfH+KDbannbbdy92nn7Z2pZ0+O+88vD3DFk980W8DXOStkU4mPjzcNmlvm8wfg9O7GL0gpFxSXWymvMAT6enEoLZ/Csgpio4IQxywducXl+Hlb8LY0X7/JYJzPX11eSTn7TuUBICIM6xaMp0fj1lZWYTvn9RoMGQVleHt6EOjrSVFZBbtP2EO9e6g/nYP86v085TYbnh5yzmstq7CxPSUbQYjvHoKH1LwvjGPWFEE4kVNESnYxAF1D/IkMtteTWVhKanYxgyKDsFTZzp6TuRSUWvH1tBDb9affd0GpleyiMkIDvAnw9mRXag4IFJVVEBXsR1SIPwBWm41tx3MIC/CmZ0Q76DQYpj1W730AICJbjTHxdbVz6Th3pdqaEqt9eiRfT4tzndVmI/FUPtYKGz0j2pHp+NhdarXh62mhqMzKvtN5eIgQFeJHl/MEWHG5laKyCkIDvM8KK4Nxbq/2x1aw71QefTq2x8tT2JmSS7dQfzoF+p7VzmqzkZRWgI+nhZjwAPadziMtv/SsuqoG3jn7oLyC03klBPl50d7XE4sjWG3GcDqvhI6BvpSU20O7T8f2hPp7Ox9bUGrlcHqBfR96WQjy8wLA00PIL7XSuZbXVlZhI6uwjIj2PmeFa0ZBKUnpBfh4etAlyI+OVV7rqdxijLG/ltySckL8vM/ZbkGpleSMAipsENs1iOyicny9LPh5eZCSXUR7X0/a+XhyLLOIsgobecXlBPl5kV9ixd/HQkGpFX9vC0VlFRxOK6CswlBeYaO43P4+OZFTjK+XhZLyCrqH+pOWX3rWp4IzeaXYjKFTkO85tTWV1hfuDfxrp9yPzWYot9nwOU8QNsSZvBKmP/sdxeUVLPnVcMb1sc+D9MhHu3k76TgBPp7kn7Di7elBmdXG40OHcGN8V576NJE3jh1lTK9wvj2YzsobLyPAx8Lr3ydzMreE07kl5JeUE+DjSeKpPIyB9+8Yw/DuIc7nXrnrFL99extvzRvJiOhQ3t+WSpcgP0bGhBLgY//vevdbW1hVfIZbo2KIiQjgL0l7sJwR3pw+krG9wwH7sMM9y3fweekpPrhzDEFRwTz90o+kZhez7neXY4C3fjzGv9Ycon/nQJ64IdbZe80tKmfXiRzuWbGDjIKfxo2D/Lx4+zej2HI0m4cT9vLoyIFkFJTyr+QkZoZ14dlZQ51tn/sskTdSjvKPawazMGEPJYU2hnYLJirEny1Hs/jxfyYCUFRmZXdqLkXlFYzrHcHNr2xg46ksIov9ePLGWEb3CKO4rIKrn/yGgGBPQvy92Hw0myUThzF1UGdyi8qZ/tjXjO0dzneHMri2SyQh/t4cOJ3PM7Pi8Pf2ZEdKDtf/ez0WD6HUauPty0bx69c3MX9cD347oRe3/+t78jOtTOgbwbvFqVg8hBs7dyXIz4sl6w5zRb8OrClL4725l/Dn93dxLLOIgV0C6dDeh2HdQ5gZF8myjcdJSisgrJ03f5zSj//7NJF1B9PYPOdK0vNLuf6FH4jpGsDbvxndqO/V82l94a7atOe+PsT249m88ut4Hk7Yy/dJGaz9w+V4eNh7eeUVNryqDRGk5ZWwbFMK4/qEM7RbSE2bdaqwGe5ZsYOisgq6hvox783N/GV6fyoM/HfDceaNjWFIVBC/X76DBRN68cb6o2w4ksnPh0by0Y4TXNm/I4uvG8KYx77mua8PkZxRyNHMQmLCA+gY6Ev3MH9yi8u5fXxPXvs+mYQdJ84K9492nADg8S8OMCgykGWb7LNtx0YF8eGdl7IjNYdVe8/g6SF8czCdlOwiOgf5EujrxR1Lt7J8/miiwwL443u7+Gz3Ke6b3IchUcEAzBsbw/y3tvLF3tPsPZnHv785zIjoEHam5DD92e9Ydfc4th/PZsGy7VTYDJHBfqy6exwpWUUczSzkX2uSeHr1IY5mFgLw8Y4TFJXZe65r96c5970xhs/3nGZsr3BuHNGV0ABv5r+1hZ8N7oynh/DJzpOcdPR0pz/7HafzSgCICQ8gOaOQ28b1YHXiGW5+bSMP/WwASWkFnMotYcX80cR2DWb2Kxu4Z8VOBkcFszMlh8KyCuaP64Ex8OG2ExQ6arr1zS387ZpB3PvODjq092HZ/NFc8eQ6Fn+xH6vNcEnPMAJ8PFnyq+H8btk23t2aysiYUEL8vfh63xnKK2yIwJr9aYT4ezG0Wwgf//ZSANr7ep31vrlvSt+zluO6BvH+tlQeeH83CTtPYrXZePLGWFfe4o1Gw121GvtO5fHMVwexGfjzh7t5d2sqxsD+0/kM6BLInhO5/Oq1jdwU35UHp/cHYNmm4zycsJcyq40v9p5m5V1jWXsgjTKrjfF9OuDnbcEYw+H0QrqF+vPnD3ez/nAmj18/hKmDOnH38h088on9ujSTB3Tkvsl98fXyoHeH9vTr1J79p/PYeCSLtfvTyCgo44b4KIL8vfjV6O689O0RAP7zPyOY0LfDOa/nWGYhn+46xfzxPVmzP42rBndm3YF0osP82X0il90ncpk3NobOQb78/bN9fLj9BP9Zn0x4O29+fUk0T60+SEpWEdcOi2LBFb24/t/r+eWrG6mwGQpLrfx5ej9+c1kP5/NN7N+R7mH+PPf1IY5lFnFNXBeevimOw+kFTH/2e/7+WSLbj+fQu0M77pnUh9ExYQT5e9G3U3sACksrePqrgwAM7BLItuP2LwxHRoey6WgWG49k4edtIa+4nNTsYu6a2BuAKwd0ZMOfJxLRzsc57v5DUgabj2aRUVDK878YSlZhGf+7ch9TB3bigWn9uHNCL257awsPJ+wFYEZsF0b1CAPg8euGMOnpb/n2YDpHMwrxtngwKDKIif078GXiGXp3aMe8sTE8+OFuJj65DoD/zhtF97AA4ruHsDE5C08Pcf5R7dupPZ/ddRmf7DzJiOhQfjySyaq9ZwB4+qZYHklIZEK/Dlg85JxQr83w7qEAfLjjBFMHduIPk/vQPSzApcc2Fg13dV4FpVYCvC1ILV9S1ceXe0+z+Iv9vD5nRK1v9LyScp5ZfYgfj2QyrFsw//j5YMA+/vvwx3sJ8rOHzTtbUgn29yKnqJzvDqXTzseTOf/ZTGGplZe+PUKXYD8KSq08seoA4/pEMKxbMM98dYgXvznMk18ewGYgxN+Lj387lvWHM3jgg934eHpQarXx+4m9uTHePsv1a7fE89+Nxyktr2DupTHOTwgDugQCcEmPMFbuPs0f3t1JlyBfxjuGcOaOjeG/G44xI65LjcEOMDMukpW7TzPl6W8pKLXy1o9HKauw8cQNsTz88V58vDz409R+WDyEFZtTuO+9nQC8cnM80eH+PLX6IKVWG+P7hBMZ7Mdb80bx26Xb6N+5PbeMiT7nU4rFQ/ifMdE88kkiXhbhD5P7IiL06tCe/xkbzUvr7H+M3v7NKMb0DD+n3jmXRvPqd0fwtAjPzhrKlU/Zg/PhGQO49sX1/PbtbeQWlwP2sfXJAzo6H9uhvX2suX/nQCKD/fjT+7uwGbhtfA+uGtLFuT8q32tBfl7839xRbDueTXg7H3qE//R+6dWhHeHtvNmcnMWJnGL6dwnEx9PClIGdWLX3DPdN7suALoGM6hHG2v1pBPp5OYerruzfkY3JWcR2Dcbf+6f487J4cO2wKMD+HQFA347tuSYukjE9w51DYq4a0CWQT383lu5h/i7/QWhsre9oGdVsDqcXMO2Z75jYvwOPXTuEIP/6vUkLSq18uP0ESWfyuSG+K7e8vonMwjLG94ngjf8Zcc4fjPIKG3Pf2Mz6w5lEhfiRklXEjw9OZPvxHB76aDcZBWUsvm4wo3uEMfvlDfxpWj9eXHuY8PbeeIiwMyWHFbddwkMf7WHrsWwAJvbrwIu/GoYxMHbxGjIKyogM9mPRzIEseHs7Y3uHs+dELsH+3gyJDGJQVBA3j+7u8mtMSivgyqfW0SMigNdvGUF0lRDKKSojyM+r1j+MpdYKRv3v13hZPLisdzgfbDtBlyBfvv/TFZRV2I8M8XQMMa3ae5rb3trKHyb14XcTe2OMYezitZzOK2H7wkkEuhgghaVWxj/xDdcOi+TPjk83YP+jeuWT64jtGswrv679QIy1+9MAmNCvAzcu+ZGc4jK+vGc8dy7dypr9adw7qQ8l5TZCArxr3Y85RWU8+/UhDp7J55Vfx58Vsq66/a2t7D6RS05RGdcPj+LRmYNcelxyRiET/vkNCyb0Omcopaol6w4zJCqoxj9yLc3Vo2U03Nuo3KJyXlyXxO8n9q71P9fiL/bz0rrDWDyEEH9v7p/Sl+uGRTl7r1UZYygotTp7KSdzirn2xfWczivB4iFU2AxeFuHG+K4s3XicF385jOmD7cdMHM0oZGHCXk5kF3E4vZDHrx9CfPcQrnhyHQsm9OKtDcfoEuzHwqsGcEnPMOfziQh//zSR135Ixhh46Gf9ufWyHuSXlLPhSBaRwX7069TeWe8La5P455cH+O+8UVzaK5ynVx/k2a8PAbX3Vl2xKTmLvp3aO48IqY/kjEICfT0J8ffmf1fuo3/nQK4bHlVj25M5xXQO8nX+sXj1uyMczypikYvBVqm4rAIfT49zfo95JeX4e1mcf1DqklFQirXCfgRIQamV0vIKwtr51KuWhnr1uyP8/bN9ADx1Y6yz1+2KtfvTGNYtpN6dlYuFHgqpzmvppmO8tO4Ig7oEcXVsl3Put9kMH20/wfg+Edw7qS9//XgP97+3i5SsIu6c0ItnvjpEp0AfxvWJwMfLwr0rdrAxOYv+nQOZ0DeCbw6kU1BqZfn80fae8qeJjO0Vzi9HdWNTchb//uYw0wd3przCxl3Lt5OcUciwbiHMuTTGOSQytFswzztOTFl66ygGRQY566sMuHF9Inj1+2Qig/24+RJ7T7G9rxeTqgwJVLrz8p5cMzTSeVTIbeN78O6WFHp2aHdBPbSRMaENfmxMlZ7+Q1cNOG/bLsFnH1p5a5Xx9Prw86756CJXe/+VwqsEeTsf+6GEzaXqPo/tGlyvx07oV/MwmbvRcG/FSsorOJlTTI+IdvV6nDH24AZ7r7My3NcdTOe9ralkFZbSr1Mgp3JL+PP0/gyOCuKDO8Zw37s7eX5tEhuOZLHp6NkTI/l6eXDr2Bh2n8jlpW+PYIzhtTkjGO34EqzqR/3ZI7ux6NNEDpzO57Pdp9iVmsu/fzmMaYPPPvr5umFRbD+ew4zYLmcFe1UjY0KJ7RrMnZf3rPOQSBFxBjuAv7cnn989Dl8vPVG7tRnQORB/bwsWDyGmmb+obC003FuxxV/sZ+nG46x/4ArC2/mQlJbPoTMFDI8OcX6BVZN9p+zzX3hZhI3JmYB9aOT2t7YS4ONJoJ8nPyRl0t7H09kD9vAQHpk5kI3J9mB/+OoBjOsTwdZj2aTnlzJ5QEd6d7QfVZFbXE5WYdlZvdKqZsZ14X9X7uPvnyWy/nAm1w6LPCfYAa4ZGsm+U3n8dkKvWl+Lr5fFeXhaQzRkKEW1PE+Lh/OL6pqGCZWGe6tVVGblva2plFltJOw4SYXN8I+V9jHImXFnn1BS3YfbU/H0EOaMieaV75JJzy/lvnd34mURPv3dWDoG+rBmfxrenh7OIwfA/rH9zbkjOXgm3zle3rOGTw1Bfl7nDc2wdj5M6NeB1YlniAz245EZA2ts187H03m0jFLV/Wt27e9xpeHean268xT5JVZC/L1Ytuk4Z/JKGNsrnAAfC2v2pVFqrWBzcjYr99jb/XFKX7qG+vPFntO8sf4okwd2ZMrATrzyXTK3/3crW49l89SNsc7Toyf2P3fMGuyHofXqUL9hoJr8anR31h1M58kbY+s91qsUaI+9LhrurcjJnGI6tPfB0+LB0o3H6NOxHbNG2MevAf48vT+n84pZtfcM72xJZdEne/H1smAMfH8onWHdQvjmYDqxUUH8v2uH4OdlwcfTg63HsrkxPoqfD41sttcyvk8Eex6Zgrenjncr1RT0f1Yr8ePhTC57fC1v/niMlKwidqbmcmN8V2bGdcHLIkwf3IkBXQLtJ1x4W3g0YS8eIqy+Zzyf/G4sUSH+HEzL55ZLovm/eaMI8vPC29OD8Y4TfBbNHNQoJyrVhwa7Uk3HpZ67iEwFngUswKvGmHNm7xKRG4FHAAPsNMb8ohHrbNPS80u5a7l9vo8v9552Ht1xed8OhLXz4cM7L6VbmH1qUV8vCxP6deDTXaeYd1mMc5jlk9+NrXHbS341HJsxLh/brJRqHeoMdxGxAC8Ak4BUYLOIJBhjEqu06Q08CFxqjMkWkbZxIGkzqJzIKq+4nEkDOrJmfxqeFqFLkC89I+xHo1Q/TPCXo7pzKreEO8b3rHP7Hh6CRy1zgSulWi9XumsjgSRjzBFjTBmwHJhZrc1vgBeMMdkAxpi0xi2z7XphbRLfJ2XwyIyBzB/Xgwqb4YekTC7rHVHrMMolPcN4/44xBPufO6+1UqptcGVYJhJIqbKcCoyq1qYPgIj8gH3o5hFjzBfVNyQi84H5AN26dWtIvW7vSHoB897cQqCfF0WlVg6lFTAzrguzRnSlwmYI8vMit7jcOce4UkrVxJWee03dw+oT0ngCvYHLgdnAqyJyzjnBxpiXjTHxxpj4iAgNJ7DPif3WhmOAfSKpu5ZvJ6uwDH8vCxHtfXh0xkAWXzcEEfskUuP6RCACYxxzrCilVE1c6bmnAl2rLEcBJ2tos8EYUw4ki8gB7GG/uVGqdFNFZVb++tEeCkqtxHcPYcXmFPacyOOlm4czZWCnGh9z3+Q+TB/UiZAAHXJRStXOlXDfDPQWkRjgBDALqH4kzEfYe+xviEg49mGaI41ZqDt6f9sJ8kqs+Hp5cMvrm0jLL2Xe2Jhagx2ge1hAs0/6r5RqfeocljHGWIEFwCpgH/COMWaviCwSkRmOZquATBFJBNYC9xtjMpuqaHdgsxne+CGZIVFB/OVnA0jLL2XKwI5nzbGtlFIN5dJx7saYlcDKausWVrltgHsdP6oam83w5w93ExMewG3je5JbVM5fP97D4fRCnr4plpmxkXQL9WdUTCgWPaVaKdUI9MyVZvDyd0dYvjmFZZuOA/D7FdtZufsU907qw8zYSDw8hPF9Is6apEsppS6Ezi3TxHan5vLEqgME+npyNLOI1Owi1idlMndsjPMCwkop1di0597Elqw7jL+3hWdmxQHw6nfJlFXYnJeLU0qppqDh3oRSs4v4fM8pfjGyG2N6huNlEZZtOo7FQxgR3fBLsymlVF003JvQ698fRUT49ZhofL0sDOwSRKnVxpCooGa93qRSqu3RcG8C+SXl3P7WVl7/IZkZsV2c1+2M7x4CwCU9dEhGKdW0NNybwPJNKXyx9zT3T+nLY9f9dJm4EY4rto/tFd5SpSml2ggdG2gCG5MziQkPOOfCzpP6d2T5/NGMitHxdqVU09KeeyOz2QybkrMYWcMXph4ewugeYc1+xSOlVNuj4d7I9p/OJ6/Eyqge2jtXSrUcDfdGtinZPqXOKP3SVCnVgnTMvZF8vvsU/1qTRKm1gshgP+cRMkop1RI03BvBjpQc7l6xgwAfT7IKy5g9smvdD1JKqSak4X6BcorKuO2tLXQI9OGjOy+lxGojxN+rpctSSrVxGu4XaNEniWQUlPHxby8lrJ1PS5ejlFKAhnuDvfhNEhuPZLHuYDp3XdGLQZFBLV2SUko56dEyDXDgdD6Pf3GAY5mFzB7ZlQVX6NS9SqmLi/bcG+DdLSl4egjv3zFGh2KUUhcl7bnXU3mFjY92nGBi/w4a7Eqpi5aGez2t3Z9GRkEZNwzXwx2VUhcvl8JdRKaKyAERSRKRB2q4f46IpIvIDsfPrY1f6sXh012nCA3wZnzfiJYuRSmlalXnmLuIWIAXgElAKrBZRBKMMYnVmq4wxixoghovGqXWCtbsT+NngzvjZdEPPUqpi5crCTUSSDLGHDHGlAHLgZlNW9bFaf3hTApKrUwZ1LGlS1FKqfNyJdwjgZQqy6mOddVdJyK7ROQ9EXHLAekv954mwNvCmJ56sQ2l1MXNlXCvafJxU235EyDaGDME+Ap4s8YNicwXkS0isiU9Pb1+lbYwm82wOvEMl/frgK+XpaXLUUqp83Il3FOBqj3xKOBk1QbGmExjTKlj8RVgeE0bMsa8bIyJN8bER0S0ri8k95zMJaOgjEn9dUhGKXXxcyXcNwO9RSRGRLyBWUBC1QYi0rnK4gxgX+OVeHH49qD9k8bY3joko5S6+NV5tIwxxioiC4BVgAV43RizV0QWAVuMMQnAXSIyA7ACWcCcJqy5RXx7KINBkYGE64lLSqlWwKXpB4wxK4GV1dYtrHL7QeDBxi3t4pFfUs62Y9n8ZlyPli5FKaVcogdr18FaYePzPaex2gzjereu7wmUUm2XThx2HoWlVib88xvS8ksJ9PVkePeQli5JKaVcouF+HpuSs0jLL+W+yX24Zmgk3p76QUcp1TpouJ/H90kZ+Hh6cOtlPfTYdqVUq6Jd0fP4ISmD+OgQDXalVKuj4V6L9PxS9p/O59Jeely7Uqr10XCvxfrDGQCM1XBXSrVCGu61WHcwnSA/LwZ20QtfK6VaHw33GpRaK1ideIZJAzpi8ahp3jSllLq4abjX4IekDPJLrPxscOe6Gyul1EVIw70Gn+46RaCvp36ZqpRqtTTcqykprxyS6aQnLSmlWi1Nr2re2ZJCfomV64bVdLEppZRqHTTcqygpr+CFtUmMjA7lkp5hLV2OUko1mIZ7Fcs3HedMXin3TOqDiB4lo5RqvTTcq0jYeZLBkUHaa1dKtXoa7g55JeXsTM3l8r46Z7tSqvXTcHf48XAmFTajhz8qpdyChrvD94cy8Pe2MKybXpBDKdX6abg7fJ+UwaiYUD22XSnlFjTJgOOZRSRnFOqQjFLKbbgU7iIyVUQOiEiSiDxwnnbXi4gRkfjGK7HpfbTjBABTB3Vq4UqUUqpx1BnuImIBXgCmAQOA2SIyoIZ27YG7gI2NXWRTMsbwwbZURvcIJSrEv6XLUUqpRuFKz30kkGSMOWKMKQOWAzNraPc34HGgpBHra3LbjmdzNLOI64ZFtXQpSinVaFwJ90ggpcpyqmOdk4gMBboaYz5txNqaxYfbT+DnZWGaTu+rlHIjroR7TefhG+edIh7A08Af6tyQyHwR2SIiW9LT012vsgltSs5iVI9Q2vl4tnQpSinVaFwJ91Sga5XlKOBkleX2wCDgGxE5CowGEmr6UtUY87IxJt4YEx8R0fJnghaUWjmUVkBc1+CWLkUppRqVK+G+GegtIjEi4g3MAhIq7zTG5Bpjwo0x0caYaGADMMMYs6VJKm5Eu1NzMQZiNdyVUm6mznA3xliBBcAqYB/wjjFmr4gsEpEZTV1gU9qZmgNAbJSGu1LKvbg00GyMWQmsrLZuYS1tL7/wsprHrtQcuoX6Exrg3dKlKKVUo2rTZ6juTMnVIRmllFtqs+Genl/KiZxiYqOCWroUpZRqdG023LcczQJgaDftuSul3E+bDfcNRzLx87IwRL9MVUq5oTYc7lnER4fgZWmzu0Ap5cbaZLJlFpRy4Ew+o3votVKVUu6pTYb7pmT7eLuGu1LKXbXJcP9pvF2PlFFKuac2Ge7bjucwtFuwjrcrpdxWm0u38gobB07nMzhSe+1KKffV5sL90JkCyipsDNRwV0q5sTYX7ntO5gIwqEtgC1eilFJNp82Fe+LJPAK8LUSHBbR0KUop1WTaXLjvOZHLgC6BeHjUdIEppZRyD20q3CtshsRTeQzsouPtSin31qbC/WhmIUVlFQzU8XallJtrU+FeeWaqzgSplHJ3bSrcv0/KoFOgLz0j2rV0KUop1aTaTLjbbIb1SRlc2iscEf0yVSnl3tpMuCeeyiO7qJyxvXWyMKWU+3Mp3EVkqogcEJEkEXmghvtvF5HdIrJDRL4XkQGNX+qF+SEpA4BLe4a3cCVKKdX06gx3EbEALwDTgAHA7BrC+21jzGBjTBzwOPBUo1d6gdYfzqR3h3Z0CPRt6VKUUqrJudJzHwkkGWOOGGPKgOXAzKoNjDF5VRYDANN4JTaOvSfziOuqR8kopdoGTxfaRAJiNM+9AAAbWElEQVQpVZZTgVHVG4nIb4F7AW/gikaprpFkFJSSUVBK307tW7oUpZRqFq703Gs6tOScnrkx5gVjTE/gT8BDNW5IZL6IbBGRLenp6fWr9AIcOJ0PQL9OevKSUqptcCXcU4GuVZajgJPnab8cuKamO4wxLxtj4o0x8REREa5XeYH2O8Jde+5KqbbClXDfDPQWkRgR8QZmAQlVG4hI7yqLPwMONV6JF+7A6TzCAryJaO/T0qUopVSzqHPM3RhjFZEFwCrAArxujNkrIouALcaYBGCBiFwJlAPZwC1NWXR9HTidT7/O2mtXSrUdrnyhijFmJbCy2rqFVW7/vpHrajQVNsOBM/n8YmT3li5FKaWajdufoXo8q4iSchv9dLxdKdWGuH24J560H4KvwzJKqbbE7cN9+/FsfDw99DBIpVSb4vbhvu14NoMjg/D2dPuXqpRSTm6deKXWCvaczNOLcyil2hy3Dvd9p/Ips9oY2i2kpUtRSqlm5dbhvu1YNgDDNNyVUm2MW4f79pQcOgf50ilIp/lVSrUtbh3ue0/kMiQqqKXLUEqpZue24V5mtXEsq4jeHfT4dqVU2+O24X48q5AKm6Fnh4CWLkUppZqd24Z7UlohAD0j2rVwJUop1fzcNtwPpxcA0EPDXSnVBrltuB9JL6RjoA/tfFya+FIppdyK24b74fQCHZJRSrVZbhnuxhgNd6VUm+aW4Z5eUEp+iZWeEXqkjFKqbXLLcD+S7jhSpoP23JVSbZNbhnvlkTI6LKOUaqvcM9zTCvHzstApUOeUUUq1Te4Z7ukF9IgIwMNDWroUpZRqES6Fu4hMFZEDIpIkIg/UcP+9IpIoIrtE5GsR6d74pbruSIYeKaOUatvqDHcRsQAvANOAAcBsERlQrdl2IN4YMwR4D3i8sQt1VUl5BanZxRruSqk2zZWe+0ggyRhzxBhTBiwHZlZtYIxZa4wpcixuAKIat0zXJWcUYgw6YZhSqk1zJdwjgZQqy6mOdbWZB3x+IUVdCD1SRimlwJWJV2r6VtLU2FDkV0A8ML6W++cD8wG6devmYon1cyS9EBGICdeeu1Kq7XKl554KdK2yHAWcrN5IRK4E/gLMMMaU1rQhY8zLxph4Y0x8REREQ+qt0+H0AiKD/fD1sjTJ9pVSqjVwJdw3A71FJEZEvIFZQELVBiIyFHgJe7CnNX6ZrrMfBqlDMkqptq3OcDfGWIEFwCpgH/COMWaviCwSkRmOZk8A7YB3RWSHiCTUsrkmVWEzJKUV0EenHVBKtXEuTXZujFkJrKy2bmGV21c2cl0NcjyriJJyG3066XVTlVJtm1udoXrgdD4A/TTclVJtnFuF+8Ez+YhALx2WUUq1cW4V7gfO5NMt1B9/b720nlKqbXOvcD+dT5+OOiSjlFJuE+6l1gqSMwp1vF0ppXCjcD+SXkiFzWjPXSmlcKNwP3jGfqSMhrtSSrlRuB9OL8RDIDrcv6VLUUqpFuc24Z6cUUhUiD8+njqnjFJKuVG4F+hMkEop5eAW4W6MITm9UMNdKaUc3CLc0/NLKSyroEeEhrtSSoGbhPuRjEJAL9ChlFKV3CLckzXclVLqLG4T7t6eHnQJ8mvpUpRS6qLgFuF+JL2QmLAAPDxqutyrUkq1PW4R7noYpFJKna3Vh3uFzZCSVUx3PTNVKaWcWn24n84roazCRvdQ7bkrpVSlVh/uxzOLAOgWqj13pZSq1PrDPct+GGT3MA13pZSq5FK4i8hUETkgIkki8kAN948TkW0iYhWR6xu/zNodyyzC00PoHOTbnE+rlFIXtTrDXUQswAvANGAAMFtEBlRrdhyYA7zd2AXW5VhWEZEhfnhaWv2HEKWUajSuXEl6JJBkjDkCICLLgZlAYmUDY8xRx322JqjxvFKyinS8XSmlqnGluxsJpFRZTnWsqzcRmS8iW0RkS3p6ekM2cY5jmUU63q6UUtW40nOv6bRP05AnM8a8DLwMEB8f36BtVJVbVE5ucbkeBqmaRXl5OampqZSUlLR0KaoN8PX1JSoqCi8vrwY93pVwTwW6VlmOAk426Nka2THHkTJddVhGNYPU1FTat29PdHQ0IjrVhWo6xhgyMzNJTU0lJiamQdtwZVhmM9BbRGJExBuYBSQ06Nka2fEs+zHuOiyjmkNJSQlhYWEa7KrJiQhhYWEX9CmxznA3xliBBcAqYB/wjjFmr4gsEpEZjkJGiEgqcAPwkojsbXBF9ZCSVQxoz101Hw121Vwu9L3m0vGDxpiVxpg+xpiexph/ONYtNMYkOG5vNsZEGWMCjDFhxpiBF1SVi1Kziwj296KdjyujS0q1fiLCzTff7Fy2Wq1ERERw1VVXAfDGG2+wYMECAB555BH8/f1JS0tztm/Xrp3ztsViIS4uzvnz2GOPXVBtO3bsYOXKlc7lhISEC95mpWeeeYaioqJG2VZb0aoPDj+RU0xUiM7hrtqOgIAA9uzZQ3Gx/VPr6tWriYys/eC18PBwnnzyyRrv8/PzY8eOHc6fBx445/zEeqke7jNmzLjgbVbScK+/Vh3uqdnFRAZruKu2Zdq0aXz22WcALFu2jNmzZ9fadu7cuaxYsYKsrKwGPdfhw4eZOnUqw4cP57LLLmP//v0AvPvuuwwaNIjY2FjGjRtHWVkZCxcuZMWKFcTFxbFixYqzPkXMmTOHO+64gwkTJtCjRw/WrVvH3Llz6d+/P3PmzHE+3x133EF8fDwDBw7k4YcfBuC5557j5MmTTJgwgQkTJgDw5ZdfcskllzBs2DBuuOEGCgoKGvT63FmrHc8wxnAiu5jxfSJauhTVBj36yV4ST+Y16jYHdAnk4avrHtGcNWsWixYt4qqrrmLXrl3MnTuX7777rsa27dq1Y+7cuTz77LM8+uijZ91XXFxMXFycc/nBBx/kpptuOqvN/PnzWbJkCb1792bjxo3ceeedrFmzhkWLFrFq1SoiIyPJycnB29ubRYsWsWXLFp5//nnAPkRUVXZ2NmvWrCEhIYGrr76aH374gVdffZURI0awY8cO4uLi+Mc//kFoaCgVFRVMnDiRXbt2cdddd/HUU0+xdu1awsPDycjI4O9//ztfffUVAQEBLF68mKeeeoqFCxe6spvbjFYb7lmFZRSXV2jPXbU5Q4YM4ejRoyxbtozp06fX2f6uu+4iLi6OP/zhD2etrxyWqU1BQQHr16/nhhtucK4rLS0F4NJLL2XOnDnceOONXHvttS7VffXVVyMiDB48mI4dOzJ48GAABg4cyNGjR4mLi+Odd97h5Zdfxmq1curUKRITExkyZMhZ29mwYQOJiYlceumlAJSVlXHJJZe4VENb0mrD/USOfcxRx9xVS3Clh92UZsyYwX333cc333xDZmbmedsGBwfzi1/8ghdffLFez2Gz2QgODq7xD8CSJUvYuHEjn332GXFxcef9I1HJx8cHAA8PD+ftymWr1UpycjL//Oc/2bx5MyEhIcyZM6fGQwGNMUyaNIlly5bV6/W0Na12zD012x7ukRruqg2aO3cuCxcudPZ+63Lvvffy0ksvYbVaXX6OwMBAYmJiePfddwF7qO7cuROwj8WPGjWKRYsWER4eTkpKCu3btyc/P7/+L8YhLy+PgIAAgoKCOHPmDJ9//rnzvqrbHj16ND/88ANJSUkAFBUVcfDgwQY/r7tqteF+Iruy567HuKu2Jyoqit///vcutw8PD+fnP/+5c1gFfhpzr/yp6ciWpUuX8tprrxEbG8vAgQP5+OOPAbj//vsZPHgwgwYNYty4ccTGxjJhwgQSExOdX6jWV2xsLEOHDmXgwIHMnTvXOewC9rH/adOmMWHCBCIiInjjjTeYPXs2Q4YMYfTo0c4vetVPxJgLnuKlQeLj482WLVsa/PhHEvby/tZUdj86pRGrUqp2+/bto3///i1dhmpDanrPichWY0x8XY9ttT331OwiHZJRSqlatOJw1xOYlFKqNq023E/oCUxKKVWrVhnuucXl5Jda9ctUpZSqRasM99Rs+xwTOuaulFI1a5Xh/tNhkBruSilVk1YZ7s4TmHTMXbUhmZmZzmPSO3XqRGRkpHO5rKys0Z7nq6++QkR48803nes2b96MiPDMM8+c97Hbtm3jiy++uKDnj4qKIicnx+X1VSUmJjqPlz969OgF1eGqI0eOsHz5cufyxo0bueeee5rluc+nVYb7iZxi/LwshAZ4t3QpSjWbsLAw5/S8t99+O/fcc49z2dvb/n/BGIPNZrvg5xo8ePBZgbV8+XJiY2PrfFxjhPuF+OCDD7j++uvZvn070dHRdbZvjP1VPdxHjRrF008/fUHbbAytM9yzi4kM8dOr4igFJCUlMWjQIG6//XaGDRtGSkoKwcHBzvuXL1/OrbfeCsCZM2e49tpriY+PZ+TIkWzYsKHGbfbo0YO8vDwyMjIwxrB69WqmTPnphMGxY8c655M5ffo0vXr1ori4mEWLFrF06VLi4uJ47733eOihh87q7ffr14/U1FTAPpHY8OHDGThwIK+++mq9X++8efMYOHAg06ZNo6SkhISEBJ5//nmWLFnClVdeCcDjjz/OoEGDGDRoEP/617/Ou7/uv/9+hg0bxpQpU9i4cSPjx4+nR48ezjnqDx8+zGWXXcbQoUMZPnw4GzduBOCBBx5g7dq1xMXF8dxzz/HVV19xzTXXAJCRkcGMGTMYMmQIY8aMYc+ePQA89NBDzJs3z/kcL7zwgsuv31WtcuKw1JwiHW9XLevzB+D07sbdZqfBMK1hVy5KTEzkP//5D0uWLDnv/DF33XUXf/zjHxk9ejRHjx7lqquucgZOdddddx3vvfce/fv3Z9SoUXh5eZ23Bj8/PxYuXMiePXucgX6+CcXefPNNQkNDKSoqIj4+nuuuu46QkBAXXi0cOHCAZcuWMXjwYK699lo++ugjZs2axaZNmwgPD+fuu+9m06ZNLF26lE2bNlFRUcHIkSMZP348/v7+5+yv3NxcJk+ezBNPPMHVV1/NI488wtdff83OnTu57bbbmD59Op07d2b16tX4+vqyf/9+brnlFjZu3Mhjjz3G888/z0cffQTYh7Uq/fWvf2XUqFEkJCTw5ZdfMmfOHCrPzD948CBff/01OTk59O/fn9tvvx2LxeLS63dFqwz3E9nFxEYF191QqTaiZ8+ejBgxos52X331FQcOHHAuZ2dnU1xcjJ/fuZ2lm266iZtvvpk+ffowe/Zs1qxZ06g1P/300yQkJACQmprK4cOHiY+v86x6AHr16uWcNG348OE1jq9/9913XHfddfj72w+Zvuaaa/j++++ZPHnyOfvLz8+PSZMmAfYhqaCgIDw9PRk8eLBz26WlpSxYsICdO3fi6enJ4cOH66zz+++/d15YZfLkycyZM4fCwkIArrrqKry9venQoQOhoaGkp6fTqVMnl16/K1pduBeWWskuKtfDIFXLamAPu6kEBAQ4b3t4eFB1zqiq0+YaY9i0aZNzjP58IiMjMcawbt06XnzxxbPC3dPT0zlWXdO0vDW1q9r2q6++4ttvv2XDhg34+fkxduzY826nuqpTBlsslho/rZxv3qyq+ws4a39UnZK4cjpigCeffJKuXbvy3//+l/Ly8rOuR1ub6jVUXXblNVwIl8bcRWSqiBwQkSQROWfqOBHxEZEVjvs3ikh0o1ZZxU/zuOsJTErVxMPDg5CQEA4dOoTNZuPDDz903nfllVeeNb5b1zzsf/vb31i8eDEeHmdHRXR0NFu3bgXgvffec66vPu1v1XabNm0iJSUFgNzcXEJDQ/Hz82Pv3r1s3ry5ga+2duPGjePDDz+kuLiYgoICPv74Yy677LIGby83N5fOnTs7jySqDOrzTXU8btw4li5dCtj/oEVFRZ3zh6Wp1BnuImIBXgCmAQOA2SIyoFqzeUC2MaYX8DSwuLELreQ8gUkPg1SqVosXL2bq1KlMnDiRqKgo5/oXXniBH374gSFDhjBgwABeeeWV825n7NixzJgx45z1999/P88++yxjxowhOzvbuf6KK65g586dDB06lPfee48bbriBM2fOMHToUF577TV69OgBwM9+9jOKioqIjY1l0aJFjBo1qpFe+U9GjhzJ7NmzGTFiBKNHj+aOO+5wef77mixYsIBXX32V0aNHc+zYMWfPe+jQoVRUVBAbG8tzzz131mMWLVrE+vXrGTJkCAsXLuQ///nPBb2m+qhzyl8RuQR4xBgzxbH8IIAx5v9VabPK0eZHEfEETgMR5jwbb+iUv2/9eJS/fryXTX+eSIdA33o/XqmG0il/VXNr6il/I4GUKsupjnU1tjHGWIFcIMyFbddbh0BfJg3oSHg7n7obK6VUG+XKF6o1HUxevUfuShtEZD4wH6Bbt24uPPW5pgzsxJSBjfeNslJKuSNXeu6pQNcqy1HAydraOIZlgoCs6hsyxrxsjIk3xsRHREQ0rGKllFJ1ciXcNwO9RSRGRLyBWUBCtTYJwC2O29cDa8433q5Ua6Vva9VcLvS9Vme4O8bQFwCrgH3AO8aYvSKySEQqv0Z/DQgTkSTgXuDcK+0q1cr5+vqSmZmpAa+anDGGzMxMfH0bftBIq71AtlLNrby8nNTU1HqdbKNUQ/n6+hIVFXXOtA+uHi3T6s5QVaqleHl5ERMT09JlKOWSVjkrpFJKqfPTcFdKKTek4a6UUm6oxb5QFZF04FgDHx4OZDRiOY1F66ofrav+LtbatK76uZC6uhtj6jxRqMXC/UKIyBZXvi1ublpX/Whd9Xex1qZ11U9z1KXDMkop5YY03JVSyg211nB/uaULqIXWVT9aV/1drLVpXfXT5HW1yjF3pZRS59dae+5KKaXOo9WFe13Xc23GOrqKyFoR2Scie0Xk9471j4jICRHZ4fiZ3gK1HRWR3Y7n3+JYFyoiq0XkkOPfkGauqW+VfbJDRPJE5O6W2F8i8rqIpInInirratw/Yvec4/22S0SGNXNdT4jIfsdzfygiwY710SJSXGW/LWnmumr9vYnIg479dUBEpjRzXSuq1HRURHY41jfn/qotG5r3PWaMaTU/gAU4DPQAvIGdwIAWqqUzMMxxuz1wEPs1Zh8B7mvh/XQUCK+27nHgAcftB4DFLfx7PA10b4n9BYwDhgF76to/wHTgc+wXpBkNbGzmuiYDno7bi6vUFV21XQvsrxp/b47/AzsBHyDG8f/V0lx1Vbv/SWBhC+yv2rKhWd9jra3nPhJIMsYcMcaUAcuBmS1RiDHmlDFmm+N2PvbpkKtffvBiMhN403H7TeCaFqxlInDYGNPQk9guiDHmW869mExt+2cm8H/GbgMQLCKdm6suY8yXxj7tNsAG7BfLaVa17K/azASWG2NKjTHJQBL2/7fNWpeICHAjsKwpnvt8zpMNzfoea23h7sr1XJudiEQDQ4GNjlULHB+vXm/u4Q8HA3wpIlvFfmlDgI7GmFNgf/MBHVqgrkqzOPs/XUvvL6h9/1xM77m52Ht4lWJEZLuIrBORy1qgnpp+bxfL/roMOGOMOVRlXbPvr2rZ0KzvsdYW7i5dq7U5iUg74H3gbmNMHvBvoCcQB5zC/tGwuV1qjBkGTAN+KyLjWqCGGon9al4zgHcdqy6G/XU+F8V7TkT+AliBpY5Vp4Buxpih2C+Q87aIBDZjSbX93i6K/QXM5uwORLPvrxqyodamNay74H3W2sLdleu5NhsR8cL+y1tqjPkAwBhzxhhTYYyxAa/QRB9Jz8cYc9LxbxrwoaOGM5Uf9Rz/pjV3XQ7TgG3GmDOOGlt8fznUtn9a/D0nIrcAVwG/NI5BWsewR6bj9lbsY9t9mqum8/zeLob95QlcC6yoXNfc+6umbKCZ32OtLdxduZ5rs3CM6b0G7DPGPFVlfdWxsp8De6o/tonrChCR9pW3sX8ht4ezr3N7C/Bxc9ZVxVk9qpbeX1XUtn8SgF87jmgYDeRWfrRuDiIyFfgTMMMYU1RlfYSIWBy3ewC9gSPNWFdtv7cEYJaI+IhIjKOuTc1Vl8OVwH5jTGrliubcX7VlA839HmuOb48b8wf7N8sHsf/l/UsL1jEW+0enXcAOx8904C1gt2N9AtC5mevqgf1ohZ3A3sp9BIQBXwOHHP+GtsA+8wcygaAq65p9f2H/43IKKMfea5pX2/7B/pH5Bcf7bTcQ38x1JWEfj618jy1xtL3O8fvdCWwDrm7mumr9vQF/ceyvA8C05qzLsf4N4PZqbZtzf9WWDc36HtMzVJVSyg21tmEZpZRSLtBwV0opN6ThrpRSbkjDXSml3JCGu1JKuSENd3XREJFgEbmzgY9dKY4ZE8/TZpGIXNmw6lyqYY6IdGmq7StVH3oopLpoOObh+NQYM6iG+yzGmIpmL6oeROQb7DMlbmnpWpTSnru6mDwG9HTMt/2EiFzumBf7bewndyAiHzkmRNtbZVK0yjnswx3zdu8TkVccbb4UET9HmzdE5Poq7R8VkW1in/u+n2N9hGOu7W0i8pKIHBOR8KpFiojFsa09jsfe49huPLDUUb+fiAx3TFK1VURWVTn1/BsReUZE1ju2MdKxfrz8NN/49sozjZVqkKY6S0t/9Ke+P1Sbcxu4HCgEYqqsqzyrzw/7Ke9hjuWjQLhjG1YgzrH+HeBXjttvANdXaf87x+07gVcdt58HHnTcnor9TMPqc+MPB1ZXWQ52/PsNjrMLAS9gPRDhWL4JeL1Ku1cct8dVvmbgE+yTvgG0wzGPu/7oT0N+tOeuLnabjH1e8Ep3ichO7HObd8U+R0h1ycaYHY7bW7EHfk0+qKHNWOzXCcAY8wWQXcPjjgA9RORfjrlfaprxry8wCFgt9qsBPcTZc7EvczzHt0Cg4/uCH4CnROQu7H8wrCjVQBru6mJXWHlDRC7HPinUJcaYWGA74FvDY0qr3K4APGvZdmkNbWqafvUsxphsIBZ7D/y3wKs1NBNgrzEmzvEz2Bgzuepmzt2seQy4Ffunkg2VQ0VKNYSGu7qY5GO/LFltgoBsY0yRI/hGN0EN32O/gg8iMhk45+IhjjF4D2PM+8BfsV/qDc6u/wAQISKXOB7jJSIDq2zmJsf6sdhnAcwVkZ7GmN3GmMXAFkDDXTVYbT0apZqdMSZTRH4Q+wWPPwc+q9bkC+B2EdmFPTw3NEEZjwLLROQmYB32WQfzq7WJBP4jIpWdowcd/74BLBGRYuAS4HrgOREJwv5/7RnsMxMCZIvIeiAQ+xWWAO4WkQnYP0kkcvZVl5SqFz0UUqkqRMQHqDDGWB297n8bY+Ia+Tm+QQ+ZVE1Me+5Kna0b8I6jV14G/KaF61GqQbTnrpRSbki/UFVKKTek4a6UUm5Iw10ppdyQhrtSSrkhDXellHJDGu5KKeWG/j8HZceiPZPUPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "H=20\n",
    "n_epochs = 200\n",
    "data_size = 20000\n",
    "\n",
    "def MINE(x_in, y_in):\n",
    "    \n",
    "    # shuffle and concatenate\n",
    "    y_shuffle = tf.random_shuffle(y_in)\n",
    "    x_conc = tf.concat([x_in, x_in], axis=0)\n",
    "    y_conc = tf.concat([y_in, y_shuffle], axis=0)\n",
    "    \n",
    "    # propagate the forward pass\n",
    "    layerx = layers.linear(x_conc, H)\n",
    "    layery = layers.linear(y_conc, H)\n",
    "    layer2 = tf.nn.relu(layerx + layery)\n",
    "    output = layers.linear(layer2, 1)\n",
    "    \n",
    "    # split in T_xy and T_x_y predictions\n",
    "    N_samples = tf.shape(x_in)[0]\n",
    "    T_xy = output[:N_samples]\n",
    "    T_x_y = output[N_samples:]\n",
    "    \n",
    "    # compute the negative loss (maximise loss == minimise -loss)\n",
    "    neg_loss = -(tf.reduce_mean(T_xy, axis=0) - tf.math.log(tf.reduce_mean(tf.math.exp(T_x_y))))\n",
    "    opt = tf.train.AdamOptimizer(learning_rate=0.01).minimize(neg_loss)\n",
    "\n",
    "    return neg_loss, opt\n",
    "\n",
    "# prepare the placeholders for inputs\n",
    "x_in = tf.placeholder(tf.float32, [None, 1], name='x_in')\n",
    "y_in = tf.placeholder(tf.float32, [None, 1], name='y_in')\n",
    "\n",
    "# make the loss and optimisation graphs\n",
    "neg_loss, opt = MINE(x_in, y_in)\n",
    "\n",
    "# start the session\n",
    "sess = tf.InteractiveSession()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "# train\n",
    "MIs = []\n",
    "for epoch in range(n_epochs):\n",
    "    \n",
    "    # generate the data\n",
    "    x_sample=gen_x()\n",
    "    y_sample=gen_y(x_sample)\n",
    "    \n",
    "    # perform the training step\n",
    "    feed_dict = {x_in:x_sample, y_in:y_sample}\n",
    "    _, neg_l = sess.run([opt, neg_loss], feed_dict=feed_dict)\n",
    "    \n",
    "    # save the loss\n",
    "    MIs.append(-neg_l)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(range(len(MIs)), MIs, label='MINE estimate')\n",
    "ax.plot([0, len(MIs)], [mi,mi], label='True Mutual Information')\n",
    "ax.set_xlabel('training steps')\n",
    "ax.legend(loc='best')\n",
    "fig.savefig('MINE.png')\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (tf)",
   "language": "python",
   "name": "tf"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}