Gaussian_MINE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/Gaussian_MINE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.gaussian import Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 400   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 6               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "g = Gaussian(sample_size=sample_size,rho=rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = g.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XvQXVd53/HvIwsL40tAsRsH7FcmhpBQj4bLW4MnF2hwsGnd0CglJWEwmRQ05DYhQ20KJkC4ZIzd8ZAplKIWCgbCpcEEikkENE08NDaDxFBFYGhtwmsb2xgQjrBxbYSf/nH2kte7tO9nn7P2Oef3mdGg85599l7nFV7Puj7L3B0REZEtuQsgIiLjoIAgIiKAAoKIiBQUEEREBFBAEBGRggKCiIgACgiSgZn9tZm9eE7P+i0z+6aZ3WNmP9ri+q+b2fnzKNsYmNm7zeyNucsh46CAIDNRVKz3FRXxN83sv5rZSR3vcZaZuZlt7VmGhwFXAc9295Pc/Tt97lNzfzezxw15T5GcFBBklv6Fu58EPAX4J8Cr5/z8HwMeDnxpzs8VWUgKCDJz7v4N4C+Ac9L3zGyLmb3azDbM7C4zu9rMfqR4+7rif+8uehrnlXx+m5m9xcxuL/68pfjZTwJfjT7/V2VlM7MXFs/+jpldlrx3rpldb2Z3m9kdZvZWMzu+eC+U7X8XZfvXZvYoM/uEmX3LzL5b/P2Mqt+Lmb3CzL5hZt8zs6+a2bOanlu872b222b2f4vPvsHMzi4+c9jMPhyV85lmdpuZvcrMvl303F5QU6aLzOyLxbP/1sx2NpVXloi764/+DP4H+DpwfvH3M5m00t9QvP5r4MXF338TuAn4CeAk4BrgvcV7ZwEObK15zuuBG4B/BJwG/G30nNrPA08E7gF+HtjGZHjpSFTupwJPB7YW97oReFn0eQceF73+UeBXgEcAJwP/Dfjzimc/AbgVeHRU1rM7PPfjwCnAPwbuB/5H8Tv8EeDLwIuKa59ZfKeriu/4DOBe4AnF++8G3lj8/SnAXcDTgOOAFxX/jtvqyqs/y/NHPQSZpT83s7uBzwJ/A/xxyTUvAK5y96+5+z3AK4Hnd5g3eAHwene/y92/BfwR8MKWn/1XwCfc/Tp3vx/4Q+DB8Ka773f3G9z9iLt/HXgHkwq1lLt/x90/4u7fd/fvAW+quf6HTCraJ5rZw9z96+5+c4fnvtndD7v7l4CDwKeK3+E/MOmNPTm5/g/d/X53/xvgWuBXS8r0EuAd7v45d/+hu7+HSbB5el15ZXkoIMgs/Ut3f6S773D333b3+0queTSwEb3eYNIy/rGWzyj7/KM7fPbW8MLd7wWOTjyb2U8Wwz53mtlhJgHt1KqbmdkjzOwdxRDUYSZDXo80s+PSa939JuBlwOuAu8zsg2b26A7P/Wb09/tKXscT+N8tvltQ9TvaAby8GC66uwjmZzLpFVSWV5aHAoLkdjuTiihYYzLE8U0mQyN9Pn97y2ffwaTCAyYVOpNhn+DtwFeAx7v7KcCrAKu538uZDK08rbj+58Otyy529z91958tyu/Am3s+t8mjzOzE6HXV7+hW4E1FEA9/HuHuH2gorywJBQTJ7QPAH5jZY4tlqX8MfMjdjwDfYjKE8xMNn3+1mZ1mZqcCrwHe1/LZfwZcZGY/W0zCvp7N/02cDBwG7jGznwJ+K/n8N5OyncykdX63mW0HXlv1YDN7gpn9gpltA/5f8bkftnxuH39kZseb2c8BFzGZ30j9Z+ClZvY0mzjRzP65mZ3cUF5ZEgoIktu7gPcyGV75eyaVze8BuPv3mYzD/69iCOPpJZ9/I7APOAD8HfCF4meNivH33wH+lElv4bvAbdEl/xb4deB7TCrLDyW3eB3wnqJsvwq8BTgB+DaTie6/rHn8NuDy4to7mUyKv6rlc7u6k8l3ux14P/BSd/9KepG772Myj/DW4vqbgN9oUV5ZEuauA3JElpWZPRN4n7tXLn8VCdRDEBERQAFBREQK2YaMzOzhTMaNtzFZZvhn7l45CSciIrOVMyAYcKK732OTJGSfBX7f3W/IUiARkRXXK4vkEHwSie4pXj6s+FMbnU499VQ/66yzZlwyEZHlsn///m+7+2lN12ULCADFDs79wOOAt7n750qu2Q3sBlhbW2Pfvn3zLaSIyIIzs43mqzJPKhf5Up4EnAGca2bHZMN09z3uvu7u66ed1hjgRESkp1GsMnL3u5lkwLwwc1FERFZWtoBQpBp4ZPH3E4DzmeRvERGRDHLOIfw4k23/xzEJTB92909kLI+IyErLucroAMfmbBcRkUxGMYcgIiL5KSCIiIzYJVdfzyVXXz+XZykgiIgIkHljmoiIlAu9ggMbhza9vvLi82b2TPUQREQym+ewUB31EERERij0BObRMwgUEEREMskxLFRHAUFEZMTmGRwUEEREMskxLFRHk8oiIgMaywRxH+ohiIhklrtnECggiMhCG8twy9gmiPvQkJGIiABgk6ONF8P6+rrrCE0RgWNb5Dt3bAfyt8jH2DMws/3uvt50nXoIIiICqIcgIgtujC3ytuZVdvUQREQyW7QlqFplJCILbSw9gy6t/bGuSFJAEBEZ2Fgr/CYKCCIiUyir/G++8zBnn35K5WfGlrIiUEAQERnY2aefwpUXnze6Cr+JAoKIyBSmae2PLVAoIIiIzMjYKvwmCggiksWiDac0WYbvoX0IIiI9LNoegzbUQxCRuVqEJZljLNM8ZAsIZnYmcDVwOvAgsMfd/yRXeURE2liEgNZXzh7CEeDl7v4FMzsZ2G9mn3b3L2csk4jM2FjX4MO4Kvscz84WENz9DuCO4u/fM7MbgccACggiksXNdx5uvGbMAW1ao5hDMLOzgCcDnyt5bzewG2BtbW2u5RKR2RljRZruLs7ZM8jRS8keEMzsJOAjwMvc/Zjw7O57gD0wSX895+KJyApIK+ETtzVXjWMMaNPKGhDM7GFMgsH73f2anGURkeXXtrVdl4do2ns3yTkklXOVkQHvBG5096tylUNEZJnnBbrI2UP4GeCFwN+Z2ReLn73K3T+ZsUwisgSqNozNYlx+VmP+OYJSzlVGnwUs1/NFRFKr2jMIsk8qi4gMZdcVewG49/4jm36+c8f2Tf87i4p/lveeF+UyEhERQD0EEVkAZePy8c/C30PPICwbDauF+rba6+YD0rmDnTu2c/Odh7nk6usXtpegHoKItLaMGT7lIeohiMho1Z1XXLaqZ6gVPm1WDsVLVUPKi3vvP8KBjUMLu3xVAUFEGs1qaWXX+7TJNST9KSCIrLCxt2TTVng4vD78LL4mvT59r+9zm+5TNa+xiBQQRKTR0Dt52/Q44rmKm+88fMxwTJNFn+DNQQFBZAWNKe9/G/GcQZ30e00bFLp+bqy/v7YUEESktaEqvLoeR9lyzvh/08+m4nmGe+8/op5CBwoIIiuoyxBQ297DvHoZdb2bsNIoDDFBt8ylY+8pzZoCgohkU1bxNgWrpjmEOCjEk9DSzNwX58yZ9fV137dvX+5iiKyEqqGbsgo6bpGHXcLXXHrBIM+vqtCnfb/s2qbv2vW+Y2Fm+919vek69RBEZJRmnT56ESv2WVMPQURqNS0JhYda1lsMHiyqlLFn/2zKj1R1fZtexNiohyAiMxVW88ziuMlZtt4XbcntPCkgiCyooSuyqpxAdUtCw3xB3FqeZwU79LPa7Ehe5gCigCAinZSt84d2PYWm1nnX1nuf3EZ9KvZlDgIxBQSRBTP0kEd6v11X7D26Yqjt8s34/Xn2DOJyzuvZyxwUFBBEVsgQFWecRuLEbVs7rfVvap2XvR/OYIivTXsG0/QU6qzafIMCgsiC6Trk0XRdVSVcl+M/98avMDwVKuo+E9tyLAUEkRUQ8vmUHTTTt0Kf5rNNn4t7BnUpKrYYnHD81pntOViFieSYAoLIgmpbCcaZP9veL14t1OfZXdRVtk1lPuH4reodDEgBQWSBDNVSLTtfYF6VfFtVQ1NpoAs/a7OZrG8yv2XvGQQKCCIjNk3FWnbaWJszBco+37cMVcrKtOuKvUcr/qHPNZB2FBBEBjLLceZdV+wF2i2zLHsvrWi7biBrWyF3aZHf98CR1quD4mGhPsdbdukZrMqKojJZA4KZvQu4CLjL3c/JWRaRMUkrpy3W/15DJIlre2xlk/R73ffAkaP5j+JgN23vZBUr8yHk7iG8G3grcHXmcoj0No+WZUgYFwJD3/HyNuPi8efa9kyqdhw3TRQ/OEVuzS6/3zbXrtqKojJZA4K7X2dmZ+Usg8iYpJVRWiHncN8D/Z5dNsxUdqLZzh3bOXjLodLlo317Bqs87DON3D2ERma2G9gNsLa2lrk0Isdqe8JXn0opjJ3Hm8TK7tl3bqCsnGUpLNp+pqwcqbLVQyHwNT1HqSlma/QBwd33AHtgch5C5uKIzERTy7apwmz7jC6b0cp6BuHz6f3Sa6o2lAVxULjk6usHy0mkYZ/pjD4giCyKtPLpsjKoSTiOsm120PD3PuPs8T0O3nKodpw/beF3mYBuc+2shoAUMMopIIiMQF3LdohspvGY/YGNQ5vW/NeVadcVe7nvgSNHg8LZp59S2gMIvYmuK4Rm1aJXRd9P7mWnHwCeCZxqZrcBr3X3d+Ysk8i00tTM4RCZISqppuygQZqzqO8zrrn0gmOGmtJnxb2IoTeQDR0wNOlcL/cqo1/L+XyRsanbXdxmKWdZXp+4Ii8b829TpjgIpGWM5xrCgTl9ewqSl4aMRAaWY2IzVPLTPDtt+TftWQhBKJ1jmEWyubH2OJaNAoLIiJQlbqs7hKZuUjlVliNoiCWx8W7qeC/BNM9ShZ2HAoLIjLTZCTyvZ8eBouy0sXiSOP1M21b1CccvTnWiQFNucf4FRWZkTK3RrkMafa5P5xK6JJmrE69AuvLi83oNz2jSNy8FBJE5mcUQSjpE1LRPId59HIJAaNmHoamgbuinbNVR2z0Idd9liMAk/SkgyMoaojU6qxZs1/tVLT+tE+8vCPqsQgqfq9qd3OW7xKk6cp3XvMoUEGSpjWnIYYghlLChLIh/DuW7oi+5+vpNE9NxS7wqV1HZCqP0c/Fpa22CSNMEeDqHoUNx5k8BQVbWEEs05z3W3TfzaJAO91zwhms3vR/2EZQJu53L3k9zLg3xe9BZyfOngCBLacyTk13zC4UW/r33T4Z4QoUcdhHH16Vj++nqobodx0FVRRwmjcvmGdp+F2g/HzKGf6tVo4AgK69PxTOvDU6XXH09B2/ZPJTS92yE0MI/sHGI57zx2sqkdTffeZhdV+zdlFAvzYUU9hzE5YzLNkQ6bpk/BQRZStOM18+ycu9z/xOO37ophUV8v9BbSHsCaeqINA3GNCeVxWXqo+r7K2Dkp4AgMoVZB480CMRCHqF0Q1jZSqGqfEZxYrpwPGfoGQRlgaVNkr2qlUfpZ8q+t4JDHgoIstTaVDpN8w1VKanrlkaWDZlMM58RegLhQPogPpw+ZFVNyxTKEK9EOrBxiJ07th+95zlr2zvtAWj6nchiUkAQ6WmoHb5l4lU79z1wZFMlG3oG8RBQOsZfFXTSoadz1rYf/XvaM6gqEzTvdyjLxZS+N3TQlOkpIMjK6NoTqNqhCw+1yu+9/8imw2bqntF1PiOuVNMKNiSRS4do0gnnrsnyysoRl7lu13IcnGQxKSCIdJTuBejaU5i29ZueaRzmEcoO5KlKBzHLTV/x8tS28wRKSz0OCgiyMpoqnarX6War9FjJeFK37hlt00uE58XCOH96n1D5Vm1YK9ufkJah7v2u8yrxM2XxKCDIQhlDCzJUwqFF3nYSddpx8niJaXyf2Bajcv9A2F+QnsIW9yCGqsjrvlPf92T2FBBk5bStdOo2W0F9moZphkBChR56IuesbW/d6n7Qj813VCdedhpvOIuDSln5676zLC4FBFkIXVrXTcMZ04jH3psq3XQZZtvUDenP0/X8Ybnozh3bjw4XVWUtDcEs3WNQleAuLGPVkM9qUkAQqVBWge+6Yu/RirQuuMTDNfEyUGiXYjpcHyrouknrsgnlVFlFH6+MOnjLoU33CdTqXy3mPuUe9jlaX1/3ffv25S6GzEnXln7aiwgt6HjMPLSuu2Y2TSvM9D5lSebiyrXq8Jm6cgObnle2cid9Vro57TlvnGQzDT2IsvvEPZk4cd5QxjDvs+rMbL+7rzddpx6CLKR5VjKhYo1X8qQJ5+pa8GG4BtoNeaXDQHElXyW+5sDGIQ5sHGLXFXuPuUeZuAxlK5yqqKJfPgoIMjp18wV1SzfbLonsU5Z0jD4sNa2aeIaH0kyEz1YFjTjFdXxiWNPwTdo72GL1gSnsl0iHseIJ66oJ9D60+3jxKCDIQpm2kum71DMVKs60pxCv2om1PdgmnQyu2lVcNiyVZiBt2sncN++QKvrllTUgmNmFwJ8AxwH/xd0vz1keGYe+Szbrrk/Hy9vu0o1b7HDsWQShpxB+fsLxW48Z7okzkjalljhx29baVn7V6p84cFQN+8RzKVtscxK8WewU7nJPBZVxyBYQzOw44G3ALwK3AZ83s4+7+5dzlUnGb9pgEefbKQsKaYUbV6CxdBlnnIQu/kwIDnEgCEM7Za3z9Gd1ZwfEqSvia+smr+OT10K5ux5mrzQTyytnD+Fc4CZ3/xqAmX0QeC6ggLDipq3s69bzB2EpZ9ueQhiOSXsK4Vnxe+nn6k44q1qlFFY1VX3HXVfsPWaPwViG01JtegYafhqHnAHhMcCt0evbgKelF5nZbmA3wNra2nxKJqPXtcIoWw4at6TLln7CQ0M+TcMxVecCpAnmQuu8bH9CrMuJZG32KsRlDAEnnWDuOqegcw+WT7Z9CGb2POACd39x8fqFwLnu/ntVn9E+hOVWVSn3aeHW7QauWs+fPrts+Wea5rpqP0RdQIjF37Fu/0LVJrn4nuFeQdXvrWxfRRhKSstU9fn0d1V3fRvqGcxW230IW2pu8EkzO2vIQiVuA86MXp8B3D7D54kcrXDS1nT4+RbbvHwzXjHUphUeL9+MewDXXHoB11x6ATt3bN/0p00G1DC0lTr79FM4+/RTNt2rqUJNd0CfuG3r0Qlmkboho3cDnzKz9wBXuPsPBn7254HHm9ljgW8Azwd+feBnyAKZdrIytLLjvD9l92pb+ZUdbB/uOYuWbNnO4XTpaJt02l3Llh620/T5Wa5IkrwqA4K7f9jMrgVeA+wzs/cCD0bvXzXNg939iJn9LrCXybLTd7n7l6a5p0iqLLVz2aHv6VBN6A2kZxi3nYyuGkZJh6nKgldcpqbgVva6jlYISZ2mSeUfAPcC24CTiQLCENz9k8Anh7ynLL5pKqmyVnPVkEsb8UqfuvMPhswOmvYUYrM86axPj0yWS2VAKDaNXQV8HHiKu39/bqUSKdGlVZtO4taN/ZdVwGWTzmEitir/UPy8dIK36bnx6/j9uuGjaagylzJ1PYTLgOdpGEcWSdW+gHvvP3LM5rKu90xTSIeVPmXnEUwrHd46eMuhTXsPhszXJBLUzSH83DwLIlKlz+alsmWcVQfRx58pe24XZfet07YSj4euFoGC02JarP+XiVSoChpxTyFeNRQyf7ZRtZs4brHDQzuGh1jCWTWcVNcz0G5fmZYCgozeECtjuozB161Gavucvs9Y9Ep8Wb/XqlBAkLlqW0FU7SiuUtWKT++RTv6G3b5tTgir2gPQpZxdxXsrqjaeaSmpDEUBQRZGmzQR6c/DjuFQ8YfhnLINZ+lzZlnBdn3GLJebDknBabEpIKywef5H23YoIZ0IDpV5l55CWU6gcCZBuEeaB2ia3dFl33OI32nblN1VZRHpSgFBFlJd4CiTLhlto0sF23TfvruMp0nZndPYyyflFBBW0Cwn/poqvjbPik8p6zM2H/YHpGkn4nuGOYNpvnv82TQz6hC/07LejpLQySwpIMjcxOP9Q1SY8TLSNE10utwUygNM2+ylVeUMLfaqQ3CmDb5Vk+Mis6CAsIKGmvgrO3+gqeKrOmimLMd+m7MHqsoUB4qDt0z+/qCzKSX1lRefN8gh8+H+QdiPMFSqiSFTVojUUUCQTrqsigmVbVqBlqVg6Kpp41YqPcqyKZNol5b9Ccdv3TRcFO8qbht8m95Xz0DmQQFhhfWtZOoq0zbDG+FQ+HAvoHRMv83ZBlXfKZQjrC5KU1ZMMxYfyrbFJpX/NZdeoGEdWQoKCNJK2RLIsmvSc3qhfBNXvA9g2sq5TDzJWza+Hyr1dO4hvW8Y2mqq5JuGdZp6BtrZK2OggCC9hEnb9LzfqlPG4s/FWULTyjjdmRvuGwJQmzmFODtpELfmp5Gu+rn3/iOd9kmIjJkCgrRSlRqi6Zq698smeVNxMCg76KZPizptlYeflSWMm3auo4l29sqYKCBIb0MtGQ3Khk/CoTShYj54y6FNa/3jyev4HqEHEz53zlrzgTVtzCuHkUgO5j7wyR4ztL6+7vv27ctdDOmgqfUfv5cGhLKNZbETt20tPaEsDhJ1lXabVnlTig0FBFkEZrbf3debrtsyj8LI4gjj9/MSDwOFbJ6hck/nAIITt21li22ejL75zsPcfOdhDmwcOpriAWa3s7fvHoYm8/79i8Q0ZCQzUTb8E69MqjqfIL0uXiaaDg+lny/L+xM+W6ZLWm2RVaCAIMD8lz9W7WVId/2GIaGm1UVVeX/anoqWajPUVbUPo+pzbZ6n5aeSkwKCzEQ6+Zqmo4ZjK/8y8VBR2SqgsiWraU+haSWTiEwoIAjQfSPWEELKh/icgrSn0HV1UFlQ6KLNOQlVS3C77qpOy93ncyJDUkBYUVUrfNocGNOl0kpb9XE6iWsuveBoBVy2NyB+nc411FW8VT2JdMObiGymgCBzk6aTiJPcHdg4dMxKIihPOzELaUCsGqqKNfUa+lCwkpyyBAQzex7wOuCngXPdXZsL5qRp8rLLZGrXlA115xmH1BLhurryxInl2uYYKsuwWvfZuv0PIssq1z6Eg8Au4LpMzx+9ZV6PvnPH9qN7DUJLPPQa4mMihxJ6JnVCcAg9kyBMSLctT9hL0dcy/7vL+GXpIbj7jQBm1nSpDKxtvqGqz4YMoaHyDvcacqgj3WsQ9Ekslya6CxV+m9TcQbr/QWRZjX4Owcx2A7sB1tbWMpdm9sa+Hr1NazuVfqfQOwiZR8PEcqh0h/quaQCB8go/KNvLUJUee2hj/3eX1TCzgGBmnwFOL3nrMnf/WNv7uPseYA9MchkNVLyV17eiSVvKbe7TZgnofQ8caWyJ90ksV3fuctUz4mWrylUkq2RmAcHdz5/VvZdZjv0AbaQt2HSsvU5ZEImXgIYTx8J7Q2kKIHXDZvNuoWsfgozB6IeMVk3ZfoBFzarZJoj0GSrp83vo+vtbtN+1yBByLTv9ZeA/AKcB15rZF919uqOsZGbqDsWpuzYVcgvFO3q79DT6KktdDeMcrx9DGWR15Vpl9FHgozmePXZlx0Z2zcczywouLlPTs9qcslY3JzGmilpkFWjISCpVtaTbnnYGzcdopq/nQeP1IuUUEEaqT0t5HkMh6TkGbSr0NnmEyr7vGId0RJaZAsKK6jJ5Gw9fNVXobSvvMVTuYyiDyJgoICyAPimU47X/s8jw2fZ+XVr7XfIqicjwFBBWTNfhmHgZbJuJbVXeIotLAWHJpBV4+vN5m0cWVREZhgLCiuk71r+Kwzer+J1ltSkgLJmqVT+5K7UxZFEVkXoKCCMyzxZp12c0zTEs04YyLXuVVaWAsKS6Vl65K71VHpoSGQsFhIFMU5HVtUjHWkGWlTkk4VvklvUi/O5FZkUBYcl0rcTGNjyiylckHwWEKQ1RoVYlgYuzguauqFN1QzxjSMLX1dgCo0gOCggZ9Klsmk4d61uhrcLY/TJ/N5EhKSBMacgKNf5sjiygfdSlva4yxtb4KgRGkSYKCHPUpyJse/7AtBXaMlaAYww8ImOmgDCQWVUy4QCZ9PyBRTbm1viYyiIyb+buucvQ2vr6uu/bty93MabWpyLskiV0Ucyr3Iv6+xEZipntd/f1puvUQ5BsVEGLjIt6CBUWoVWZjpHv3LEdGHeZRWT+2vYQtsyjMFIu7DUQERkDDRklFmFlik4WE5FZUEDIYBGCjoisHgWExJhb3YseSBatvCKrRgEhgzEHHRFZXVpllFHfU80WLZBoNZRIXqNeZWRmV5rZV8zsgJl91MwemaMcuV158XmqFEVkNLL0EMzs2cBfufsRM3szgLu/oulzy9pDmEXLeYy9iDGWSWQVjLqH4O6fcvcjxcsbgDNylENERB6SfQ7BzP478CF3f1/F+7uB3QBra2tP3djYmGfx5mLIlrPG60UklT2XkZl9Bji95K3L3P1jxTWXAUeA91fdx933AHtgMmQ0VPk0fCEistnMAoK7n1/3vpm9CLgIeJbn7qZkNmRQqlrSqgAoIk2y7EMwswuBVwDPcPfvz/PZi765S0RkVnJtTHsrsA34tJkB3ODuL81UlsGMKbikPQMFQBFpkiUguPvjcjwXNIQiIlJFqSsGMOZWuAKgiLS1sgFBFaOIyGbZ9yF0MfadymqFi8gYjXqnsoiIjM/KDhnNgnoGIrLI1EMQERFAAWHhXXL19cecqyAi0ocCgoiIAJpDyGbaFUlj3vsgIotJPQQREQG0D2Hu2pxX0KW1r56BiDTRPgQREelEPYRMylr2Ou1MRGZBPYSElmeKiNTTKqNMylr9ykwqIjktfUDQ8kwRkXaWPiAsIgUrEclh6QOChmFERNpZmUllERGpt/Q9hEA9AxGReuohiIgIoIAgIiIFBQQREQEUEEZDO6lFJDcFBBERATKtMjKzNwDPBR4E7gJ+w91vz1GW3LSTWkTGIlcP4Up33+nuTwI+AbwmUzlERKSQpYfg7oejlycCi5ODe2DaSS0iY5FtY5qZvQm4GPgH4J/mKoeIiEzM7IAcM/sMcHrJW5e5+8ei614JPNzdX1txn93AboC1tbWnbmxszKK4IiJLq+0BOdlPTDOzHcC17n5O07XLdGKaiMi8jPrENDN7fPTyl4Cv5CiHiIg8JNccwuVm9gQmy043gJdmKoeIiBRyrTL6lRzPFRGRatqpLCIigAIqIniCAAADRklEQVSCiIgUFBBERAQYwbLTLszsW0wmoZfFqcC3cxdijvR9l5u+73jtcPfTmi5aqICwbMxsX5u1wctC33e56fsuPg0ZiYgIoIAgIiIFBYS89uQuwJzp+y43fd8FpzkEEREB1EMQEZGCAoKIiAAKCNmZ2ZVm9hUzO2BmHzWzR+Yu0yyZ2fPM7Etm9qCZLdWSvcDMLjSzr5rZTWb273KXZ9bM7F1mdpeZHcxdllkzszPN7H+a2Y3F/49/P3eZhqSAkN+ngXPcfSfwf4BXZi7PrB0EdgHX5S7ILJjZccDbgOcATwR+zcyemLdUM/du4MLchZiTI8DL3f2ngacDv7NM/74KCJm5+6fc/Ujx8gbgjJzlmTV3v9Hdv5q7HDN0LnCTu3/N3R8APgg8N3OZZsrdrwMO5S7HPLj7He7+heLv3wNuBB6Tt1TDUUAYl98E/iJ3IWQqjwFujV7fxhJVGPIQMzsLeDLwubwlGU6uA3JWSpvzpc3sMibd0ffPs2yz0PY87SVlJT/T2u4lY2YnAR8BXubuh3OXZygKCHPg7ufXvW9mLwIuAp7lS7AxpOn7LrnbgDOj12cAt2cqi8yAmT2MSTB4v7tfk7s8Q9KQUWZmdiHwCuCX3P37ucsjU/s88Hgze6yZHQ88H/h45jLJQMzMgHcCN7r7VbnLMzQFhPzeCpwMfNrMvmhm/yl3gWbJzH7ZzG4DzgOuNbO9ucs0pGKBwO8Ce5lMOH7Y3b+Ut1SzZWYfAK4HnmBmt5nZv8ldphn6GeCFwC8U/71+0cz+We5CDUWpK0REBFAPQURECgoIIiICKCCIiEhBAUFERAAFBBERKSggiPRUZL78ezPbXrx+VPF6R+6yifShgCDSk7vfCrwduLz40eXAHnffyFcqkf60D0FkCkUag/3Au4CXAE8uspyKLBzlMhKZgrv/wMwuAf4SeLaCgSwyDRmJTO85wB3AObkLIjINBQSRKZjZk4BfZHJ61h+Y2Y9nLpJIbwoIIj0VmS/fziQn/i3AlcC/z1sqkf4UEET6ewlwi7t/unj9H4GfMrNnZCyTSG9aZSQiIoB6CCIiUlBAEBERQAFBREQKCggiIgIoIIiISEEBQUREAAUEEREp/H+CJu4anYEAdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.mine import MINE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "lr = 1e-5              # learning rate\n",
    "ma_rate = 0.1          # rate of moving average in the gradient estimate \n",
    "\n",
    "mine_list = []\n",
    "for i in range(rep):\n",
    "    mine_list.append(MINE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                          batch_size=batch_size,lr=lr,ma_rate=ma_rate))\n",
    "dXY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Previous results loaded.\n"
     ]
    }
   ],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    mine_state_list = checkpoint['mine_state_list']\n",
    "    for i in range(rep):\n",
    "        mine_list[i].load_state_dict(mine_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEICAYAAAB/Dx7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4FNX6B/Dvm5AQIAEChF5Cr0qLFJEqIEVUELsieBW7XuWq2L12r8pVr3IVuaBYEFQsP0GK0pHeOwQIhJYChCRAEpKc3x8zu8zuzu7O7k7b3ffzPHmyOzs7c2Zn950zp5IQAowxxsJHjNUJYIwxFhgO3IwxFmY4cDPGWJjhwM0YY2GGAzdjjIUZDtyMMRZmOHCHESJaSkT3WpyG1kS0mYgKiOgxDeu/QkRfy48bE1EhEcUan9LIRUR3ENFCq9PBrMOB22aIKIOILsgBLouIphNRYoDbSCUiQUQVDEji0wCWCiGShBAfBfJGIcQRIUSiEKLMgHRFJLVzKYT4Rggx2KD9WZ45YP5x4LanEUKIRABdAFwB4AWL06PUBMBOqxOhZNAFijHb4sBtY0KIYwB+B9DB/TUiiiGiF4joMBFlE9EMIqomv7xc/p8n59x7ElELIlpGRGeJKJeIZnnbLxFdR0Q7iShPzoG1lZcvBtAfwMfydlupvLepvJ8CIloEoJbiNWfukYhuJaINbu99goh+lR9XJKL3iOiIfOfxKRFVkl/rR0RHiegZIjoJYLq8/GkiOkFEx4noXnlfLQLY3gT5szxBROMU6apERO/Ln/VZIlqpeG8PIvpL/qy2ElE/H59rfSL6kYhyiOiQsqiJiLoR0QYiypfTN8nHuRxLRCsV7xVE9BAR7Zc/99eIqDkRrZa3N5uI4uV1k4noNzkNZ+THDeXX3gDQW3F+P5aXtyGiRUR0moj2EtHNin0PI6Jd8n6PEdE/vB0/05EQgv9s9AcgA8BA+XEjSLnb1+TnSwHcKz++B0A6gGYAEgHMAfCV/FoqAAGggmK7MwE8D+linQDgKi/7bwXgHIBBAOIgFY2kA4h3T4OX968GMAlARQB9ABQA+No9XQAqy6+1VLx3PYBb5ccfAPgVQA0ASQD+D8Bb8mv9AJQCeEfeTyUAQwCcBNBe3vZX8r5aBLC9V+VjHgbgPIBk+fVP5ONuACAWwJXyfhsAOCWvHyN/ZqcApKh8LjEANgJ4CUC8fN4OArhG8bndJT9OBNDDx7kcC2Cl4rmQj62qfPzFAP6U91ENwC4Ad8vr1gRwo/wZJQH4HsDPim25nF8AVQBkAhgnn7cuAHIBtJdfPwGgt/w4GUAXq39D0fBneQL4z+2ESIG7EEAegMMAJgOoJL/m/FHJP8yHFO9rDeCi/ONS+7HPADAFQEM/+38RwGzF8xgAxwD0c0+DynsbywGwimLZt1AJ3PLzrwG8JD9uCSmQVwZAkC4ezRXb6QngkPy4H4ASAAmK16dBDsTy8xbyvlpo3N4Ft88rG0AP+fgvAOiocrzPQL5YKpYtcARJt+XdARxxW/YsgOny4+UA/gmglts6audyLDwDdy/F840AnlE8fx/AB17OWScAZxTPXc4vgFsArHB7z2cAXpYfHwFwP4CqVv92oumPi0rs6QYhRHUhRBMhxENCiAsq69SHFNgdDkMK2nW8bPNpSAFsnVwMco+X9Vy2K4Qoh5TjaqAh3fUhBYFzbuny5lsAt8mPb4eU8zsPIAVSAN8oF0HkAZgvL3fIEUIUue07U/Fc+VjL9k4JIUoVz89DyvnWgnSHckAl/U0A3OTYprzdqwDU87Jufbd1n8Ol8/U3SHc7e4hoPRFdq7INX7IUjy+oPE8EACKqTESfycU++ZAuGNXJe0ufJgC6u6X7DgB15ddvhHTHcVguIusZYLpZELhSJ3wdh/SjcnDkdrOgEmSFECcB3AcARHQVgD+IaLkQIl1lu5c5nhARQSqyOaYhTScAJBNRFUXwbgwpR6hmIYBaRNQJUgB/Ql6eCynYtBdSOb8a922eANBQ8byR4rGW7XmTC6AIQHMAW91ey4SU475Pw3YyIeXwW6q9KITYD+A2IooBMArAD0RUE94/u2BNgHR31l0IcVL+7DdDuqhDZX+ZAJYJIQZ5Sfd6ANcTURyARwDMhutnzwzAOe7wNRPAEyRVBiYCeBPALDnXmAOgHFIZJwCAiG5yVEIBOAPpB6rWLG82gOFEdLX8Y5wAqcz0L38JEkIcBrABwD+JKF6+QIzwsX4pgB8AvAup7HmRvLwcwOcA/k1EteX0NyCia3zsfjaAcUTUlogqQypLduwnmO0p3zsNwCS5cjFWriCsCKmoZwQRXSMvT5ArOhuqbGodgHySKlQryet3IKIr5PTcSUQp8v7y5PeUQeVchigJ0kUsj4hqAHjZ7fUst339BqAVEd1FRHHy3xXy5xxPUpvyakKIiwDyof6dYjrjwB2+pkGqgFsO4BCkXOGjACAXN7wBYJV8e9sDUrPCtURUCKki63EhxCH3jQoh9gK4E8B/IOU2R0BqnliiMV23QyrPPQ0pKMzws/63AAYC+N6tqOIZSJWia+Rb+j8g5RRVCSF+B/ARgCXy+1bLLxUHsz03/wCwHVLl6WlIlaIxQohMANdDKvLIgZQ7fQoqvyshtV0fAalM+RCkz3YqpMpDQKpc3Smfnw8hVdIWeTmXofgAUmVuLoA1kIqMlD4EMFpucfKREKIAwGAAt0K6GzuJS5XCAHAXgAz5M30A0neHGYyE4IkUWOQhqQnjDgAV3S4IjIU9znGziEFEI+Xb92RIucL/46DNIhEHbhZJ7odUZHEAUlnrg9YmhzFjcFEJY4yFGc5xM8ZYmDGkHXetWrVEamqqEZtmjLGItHHjxlwhRIr/NQ0K3KmpqdiwYYP/FRljjAEAiMhXL2MXXFTCGGNhhgM3Y4yFGQ7cjDEWZjQFbiKqTkQ/ENEeItrNI4Axxph1tFZOfghgvhBitDyTRmUD08QYY8wHv4GbiKpCmslkLADIgw1pHXCIMcaYzrQUlTSD1I14OhFtJqKpRFTFfSUiGk/SnHkbcnJydE8oY4wxiZbA7Zhn7r9CiM6QpoCa6L6SEGKKECJNCJGWkqKpDTljLIJknj6PZfs402YGLYH7KICjQoi18vMfIAVyxhhz6vfeUtw9bZ3VyYgKfgO3POVVJhE5Bp2/GtKs0Ywx5lRWzgPWmUVrq5JHAXwjtyg5CGCccUlijDHmi6bALYTYAiDN4LQwxphmeedLULFCLCrFe5ugPnJxz0nGmF+3f74G13280upkuOj06iIM+XC51cmwhCGjAzLGzHPybBGKLpYhtZZHK13d/HXglGHbDsXhU+etToIlOHAzFuZ6vPUnACDj7eEWp4SZhQM3Yyys7D6Rj/wLFy1NQ37RRZwvLkPdagmW7J/LuBljpsg8fR6rgyhy2XU8H93f/AOnz0kjbQz9cAVumbJG7+QFZPCk5c47HStw4GaMmaL3v5bgts8DD7iTl6YjK78YK9NzdUlHRu45HMwpxPsL96LoYllQ2ziZX6RLWoLFRSWMRZBJC/diSId6aFe/akjbmb0+E2/M243NLw5CTAwFtY3s/CIUXSxH45r2Gky033tLnY+rVKyAB/o2ty4xQeIcN2MRori0DB8tTseN//1L0/o/bT6KgiL1suIXft6Bsxcu4mJ5uaZtLdmb7bGs25t/os+7SzyWnyos1rRNh70nCwAAx85cCOh9DkfPnMeb83ajXKVnZ0mp6/HNWn8E2T5y06cKi3GhxDWXnl1QhBd+3o7ZGzIhhDm9RzlwMxYhHDGjXEPw2HY0D0/M2ornftqhy77HTV+ved0bJq/Cn7uz8OVfGZrW359dCABYtOukz/XmbjuBsyqVlo98uxlTlh/EzuP5Pt+fnV+EZ37cjnFfeB5LaVk5hBDo+vofHu3Zn5i1BV+vOYKnf9iGBTuz/B2OLjhwMxYhHAE7hvwXbZyXc41ZQZbV3j1tHYZ9uEL1tT93+w5emacv4G9fbsDLv+4MaJ+bjuTh7d/3qL52KPccHv52E56ctcXjtVKNdw2lco7cUQnqkF1QhBbP/44Zq6VJ2B0XEodV6ZcqXPO93MHojQM3Y2Gm1fO/4/6vNngsd5QEaIjbIVu2Lwe7TqjnYGetzzRsv58uO+Cx7LXfduF8SSkA4FhecMUp7oQQOHteCsKZp6Vt/rzlmMd6F8u0XRT0xoGbsTBTUlaueksuAshxqzmWdwFjp69DYXGpYpue63krF3dw331OQTEmLdoXcHp6vvWnpvf9b+UhZJ6+1IMyK7/IWS7uj9rxCQF8v+EoOr66EHtOXro4HVe5KExZftDl+bztJzTtN1TcqoSxCOHIcSsDr1+KwPX+wr1YujcH83f4Lks+Veh75kKCa+QeM20ddnvJnW8+cgbFpeXo0aymx2snzhbhoz/3+9yXw3lFhWH3N6X21e/f1BE3dm3oDM6nzvmuFFVecJbtlyaEWLo3xzlcbVa+5/vdP4s9J7RdMELFOW7GbCo9u0BzBR4AlxYNv2077vHarPVHnLllrXlytcy7v6pP9/coc63uRk7+C7eqdKZZczCwjjpPzt7qsWzC91vx8eL9zkrJh7/Z5LFO3oUSfLv2CO6cuhYr9kntxE/mF2HuNinn/Pbve/Dugr1e93sw17W824xiKoADNwtzZ89fROrEuVi8x5za/GDNWJ2BbUfznM+FENiQcdpn87Fr/7PSowJvtI+mfspNbTmS5/Lalsw8PPPjdjw7Z3tgCQ+Ce/DS0kLuzLkSl+Z6asE8GO8tvFTUcq7Es7PN9FUZeO6n7ViZnounf9wW8PaX7nWdqs2k1oAcuFl425sl3Zr+d6lnpVUghBABty8OxEu/7MR1H69yPp+/4yRGf7oaM9dJFXmHcs+huNQ1sBRd9Kz42nD4jNd9+GoGeEHuIZhbWIzl+3JcymK3ZOa5XFSy8otQEkKlm3tRiRadX1uEZs/Nw6Yj3o9PD4MmLTN0+2b1qOTAzaLC7hP5PnO3n684iK6v/4HDp86Zkp4jcmXapiNncPb8RfR/b2nIuWFl/xJvt+z5F0oxZto6fCk3bQOAGz5Zhes+XuWsfPNVNAC4FslcKCnzuNuZu/0EJi9NDzD1kpX79enW7o17U75wxYGbRby/DuRi6Icr8NWaw17XWbJHuuUNtneeVrluPe9+2HgU5+SmbIEOwPTegr0uQVSolD7/uvU4rv9kFfLkpm3emvABwJqDpz2WuV/rrn5/qbO9MyD1sLznC8+miYGUzetpj8bWJOGOAzeLeI7B9nf56DmnFvSMkPb6H7hlympd9vbxknTkF6k33cu/UIqs/CI8NnMztmbm4SGVijkAWJfhGax9OZBzDumKXGuGSXcozBUHbhYRQq0UcuY2TWgVsO3o2YDfsy+rAKkT5/pcR1nGPWtDprNZXCh+3HQUK/a7VsBpqU9QazrH9MOBm0WN79ZnqnYe+dpHEYoZzvlpd710bzZ+2uzZaw+AS9s8I1o0PP/TDmS4TQ+mZSwUZiwO3MxSh0+d06XbsNb2s2/O2+3yfPeJfLzw86WBltRaRFz/8Urc/Olqj+UFRRdx9Iz3OQ+3Zubh163Hvb7uMOjf0oS3Qkgj/K075Fp8MXb6er8DJAFw6T2opxd/dh2ISktamLG45ySzTG5hMfq+uxR3dG+MN0ZeFvD7X/l1p6bAqKTsYZeeXeDRvE7tArDVrWhDCIGNh8/grv+tw4WLZV7nerz+E6n533Ud63u85i3T+tpvu/D1miNY+EQfl+UZueplycqyeStmhdnoo3lisMxq2RPONAVuIsoAUACgDECpECLNyESx6OCYNzDYGcS/0NhywX3MZYeBk5Z7LFPG7bnbTuDhbz0r9e7/aiMW7rrUBE4Igf8sTset3RqhdlJocxA6xthwtAJxOGJQbtpuvll7OKhxTaJNIEUl/YUQnThos2A9O0cabN6dHoPP+9pEsBPLTl15UHW5MmgDUmXjpEX78PfvtiA7vwhfrc5wacGitSjoZH5RwOXU78z33eY63HClpjZcVMJMM3PdEcxcdwQ3pzUCAJDOAzsUFpfivQV7MXFoGyTExaqu4wiM3sahJiKsPXgK1SrHQWXCFCza5dm13tGu+cz5i+gmt+SomnDpp6U2qpy35oeOopubP/MsU1czc90RTeuxyKI1cAsAC4lIAPhMCDHFfQUiGg9gPAA0btxYvxSyiFReLpxdvPVqozB5STq++CsDDapXwn19mvlc11tgLBfC5zgZ78xXH8gfgMsIeMr21W/N8/4exoKhtaiklxCiC4ChAB4moj7uKwghpggh0oQQaSkpKbomkkWet+fvwZAP1GdQ0WLZPte2xURwDr9ZWi5wIKcQ54pLIYTA+4oyU0exxeFT6mXG933p2QvQH383DvN3eg6TWlgUwNCrjLnRFLiFEMfl/9kAfgLQzchEscinLOv2Vq5bXi5w/1cbVIf4fPoHz2E8HTWLAgJXv78MY6evcwZzh9/9jDVd4KNN9cu/7ECpSnn1qMnaJudVmhzioFgsuvkN3ERUhYiSHI8BDAagzwyjLOIdy7vgt/LRW3lvQXEpFuzMwn0z/OeC12eccc784tjd+owz6PzqIo91X/9tl9/tqfly9WGPziiMWUFLjrsOgJVEtBXAOgBzhRDzjU0WiwS7juej19uLMcrHGNKAjxYhARZ+O7piT11xqTWIWg566spDgW2YMZvxWzkphDgIoKMJaWERxtGRYvORPOSd9z3dlUPe+RJUjq+A+AqX8hSOIuRzxaUoKCrFqMmrfDYbO3PenJm2GbMKNwdkpnAva/am06uL0L5+VVzTvi66Na3h8t7Rn672OnchY9GEAzczha822+5FJTuP57uMh3GupAxl5YKDNmMyHmSKmeKsovdiVn6RR5duf5o/N0/vJDEWtjhwM8OsVYxy1/+9pc7HauNEl5SW+x1vmjEm4cDNdFVeLnDzZ6uxdG+25kGghBB4Sq1dNmNMFQdu5tPiPVn4ZEk6ysoFzpdcalqXU1DsHItaCIEnZ23BqvRcFBSVYt2h03hs5mbN+xAATp/T1uqEMcaVk8wPx0Sw+7IK8MuW486xp6944w8AQMbbw1FaLjBn8zHM2XwMW18aDAAoCWByhBNni9CidqLOKWcscnGOmwEAlu/LwT1frPfay/GXLdomLNh5XJp0oOhiYLParNifG9D6jEUzDtxRIj27EAMnLcO2o3lYqzL2x71fbsDiPdmacsqOUf0AadCmKcsv9VS8fepafRLMGPOKi0qixOSl6UjPLsR1H0vTaWW8PRz7sgqQXDkeKUkVNW/HveVHy+d/1zWdjDH/OMcdoU6fK8HUFQchhMDag6cwZ5PnLOGD/73c2UxPbaCnco29HRlj5uLAHaEmzN6C1+fuxtajZ31OIlsoD8J0sUwK0gt3Xprh5R/fcxM9xuyIA3eEcszAojZ+tDvleNf/p5g1fc5mz1w6Y8x6HLgj3OhP/c9dqJyqS+dpIBljBuDAHaUe+XaT6nKSB1FdsifbzOQwxgLAgTvMFRaXah4yVem3bSdUl586V4zl+3Iw7ov1oSaNMWYQDtxhoLC4FI9/txln3LqFl5ULdHh5AV74ebtz2anCYnyvmM8xUOszzmDMtHVBv58xZjxuxx0GvllzGL9sOY46VRPw3LC2zuWl5VLF48x1mQAI9/RKxaB/L7colYwxs3DgDkPL9+WgSc3KmLH6sHPZzHVHMHPdEQtTxRgzCwfuMMRFGYxFNy7jDiPeBoBijEUXDtw2I4TAodxzLst+3HQUAPD5ikNWJIkxZjMcuC02f8dJXP/xSgghUF4u8O26I+j/3lKsU0z7tS+r0PlYOXcjYyw6ceA2WWlZOZ6dsw3H8i4AAB6buRlbj55FSVk5mj03D8//tAMAcPNnqzFtpWcO+/ft6u2vGWPRQ3PgJqJYItpMRL8ZmaBIt/rgKcxcl4mn3eZY3JhxxmPdV3/bhSV7XXswTpyz3WM9xlh0CSTH/TiA3UYlJNo4upY7eJuAYNx07sHIGHOlKXATUUMAwwFMNTY5kc+9YYjaONiMMeaL1hz3BwCeBuB1jFAiGk9EG4hoQ05Oji6JizTZBUUokyO3YxQ+xzjYjDGmld/ATUTXAsgWQmz0tZ4QYooQIk0IkZaSkqJbAiNFbmExur3xJ96dv9fqpDDGwpyWHHcvANcRUQaA7wAMIKKvDU1VmCovFzhfUqr6Wm5hMQBg14l8M5PEGItAfgO3EOJZIURDIUQqgFsBLBZC3Gl4ysLQa3N3od1LC1xmQQeADi8vwPgZPm9YGGNMMx6rREez10vDqZaUlqNihVjn8sLiUufcjkrchZ0xFoyAArcQYimApYakJIKQhvm/VuzPRdNn55mQGsZYpOEct46U+eecgmIkxMWgQgx3TmWM6YsDt44cJR9rDpzCvTM2AACqVYqzMEWMsUjE2cEQlZaVY9LCvSgsLkW5HLl3HD/rfJ0HhWKM6Y1z3CH6ZctxfLQ4HflFpc6ikj0nCixNE2MssnGOO0TFpVJn0qKLl5oAzt950qrkMMaiAAdunazPOA0edoQxZgYuKtHJgZxz/ldijDEdcOAO0r6sAvy5O5tbjTDGTMeBO0Bl5QKLdmXhydlbcL6kzP8bGGNMZxy4A/T+wr2YvPSA1clgjEUxrpwMwNkLFzloM8Ysx4E7AJ8t46DNGLMeB+4AcG6bMWYHHLg1Ki/nRtqMMXvgykkNPvxjP7YfO+t/RcYYMwEHbg3+/cc+q5PAGGNOXFTCGGNhhgM3Y4yFGQ7cfmQXFFmdBMYYc8GB24/0rEKrk8AYYy44cPvw14Fc3D51rdXJYIwxFxy4vcjKL8Ltn3PQZozZDwduFUIIdH/zT6uTwRhjqvwGbiJKIKJ1RLSViHYS0T/NSBhjjDF1WjrgFAMYIIQoJKI4ACuJ6HchxBqD02YZIrI6CYwx5pXfHLeQOJpWxMl/ETtwx7J9OUidONfqZDDGwtDnY9JM2Y+mMm4iiiWiLQCyASwSQnjU2hHReCLaQEQbcnJy9E6naX7ZfMzqJDDGwtSgdnVM2Y+mwC2EKBNCdALQEEA3Iuqgss4UIUSaECItJSVF73SaJuMUT/rLGAtcu3pVTdtXQK1KhBB5AJYCGGJIamxg05E8q5PAmIsrm9e0OglMAzMnDtfSqiSFiKrLjysBGAhgj9EJY4xJHurXwuokMJvR0qqkHoAviSgWUqCfLYT4zdhkWaOMJ0tgNsSNnJg7La1KtgkhOgshLhdCdBBCvGpGwqzw5OwtVieBMQ+C8xO29eGtnSzZL/ecVPhly3Grk8AYY35x4GaMsSBZdTfEU5cB+GXLMUxewjO4M8YCU66I3GbWRXDgBvD4d1y2zeyLKyfty6ocNxeVMFupnVTR6iTYDldOhgczL7BRH7i/XnPY6iQwhXdGX251Ehhz4asDlFXX1KgP3C/9ssPqJNjWiI71Td9nxQpR/5VkOrqqRa2QtzH0sno6pERfUf8r4T433jVPqWJ1EhiAuFh7FnL//HAvq5Pg19f3dte03o8PXhnU9oVF5VhRH7iZzfCF1EPNRHuW+ydXNm9sDqOl1qysed3K8bHOx1xUwhg4brtrVqsK+FMxnq+Lo8/7HW5VYr6feextxliA7NDKJ6oD98JdJ61Ogq3Z4QvKmN0IRTa7b+tLcw/Uq1bJtDREdeCetz16A/eL17azOgmq+GLh3xsjPeYxYTrr3dJ7a5TK8Zf6LSbEXSrvvqlrQ0PTpBTVgTua/e2qpn7XCcceexlvDw/4PS1rJxqQEuOMNjFARKtbr2isab0qiopKMycZj9rAnXe+xOokMJ21CDIA+/q9mTWHoDcC9r0LId/VdlGhQmwMXrtBugNqXEN7y5RQRW3gfuDrjVYngemsfvXgyhhjLLy16Na0RsDv8RYw7+vt/y5KL48OMH5WnvrVEtA8pQqmjTVn5nSt3Ntu39m9Mfa/MRR1qyWYloaoDNx/HcjFmoOnrU6GbdyS1sjqJAQkJakiHujb3Pl8wqBWIW2vc+PqoSbJFszMmT/Yr7nHsv6t9Z0k/MauDfHnhH4Y0Mb8u54UxZg5/q7rRIS4WHNDaVQG7ts/X2t1Emyldyv/3YJn3NMND/f3/LEGo2aVeK+vKWvsP7qts+o6cx68EhOHtnE+79CwGgA/7W196N1S34CjJ7VjskPdg7KCzmHKGH1zxo9d3VLX7QXC152QHTpERWXgtpOfHlLvatu2XlXT0qClrLJPqxQkV/YecAMRr3E8kngvXb0b6VyWGOPj8G0QI1WpVVBaGdBHdWngkuvUY5RHs3OxDp/c3sXn68qek1bhwG2hN0Z2QOfGyc7ndaoqbs8C2E6gXY9TNP6o3AO6XrXmvraivN3XevG6rIGU4x7XK9Xvun8fqJaLsy7i3dmjiWX71pPUw/OSu8L4uIZfbr9Bpdxx4DZBz2bqw0I2ren6ZV/73MCgtt+/Te2A1vdVVGEGrRcAZRtZX2olVkTG28PRr7W2z6F1nSRN6wHG5WIbVK+ERwe0wHUd62Pzi4O8rqdWbE2Q7oCspudnM+Wurs7WGYFa9lS/gNbv2iTZ/0pe3NCpvi1a+XDgNkElb7dWQX7xlz3VDyM7N1BsxnVD4Tz4j/I3YVY+2FcAMqrJ26qJAzBhcGsAQHIQF9LrOtb3aNlhh4DiEGhQb1uvalC59CtSk9GkZmCjWAZ6RuspWosEc66MEHWBO6eg2PR9Ptxf36ZTNarEo4KiYNb9R/L+zR19vt/xA58+7grc2UNbRwPAvuW90cZxx5KUwDMPatGrhesdb6DXNytatfgTdYH727VHTN9nQlzgH3MgOZbYIO9Z+7eujddvuAyt66oXHbhvtlOYN5trmFwJN6k0fQzXC1IwOexwLntWo6Xt/jf39tC0rXCaxMNvSomoEREtIaLdRLSTiB43I2FG8TWw1OVyszIAeO+mjhjVpYHXdQPh7XY72NtwInLJNcS4NYvwt13hlufQ2uOwS+NkrHymv6Z17WjlMwPQoHolj+O3sgOON00U40MrUztrfA/E+moGE6RaNmjiFow3R16m27bmPnYVXru+vW7bM5KWS0wpgAlCiLYAegB4mIjsOUKQ+rRFAAAWeklEQVSRBsfyLqguH9GxPgYrujeP7toQk27uhHmP9TYraQCAz8ek4ekhrW3RVldNw+TKmlpvBMvfjCID24Z+21o1wfw6gGvaX0r32CtT/a4/fewVqsu7e6no1ooIuFdlnBr34gS9BToA05f3dNO0XpWKnsVF08epf3b+tKidhLt6pvpdz/3CbwW/gVsIcUIIsUl+XABgNwB9sqIWyDt/UXX5I/1b4KF+nmXRXisWAxBIEB7Uro5qOpTiY2NCur3Xeoutx7WjVR3P3LzWdtxqpt4deiePV64zP1elbDGRluq7VcOnd/puR+zN/X21dZCKNWEqNPfv2Ls3+a53cdc3hFYzXRoH32okXAT0CyKiVACdAXh0PSSi8US0gYg25OTk6JM6nZ0rLlVd/sLwtmhdNwkxMYReLWriqWtaG5qO+hrGNPBW3LHpxUEegc/9wmDErXSwklRyt9O85CY9GHQY7rk0n61KvLw2qnMDdFQUrelJWWwRyJyGWtvnqwn1LiTQNv7BDgimtOiJPh7L/ntHF1SrZNwdlfJ0vDVKv2KaQGkO3ESUCOBHAH8XQuS7vy6EmCKESBNCpKWkWN/GVE1uoWeLkoNvDsO9vZs5n39zbw/dW4EoZbw9HJc3lCr5KvqotHT/HSTJwcYRlH39nPWY2VqpblV9B89pWkufSYj/MbgVBqi0YZ/kp1WNu2CKpcb2SsUvj1wV+Bs1MnOIUAC4W0PxjV7S3xiKBX/3DLqBaqnSHl/rjOzTx16Bj29XH1JBqw71jblwa6GpPRERxUEK2t8IIeYYmyTjbDx8xmOZe8Weu5qJ+rfbfOfGy3Fli5ro3CjwVhpafs8xMYQWtRORnl0YROpcbX5xkM8LTLCa1qqCQ7nnPJYHUozyyAD1sSy0dtxx8FWZa0T89Fd57L7PFin6jReuloGvUSU+qLs0ZToDuTOooENX9p3/vCak9zs6rT3y7eaQ02IFLa1KCMD/AOwWQkwyPknGeXL21oDfUzUhDv+9I7gyR2+qVY7DmJ6pPnNVgfyMHOuO7toQiyf09bt+IFUryVXiVQcUMsKL17bz6GUazMQI/ngEmSCCs5mdXfxlLgLlfuEIJOjahVqlpC969jR9eUR7tKydiJYq9Tdm0XLp6wXgLgADiGiL/DfM4HTp7sip80G/N9QyYyNybcoyScf2OzashmZy7szIG+1AfueBBIW/XdXUlCKCQMKUET0n43SqHBzZpUHAs/cQAQ+FOMqjsoVMIEZ0rB9Sd/NQpOjY3PGK1BpY9GTfgO/s9KSlVclKIQQJIS4XQnSS/+aZkTg9zd95wmPZr4/00nUfwzWWr2niFsDcg42/CtSGyaFPXGpEDHW0lTcrl6elWVkFk4YH7CJ3YNLrDqZ2UgIWPen/DstdqBWRIzrW97uO2tn9z22d8eOD6qNh6u25YW1CHqfdzsKnq1CILpZ5fpUclYT+hJoL1CPX5gh0ji35a6aoR5O3Ljrkjtw/u/94GWNb9b1Bfm7Ka4KWZmX1TJi5hEABNS11HLmjmCTc5sW02vg+zfGoheN5Gy1qBjtYlZ7r8vzGLgZMuOojzrw16jLn8KPuWtdJwt6sAk2b0noRqeOjJYjW3PiVzfVtneLNS9e2Qw0dB+9x7yDRsnYi9isqah2BPbVmZXw+Jg0tamsfLdAsTWpWxvg+zXBbN2ksmXmP90ZZufpdSvv6VT0uPvGxMSgpK9e0r0DvfcKwSDziRE3g/uvAKZfn/gZi0mJ014b4YeNR5/M6SerBkgjOH6Ca7x/siVOFoU9erPw9JcTF4o2RHfD8Tztc1vnsrq7o4aP33d09m+DL1YdDTos3arnoe9x68lWOj8X5kjLd9ukt6MXEkGqTMi0CDnZ+3jH88nqYu+1ScR4R4blhbZ3P42Jj4K1Ida7JvXvtOERAtImaohKl7gFO0Kr8msYrmjIpO0p8cEsnPD0kuI47VRPiPNo2V3cbmtXXzz6QIoVr2tc1tIOCHkKdYcT984iLjQm6Isn9k21f35iZiT65vYti7svgA+OPD16JeY/1drlQVIqLdc5AHkrIHdqhrstzs9uah8IO3dT1FJWBO9Bu045ih7QmyVjiZdD2Gzo30LWWedLNnVBFJYApfyrf3Nsd3z/QU7d9AqF1Rw9EsDldPSTKTcku91J05YvPikwfjBrX213XJslo53ZxGdSuDu5xG1+meUrgnaDsEqeDGW0z0kTFJ+DegkGtG7YvlzWshp8f7oXvxvdAA8Uwklq/yMF832tUifcoQnDXq0UtXJHq/e4hmGAxuH1d/yt58ceTfdCohrby8w9u6eT1Ncf5CTZQ+Mtd1a2WgJ8f7oW3b7w8uB3A/m2fHef+hwd64l+jgz9O3/tgVomKMu4v/8oIeRudgujlaDY9Ykmw5ZdD2tdFi9pJSK1ZBZmn1UdgVKpSsQK+vKcbSko9K9Bm3NMNC3aeNHSo0aDPJ3kfcqBtvarYfcJjNAhLXdawGipWULlzU5xnLd+bNc9ejV+3HsO87Sc1VW4HWhwZCLUMydgrUw3dp91ERY77lf/b5Xys1xjbgQg65yjUH2vdvlG3tr4HZdK+076tUjConWdnjkY1KruMH2Ml9+PxdXR1q3q/0BgxbEAolMeh5e6hbrUE3Ne7Gf54sg86arjodW9WE3tfH+JznVdG6Dc69CvXtdc8Tkko7DDDOxAlgVtJLfcRTuxQzhhKzt4O6Vez/CnjJohoV68qqibEob08KFGtJPWmjw2TpQpEI4KDo8lj+yDK9R2IKKCmk/5+a2N7+S4KNMo7NwY/qt9jNmkbHhVFJUp9WprTNtmVTaOVihYpiXjqmtYuswFFg8aKGWeUfOWitWpTTwp2T13TGsMvq4c2ddVbprw16jIM7VAXbevp13LFcZG9qmUt/PFk36AqJSPJowNahNRu38pu7koRH7gL3cbgHtIh+Mo3d3YPx8Gkr1rlOEOHtTVDsHNwqpkwuDXa1a+KJ2a5DlAWzF1HXGyMz2KGxIoVMEyv232Vj0CPMbBddkHKx3b/NUgCnRHeriK+qGTvSdfKIr2/YL8+0gtPmjAmgpZ2qFa2dLBTO9lQZi3q0cy1gishLhYjO1/qZRsm8ckvLceRXNne7f2jWcTnuG/5bI2h27+8YXXNY54ESi0YqtWoR0gssdyqiQNQo3I82r403+qkmMrbJVePTM6Kp/vj7AX16QKD0aVxdTzoZ2o/NUM71MOcTccCmrXo1evbex2mwmoRn+MuVXR1DnRmFH/sk8dU5/jdjeps/RShVzaXutlXs3EurkH1SiHl1pUTO4zs3AB3dPc+zIHdVDKo7LZRjcrooCH4aR2rZs5DvVRbIvkzqF0dHHprWEAdv8b0TEVnm85fGfGBW2mUEQNLaWD17XUFEyaHdfC2pycGtcLuV4cENaSo3mWz+nC9bL92Qwd0bZKMxRP6YvWzA/DvWzrZ9kevpm61BHwx7gqPbu1mCXRihGCESzm8FlEVuPVm5tfA0WEkxscZc78DqCCvbGQ39vF9tLW3JgRf9qz3nJeh8HbO7+rRBADQLCUR9aqFPhZ6qPx9N+/o3sRjWb/WtXGr22Bo9/RqqmsrF6aPiC/jNpKRXd4B15YLn49Jw8Gcc357wSld16k+9mcX4qH+zfH1miM+91UrMV51zHJfXrq2HepXNz6oBlrxaZeclZWVxTd0aoBZGzK9zt5095WpuKFzA3T850KX5X1bpeDgm8Pw+KwtGHtlE3RtEj29EcNJRAfu4tJLQ4PqPG2f6ZIS4jT1WFOKi43BxKFtNK279rmBwSTLyeZDd1jKrAGmlN4Y2QHPDWuLOB8T8zp+E1XcZuSJiaGAJrxg5ovoopKpKw45H694ZoBl6bBLDtCX2BgKeG5N6bB8v+eNkR3Qt1UKLouQDj2Jchm9Y0yXJwZKTUHtNuhUhdgYvxXBSQlxeGF4W3w3vkdQ+7DigsQkEZ3jfnfBXufjOknGDVjELmmWUgXL9uU4n7epW1XTvI/hYtLNHTFrfaazzuHxgS3x+EB7dIMOhl3GhGGBiegct1IFH7eMkcIOmb5nh7bFV3+LnEDtrlZiRTzcv0VY3EWxyBXROW67CPYn3iWMmpM5xFeIQe+WKXh0QAsczyvSZZt2uCAxY7SrVxW7FEPhNqheCf3b+J/gOdpx4LaxgUF0NDCbo1jcva34hMHBTePGwoceNx3f3d8DJ89eusCvmmhdXVQ48Vt+QETTiCibiHb4W9dOzikGl7JiKMafH+7lnC0nku+qB7SpjXG9UvHP6zpYnRTbicYbhR8e6ImpY7RPDVg1IQ6tLJzGLlxpKfj9AoDvEdFtaNGuLOfjB/s2N2QfvmrVOzWq7uzmHUqvMLsPr1ohNgYvj2iPFAMrf4O98PW2ZAhfT5F84XaXllojLO4Uw53fiCKEWE5EqcYnRV+vz700600o40/44u8H+frIDri/b7OQpuCafX9Pj6FpvfGXw+PR3pieouh6ZDu6lXET0XgA4wGgcWPrB9fJLSwxfB/+Ks0qVogNadB2QBpW1N/g7VpydNteGRz0DOXuzG6zzJWTjLnSrY2cEGKKECJNCJGWksK1wnZTNSEOleP1rYvmHBcLdzenWTPwXKgiv3EzizrVKklFQo1rqE9HxpjDv0Z3RMbbw61ORsAivjng0n/0M3V/H9/eGafPGV9Mo9SmrlQck+pl3sRo06lRdUwdk4argqycnDW+B26ZYuwEHIyFwm/gJqKZAPoBqEVERwG8LIT4n9EJC0VJabnzcWot4+aYUytbvvby+obtz5ub0xqhff1qmgasjxahtGzo3qymjilhTH9aWpXcZkZC9JR33twcr9WIiIO2DTl6vg7XawJgxmQRWVSyQNGGm4W/YZfVw18HTplaFPTHk31QVu5/PV9a1E4My/JTzbh22jIRGbhf/DmsOnkyP+pVkyZraJZi3hRmoTbjZMxI3KqEMcbCTEQHbscg90bhO0XGmBUiOnCP6ek5ISoL3ZgrU1ErsSKGmlTpxj0nGXMVkWXcDslV4q1OQkRqnpKIDS+ENkdlMPgOhzFJxOW4o60pYDTgDLc9xUTTsIc2E3GB+0DOOauTwAzCccJeaiVW5HNikYgL3MrhXI0SJ8/2cnVbHneYRbcHDBrrnvkWcYG7rFy6sXbMwm2ERvLgRaFMkMBYJLimfV0AQL/WPCKomSIu8mw7ehYAMMOUmca59JVFt06Nqkd271Cbirgct0OSgblhLtZjjFkpYgM3GVhrwvlsxpiVIipwm98UkPPeZjB7qjTG7C6iAvdXqw9bnQRmKGsvlFwBx+wioion31+0DwBQNSGiDovZwJpnr0b1ynFWJyOq8e/6koj8JP41uqPVSWARpq48tCyzxp8T+iK5Mg9h4RCRgXtIh7rG7oCLXBkzVXMTx2IPBxFTxn3G5Al6Ae6CzRizRsQE7s6vLbI6CYwxZoqICdwOiyf0tToJjDFmqIgI3LuO5zsfmzkvITMHVykw5ioiAvewj1aYur8nBklTotWvVsnU/UY7rlNgTBL2rUqO511wPt73+lBT9jmiY32M6FjflH0xxpg7TTluIhpCRHuJKJ2IJhqdKK1Oni3ClW8vBgBUqxSH+AoRcQPBGGM++c1xE1EsgE8ADAJwFMB6IvpVCGH8jAWy0rJyFBaXoqCoFIXFpTh74SLumLrWOfb2pJs7YlSXhmYlhzHGLKWlqKQbgHQhxEEAIKLvAFwPQPfAPWryKly4WA4hBEpKy1FQXIrColJcuFimun6bukl4eUR79GxeU++kMMaYbWkJ3A0AZCqeHwXQ3X0lIhoPYDwANG7cOKjE1K2WgItlAjEExMXGICkhDkkJFZBYUf5LqIAk+X+D6pW4BQljLCppCdxqdfkeLbSEEFMATAGAtLS0oFpwTb6jazBvY4yxqKKlNu8ogEaK5w0BHDcmOYwxxvzRErjXA2hJRE2JKB7ArQB+NTZZjF3StFYVAEDvlrUsTglj9uC3qEQIUUpEjwBYACAWwDQhxE7DU8aYrFWdJKx/fiBqJfKwnowBGjvgCCHmAZhncFoY8yolqaLVSWBR7Jt7uyO3sNjqZDiFfc9JxhgzWq8W9iqm466GjDEWZjhwM8ZYmOHAzRhjYYYDN2OMhRkO3IwxFmY4cDPGWJjhwM0YY2GGAzdjjIUZEkL/qViJKAfA4SDfXgtAro7JsQs+rvDCxxV+wv3YmgghUrSsaEjgDgURbRBCpFmdDr3xcYUXPq7wE8nH5o6LShhjLMxw4GaMsTBjx8A9xeoEGISPK7zwcYWfSD42F7Yr42aMMeabHXPcjDHGfODAzRhjYcY2gZuIhhDRXiJKJ6KJVqdHCyLKIKLtRLSFiDbIy2oQ0SIi2i//T5aXExF9JB/fNiLqotjO3fL6+4nobouOZRoRZRPRDsUy3Y6FiLrKn1W6/F6y8LheIaJj8nnbQkTDFK89K6dxLxFdo1iu+v2U52JdKx/vLHleVqOPqRERLSGi3US0k4gel5dHwvnydmxhfc50J4Sw/A/SXJYHADQDEA9gK4B2VqdLQ7ozANRyW/YvABPlxxMBvCM/HgbgdwAEoAeAtfLyGgAOyv+T5cfJFhxLHwBdAOww4lgArAPQU37P7wCGWnhcrwD4h8q67eTvXkUATeXvZKyv7yeA2QBulR9/CuBBE46pHoAu8uMkAPvktEfC+fJ2bGF9zvT+s0uOuxuAdCHEQSFECYDvAFxvcZqCdT2AL+XHXwK4QbF8hpCsAVCdiOoBuAbAIiHEaSHEGQCLAAwxO9FCiOUATrst1uVY5NeqCiFWC+nXMkOxLUN5OS5vrgfwnRCiWAhxCEA6pO+m6vdTzoUOAPCD/H7lZ2QYIcQJIcQm+XEBgN0AGiAyzpe3Y/MmLM6Z3uwSuBsAyFQ8PwrfJ8suBICFRLSRiMbLy+oIIU4A0pcQQG15ubdjtPOx63UsDeTH7sut9IhcbDDNUaSAwI+rJoA8IUSp23LTEFEqgM4A1iLCzpfbsQERcs70YJfArVZ+Fg7tFHsJIboAGArgYSLq42Ndb8cYjsce6LHY7Rj/C6A5gE4ATgB4X14eVsdFRIkAfgTwdyFEvq9VVZbZ9rgA1WOLiHOmF7sE7qMAGimeNwRw3KK0aCaEOC7/zwbwE6Tbsyz5VhPy/2x5dW/HaOdj1+tYjsqP3ZdbQgiRJYQoE0KUA/gc0nkDAj+uXEjFDhXclhuOiOIgBbZvhBBz5MURcb7Uji0Szpme7BK41wNoKdf2xgO4FcCvFqfJJyKqQkRJjscABgPYASndjtr5uwH8Ij/+FcAYuYa/B4Cz8u3sAgCDiShZvv0bLC+zA12ORX6tgIh6yGWMYxTbMp0juMlGQjpvgHRctxJRRSJqCqAlpEo61e+nXP67BMBo+f3Kz8jI9BOA/wHYLYSYpHgp7M+Xt2ML93OmO6trRx1/kGq+90GqCX7e6vRoSG8zSDXVWwHsdKQZUhnanwD2y/9ryMsJwCfy8W0HkKbY1j2QKlXSAYyz6HhmQroFvQgpt/I3PY8FQBqkH9sBAB9D7rVr0XF9Jad7G6Qffj3F+s/LadwLRUsKb99P+XuwTj7e7wFUNOGYroJ0e78NwBb5b1iEnC9vxxbW50zvP+7yzhhjYcYuRSWMMcY04sDNGGNhhgM3Y4yFGQ7cjDEWZjhwM8ZYmOHAzRhjYYYDN2OMhZn/ByZy1mYbU64IAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 10     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                mine_list[i].step()\n",
    "                dXY_list[i,-1] = mine_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    mine_state_list = [mine_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'mine_state_list': mine_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 4.982193620464953 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = g.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = dXY_list.copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX5x/HPk4UEwk5YBWRxQxQR4oZocQUBtVVb97pUsHVf8Ne6tVbc2tqq1WpF61LFHVtUELQqKIsgKLKqbAHCmgAhCRCynd8f92aYJJNkEjOZZPJ9v17zyp27PmdmMs/cc84915xziIiIlBcX7QBERKRhUoIQEZGQlCBERCQkJQgREQlJCUJEREJSghARkZCUIEREJCQlCBERCUkJQkREQkqIdgDBUlNTXa9evaIdhohIo7Fw4cIs51zHSOy7QSWIXr16sWDBgmiHISLSaJjZukjtW1VMIiISkhKEiIiEpAQhIiIhKUGIiEhIShAiIhKSEoSIiISkBCEiIiEpQYiISEhKECIiEpIShIiIhKQEISIiISlBiIhISEoQIiISkhKEiIiEFNHhvs0sHcgFioEi51xaJI8nIiJ1pz7uB3GKcy6rHo4jIiJ1SFVMIiISUqTPIBzwkZk54Fnn3ITyK5jZWGAswICkJBg2rOwKv/gFXHcd7NkDI0dWPMKVV3qPrCy44IKKy3/zG7jwQtiwAS6/vOLy22+Hs8+G77+Ha6+tuPyee+D002HRIrjllorLH3oIhgyBOXPgrrsqLn/8cRg4EP73P3jggYrLn30WDj0U3n8f/vrXistfeQV69IA334Rnnqm4/J13IDUVXnrJe5Q3dSq0aAFPPw1vvVVx+YwZ3t9HH4UPPii7rHlz+PBDb3r8ePjkk7LLO3SASZO86TvvhLlzyy7v3h1efdWbvuUW7zUMdsghMMH/SIwdCz/8UHb5wIHe6wdw2WWQkVF2+QknwMMPe9Pnnw/bt5ddftppcO+93vRZZ8HevWWXjx4N48Z50+U/d6DPnj573nRD/+xFUKQTxInOuU1m1gn42My+c859HryCnzQmAKS1auUiHI+IiITJnKuf72Qzuw/Ic849Wtk6aWlpTvekFhEJn5ktjFQHoIi1QZhZipm1Kp0GzgSWRup4IiJStyJZxdQZ+I+ZlR7nNefctAgeT0RE6lDEEoRzbg1wVKT2LyIikaVuriIiEpIShIiIhKQEISIiISlBiIhISEoQIiISkhKEiIiEpAQhIiIhKUGIiEhIShAiIhKSEoSIiISkBCEiIiEpQYiISEhKECIiEpIShIiIhKQEISIiISlBiIhISEoQIiISkhKEiIiEpAQhIiIhKUGIiEhIShAiErOycvJZszWnynWcc3y/KRvnXD1F1XgoQYhIzLr0iU/4zYQvKCgq5p8fLWdDVl6FdT5enMFN/5rNjGWbohBhw6YEISIx7+yHp/GfeWu55pmZFZb99b3FAGzI2k1RcQmTv0qnqLikyv0559i4Y3dEYm1IlCBEJCarWfILikLODy7j8PFTyix7fdYqnp62jFEPfVjlvv/41kKu/scM0rflVrne/JXbWJGxM8yIG56EaAcgItG3YHUm97z+FePOOYozjuoe0WMVlziWrNvOwN6p1a6bu7eQxHgjuVnZr6rSL/bp946ixDle/PR71mXmct+FacSZVfjiD7Yley9d27WgpFwy/OjbDfRMbVlh/esmfMGxB3fiylMO5bevfMmi9O2BZdc++znT7hmJmZWJ66IT+3Lx0IO4942vABjcJ5WRg3pydO9UUpITqy13Q6EEIdIIfbcxmxZJCRW+0GZ/t4WUpISwvnzB+zU94oGpgedzv99SowThnGPD9t1l4vhq1TZen7WKS046mLS+HSts8+bsVbw84wceuew4ju6dSlFxCaMe+pALh/Tl6tMOK7PuBY9+BMDvfjaQU444AIDs3fsCyy//+6eMPaMfb81ZDcArM37gsp8cUmXMVz71Wcj5mTn5ZObkl5k37Zv1rN6aw+qtObw+a1XI7UY8MJWHLz2OQX32v+ZvzF7NG7NXB54vXJPFwjVZgJfUGgtVMYk0IIvSs/jHtKXVrnfzC7MZE1SffvnfP+Xaf37O/W8v5LevzuO/89cGlm3ZuYct2XtC7mfVlrI9fMr/Uv/yh62Vbgvw9PRljHlmJqu37AJg9ZYc7nn9K5Zt2Mmf/7sI8JLW8PFTWLs1h2c/Xs7LM34A4JWZPzB8/JRAdc6bc1ZTWFzCvz75jrz8wjJVQY9O/jYw/dWqzMD0tl17eeCdrwPPO7ROZuSD+xPej/H1miwe+2BJWOveOXEei9dtr37FRkYJQqQB+e0r83jvq3UUVtFIGlx98sOmbB54ZyHbdu0lPXN/ffgz05fzwcJ1lDjHFU99xhVPfsaYZ2ZWaGN46bPvyzz/ZMnGwPQTU5bwhzcXcMWTnzF8/BT+MnlRhVje+2odANc9N8v/+0Vg2a49BQDc//ZCAH494Qve/XJ/4lq2oWLd/OiHPuStOat57YuV7Mjbf6YQHx/H8PFTWLA6k0ff+7bCdqWenFoxuT495iQAXrx+WMht+nZuHXL+nRPnVXqcUO7495dhrVe+aqshU4IQiZKt2Xu4a+I8CotLSN+Wy9zvtwaW3fDcLIaPn8Kmcj1ldgZ9aQLc+K/ZfLFiS8j9Pzl1KWcFVR+tz8rjs6WbKCgq5uNvMwJfuOXN+c7b39Sv15eZ/7/FG/li+WYWrslk+PgpZOXk07r5/vr0UA3co2r5a37Sl2u55PFPAs/3FRYDcPdr82u8rx6pKUy/dxTd2qeEXP702JNqvM/p947i7dvPqHa9P11+XJnnH94zkji/vaIxsIbUayEtLc0tWLAg2mGI1IuqGlKDBddZh7tNZcZfdEyg4TTYwV3bsHKzV030k8O7csxBnar8pR7KpDvO5Py/fESHVklsz91X/QY/Uvf2KWRU0tV0wIHtWbxuB1D29SsqLqGouIQrn5rBzt37ePO202mbksR3G7O5+YXZAKQkJbB7X9keUKmtk8ny2yf+duUJ9O/RHqjYhgPw1DVDueH5WYFj5xcWkxBnJMRH5ve4mS10zqVFYt8RP4Mws3gz+8bMPoj0sUTCVVBUzFtzVlNQVMwPm7IpKCquk/1++cPWCr/6f6xf/PVjvtuYHda6N448gguH9K10+fOfrAg5/6lrhgamN+7YXePkAHD+X7wG5V+deliV691z/iBeu+W0wPPuHVK4ceQRNTrWtWf047GrhgBemf994ylllv/58uNDbpcQH0dyswTeuO10pt87irYpSQAcdkBbpt87iun3juKlG8rua/q9o3jiqhMDz0uTA4CZ8Va5M4mDurTm/ON787crTwAgOTE+Yskh0uqjF9PNwAogdEWfSBRc8JeP2FfkNYiWKv2luWtPAS2TE9iSvZev12TRoWUSJxzaOdCVsbwS5/jPvLXszi9i4hcry+wrWEFRMTe/MIcrhh1CfmH4CWnXngJufmE2b9x6OgAJccZ7d55VoTE2+JhvzllNKOsyK15J/OszDwfgoUuO5a7X5nPoAW0rNF7XRJuUJP58+fH83ysV6+T/OfYkevt1/lPvPgsw4uO81/WrVZl8+cPWCts8e+3JXPvs5xzUpXUgriGHdaF1i2Zlyjz93lFs2bmH9q2SMLNa9xZq3aJZYPq2swcA3hnE4L4dOblfl4rlDVofvKQx9ozDa3XshiaiCcLMugOjgAeB2yJ5LJFwlTjHvqKKjcCbd+7BOcdV/5gRcrvp944KVPHcfs4AzjyqB0CZev5Q3p67muf/tz8R/eHNqqtRjz+kc8gvyose+x/gfTnGxxk/P6EPb89dA8DNo44ss+77d47g7IenAfCv637Co5O/ZUUlZyE/PbYXAIP9LqlTFu5ve3j22pO5c+I87jh3YJlG27S+HUO2X4BXjXNUrw5l5j3/m5+wYHVmIDkAxMeV/VV97wWDuO2luYw5/TDG+Q2+v/vZQHp1ahX4sn9y6hI+WLieDq2SQx67S7sWIefX1hE9958tPHTJsdWuX3rWECsi2gZhZu8ADwOtgHHOudFVra82CKlLxSUljHzQ60IZ6mKmmurbuTWrgwZ+e/3W02iXklShDrr0eMUljosf+x85ewvDPsYto4/k428zQvbwKXXMQR154GLvy2rjjt08/O43/OWXx9O8XBfVj7/N4Mie7enSrkWFuvJHf3k8B3ZqRd7ewjKNt6Mf+jDQg2pY/27ced7RgWXOOdK35dK7c2vWbs3h1xP291gqX3Yzo7C4hPe+SmdQ79QyiaE6Jc4Fku5//m84LZL2l6u4pITd+UVlfuVHgnOO4hIXdtVQ9u59xMUZrZtHNq5QItkGEbEzCDMbDWxzzi00s2FVrDcWGAvQs2fPSIUjjdS+wmKSEuMBWLs1h5bNE+nYunml63/4zXqWrd/Jx4szysx/ZeZKfjnsELZW0ae/OqvLjQp6w/OzSE4M/S+0r6iEcx+ZVqP9P3XNUA7u2ob1WXks27CTP16YRlJiPL97tWx3yz6d9n/ZHtA+pUz7QbDgC97MjLfHncHPH/2Y288ZwJEHer/wy3+hBXev/dlxvcssM7PAF32zhPjA/F8M6ctbc1ZzxlHdue3sAYFEnBgfx/nH9wm7/KWCe/kEJwfwzjoinRzAK2tCfPi9jUrbMmJNJKuYTgTOMbORQDLQ2sxedc5dFrySc24CMAG8M4gIxiONzMbtu7n66RlcPPQgRhzdI/CLtXzd8vbcfC55/BNuHHlEyH7wABO/WEl8nPHvmT8E5k29eyTxccYz05fx3/npNY5vWP9uTArq1w8wqE8qX6/JqrSh+sXrh3HVP2aQkpTApDvOZE9BEV8s38yQQ7sEvviuOuVQBvVO5ZiDOgFwydCDeC3oKt5LTjqoxrGClwxqUi/ft0vlv/oP6LD/rONXpx3Gz4f0qdNfz2/edjoFIaoBpX5FrGndOXenc667c64XcBHwafnkIE1XcUkJHyxcV+mAalk5+Vz99AzAG0Dt+qALsIK3KS5xgf7ylSWHUsHJ4ey0AwONo2PP6Mfjfo8YgOtG9Afg2jMP55WbTi3T4yZYcHL41WmHMf3eURzj1+P/ppLql27tvT757/7fcMyMlKRERhzds8yv4mYJ8YHkAATOoMDrnVP+audIuH5EfxJr0POmrqtW2qYk0alN5WeKUj80FpNExefLN/Pk1KVs3LGba8v1+Fi2YQe3vTS3zLy8/P1J4U//XcQffuFVub5dSW8dgH4HtOUvV5zALx79mD3lEtENZ+3vVhkfF0e/7u144bphJCbE0alNc0YN6llp/fPTY4YGrhwu9Qu/a+ms70JftBbcBlJTefn72zDKN/7WtaeuGcqbs1cxcpCqe6WerqR2zs2oroFamobd+wpZvSWHR/7jDdvw7pdrGT5+Ck9M2T/mTfnkUN4c/4rjvQVFvFhuqIhgj101hMT4OO746VFl5p93fO+Q6x/QISXwq7Wy5PDYVUPo26VNmXmPXLb/atnLTq44UNzk342odXIAyNlbUOtta+rgrm2454LBYTXOvnbLabz7f2fWQ1QSLTqDkHp13p8/Cjl/6tfrA0M7JMbHUVhcwi9/cggfL85g806vYfnAji0D/fi/27iTm1+YE9j+yJ7tWbJ+B0+POYmv12bSsVXzwJdyUlCDKsA1p/WrcdxXn3ooL3z6PYd285LDw5ceF+j2eXTQyKnBI3oC/PWKE0hOLHv8mjr58G5MX5TR4LpQVtbVVGKHhtqQiHPOMeXr9Zx8eFd+/ujHYW9X2qBa4hxrt+ZyQPsWnPun6VWuG8qSddsD/erfHndGndWXb9u1l6TE+AoXSoE3qmmXds1JSWo8Y/9L49Soh9qQ2LV7XyFTv15f7V3I3p67hienLq2QHB669FjevaP6Koo4M/p2aU1ys4SQDZd/vDD8/426bEzt1KZ5yOQAXg8gJQdp7JQgpFolzrF7XyFZOfnc/MJsho+fQolznPfnj3hiyhLmfr81kCR27Slg+PgpfLpkI8PHT6G4pOxwFqWuOe0wBvfpSEpyIpeedHBgfudqeq5cdGLFcYYGh7gpTbA2MdpHXSTS1AYh1Qo1lETwvD/64/1fc9phPO8ngz/5N4spbYwu7+dBA8pl5uwNTD921RCemLIkMD5QeeXr8z+466xqu2P2TG3JKUd046fHhm6cFpHQlCAkpDdnr2Lu91t5/OoTq1/Z93yIM4XPl2+uMO/De0aWeR7cY6ZDq2Tuv+iYSo8RPJzE5ScfHHZf/d/97OjqVxKRMlTF1MR9szaL4eOnsCFr/yifw8dP4YVPv2fFxuzAGPg/1kH+VbnXDT+8wg1TrvcvTLvmtKqHiYayt8Qc2q9rncQmIqHpDKKJyssv5J7X5gdG+Fy4JpNrgu5xXOrSJz4p83xg7w4sWuvde/dnx/XmylMOrTDm0G1nD+Bv7y8uM+/xq0+s9Nd+Qnxc2ENAJDfbX8XUq1OrsLYRkdpRgmhCiktcYHiJ0pu7lHpm+vIqt733gkH07dKGLm2bM+nLtTz3vxVce0Y/zIyWyQllrnQePrCH9/BHTf0xVxGXV9oGEdd47too0mgpQTQRu/cVct6fP6JtSjN+FUZVzs2jjixzdXNwdc4FJ/ThghP2j9I56Y7hOOe46V+zy+x72j0jcVBnyQEInIWUNJzLd0RiltogmojJ/mil2bsL+Ot7+6t/7jovdOPtkUE3SgmHmfHkNUMZGHRVsZnV+Q3ae6S25Kyje/DKTafW6X5FpCKdQcSQEufKfCEXFpcw+qEPK13/wI4t+Un/bjz07jeBeUf3TmXjjt10Dboz1zNjT4pMwLV0y+gB0Q5BpEnQUBuNWM7eArJ3ewO5pW/L5cFJX/PyjafQpW0Llmfs5NYX51S5/ft3jihz45fy9hUWk5gQV+dnASJSdxrlHeUksnbk5XPxY59UmH/Fk5+FtX04DcdJP3KQORFp3JQgGrjtufl89G0GFw89iMycvUxflMFbs1ex70febasuG45FJDYpQTRgd702n4WrMwGYvmhDYNjrcJ17TC+6tW8R6MJ6/0Vp/P4NVeGJSHiUIBoY5xyFxSVc/9ws1gdd3VzT5NC3c2uuG9GfEufYkr2XYw7qyMBeHTjhkM5cfeqhdR22iMQgNVI3IAvXZHLXxPlhrfvM2JPo3akV73y5huf/9x0n9evKWUf34K7X5nP3+YM4/pBOVTZAi0hsiGojtZkdAjwDdHbOHWFmA4BznHMPRCKgpqyy5PCny4/jt6/MKzOvT2dvbKPzjuvDod3aMuBA717F4Q5ZISJSnXAulHsOuBMoBHDOLQYuimRQTdGnSzZWmDf17rOY/NvhDOyVysOXHscZR3UH4Oy0AwPrxMdZIDmIiNSlcNogWjjn5pfr9VJU2coSHuccs7/bQsvmiQzslRq4fwLApDvOJGdPAfFxccQ383L4oD6pDOqTyrhzjopWyCLSxISTILLMrC/gAMzsAqDiIP8Slj+88RVfrtxWZt5RvfafAZRWEbVM1u0qRSS6wkkQ1wMTgMPMbCOwFrg0olHFqO827qyQHAC+TfeGz75l9JH1HZKISKXCSRDOOXe6maUAcc65XDPTvRtr4eYXqh764qyje9ZTJCIi1QsnQUwCBjnndgfNewcYHJmQYlPpvRFgfzXStG/W89gH3pDa/Q5oG5W4REQqU2mCMLPDgP5AGzM7L2hRayA50oHFkvyC0G36wVegPHbVkPoJRkQkTFWdQRwKjAbaAmcHzc8FxkQyqFjyw6Zs0jNzA8+Db7RzSv9uPP7BEv7wi8EaG0lEGpxKE4RzbjIw2cxOcM7NrceYGr0FqzNZvG47fTq15uH/7L/XQuc2zbkm6I5ryc0SdGGbiDRY4bRBfGNm1+NVNwWqlpxzV0csqkbGOcfohz6kqMTRLCGOghAjrXZr34IXrz8lCtGJiNROOFdSvwJ0AYYDM4HueNVM4hvxwFSK/Jskh0oOAId3b1efIYmI/GjhJIiDnHP3Arudcy8DowB12K+h1Vtyoh2CiEiNhFPFVOj/zTazI4AtQK/qNjKzZOBzIMk/zjvOuT/UMs5GYfq9o1izNYf1mXll2h4AxpzRL0pRiYjUTjhnEBPMrB1wL/AesBz4cxjb7QNOdc4dBQwERpjZ8bWOtIGauWxTmed9Ordm2BHdKqx3aDdd5yAijUu1ZxDOuef9yZlAn6rWLbedA0rveJPoPxrOzSfqyEPvemcKV596WJn57/7fmWTnFXD10zMAja0kIo1POPeDaAv8Eq9aKbC+c+6mMLaNBxYCBwH/cM7Nq2aTRqH0quinrhkamHfW0T3KrJOSlEhKUiJ3nHuUkoOINErhtEFMBb4ElgChu+hUwjlXDAz0k8x/zOwI59zS4HXMbCwwFqBnz4Y/FtHv3/gqMH3D87MA6NWxFa1bNAu5/ukDutdLXCIidS2cBJHsnLvtxxzEOZdtZjOAEcDScssm4I0WS1paWoOtgsovLObcR6aFXNYjNaWeoxERibywroMwszFm1tXM2pc+qtvIzDr6Zw6YWXPgdOC7Hxlv1HyXsTPk/KGHdeHu8wfVczQiIpEXzhlEAfAX4G72NzI7qm+w7gq87LdDxAFvOec+qG2g0fTq5yt5ZeYPIZe1SWmmcZREJCaFkyBuw7tYLqsmO/bvXX10raJqYP47f22Z56/cdCqX//1TAI7S/aBFJEaFkyCWAXsiHUhDlpQQTy6FPPfrk8nNLyzTIH1QlzZRjExEJHLCSRDFwCIz+wzv4jcgvG6usWDB6kyycvMB6NmxFeANzlcqtbVujSEisSmcBPFf/9Ek3f3afACOP7hTYF5wm0NSYny9xyQiUh/CuZL65foIpCHK2J4XmP7jRcdEMRIRkfpX1S1H33LO/cLMlhBiiAzn3ICIRtYA/OrpmQAc1atiQ3S7lCQ6tEqq75BEROpNVWcQN/t/R9dHIA1NSVA7w8OXHldh+Ru3nV6f4YiI1LtKL5Rzzm32J69zzq0LfgDX1U940fP+V+kAXDf8cOLjdJ2DiDQ94VxJfUaIeWfVdSANyeT5a3l6+nIABvXpGOVoRESio6o2iN/gnSn0NbPFQYtaAbMjHVg0lSYHgB6pLaMYiYhI9FTVBvEa8CHwMPC7oPm5zrkdEY0qivbsKwpMv3/niChGIiISXVW1QexyzqUD9wBb/LaH3sBlpYPwxaKJX6wE4E+XH0ezBF3jICJNVzhtEJOAYjM7CPgXXpJ4LaJRRdE7c9cAcGTPagesFRGJaeEkiBLnXBFwHvC4c+5WvJFaY87eAq96qXmzeOLjwnlpRERiVzjfgoVmdjHebUdLh+uOyXtoLl3vNa38/udpUY5ERCT6wkkQVwEnAA8659aaWW/g1ciGFR3frM0iMT6O/j3aRTsUEZGoC2cspuVm9lugp/98LfBIpAOLhm/Tt9Ove1sNwCciQhhnEGZ2NrAImOY/H2hm70U6sPq2cftuVm3JoXObFtEORUSkQQiniuk+4FggG8A5twivJ1NMmbdyKwCD+6ZGORIRkYYhnARR5JzbVW5ehdFdG7tVW3Jo06IZP+nfLdqhiIg0COHcMGipmV0CxJvZwcBNwJzIhlW/nHN8szaLo3unEmcamE9EBMI7g7gR6I93u9HXgF3ALZEMqr5N+HgFO/L2MeBAXRwnIlIqnF5Me4C7/UdMenfeWgCOP6RzlCMREWk4dLmwL86MDq2Sox2GiEiD0eQTRFZOPgBjz+gX5UhERBqWJp8gPlmyEYBDD4jZAWpFRGqlqhsGPUkV3VmdczdFJKJ69sKn3wFwSNc2UY5ERKRhqaqRekG9RRElmTl7A9MJ8U3+ZEpEpIxKE4Rz7uX6DCQavlqVCcCz154c5UhERBqeqqqYqhxvyTl3Tt2HU7+WrNtO+5ZJHNhR950WESmvqiqmE4ANwOvAPCDmLjFelrGT/j3aYbp6WkSkgqoq3rsAdwFHAE8AZwBZzrmZzrmZ9RFcJG3PzWdr9l4O7657P4iIhFJpgnDOFTvnpjnnrgCOB1YBM8zsxnB2bGY9zOwzM1thZsvM7OY6irlOfJu+HYAjdO9pEZGQqhxqw8ySgFHAxUAv4O/Au2Huuwi43Tn3tZm1Ahaa2cfOueU/It46s2T9DlKSEujbRd1bRURCqaqR+mW86qUPgT8655bWZMfOuc3AZn8618xWAAcADSJBrN2WQ98urYmPU/uDiEgoVbVBXA4cAtwMzDGzHP+Ra2Y5NTmImfUCjsZr7C6/bKyZLTCzBZmZmTXZba0VFZewZksOfTq3rpfjiYg0RlVdB1EnV46ZWUtgEnCLc65CYnHOTQAmAKSlpdXLjYh+2LyLfUUl9FMDtYhIpSJ6+bCZJeIlh4nOuXDbLiJuRcZOAPpp/CURkUpFLEGYd3HBv4AVzrm/Reo4tfH9xmw6tWlO57Ytoh2KiEiDFckziBPx2jFONbNF/mNkBI8Xtu83ZXNoN/VeEhGpSjj3pK4V59wsGuDV19m797Eley+jBx8Y7VBERBq0JjeE6aotXjv5wTqDEBGpUpNLEDOXbQKgTyd1cRURqUqTSxB5+YUAtG7RLMqRiIg0bE0uQWRs383xB3eKdhgiIg1ek0oQBUXFZGzfrSuoRUTC0KQSxPrMPEqco7cShIhItZpUglizzevB1KdzqyhHIiLS8DWpBLF2ay5JCXF0bZcS7VBERBq8JpUg1mzL4cBOrTTEt4hIGJpMgnDOsXZrrq5/EBEJU5NJEDvy9rFrTwG91f4gIhKWJpMg0jNzAejVSQlCRCQcTSZBZGTlAdAztWWUIxERaRyaTILYsH03LZISaJeSFO1QREQahSaUIPLo0aEl3n2MRESkOk0mQazdmsuBHVW9JCISriaRILJ3+z2Y1EAtIhK2JpEgNvgN1D3UQC0iErYmkSDWqQeTiEiNNYkEsSErj+TEeDq2aR7tUEREGo0mkSDWZebRI7UlcerBJCIStiaRINZn5dKroxqoRURqIuYTRH5hMdtz99GtfYtohyIi0qjEfILYsnMPAN10DwgRkRqJ+QSxaeduALrqDEJEpEZiPkGUnkF0bacEISJSEzGfIDbt3ENKUgKtkhOjHYqISKMS8wli885/53AZAAARCUlEQVQ9dG3XQoP0iYjUUBNJEGqgFhGpqZhOEMUljq3Ze+im9gcRkRqLWIIwsxfMbJuZLY3UMaqTlbOXohJHFyUIEZEai+QZxEvAiAjuv1qbA9dAKEGIiNRUxBKEc+5zYEek9h+OTeriKiJSazHdBrFl5x4S4ozU1hrFVUSkpqKeIMxsrJktMLMFmZmZdbrvTTv30KVtC+Lj1MVVRKSmop4gnHMTnHNpzrm0jh071um+N+/crQZqEZFainqCiBTnXOAiORERqblIdnN9HZgLHGpmGWb2q0gdK5TcvYXs3lekHkwiIrWUEKkdO+cujtS+w7E5u7QHk66iFhGpjZitYtq8Q11cRUR+jJhNEKX3gVAjtYhI7cRsgti8cw/tWyaRnBgf7VBERBqlmE0QW3ftpUtbnT2IiNRWzCaILdl76NJWV1CLiNRWTCaI4pISMnfl01lnECIitRaTCSIrJ58S5+isMwgRkVqLyQSxdddeADq30RmEiEhtxWaCyC5NEDqDEBGprZhMEJt37sGATqpiEhGptZhMEFuy95DaOpnE+JgsnohIvYjJb9At2RrFVUTkx4rYYH3RtCV7D4P61O29JUSaqsLCQjIyMsjPz492KE1acnIy3bt3JzExsd6OGXMJoqComO25+3QVtUgdycjIoFWrVvTq1Qsz3Z0xGpxzbN++nYyMDHr37l1vx425KqYtfg+mrmqgFqkT+fn5dOjQQckhisyMDh061PtZXMwliM0axVWkzik5RF803oOYSxDrs/IA6JHaMsqRiEismDFjBqNHj64wf9GiRUydOrVW+3zooYcC0+np6RxxxBG1ji9SYi5BrMvMo33LJFo3bxbtUESkHhUVFdX7MatKENXFE5wgGqqYa6Ren5lHT509iMSU8ePHM3HiRHr06EFqaiqDBw9m3LhxDBs2jCFDhjB79mzOOeccLrjgAq6++moyMzPp2LEjL774Ij179uTKK69k9OjRXHDBBQC0bNmSvLw8ZsyYwX333UdqaipLly5l8ODBvPrqq5gZ06ZN45ZbbiE1NZVBgwZViKmgoIDf//737N27l1mzZnHnnXeyYsUKNm3aRHp6OqmpqZx55pksWLCAp556CoDRo0czbtw4pk2bxt69exk4cCD9+/fnwQcfpLi4mDFjxjBnzhwOOOAAJk+eTPPm0W1LjakzCOcc67NyObBjq2iHIhK7hg2r+Hj6aW/Znj2hl7/0krc8K6vismosWLCASZMm8c033/Duu++yYMGCMsuzs7OZOXMmt99+OzfccAO//OUvWbx4MZdeeik33XRTtfv/5ptvePzxx1m+fDlr1qxh9uzZ5OfnM2bMGN5//32++OILtmzZUmG7Zs2acf/993PhhReyaNEiLrzwQgAWLlzI5MmTee211yo95iOPPELz5s1ZtGgREydOBGDlypVcf/31LFu2jLZt2zJp0qRqY4+0mEoQmTn57C0opmdHnUGIxIpZs2Zx7rnn0rx5c1q1asXZZ59dZnnpFzPA3LlzueSSSwC4/PLLmTVrVrX7P/bYY+nevTtxcXEMHDiQ9PR0vvvuO3r37s3BBx+MmXHZZZeFHe8555xTq1/+vXv3ZuDAgQAMHjyY9PT0Gu+jrsVUFdO6zFwAnUGIRNKMGZUva9Gi6uWpqVUvD8E5V+XylJSUSpeV9vxJSEigpKQksL+CgoLAOklJSYHp+Pj4QNtBbXsNBccTfFygym6q5ePYu3dvrY5fl2LqDGJdpteD6UC1QYjEjKFDh/L++++Tn59PXl4eU6ZMqXTdIUOG8MYbbwAwceJEhg4dCkCvXr1YuHAhAJMnT6awsLDKYx522GGsXbuW1atXA/D666+HXK9Vq1bk5uZWup9evXqxaNEiSkpK2LBhA/Pnzw8sS0xMrDaOaIupBJGemUu7lCRat1APJpFYccwxx3DOOedw1FFHcd5555GWlkabNm1Crvv3v/+dF198kQEDBvDKK6/wxBNPADBmzBhmzpzJsccey7x586o86wBvWIsJEyYwatQohg4dyoEHHhhyvVNOOYXly5czcOBA3nzzzQrLTzzxRHr37s2RRx7JuHHjyjR2jx07lgEDBnDppZeG+1LUO6vu9K0+paWlufINUDVxzdMz6NY+hfsvOqYOoxJp2lasWEG/fv2iGkNeXh4tW7Zkz549nHzyyUyYMCFkz6JYF+q9MLOFzrm0SBwvZtogcvYWsGH7bk4f0D3aoYhIHRs7dizLly8nPz+fK664okkmh2iImQSxOH07AP17to9yJCJS16rqMiqREzNtEPNWbqNlcgKHd28b7VBERGJCTCSIEueYv2obaX07ER8XE0USaVAaUltlUxWN9yAmvk2//GEr2bsLGHpYl2iHIhJzkpOT2b59u5JEFJXeDyI5Oblej9vo2yCcc7wzdw2d2zRnyGGdox2OSMzp3r07GRkZZGZmRjuUJq30jnL1KaIJwsxGAE8A8cDzzrlH6voYny3dxLINO7luRH9VL4lEQGJiYr3exUwajoh9o5pZPPAP4CzgcOBiMzu8rvbvnGPGsk387f3FHN69HaMH96yrXYuICJE9gzgWWOWcWwNgZm8A5wLLw9m4xDmKS/xHcQm5+YXk7ClgR94+Vm3exZcrt7Fy8y76HdCWP16YprMHEZE6FskEcQCwIeh5BnBcVRus2rKLUQ9OpbjEUVVzmAF9OrfmxpFHMGJgDxLilRxEROpaJBNEqKEQK3zvm9lYYKz/dN/Ue0YtDfcAz9QysChJBbKiHUQEqXyNm8rXeB0aqR1HMkFkAD2CnncHNpVfyTk3AZgAYGYLIjWmSLTFctlA5WvsVL7Gy8xqP4BdNSJZN/MVcLCZ9TazZsBFwHsRPJ6IiNShiJ1BOOeKzOwGYDpeN9cXnHPLInU8ERGpWxG9DsI5NxWYWoNNJkQqlgYglssGKl9jp/I1XhErW4O6H4SIiDQc6h8qIiIhNYgEYWYjzOx7M1tlZr+Ldjw1YWbpZrbEzBaV9iYws/Zm9rGZrfT/tvPnm5n93S/nYjMbFLSfK/z1V5rZFVEszwtmts3MlgbNq7PymNlg//Va5W9buzvD113Z7jOzjf77t8jMRgYtu9OP83szGx40P+Tn1e+QMc8v85t+54x6Y2Y9zOwzM1thZsvM7GZ/fqy8f5WVr9G/h2aWbGbzzexbv2x/rCoeM0vyn6/yl/eqbZmr5JyL6gOvAXs10AdoBnwLHB7tuGoQfzqQWm7en4Hf+dO/A/7kT48EPsS7RuR4YJ4/vz2wxv/bzp9uF6XynAwMApZGojzAfOAEf5sPgbOiXLb7gHEh1j3c/ywmAb39z2h8VZ9X4C3gIn/6n8Bv6vm96woM8qdbAT/45YiV96+y8jX699B/PVv604nAPP89CRkPcB3wT3/6IuDN2pa5qkdDOIMIDMnhnCsASofkaMzOBV72p18Gfho0/9/O8yXQ1sy6AsOBj51zO5xzO4GPgRH1HTSAc+5zYEe52XVSHn9Za+fcXOd9mv8dtK+Iq6RslTkXeMM5t885txZYhfdZDfl59X9Jnwq8428f/DrVC+fcZufc1/50LrACb0SDWHn/KitfZRrNe+i/B3n+00T/4aqIJ/g9fQc4zY+/RmWuLq6GkCBCDclR1Zve0DjgIzNbaN5V4QCdnXObwftQA538+ZWVtaG/BnVVngP86fLzo+0Gv4rlhdLqF2petg5AtnOuqNz8qPCrHI7G+yUac+9fufJBDLyHZhZvZouAbXhJeXUV8QTK4C/fhRd/nX7HNIQEEdaQHA3Yic65QXij1l5vZidXsW5lZW2sr0FNy9MQy/kM0BcYCGwG/urPb7RlM7OWwCTgFudcTlWrhpjX4MsYonwx8R4654qdcwPxRp04FuhXRTz1UraGkCDCGpKjoXLObfL/bgP+g/fGbvVPx/H/bvNXr6ysDf01qKvyZPjT5edHjXNuq/+PWQI8h/f+Qc3LloVXRZNQbn69MrNEvC/Pic65d/3ZMfP+hSpfrL2HzrlsYAZeG0Rl8QTK4C9vg1d9WrffMfXRAFNN40wCXiNYb/Y3nvSPdlxhxp4CtAqanoPXdvAXyjYK/tmfHkXZRsH5/vz2wFq8BsF2/nT7KJarF2UbcuusPHhDsBzP/kbOkVEuW9eg6Vvx6m8B+lO2sW8NXkNfpZ9X4G3KNiheV89lM7x2gcfLzY+J96+K8jX69xDoCLT1p5sDXwCjK4sHuJ6yjdRv1bbMVcZVnx/gKl6ckXg9ElYDd0c7nhrE3cd/ob8FlpXGjlcX+Amw0v9b+s9leDdRWg0sAdKC9nU1XoPSKuCqKJbpdbzT9EK8Xx2/qsvyAGnAUn+bp/Av1oxi2V7xY1+MN1ZY8JfN3X6c3xPUW6eyz6v/eZjvl/ltIKme37uheNUGi4FF/mNkDL1/lZWv0b+HwADgG78MS4HfVxUPkOw/X+Uv71PbMlf10JXUIiISUkNogxARkQZICUJEREJSghARkZCUIEREJCQlCBERCUkJQholM5thZhG/x7CZ3eSPHjqx3Pw0M/u7Pz3MzIbU4TF7mdkloY4lUp8iekc5kYbIzBLc/vFtqnMdXl/ytcEznXMLgNKbxQ8D8vAulKyLGHoBlwCvhTiWSL3RGYREjP9LeIWZPeePcf+RmTX3lwXOAMws1czS/ekrzey/Zva+ma01sxvM7DYz+8bMvjSz9kGHuMzM5pjZUjM71t8+xR+w7St/m3OD9vu2mb0PfBQi1tv8/Sw1s1v8ef/Eu1DpPTO7tdz6w8zsA3/QuF8Dt5p3L4KTzKyjmU3yY/jKzE70t7nPzCaY2UfAv/3X5wsz+9p/lJ6FPAKc5O/v1tJj+fto778+i/3XY0DQvl/wX9c1ZnZT0Osxxbz7DCw1swt/3LsqTUp9XumpR9N64P0SLgIG+s/fAi7zp2fgX7kLpALp/vSVeFeHtsIbfmAX8Gt/2WN4A7SVbv+cP30y/vAZwENBx2iLd+Voir/fDEIMYQIMxrsSNwVoiXdV/NH+snTK3e/Dnz8M+MCfvo+g+xHg/fIf6k/3BFYErbcQaO4/bwEk+9MHAwvK7zvEsZ4E/uBPnwosCtr3HLwhFlKB7XhDRp9f+jr567WJ9udCj8bzUBWTRNpa59wif3ohXtKozmfOG+8/18x2Ae/785fgDUlQ6nXw7vNgZq3NrC1wJnCOmY3z10nG+5IG/x4HIY43FPiPc243gJm9C5yEN/RBbZwOHG77b7bW2sxa+dPvOef2+tOJwFNmNhAoBg4JY99D8b70cc59amYdzKyNv2yKc24fsM/MtgGd8V6zR83sT3hJ5otalkmaICUIibR9QdPFeAORgXdmUVrFmVzFNiVBz0so+5ktP05M6bDG5zvnvg9eYGbHAbsribGub5sZB5wQlAhKY6BcDLcCW4Gj/G3yw9h3VcM2l3+tE5xzP5jZYLxxeB42s4+cc/eHVQpp8tQGIdGSjle1A3BBLfdxIYCZDQV2Oed2AdOBG83/Njazo8PYz+fAT82shZmlAD/DG00zXLl4VWKlPgJuKH3inyGE0gbY7Lxhqi/HG3Uz1P7Kx3qpv99hQJar4p4PZtYN2OOcexV4FO+WqyJhUYKQaHkU+I2ZzcGrM6+Nnf72/8QbmRVgPF7VzWIzW+o/r5LzbmP5Et6omPOA551zNaleeh/4WWkjNXATkOY3JC/Ha8QO5WngCjP7Eq96qfTsYjFQ5Dcs31pum/tK943XmH1FNbEdCcw3705ldwMP1KBc0sRpNFcREQlJZxAiIhKSEoSIiISkBCEiIiEpQYiISEhKECIiEpIShIiIhKQEISIiISlBiIhISP8PzuaomW7jz7EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>.9*mi):\n",
    "            plt.axvline(t,label='90% reached',linestyle=':',color='green')\n",
    "            break\n",
    "plt.xlim((0,mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}