MINE_implementation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "assert sys.version_info >= (3, 5)\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "import pandas as pd\n",
    "import tensorflow_probability as tfp\n",
    "np.random.seed(42)\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#leaky_relu = tf.keras.layers.LeakyReLU()\n",
    "randN_05 = keras.initializers.RandomNormal(mean=0.0, stddev=0.05, seed=None)\n",
    "dim_n = 3 # Dim 4+ needs considerable fine-tuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_A = keras.layers.Input(shape=[2*dim_n])\n",
    "input_B = keras.layers.Input(shape=[2*dim_n])\n",
    "    \n",
    "transform = keras.models.Sequential([\n",
    "layers.Dense(30, kernel_initializer=randN_05, activation=\"relu\"),\n",
    "keras.layers.Dropout(rate=0.3), # To regularize higher dimensionality\n",
    "layers.Dense(30, kernel_initializer=randN_05, activation=\"relu\"),\n",
    "keras.layers.Dropout(rate=0.3), # To regularize higher dimensionality\n",
    "layers.Dense(1, kernel_initializer=randN_05, activation=None)])\n",
    "\n",
    "output_A = transform(input_A)\n",
    "output_B = transform(input_B)\n",
    "output_C = tf.reduce_mean(output_A) - tf.math.log(tf.reduce_mean(tf.exp(output_B))) # MINE\n",
    "#output_C = tf.reduce_mean(output_A) - tf.reduce_mean(tf.exp(output_B))+1 # MINE-f\n",
    "MI_mod = keras.models.Model(inputs=[input_A, input_B], outputs=output_C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_func(inp, outp):\n",
    "    '''Calculate the loss: scaled negative estimated mutual information'''\n",
    "    return -outp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "MI_mod.compile(loss=loss_func, optimizer=keras.optimizers.Nadam(lr=0.001))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIGNAL_NOISE = 0.2\n",
    "SAMPLE_SIZE = dim_n*6000\n",
    "SIGNAL_POWER = 3\n",
    "\n",
    "x_sample = np.random.normal(0., np.sqrt(SIGNAL_POWER), [SAMPLE_SIZE, dim_n])\n",
    "y_sample = x_sample + np.random.normal(0., np.sqrt(SIGNAL_NOISE), [SAMPLE_SIZE, dim_n])\n",
    "\n",
    "x_sample1, x_sample2 = tf.split(x_sample, num_or_size_splits=2)\n",
    "y_sample1, y_sample2 = tf.split(y_sample, num_or_size_splits=2)\n",
    "joint_sample = tf.concat([x_sample1, y_sample1], axis=1)\n",
    "marg_sample = tf.concat([x_sample2, y_sample1], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 9000 samples\n",
      "Epoch 1/30\n",
      "9000/9000 [==============================] - 0s 44us/sample - loss: -0.0505\n",
      "Epoch 2/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -0.7593\n",
      "Epoch 3/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -2.0280\n",
      "Epoch 4/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -2.6971\n",
      "Epoch 5/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.0710\n",
      "Epoch 6/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.2430\n",
      "Epoch 7/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.4708\n",
      "Epoch 8/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.5317\n",
      "Epoch 9/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.5812\n",
      "Epoch 10/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.5552\n",
      "Epoch 11/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.6636\n",
      "Epoch 12/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.7999\n",
      "Epoch 13/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.5859\n",
      "Epoch 14/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.8038\n",
      "Epoch 15/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -4.0753\n",
      "Epoch 16/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.9197\n",
      "Epoch 17/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.8644\n",
      "Epoch 18/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.6294\n",
      "Epoch 19/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.7859\n",
      "Epoch 20/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.8225\n",
      "Epoch 21/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.8013\n",
      "Epoch 22/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -4.0327\n",
      "Epoch 23/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.9841\n",
      "Epoch 24/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.9979\n",
      "Epoch 25/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.9641\n",
      "Epoch 26/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -4.0970\n",
      "Epoch 27/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -4.0456\n",
      "Epoch 28/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -3.8338\n",
      "Epoch 29/30\n",
      "9000/9000 [==============================] - 0s 12us/sample - loss: -3.9486\n",
      "Epoch 30/30\n",
      "9000/9000 [==============================] - 0s 11us/sample - loss: -4.1012\n"
     ]
    }
   ],
   "source": [
    "history_mi = MI_mod.fit((joint_sample, marg_sample), x_sample, epochs=30, batch_size=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def theoretic_mutual_information(power,noise):\n",
    "    return dim_n*0.5*np.log2(1 + power/noise)\n",
    "\n",
    "\n",
    "df_mi = np.log2(np.exp(1))*np.abs(pd.DataFrame(history_mi.history)) # scaling from Nats to Bits\n",
    "df_mi = df_mi.rename(columns={'loss': 'Approximated mutual information'})\n",
    "history_range = range(len(df_mi))\n",
    "\n",
    "true_mi = [theoretic_mutual_information(SIGNAL_POWER, SIGNAL_NOISE) for x in history_range]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAE3CAYAAADrF0u3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4lFX6xvHvSSGdACGEJgQQAiSBQADpJNKsKKKrro1lV1bX+nMt6NrXbera11V3FdDFggVBV9cGodgQEJSWQAwdQgqB9Dbn90dCDCVkBpKZJHN/ritXkpl33nnmySR3zjvvnGOstYiIiIj7+Hi6ABEREW+j8BUREXEzha+IiIibKXxFRETcTOErIiLiZgpfERERN1P4ioiIuJlT4WuM6WeMWWyMOWiM2WqMmdrYhYmIiLRU9YavMcYPWAh8CLQDZgL/Mcb0aeTaREREWiRT3wxXxpg44BsgzFZvbIz5FPjWWntf45coIiLSsjhz2NnUcVlcA9ciIiLiFfyc2GYzsB+4wxjzJJAMjAOWHL2hMWYmVYelCQwMTOzWrVsDltryORwOfHx0Dpwr1DPXqWeuU89c5409S0tLy7bWRjqzbb2HnQGMMQOAZ6ka7a4CsoBSa+2v67pNTEyMTU1Nda5iASAlJYWkpCRPl9GsqGeuU89cp565zht7ZoxZba0d4sy2zox8sdb+QNVo9/AdfAXMPbnyREREvJuzbzUaYIwJNMYEG2NuBzoBcxq1MhERkRbK2QPyVwF7qXrtdzww0Vpb2mhViYiItGDOHna+A7ijkWsREQ8pLy9n165dlJSUeLqUJik8PJxNmzZ5uoxmpSX3LDAwkK5du+Lv73/S+3AqfEWkZdu1axdhYWFER0djzPHeXejd8vPzCQsL83QZzUpL7Zm1lpycHHbt2kWPHj1Oej/edR64iBxXSUkJERERCl6RehhjiIiIOOWjRApfEQFQ8Io4qSF+VxS+ItJkLFiwAGMMmzdv9mgd999/P59//vkp7ycvL4/nn3/e5ds9+OCDPP7446d8/0d7//332bhx40nfftu2bcTFHTu54Z49e7j44ovrvf3bb79Nv379SE5OPukaXJWSksJXX31V8/0LL7zAq6++6rb7r4vCV0SajDfeeIPRo0fz5ptvNtg+KyoqXL7Nww8/zIQJE075vk82fBvLqYZvXTp37sw777xT73Yvv/wyzz//PEuWHDNB4nGdzM/uaEeH73XXXcfVV199yvs9VQpfEWkSCgoK+PLLL3n55ZePCN+UlBTGjh3L1KlT6d+/P9dddx0OhwOA0NBQfv/73zN48GDGjx9PVlYWAElJSdxzzz2MGzeOp59+mu3btzN+/HgGDBjA+PHj2bFjBwAXXHBBzSjoxRdf5IorrgBg+vTpNWESHR3NQw89xIgRIxgyZAhr1qxh8uTJ9OrVixdeeKGm9vHjxzN48GDi4+NZuHAhALNmzSI9PZ2EhATuuKPqDSOPPfYYQ4cOZcCAATzwwAM1j/NPf/oTMTExTJgwgbpmB5w+fTrXX389ycnJ9OzZk6VLlzJjxgz69evH9OnTa7YLDQ2t+fqdd95h+vTpfPXVVyxatIg77riDhIQE0tPTSUpKYtWqVQBkZ2cTHR0NVI1wx4wZw+DBgxk8ePAR4XU8tUfEc+bM4aKLLmLq1Kn07t2bO++8E6j6h2bFihVcd9113HHHHZSUlPCrX/2K+Ph4Bg0aVBPIc+bM4ZJLLuH8889n0qRJpKSkMG7cOH7xi1/Qp08fZs2axbx58xg2bBjx8fGkp6cD8MEHH3DGGWcwaNAgJkyYQGZmJtu2beOFF17gySefJCEhgeXLlx9xVGHt2rUMHz6cAQMGMHXqVA4cOFDz/LnrrrsYNmwYffr0Yfny5Sd8/CfFWtsoH3369LHimiVLlni6hGZHPXPd8Xq2ceNG9xdylNdee83OmDHDWmvtiBEj7OrVq621VfUGBATY9PR0W1FRYSdMmGDffvtta621gP3Pf/5jrbX2oYcesjfccIO11tpx48bZ66+/vmbf5513np0zZ4611tqXX37ZXnDBBdZaa/ft22d79eplly1bZnv37m1zcnKstdZec801NffRvXt3+8QTT1hrrb311lttfHy8PXTokN2/f7+NjIy01lpbXl5uDx48aK21Nisry/bq1cs6HA6bkZFhY2Nja+r45JNP7LXXXmsdDoetrKy05557rl26dKldtWqVjYuLs4WFhfbgwYO2V69e9rHHHjumR9dcc4299NJLrcPhsO+//74NCwuzP/zwg62srLSDBw+233//vbXW2pCQkJrbvP322/aaa6455nEd7tN3331XU3f37t2ttdYWFhba4uJia621aWlpNjEx0Vprj3k8h9W+fPbs2bZHjx52586dtri42Hbr1s3u2LHjmPt7/PHH7fTp06211m7atMmedtpptri42M6ePdt26dKl5mexZMkSGx4ebvfs2WNLSkps586d7f3332+ttfapp56yt9xyi7XW2tzcXOtwOKy11v7rX/+yt912m7XW2gceeOCIXtb+Pj4+3qakpFhrrb3vvvtq9jVu3Lia2//3v/+148ePP+YxH+93BlhlncxIvdVIRI7w0Acb2LjnUIPus3/n1jxwfuwJt3njjTe49dZbAbjssst44403GDx4MADDhg2jZ8+eAFx++eWsWLGCiy++GB8fHy699FIArrzySi666KKa/R2+HODrr7/mvffeA+Cqq66qGY1FRUXx8MMPk5yczIIFC2jXrt1xazvnnHMAiI+Pp6CggLCwMMLCwggMDCQvL4+QkBDuueceli1bho+PD7t37yYzM/OY/Xz66ad8+umnDBo0CKgaMW/ZsoX8/HymTp1KcHAwAFOmTKmzT+effz7GGOLj44mKiiI+Ph6A2NhYtm3bRkJCwgn77Izy8nJuvPFG1q5di6+vL2lpaS7dfvz48YSHhxMYGEj//v3Zvn07p5122hHbrFixgptuugmAvn370r1795r7mThx4hE/i6FDh9KpUycAevXqxaRJk4Cqn8fhEfOuXbu49NJL2bt3L2VlZfW+DejgwYPk5eUxblzVzMnXXHMNl1xySc31h59LiYmJbNu2zaXH7wyFr4h4XE5ODosXL2b9+vUYY6isrMQYw6OPPgoce3ZpXWeb1r48JCSkzvurvd2PP/5IREQEe/bsqXP7gIAAAHx8fGq+Pvx9RUUF8+bNIysri9WrV+Pv7090dPRx34pireXuu+/mt7/97RGXP/XUU06fQVtfLUc/vhO9JcbPz6/mEH7t7Z588kmioqJYt24dDoeDwMBAp2o7ukYAX1/f4752a0+wqM/RP7ujH2ftHhze90033cRtt93GlClTSElJ4cEHH3Sp5qMdvo+66j9VCl8ROUJ9I9TG8M4773D11Vfz4osv1lw2btw4VqxYAcDKlSvJyMige/fuvPXWW8ycOROoWrbunXfe4bLLLuP1119n9OjRx93/yJEjefPNN7nqqquYN29ezXYrV67k448/5vvvv2fcuHFMmjTppCZOOHjwIB06dMDf358lS5awfft2AMLCwsjPz6/ZbvLkydx3331cccUVhIaGsnv3bvz9/Rk7dizTp09n1qxZVFRU8MEHHxwT0K6Iiopi06ZNxMTEsGDBgprJLo6uJzo6mtWrVzNs2LAjTpg6ePAgXbt2xcfHh7lz51JZWXnStdRl7NixzJs3jzPPPJO0tDR27NhBTEwMa9asOan9HTx4kC5dugAwd+7P6/6EhYVx6NCxR3LCw8Np27Yty5cvZ8yYMbz22ms1o2B30AlXIuJxb7zxBlOnTj3ismnTpvH6668DMGLECGbNmkVcXBw9evSo2TYkJIQNGzaQmJjI4sWLuf/++4+7/2eeeYbZs2czYMAAXnvtNZ5++mlKS0u59tpreeWVV+jcuTN///vfmTFjxglHZHW54oorWLVqFUOGDGHevHn07dsXgIiICEaNGkVcXBx33HEHkyZN4pe//CUjRowgPj6eiy++mPz8fAYPHsyll15KQkIC06ZNY8yYMS7XUNtf//pXzjvvPM4888yaw7VQdTj/scceY9CgQaSnp3P77bfzz3/+k5EjR5KdnV2z3e9+9zvmzp3L8OHDSUtLO+FRhJP1u9/9jsrKSuLj47n00kuZM2fOESNcVz344INccskljBkzhvbt29dcfv7557NgwYKaE65qmzt3LnfccQcDBgxg7dq1dT5/GoNT6/meDK3n6zpvXP/yVKlnrjtezzZt2kS/fv08U1A9UlJSePzxx/nwww+PuS40NJSCgoJGr6GlTpXYmFp6z473O+PKer4a+YqIiLiZXvMVkSYtKSmpzqMb7hj1ijQGjXxFRETcTOErIiLiZgpfERERN1P4ioiIuJnCV0Q8Licnh4SEBBISEujYsSNdunQhISGBNm3a0L9/f7fWsnbtWj766KOa7xctWsQTTzxxUvuKjo4+4v2zDeXPf/7zKd2+riULnV1u7/LLL2fAgAE8+eSTp1SHK5566imKiopqvj/nnHPIy8tz2/03NJ3tLCIeFxERwdq1a4GqYAgNDeX2229n27ZtnHfeeQ1+fxUVFfj5Hf/P39q1a1m1alXNfM5Tpkxx6/qzzvjzn//MPffc0+D7ve666+rdZt++fXz11Vc1s3g540T9dtZTTz3FlVdeWTP/de1/kJojp0a+xphoY8xHxpgDxph9xpjnjDEKbhFpdJWVlVx77bXExsYyadIkiouLAUhPT+ess84iMTGRMWPGsHnzZoA6lw+cPn06t912G8nJydx1110UFhYyY8YMhg4dyqBBg1i4cCFlZWXcf//9vPXWWyQkJPDWW28xZ84cfv/73wOQmZnJ1KlTGThwIAMHDqxZau/CCy8kMTGR2NhYXnrppXofU2hoKHfddReJiYlMmDCBlStXkpSURM+ePVm0aBFQtbTejTfeWHOb8847j5SUFGbNmkVxcTEJCQlcccUVxyxw//jjj9fMa/yvf/2LoUOHMnDgQKZNm3bEyPF4ao+I61pWb9KkSezfv79mxqi6luU755xzjljW0dnlEK+//nqGDBlCbGxszZKLzzzzDHv27CE5ObnmH6HaRxWeeOIJ4uLiiIuL46mnngKqljns16/fcZ87TYIzSx8BHwFzgECgI/AjcPOJbqMlBV2n5fFcp565rqkuKXhY7SXfMjIyrK+vb81SeZdccol97bXXrLXWnnnmmTYtLc1aa+0333xjk5OTrbV1Lx94zTXX2HPPPddWVFRYa629++67a/Z14MAB27t3b1tQUGBnz55dszShtVVL5F177bXWWmt/8Ytf2CeffNJaa21FRYXNy8uz1tqa5e+KiopsbGyszc7OttZWLUeYlZV1zGME7EcffWSttfbCCy+0EydOtGVlZXbt2rV24MCBNfdbu45zzz235mdXe8nAo5f5e+yxx+wDDzxgrbU1dVhr7R/+8Af7zDPPHNPjunpf17J6R99fXcvyjR49+ohlHZ1dDvFwLysqKuy4cePsunXrjtvLw98fXo6xoKDA5ufn2/79+9s1a9ac8LnTENy1pGAP4DlrbQmwzxjzP8D9s6+LSOP7eBbs+7Fh99kxHs7+60ndtEePHjXL5B1e3q2goICvvvrqiCXgSktLgbqXDwS45JJL8PX1BaqW91u0aFHNSK+kpKRmlFyXxYsX17wm6uvrS3h4OFA1MluwYAEAO3fuZMuWLURERNS5n1atWnHWWWcBVcviBQQE4O/vT3x8fIMuX7d+/Xruvfde8vLyKCgoYPLkyS7dvr5l9epblq/2so7g3HKI8+fP56WXXqKiooK9e/eyceNGBgwYUGeNK1asYOrUqTXzT1900UUsX76cKVOmHPe501Q4G75PA5cZY1KAtsDZwH2NVZSIyGFHL09XXFyMw+GgTZs2Na8Tn0hdywxaa3n33XeJiYk5Yvtvv/3WpfpSUlL4/PPP+frrrwkODiYpKemEy/gB+Pv719RV1xJ5tZf7g7qXBjzRdtOnT+f9999n4MCBzJkzh5SUFJce26kuq1fX0oB1LYeYkZHB448/znfffUfbtm2ZPn16vb20J1if4HjPnabC2fBdClwLHAJ8gbnA+0dvZIyZCcwEiIyMdPkH7e0KCgrUMxepZ647Xs/Cw8N/Xmpu9B8a545rLWV3IqWlpfj7+5Ofn09BQQEOh6OmttLSUkpLSzHG0K1bN1599VWmTp2KtZb169cTHx/PsGHDmD17Npdffjnz5s1j+PDh5OfnU15eTnFxcc2+kpOT+fvf/87jjz+OMYZ169YxcOBA/Pz8yM3NrdmupKQEay35+fmMHTuWJ598khtuuIHKykoKCwvZt28fYWFhVFZWsnr1ar755huKiorIz8/HWktBQcFxV+up/ZgOP97a13Xo0IHVq1dz8OBB9uzZw8qVK2v26+/vT25uLv7+/gQHB5OZmcm2bdsIDQ1l4cKFTJgwgfz8fA4dOkRYWBi5ubm8+uqrdOrUifz8/OPe59G1HH58h38Oh3tQ+2fi4+NDeHg4n3zyCSNHjuTf//43I0aMqHnsh28PHNH/o3+uh6/bu3cvQUFB+Pj4kJ6ezkcffVTz8wsJCWHv3r01vTzc28TERK6//npuuOGGmn+oXnrppTqfO0c/5pNVUlJySn976g1fY4wP8AnwIjASCAVeAf4G3Fl7W2vtS8BLULWqkVabcY1W6HGdeua6ulY1aior0AQEBBAQEEBYWBihoaH4+PjU1BYQEEB5eTlhYWG8+eabXH/99fz973+nvLycyy67jJEjR/L8888zY8YMnnvuOSIjI5k9ezZhYWH4+/sTFBRUs68//vGP3HrrrYwaNQprLdHR0Xz44Yecc845PP3004wZM4a7776bwMBAjDGEhYXx/PPPM3PmTObNm4evry///Oc/mTp1KnPnzmXUqFHExMQwfPhwgoODCQsLwxhDaGjocXtb+zEdfry1r5s4cSJz585l5MiRxMXFMXjw4Jr9zpw5k1GjRjF48GDmzZvHAw88wIQJE+jRowexsbE1+3vkkUcYP3483bt3Jz4+vmaloePd59G1+Pr6EhISQlhYWM0/PMf7mbz22mtcd911FBUV0bNnz5p+G2Nqbg8c0f+j93H4upEjR5KYmMjw4cPp2bMno0ePJjAwkLCwMK677jouueQSOnXqxJIlS2p6O2bMGGbMmMH48eMBmDlzJqNHj2bbtm11PncaQmBgIIMGDTrp29e7pKAxpj2QBbSx1h6svuxC4BFrbVxdt9OSgq5TkLhOPXNdc1tSsClo6cvjNYaW3rNGX1LQWpsNZADXG2P8jDFtgGuAdSdRr4iIiNdzdoari4CzqBoBbwUqgP9rrKJERERaMqdOuLLWrgWSGrcUERER76C5nUUEOPFbNkTkZw3xu6LwFRECAwPJyclRAIvUw1pLTk4OgYGBp7Qfzc8sInTt2pVdu3aRlZXl6VKapJKSklP+Y+ttWnLPAgMD6dq16yntQ+ErIvj7+9OjRw9Pl9FkpaSknNJ7Or2RenZiOuwsIiLiZgpfERERN1P4ioiIuJnCV0RExM0UviIiIm6m8BUREXEzha+IiIibKXxFRETcTOErIiLiZgpfERERN1P4ioiIuJnCV0RExM0UviIiIm6m8BUREXEzha+IiIibKXxFRETcTOErIiLiZvWGrzGm4KiPSmPMs+4oTkREpCXyq28Da23o4a+NMSFAJvB2YxYlIiLSkrl62PliYD+wvBFqERER8QrGWuv8xsYsBpZZax+s4/qZwEyAyMjIxPnz5zdEjV6joKCA0NDQ+jeUGuqZ69Qz16lnrvPGniUnJ6+21g5xZlunw9cY0w3IAE631mbUt31MTIxNTU11at9SJSUlhaSkJE+X0ayoZ65Tz1ynnrnOG3tmjHE6fF057Hw1sMKZ4BUREZG6uRq+cxurEBEREW/hVPgaY0YCXdBZziIiIqfM2ZHvNcB71tr8xixGRETEG9T7Pl8Aa+1vG7sQERERb6HpJUVERNxM4SsiIuJmCl8RERE3U/iKiIi4mcJXRETEzRS+IiIibqbwFRERcTOFr4iIiJspfEVERNxM4SsiIuJmCl8RERE3U/iKiIi4mcJXRETEzRS+IiIibqbwFRERcTOFr4iIiJspfEVERNxM4SsiIuJmCl8RERE3U/iKiIi4mdPha4y5zBizyRhTaIxJN8aMaczCREREWio/ZzYyxkwE/gZcCqwEOjVmUSIiIi2ZU+ELPAQ8bK39pvr73Y1Uj4iISItXb/gaY3yBIcAiY8xWIBB4H7jDWltc1+2Ci3bD7HMbrFBvkJCXBxltPF1Gs6KeuU49c5165jr17MScGflGAf7AxcAYoBxYCNwL/KH2hsaYmcBMgLiOAeTl5TVosS1dZWWleuYi9cx16pnr1DPXqWcnZqy1J97AmLZALjDdWju3+rJpwL3W2kF13S4mJsampqY2ZK0tXkpKCklJSZ4uo1lRz1ynnrlOPXPeweJyPt+YSWrqZs4enUjP9qGEB/t7uqxGVV7p4JMN+zh/YJfV1tohztym3pGvtfaAMWYXcOKUFpFGl11Qyicb9jHm9Ei6RQR7uhwRACodli+3ZvPO6l18smEfpRUOAF764SsA2gb7E90+hB7tQ+gREVLzdXT7EEIDnD31qOnJLSzjjZU7eO3r7ew7VOLSbZ191LOBm4wx/6PqsPOtwIeulSkip2L19lxumPc9+w6V4GNgYv8ofj26J0Oj22KM8XR54oV+yirg3TW7eG/NbvYeLKFNsD+XDT2NiwZ3ZcO6NUT2jGVbdiE/ZReyLbuQr9NzeG/NkefrRoYF0CPi5zDu0T6Y+K5t6NImyEOPqn6b9h5i9pcZvL92D2UVDsb0bs+fpsYx4W/O78PZ8P0j0B5IA0qA+cCfXC1YRFxnreWVL7fxl4820blNEK/OGMY3P+Xw+sodfLIhk/gu4fx6dA/Oie9EKz/Nm9NUWGv57497WZ6WTZe2QXSPCKZ7RAjd2wXTJti/2f7DdKiknP/+sJd3Vu9i9fYD+BhIiunAfef1Z3y/DgT4+QJwIN2HpP5Rx9y+uKySbTlVYZxx+HN2IV9s3k92QSkAfj6GGaN7cPP43k1mZFzpsHy2MZPZX2bwbUYuQf6+XJLYlekjo+kdFeby/px6VNbacuB31R8i4ib5JeXc9e4PfPTjPib2j+LxSwYSHuTP2D6R3HRmb95ds4tXvszg1rfW8pePN3H1iGh+OawbbUNaebp0r5Z5qIR731/PZxszCQvwI7+04ojrwwL9iI4IoVtEMN3bBf/8dUQwUWGB+PicXDA7HJbSCgfF5ZUUl1fi72uICAnA9yT3V3u/X6Xn8M7qnfxvwz5Kyh2c3iGUu8/uy9RBXejQOtDpfQW18qVfp9b069T6mOvyS8rZll3Ef77ZzkvLfmLh2t384dz+nD+gk8f+WTlYVM5bq3Yw96vt7M4rpkubIO45py+XDul2Sq9lN41/KUTkGJv3HeJ3/1nD9twi7j67LzPH9jziD1BQK1+uHN6dXw7rxtItWbyyIoPHPknl2cVbuGhwV2aMiub0Dq7/Ry4nz1rL26t38ccPN1JW4eCec/oyY1QPyistO3KL2J5TWP25iG05hazffZD/rd9HpePnU2oC/Hzo1q4qiLu2rXpdv7issiZQS8orj/z+iOscx9TkYyAiNIAOYQFEhtX+HHjM90GtfI+47bbsQt5ds4t3V+9iz8ESWgf6cXFiVy5OPI2BXcMbPBDDAv2J7xrO3y4ewKXDTuP+heu5+Y3veePbHTx0QSx9TmKEebK27s9n9pfbeG/NborLKzmjRzvuO68fE/pF4ed76keYFL4iTdB7a3Zxz4IfCQv05/XfnMEZPSPq3NbHx5Ac04HkmA6k7stn9pcZvLN6F69/u4NxfSL59egejOndvtke5mwudh0o4u73fmT5lmyGRbfjr9Pi6RkZCoCfL8R0DCOm47HhUVHpYE9eCdtyCtmeW8SOnEK251QF9Dc/5eLrYwjy9yWolS+B/r4E+fsQ1MqXNsH+1d9XXRfkX319zdc+lFY4yMovZf+hUrIKStmfX8KmvYfILig7IvAPCw3wo0NYAO3DAiircLB2Zx4+Bsb0juSec6uCJ9Df95jbNYbB3dqy8IbRvLFyB499kso5Ty/nV6OiuWVCn0Y7FO1wWFLS9jP7y20s35JNKz8fLkzozDUjo4ntHN6g96XwFWlCSsorefjDjbz+7Q7O6NGOZ385iA5hzh/Si+kYxl+nDeCOyTHM+3YHr369natfWUmfqFBmjOrBhYO6NGL13snhsMz7djt//XgzFnj4gliuPKO704eO/Xx96BYR7Naz1x0OS25RWVUw55dWfy4hq+brUixw51kxXDSoKx3DnX8ONiRfH8OVw7tzTnwnHv3fZv61PINF6/Zwzzn9mDKwc4P8Q+lwWFbvOMDHP+7jkw372J1XTFTrAO6YHMNlQ08jIjSgAR7JsRS+Ik3Eztwirp+3mvW7D3HduF7cPqnPSR/eiggN4ObxvfntuJ58uG4vL6/IYNZ7P/LoJ6nEhFey3rGFXpGh9OoQSveI4JqTZMQ127ILufPdH1iZkcuY3u3589R4TmvX9N8C5uNjaB8aQPvQAPo1g5n624W04q/TBnDp0NO4f+EGbnlzLW+s3MHDF8Sd1KHoikoH32bk8vH6vXyyIZOs/FJa+fkwtnd7Zp3dl7PiOuLfAIeWT0ThK9IEfLEpk/97ay0W+NfVQ5h4nLNET0aAny/TErty0eAufJuRy6tfb+PbLZl8/WlazTa+PoZu7YLpFRlSFciRofTqUPV1m2CduHU8lQ7LKysy+Ptnqfj7+vDotAFcMqSrDu03skHd2vL+DaN487sdPPq/VM5+ejm/GhnNLRN6ExZ44pOfyiocfJmezf9+3MenG/dxoKicIH9fkmIiOTu+E2f27eDWM6sVviIeVFHp4InP0ng+JZ3+nVrzzysH0z0ipMHvxxjD8J4RDO8ZQUpKCkNHjCYju5D0rALS9xeQnlX19bIt2ZRV/HzSTkRIqyPCeFC3tiR2b9vg9TUnaZn53PnOD6zdmceEfh3409R4olw421dOja+P4YozunMn81XdAAAeXUlEQVR2XCce+2QzL39ZdSj6D+ceeyi6pLySpWlZ/G/9Pj7flEl+SQWhAX6M79eBs+M6Mq5Ph2NOMnMXha+Ih2Tll3LzG9/z9U85XD7sNB44P9ZtJ7OEBPgR1yWcuC5HnkRS6bDsPlDM1qx80vdXh3NWAZ9syCS3cCcAFw3qwgPnx7b4KQOPVl7p4IWUdJ5dvJWQAF+eviyhwV53FNe1C2nFXy4awKVDu3H/wvXc8uZaXv92B3ef049dB4r4eP0+lmzeT1FZJeFB/kyO7cjZcR0ZdXp7t/2enYjCV8QDVmbkcuPrazhYXM7jlwzk4sSuni4JqD4EXX3yz5l9j7wut7CMOV9m8I+UdFZszebPU+OZ0ECHxxtSpcOSU1hK5sFSMg+VkJlfQl71IcawQD9aB/lXfQ70p3Vg1ddhgX4nfH19+6FKHn3uSzbuPcR5Azrx4JRY2jfSiTjimoTT2rDgd6N467udPPrJZi78x5dA1VGbCxK6cE58R4b3jGj013BdpfAVr+NwWL5Mz2b3gWIuHNTFrf8Fl1U4eGlZOk9+voXT2gYxd8aw40420BS1C2nFbZNimBTbkdvfXsdvXl3F1EFdeOD8/m55bdhay6GSCvYfKmHfoRIyD1WHa/XHvkOl7D9Uwv780uO+jaY+wa18q4PYn9bVn8MC/fAxhg9/KCEi1PLiVYlMju3YCI9OToWvj+GXZ3Tj7LiOLFy7m76dWjM0ut0pTy7SmBS+4jV2HSji7VW7eGf1LnbnVS1F/Y+UrTw0JZYz+zb+CO7bn3K49/31bNlfwLnxnfjLtHha13OSSFMU1yWcRTeO5rklW3l+ydaaUXBDnSRWm8NhWbx5P7O/ymDN9jyKyyuP2aZ1oB8dwwOJah3I6ZHt6RgeQFTrwJqPjq0DaRPsT3FZJfklFRwqKedQSTn5JRVV3xcf/vrIy/OKytiRW0RBaQWjOvvx3K/Hed2h9uambUgrpo/q4ekynKLwlRatpLySzzZmMn/VTlZszQZg9OlVbycID/LnoQ82MGPOKib2j+KB8/vXzCjUkHILy/jLR5t4e/UuurQJ4t9XD2mSh2td0crPh9sm9mFS/yjueOcHrn11FRcmdOaB82MbZGrLorIK3l2zm9krMvgpu5BO4YFcNuw0OocH0aF1AB1rhauzJ8wE+vuedG0pKSkKXmlQCl9pkTbuOcT8VTt5f+1u8orK6dImiFvG9+bixK5HBOzHt4zl5RUZPPPFFiY8sZSbzuzNb8b0aJD3vToclrdX7+QvH2+moKSC65N6cdOZpxPcquX82sV1CWfhDaN4PmUrzy3eyoqtOfxpatxJH5rdd7CEuV9v4/Vvd3CwuJyBXcN55vJBnO2G912KuFPL+SsgXu9gcTmL1u5m/qpd/Lj7IK18fZgUG8WlQ09jVK/2x51xqJWfD9cn9WJKQmf++MFGHvsklXfX7OLhKXGM7t3+pGtJ3ZfPHxb8yKrtBxga3ZY/TY1367y07tTKz4dbJ/RhUv+q14J/+9pqpgzszENTnB8F/7jrIC+v+IkPf9iLw1omx3bk16N7kNhdyyVKy6TwlWbNYS1fpWcz/7udfLy+ahHvvh3DeOD8/lyY0MXpP/5d2gTxwlWJLEndz4OLNnDly99y7oBO3Hduf5em1isqq+DpL7bw8vIMwgL9ePTiAVw8uOtJr1LTnPTv3JqFN47inynpPLt4C1+lZ/PIhfGcFXf8UXClw/L5pkxeXpHByoxcQlr5cvWIaH41KrpZzBIlcioUvtLsWGvZuPcQn23MZN6XxWQVf0tYoB+XDOnKpUO6Edel9UmPlpJjOjDi1gheWJrO8ynppGzez/9N7MM1I6PrPez5+cZMHli0gd15xfxiSFdmnd2Pdl62tJ+/rw83j+/NxP5R3P72Oq77z2rOrx4FH+5FYWkFb6/ayeyvtrE9p4gubYK499x+/GLoac3yBDSRk6HwlWahvNLBtz/l8tnGfXy+aT+784oxBmLa+nDPlAGcFdupwWaqCfT35dYJfZg6qAsPLtrAI//dxNurdvHHC+MY1qPdMdvvzivmoUUb+HRjJn2iQpn/2xHH3c6b9OvUmvdvGMULKek8s3gLX6dnM+vsfmzJzOf1lTvIL6lgcLc23Dm5L5NjG2aJNpHmROErTdbB4nKWpmXx2cZMUlL3k19SQYCfD2N6t+fm8adzZt8oNqz+mqRBjTNBRfeIEF6ZPpRPN2by8Acb+cWLX3PR4C7cfXY/IsMCKK90MPvLDJ78bAsWy11n9eXXo3vQyk9BAlWj4JvG92ZibNUo+Pa31+Fj4Oz4Tvx6dA8Gd/PuaSrFuyl8pUnZdaCIzzdm8vmm/XzzUw4VDktESCvOjuvIhH5RjOkd6da5WI0xTI7tyJje7Xlu8Vb+tfwnPtuYybVjevLRj3vZvC+f8X078OCUWL1OWYe+HVuz4HejSEnNol+nsEZ5O5dIc6PwFaccLCpne24h23KK2J5d/Tmn6nNRWQVtg1sRHuRP2xB/2gS3om2w/8+XBbeqdXnVda0D/fHxMVhr2bDnEJ9uzOTzjZls3HsIgF6RIfxmTE8m9u9AwmltPT5TTXArP+48qy8XDe7K/QvX88RnaXQKD+TFqxKZ1D9KZ+TWw9/Xp1Em4RBprhS+AlSdxJRbWHZEqNb+nFdUfsT2HVsH0j0imPF9OxAa6EdeUTl5RWXkFZezae+hmu/rmuXPGAgP8sfHGHILy/AxMKR7O+45py8T+kXRMzLUDY/adad3CGXeb85gzY4D9O3YmhA3LkEmIi2HU385jDEpwHCgovqi3dbamMYqShqPtZa9B0tIzcwnbV8+qZn5bMksYFt2IfmlFTXb+Rjo3CaI6IgQzo3vRHRECN0jgoluH0K3dsFOzYfscFjySyo4UB3KB4rKyCsq40Dhz0FdXFbJGT0jSI6JJKKZTFRvjCGxu3efUCUip8aVf9tvtNb+u9EqkQaXW1hG6r580jLzjwjb/JKfQzaqdQB9osKYlti1KlyrQ7Zr2+BTPnHIx8cQHuyvaflERI6iY2YtQKXD8sOuPFKrwzUtM5/UfQVkF5TWbBMe5E9MxzAuTOhCn45hxESF0Scq1C2r0YiIyJFcCd+/GGP+CqQCf7DWpjROSeKKzfsOcde7P7JuZx4AQf6+9IkKJTkmkpiOYfSJCqNvxzAiwwJ0UpCISBNhrK1/3UtjzBnARqAMuAx4Dkiw1qYftd1MYCZAZGRk4vz58xu84JasoKCA0FDnTjQqd1g+SC/nvz+VE+wHF/dpRb8IX9oHGXy8KGRd6ZlUUc9cp565zht7lpycvNpaO8SZbZ0K32NuZMz/gP9aa5+ta5uYmBibmprq8r69WUpKCklJSfVut3p7Lne9+yNb9xcwdVAX7juvv9dNY3iYsz2Tn6lnrlPPXOeNPTPGOB2+J/uarwW8Z3jVRBSUVvDY/zbz6jfb6RwexOxfDSU5poOnyxIRERfVG77GmDbAGcBSqt5qdCkwFri1cUuT2pak7ufeBevZc7CYa0ZEc/vkGEL1HlMRkWbJmb/e/sAjQF+gEtgMXGit1TFlN8gtLOOPH25kwfe7Ob1DKO9cN5LE7poTV0SkOas3fK21WcBQN9QitVhrWbRuDw99sJH8knJuHt+bG5J7EeDnvnmNRUSkcei4ZRO0J6+YPyz4kSWpWQw8rQ2PThtATMcwT5clIiINROHbhDis5dWvt/G3jzfjsHDfef2ZPjLa44sKiIhIw1L4NhFb9+fz15UlpB3YwJje7fnz1HgtUSci0kIpfD0sLTOf5xZv5cMf9hDkB49fMpBpg7toNioRkRZM4esh63cf5B9LtvLx+n0Et/Ll2jE9ifXdx5TErp4uTUREGpnC183W7szj2S+28MXm/YQF+HHTmafzq1E9aBfSipSUTE+XJyIibqDwdZOVGbk8u3gLy7dk0ybYn99P7MPVI6MJD9JyeyIi3kbh24istXyVnsMzX2zh24xc2oe2YtbZfblyeHfNTiUi4sWUAI3AWktKahbPLN7C9zvyiGodwP3n9efyYd0IaqVJMkREvJ3CtwE5HJZPN2by3JItrN99iC5tgnjkwjguTuxKoL9CV0REqih8G8i27EKu+89qNu/Lp3tEMI9OG8DUwV3w9/XxdGkiItLEKHwbwJ68Yq7497cUlVXw1KUJnDegE34KXRERqYPC9xRlF5Ry5b+/5VBxOW/MHE5cl3BPlyQiIk2cwvcUHCwq56qXV7L3YAmv/XqYgldERJyiY6MnqbC0gulzVpK+v4CXrk5kSHQ7T5ckIiLNhEa+J6GkvJKZr63ih10H+ccvBzOmd6SnSxIRkWZEI18XlVc6uPH17/lyaw6PXTyAs+I6erokERFpZhS+Lqh0WG5/ex2fb8rkjxfEctFgLYIgIiKuU/g6yVrLve+vZ+HaPdx1Vl+uGhHt6ZJERKSZUvg6wVrLnz/axBsrd3BDci+uT+rl6ZJERKQZcyl8jTG9jTElxpj/NFZBTdGzi7fyr+UZXDOiO7dPivF0OSIi0sy5OvL9B/BdYxTSVL2yIoMnPktj2uCuPHB+LMYYT5ckIiLNnNPha4y5DMgDvmi8cpqW+d/t5OEPN3J2XEf+Ni0eHx8Fr4iInDqnwtcY0xp4GPh945bTdHz4wx5mvfcDY/tE8tRlCZqrWUREGoyx1ta/kTFPA3ustX8zxjwInG6tvfI4280EZgJERkYmzp8/v4HLdY91WRU8s6aUXm18+P2QQAJ83TPiLSgoIDQ01C331VKoZ65Tz1ynnrnOG3uWnJy82lo7xJlt653hyhiTAEwABtW3rbX2JeAlgJiYGJuUlORMDU3K1+k5PP/5Svp3DmfetWfQOtDfbfedkpJCc+yZJ6lnrlPPXKeeuU49OzFnppdMAqKBHdUnG4UCvsaY/tbawY1Xmvv9sCuP38z9jm7tgpk7Y5hbg1dERLyHM+H7EvBmre9vpyqMr2+Mgjyl0mG54+0faBPciv/85gzahbTydEkiItJC1Ru+1toioOjw98aYAqDEWpvVmIW528K1u0nNzOfZywcR1TrQ0+WIiEgL5vKqRtbaBxuhDo8qq3Dw5Odp9O/UmnPjO3m6HBERaeH0/hngze92sDO3mDvOitF7eUVEpNF5ffgWlVXwzBdbGdajHUl9tC6viIg0Pq8P39lfbiO7oJS7zorR1JEiIuIWXh2+B4vKeXFpOuP7diCxeztPlyMiIl7Cq8P3hWXp5JdWcPtkrVQkIiLu47Xhu/9QCbO/zGDKwM7069Ta0+WIiIgX8drwfWbxFioqLbdN7OPpUkRExMt4ZfhuzynkzZU7uWzYaXSPCPF0OSIi4mW8Mnyf/CwNP1/DzWf29nQpIiLihbwufDfvO8TCdXuYPrIHHTSNpIiIeIDXhe/jn6QSGuDH9eN6eboUERHxUl4Vvqu25fL5pv1cN64X4cFaLlBERDzDa8LXWsujn6TSPjSAX42K9nQ5IiLixbwmfJemZbEyI5ebx59OcCuXF3MSERFpMF4Rvg6H5bFPUjmtXRCXDe3m6XJERMTLeUX4frR+Lxv2HOL/JvShlZ9XPGQREWnCWnwSVVQ6eOLTNPpEhXJBQhdPlyMiItLyw/ed1bv4KbuQ2yfF4OujJQNFRMTzWnT4lpRX8tTnWxjUrQ0T+0d5uhwRERGghYfva19vZ9+hEu6c3BdjNOoVEZGmwanwNcb8xxiz1xhzyBiTZoz5TWMXdqryS8p5PmUrY3q3Z0SvCE+XIyIiUsPZke9fgGhrbWtgCvCIMSax8co6df9ansGBonLunNzX06WIiIgcwanwtdZusNaWHv62+qPJTo6cU1DKy8t/4pz4jsR3Dfd0OSIiIkdw+jVfY8zzxpgiYDOwF/io0ao6Rf9Ykk5xeSW3TYzxdCkiIiLHMNZa5zc2xhcYASQBf7PWlh91/UxgJkBkZGTi/PnzG65SJ2UXO5i1rJiRXfyYERfg9vs/FQUFBYSGhnq6jGZFPXOdeuY69cx13tiz5OTk1dbaIc5s61L41tzImBeAjdbaZ+raJiYmxqamprq871N1x9vrWLhuDym3J9G5TZDb7/9UpKSkkJSU5OkymhX1zHXqmevUM9d5Y8+MMU6H78m+1ciPJviab15RGe99v5tfDuvW7IJXRES8R73ha4zpYIy5zBgTaozxNcZMBi4HFjd+ea5ZviWbSodlSkJnT5ciIiJSJ2fW1rPA9cALVIX1duBWa+3CxizsZCxLyyI8yJ+BXdt4uhQREZE61Ru+1tosYJwbajkl1lqWbcli9OntNYeziIg0aS1mesm0zAIyD5Uytk97T5ciIiJyQi0mfJem7QdgbJ9ID1ciIiJyYi0mfJelZdMnKpRO4TrLWUREmrYWEb5FZRWszMhlbG+NekVEpOlrEeH77U+5lFU6GBej8BURkaavRYTv0rQsAv19GBrdztOliIiI1KtFhO+yLVmc0SOCQH9fT5ciIiJSr2Yfvjtzi/gpq1BnOYuISLPR7MN32ZYsAMYpfEVEpJlo/uGblkWXNkH0igzxdCkiIiJOadbhW17p4KutOYzt0x5jNKWkiIg0D806fL/fkUd+aYUOOYuISLPSrMN3WVoWvj6GkadrPmcREWk+mnf4bsli0GltaB3o7+lSREREnNZswzenoJQfdx/UIWcREWl2mm34rtiajbVaxUhERJqfZhu+S9OyaBvsT1yXcE+XIiIi4pJmGb4Oh2VZWjaje0fi66O3GImISPPSLMN3075DZBeU6vVeERFplppl+C5LywZgbG+9xUhERJqfZhq+WfTtGEaH1oGeLkVERMRl9YavMSbAGPOyMWa7MSbfGPO9MeZsdxR3PIWlFazansu4GB1yFhGR5smZka8fsBMYB4QD9wHzjTHRjVdW3b5Oz6G80jKut8JXRESaJ7/6NrDWFgIP1rroQ2NMBpAIbGucsuq2bEsWQf6+JEa3dfddi4iINAhjrXXtBsZEAduBBGvt5qOumwnMBIiMjEycP39+Q9VZ485lRXQK8eH/Elve670FBQWEhoZ6uoxmRT1znXrmOvXMdd7Ys+Tk5NXW2iHObOtS+Bpj/IGPgXRr7W9PtG1MTIxNTU11et/O2J5TyLjHUnhoSizXjIxu0H03BSkpKSQlJXm6jGZFPXOdeuY69cx13tgzY4zT4ev02c7GGB/gNaAMuPEkazsly9KyAE0pKSIizVu9r/kCmKqV6l8GooBzrLXljVpVHZamZXFauyCiI4I9cfciIiINwtmR7z+BfsD51triRqynTmUVDr5Oz2Fcn0iq/hcQERFpnpx5n2934LdAArDPGFNQ/XFFo1dXy+rtBygsq2Ss3mIkIiLNnDNvNdoOeHyouWxLFn4+hhG9IjxdioiIyClpNtNLLk3NIrF7W8IC/T1dioiIyClpFuGblV/Kxr2HdJaziIi0CM0ifJdvqXqLkZYQFBGRlqBZhO/StCwiQlrRv1NrT5ciIiJyypp8+DocluVbshnbJxIfH4+f9yUiInLKmnz4bthziNzCMsb2ae/pUkRERBpEkw/fpWn7ARij9/eKiEgL0eTDd1laNnFdWtM+NMDTpYiIiDSIJh2++SXlrNlxQLNaiYhIi9Kkw/er9BwqHFbv7xURkRalSYfv0rQsQlr5MrhbW0+XIiIi0mCabPhaa1mWlsXI09vTyq/JlikiIuKyJptqGdmF7DpQrEPOIiLS4jTZ8F2aVj2lpE62EhGRFqbJhu+ytCx6tA+hW0Swp0sRERFpUE0yfEvKK/nmp1zG9tasViIi0vI0yfBdte0AxeWVer1XRERapCYZvsu2ZOHvaxjeM8LTpYiIiDS4phm+aVkMjW5HSICfp0sRERFpcE0ufDMPlbB5X74OOYuISIvlVPgaY240xqwyxpQaY+Y0ZkGfbswE0HzOIiLSYjl7XHcP8AgwGQhqrGLKKx28uDSdgV3D6dcprLHuRkRExKOcGvlaa9+z1r4P5DRmMe+t2cWuA8XcMqE3xpjGvCsRERGPaTKv+ZZXOnhuyVYGdA0nOaaDp8sRERFpNMZa6/zGxjwCdLXWTq/j+pnATIDIyMjE+fPnO73vZbvKeWV9GbcODiChg3ee5VxQUEBoaKiny2hW1DPXqWeuU89c5409S05OXm2tHeLMtg0avrXFxMTY1NRUp/ZbXungzL+n0CaoFYtuHOW1h5xTUlJISkrydBnNinrmOvXMdeqZ67yxZ8YYp8O3SRx2XvD9bnbmFnOrXusVEREv4NTxXWOMX/W2voCvMSYQqLDWVpxqAeWVDp5bvJX4LuGc2Vev9YqISMvn7Mj3XqAYmAVcWf31vQ1RwILvd7Mjt4hbxmvUKyIi3sGpka+19kHgwYa+84pKB/9YspW4Lq0Z30+jXhER8Q4efc13wfe72Z5TxC3j+2jUKyIiXsNj4VtR/b7e2M6tmaBRr4iIeBGPhe/7a/dUj3r1Wq+IiHgXj4RvRaWDZxdvoX+n1kzsH+WJEkRERDzGI+G7sHrUq/f1ioiIN3J7+GrUKyIi3s7t4bto3R625RRp5SIREfFabg3fqlHvVvp3as0kjXpFRMRLuTV8F63bQ0Z2ITfrDGcREfFibgvfiuo5nPtp1CsiIl7ObeH7wQ97+Cm7kFvGn46Pj0a9IiLivdwSvpUOy7NfbKVvxzAm9e/ojrsUERFpstwSvh+sqxr13jqht0a9IiLi9Ro9fCsdlmcWb9GoV0REpFqjh++HP+zhp6xCbhmvUa+IiAg0cvhWOixPf1E16p0cq1GviIgINHL4Hh713qxRr4iISI1GDd9nvthCTFQYZ2nUKyIiUqPRwrew3JKuUa+IiMgxGi1880otfaJCOTtOo14REZHanApfY0w7Y8wCY0yhMWa7MeaX9d2m3AG3jO+jUa+IiMhR/Jzc7h9AGRAFJAD/Ncass9ZuqOsG/j5o1CsiInIc9Y58jTEhwDTgPmttgbV2BbAIuOpEt4sI8tGoV0RE5DicOezcB6i01qbVumwdEHuiGwX6nkpZIiIiLZczh51DgYNHXXYQCDt6Q2PMTGBm9belxpj1p1ae12kPZHu6iGZGPXOdeuY69cx13tiz7s5u6Ez4FgCtj7qsNZB/9IbW2peAlwCMMaustUOcLUTUs5OhnrlOPXOdeuY69ezEnDnsnAb4GWN617psIFDnyVYiIiJSt3rD11pbCLwHPGyMCTHGjAIuAF5r7OJERERaImcn2fgdEATsB94Arj/R24yqvXQqhXkp9cx16pnr1DPXqWeuU89OwFhrPV2DiIiIV2n09XxFRETkSApfERERN2vw8D2ZeaC9nTEmxRhTYowpqP5I9XRNTY0x5kZjzCpjTKkxZs5R1403xmw2xhQZY5YYY5x+r11LVlfPjDHRxhhb6/lWYIy5z4OlNgnGmABjzMvVf7fyjTHfG2POrnW9nmdHOVHP9Dw7MWfndnaFy/NACwA3Wmv/7ekimrA9wCPAZKpO/gPAGNOeqrPxfwN8APwReAsY7oEam5rj9qyWNtbaCveW1KT5ATuBccAO4BxgvjEmnqr5DvQ8O9aJenaYnmfH0aDhW2se6DhrbQGwwhhzeB7oWQ15X+JdrLXvARhjhgBda111EbDBWvt29fUPAtnGmL7W2s1uL7QJOUHP5Diq31b5YK2LPjTGZACJQAR6nh2jnp6t9khRzURDH3Y+qXmgBYC/GGOyjTFfGmOSPF1MMxJL1XMMqPljkI6ec87YbozZZYyZXX0EQWoxxkRR9TdtA3qeOeWonh2m59lxNHT4Oj0PtBzhLqAn0IWq98Z9YIzp5dmSmg0951yXDQylah7aRKp6Nc+jFTUxxhh/qnoyt3pkq+dZPY7TMz3PTqChw9fpeaDlZ9bab621+dbaUmvtXOBLql47kfrpOeei6qVBV1lrK6y1mcCNwCRjzNF99ErGGB+qZvAro6o3oOfZCR2vZ3qenVhDh6/mgW4YFtBiyM7ZQNVzDKg576AXes654vBMO17/nDPGGOBlqk4YnWatLa++Ss+zOpygZ0fT86yWBg1fzQPtOmNMG2PMZGNMoDHGzxhzBTAW+MTTtTUl1b0JBHwB38P9AhYAccaYadXX3w/84M0nwRxWV8+MMWcYY2KMMT7GmAjgGSDFWnv0YVVv9E+gH3C+tba41uV6ntXtuD3T86we1toG/QDaAe8DhVSdev7Lhr6PlvQBRALfUXX4Kg/4Bpjo6bqa2gdVZ1Taoz4erL5uArAZKAZSgGhP19sUPurqGXA5kFH9O7oXeBXo6Ol6Pf1B1WuTFiih6jDz4Y8rqq/X88yFnul5duIPze0sIiLiZppeUkRExM0UviIiIm6m8BUREXEzha+IiIibKXxFRETcTOErIiLiZgpfERERN1P4ioiIuJnCV0RExM3+H1ighzAH+AgAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_mi.plot(figsize=(8, 5))\n",
    "plt.plot(history_range, true_mi, label='Theoretical mutual information')\n",
    "plt.grid(True)\n",
    "plt.gca().set_ylim(0, dim_n*3)\n",
    "plt.legend()\n",
    "#plt.tick_params(axis='x', colors='white')\n",
    "#plt.tick_params(axis='y', colors='white')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}