MixedGaussian_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/MixedGaussian_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 mixed gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.mix_gaussian import MixedGaussian"
   ]
  },
  {
   "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 = 1               # 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",
    "mg = MixedGaussian(sample_size=sample_size,rho1=rho,rho2=-rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = mg.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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X/wZXV93/HnG5AN8iNKoCEoXzBoTCzDqN36Y2qijSRCJElDY2riiK0Td4hJJ2bsahXiT8wodBgz0TpuK9E1/mzF6AgGsG1CbRbr4hBcRFvWuAvCgrqSBaTg1333j3s+y/mePb9/3M85574eMwx7v/fccz7n3Hs+7/P5be6OiIjIEbETICIi46CAICIigAKCiIgkFBBERARQQBARkYQCgoiIAAoIEoGZ/bWZ/e6SjvV7ZnaPmT1gZj9RY/tvmdk5y0jbGJjZB83s0tjpkHFQQJBBJBnrQ0lGfI+Z/bmZHddwH2eYmZvZUS3T8BjgCuCX3f04d/9em/2U7N/N7Ml97lMkJgUEGdKvuvtxwDOBfwpcsuTj/yTwY8CtSz6uyCQpIMjg3P3bwOeBs7LvmdkRZnaJme0xs3vNbLuZ/Xjy9g3J/+9LShrPzfn8JjN7t5ndlfz37uRvPwN8I/X5/56XNjN7eXLs75nZxZn3nmVmO8zsPjO728zeY2ZHJ++FtP1dkrZ/ZWaPN7PPmdl3zOz7yb+fWHRdzOz1ZvZtM7vfzL5hZi+sOm7yvpvZq83s/yaffbuZnZl85oCZfTKVzheY2Z1m9kYz+25ScntZSZrON7Obk2P/rZmdXZVemRF313/6r/f/gG8B5yT/Po3FU/rbk9d/Dfxu8u9XArcDPw0cB1wFfDh57wzAgaNKjvM24EbgHwEnA3+bOk7p54GnAQ8AvwBsYlG9tJ5K9z8BngMclezrNuA1qc878OTU658A/iXwWOB44L8Af1lw7KcCdwCnptJ6ZoPjfhY4AfjHwMPAf0uu4Y8DXwNekWz7guScrkjO8fnAg8BTk/c/CFya/PuZwL3As4EjgVck3+OmsvTqv/n8pxKCDOkvzew+4IvA3wB/krPNy4Ar3P2b7v4A8AbgpQ3aDV4GvM3d73X37wBvBV5e87O/CXzO3W9w94eBPwYOhjfd/SZ3v9Hd1939W8D7WWSoudz9e+7+KXf/gbvfD7yjZPsfschon2Zmj3H3b7n77gbHfZe7H3D3W4FdwHXJNfwHFqWxZ2S2/2N3f9jd/wa4GvitnDS9Cni/u3/J3X/k7h9iEWyeU5ZemQ8FBBnSv3D3x7n76e7+and/KGebU4E9qdd7WDwZ/2TNY+R9/tQGn70jvHD3B4FDDc9m9jNJtc8+MzvAIqCdVLQzM3usmb0/qYI6wKLK63FmdmR2W3e/HXgN8BbgXjP7uJmd2uC496T+/VDO63QD/veTcwuKrtHpwGuT6qL7kmB+GotSQWF6ZT4UECS2u1hkRMEaiyqOe1hUjbT5/F01j303iwwPWGToLKp9gvcBXwee4u4nAG8ErGR/r2VRtfLsZPtfCLvO29jdP+ruz0vS78C7Wh63yuPN7NjU66JrdAfwjiSIh/8e6+4fq0ivzIQCgsT2MeCPzOxJSbfUPwE+4e7rwHdYVOH8dMXnLzGzk83sJOBNwF/UPPZ/Bc43s+cljbBvY+M9cTxwAHjAzH4W+L3M5+/JpO14Fk/n95nZicCbiw5sZk81s180s03A/0s+96Oax23jrWZ2tJn9PHA+i/aNrP8EXGRmz7aFY83sxWZ2fEV6ZSYUECS2K4EPs6he+XsWmc2/BXD3H7Coh/9fSRXGc3I+fymwE7gF+CrwleRvlZL6998HPsqitPB94M7UJv8O+B3gfhaZ5Scyu3gL8KEkbb8FvBs4Bvgui4buvyo5/Cbgncm2+1g0ir+x5nGb2sfi3O4CPgJc5O5fz27k7jtZtCO8J9n+duBf10ivzIS5a4EckbkysxcAf+Huhd1fRQKVEEREBIgYEMzsx8zsf5vZ35nZrWb21lhpERGRiFVGZmbAse7+gC3mnPki8IfufmOUBImIrLhWk4b1wReR6IHk5WOS/9SgISISSbSAAJAM2LkJeDLwXnf/Utn2J510kp9xxhnLSJqIyGzcdNNN33X3k6u2ixoQ3P1HwNPN7HHAp83sLHffld7GzLYAWwDW1tbYuXNnhJSKiEyXme2p3mokvYzc/T4WE56dm/PeNnff7O6bTz65MsCJiEhLMXsZnZyUDDCzY4BzWAzXFxGRCGJWGf0Ui1GeR7IITJ90989FTI+IyEqL2cvoFg6foldERCIZRRuCiIjEp4AgIiKAAoKIyOhs3b6Drdt3LP24CggiIgJEHpgmIiKPCqWCW/bs3/D68gufu5Tjq4QgIhJJrKqhIioh1LTsSC0iqyfkL7HyGwUEEZEli101VEQBocJYvzgRma9Y+YsCgojIksWuGiqigFBhrF+ciEjfFBBERCIZ2wOmAkJNY/viRORRKsH3Q+MQREQEUAlBRCZMvQD7pRKCiIgAKiGIyISpF2C/VEIQERFAJQQRmQGVDPqhEoKIiAAKCCIiklBAEBERQAFBREQSCggiIgIoIIiISEIBQUSiG9vawqtKAUFERAANTBORiDQ53biohCAiIkDEEoKZnQZsB04BDgLb3P1PY6VHRJZPk9ONS8wqo3Xgte7+FTM7HrjJzK53969FTJOISFQxg2O0gODudwN3J/++38xuA54ADBIQ5voEMtfzktWi3+84jKJR2czOAJ4BfCnnvS3AFoC1tbWlpks2UvCRuYv5Gx9DA3v0gGBmxwGfAl7j7gey77v7NmAbwObNm73p/sdwkYcw1/MSkXiiBgQzewyLYPARd78qZlqk2JDBR4FMxmAMD1hjaGCP2cvIgA8At7n7FUMdZwwXeQjLOq+t23ewe98BzjzlhEH2LyLjEbOE8M+AlwNfNbObk7+90d2viZgmyXHmKSdw+YXPHaRkoCovGUqT39SYHhxjHjtmL6MvAras4801oxmyZAAbM2yVFETmLXqjsgxjiCedUFLow5ieyGReupQ+V/13qIAguZRhy1Tt3ndYZ0WpSQFhILEy0qnVzY81XTINeb/vbLWmfmP1KSCsmKYBQjeTjEnZ7zf7MHTspu7Z29gfqPqmgNCz2E/oquqRZRpTSTjb6UEdIJpTQFgRsQOVSBd1nv777B69qveLAkLPxvKEHuuJbe43jCzEzjDVTjAMBYQVsYxAtXvfAbZu36Gbc+LGGNyLfr956zD3ke6xPNgtmwLCQFblB5R9UlRQWA3LzDDLjqHfWb8UEFbMUCWD4MGH1xUUJip2NVAd2bRUpbHrOYzp3JdBAUE6CY14u/cd4MGH1wH17lglyygZjDlAzY0CgnSWDgoKBtPVZzXQMmbhheJgoWDSjgKC9CLd3U+kD6vasBuTAoL0Tk9l0zaF/vtVwULBpB0FBBEZTB8TzSkzXx5zb7xMcTSbN2/2nTt3xk6GVIjxVKYnwXG54LJrgUc7GOh7icvMbnL3zVXbqYQgIp1l249Cj7OqqiMF8nFRQJDexSgZqN2if02vZageUk+z6VJAkF7VmZ5YmfVwhrrGVU/42Unnzj79xA3bVX2uS7r1u+qPAoJMmnqT9K/petp5I9VBJYUpUqOy9CKbiYQnxPT4hPRTZNX6zE0z+LrbzzVwlF3/rvs6dtNRh9oE8vabHane5Nh9lAz6OOe5U6OyrBRlAv1Jl7pCySBkukXbp7fVdzFdKiFIr8rqmus8RQ711FdU133V617Uab9j02cJKL2von8Pdeyq9CzzuHOgEoKI1FI12rcuZcjTpxKCLFXsuv4wYCpbUglWMVNr+p0su85ebQXdqYQgvVPRfF6WNYZDv5vpUECQpaqbKQyVeYQ2g+zI2qkPbFtGuofo4ltnX+pavDwKCFJpCqOBx5imIfVxvkNntG1+N6v2PY5N1IBgZlcC5wP3uvtZMdMi/YndTlDHXHqqZNtElllS6CIbLC647NrKgWxNjzvV7zSm2CWEDwLvAbZHToeUWHaRvclxupZepraoT/Z8j7Du+xzq+2zyu3nokfUN3ZKVmccRNSC4+w1mdkbMNEh/6mbOY6yCmmrGczDpJBjGVQxxHkN8P2E8QxgFfbDHzo5j/H1NRewSQiUz2wJsAVhbW4ucmtXQV7/0tsdtciPXeQrNey87UC5WptH0uGG7bFXRmNU9tzNPOYFde/dzzNFHKfOOZPQBwd23AdtgMQ4hcnKkRN0qgqLtyqpvqvbZx8pcU5Ktby+6Ln3OFVRntHITeb+DEOi6UK+k9kYfEGR5ll3UzmbiXW7kvAbJovMJwtP1kNUtRbpe6zpBtGlaYmacu/cdYOv2HdFLbKtOAUF616SKoEze5GpFpYqp1Rf3VaKpKhl0uS55AXrr9h1s3b6j15XQ0iWOPo39NzBGsbudfgx4AXCSmd0JvNndPxAzTatsWUXtqsyq7Om3SUaabrjMm4VzmcEje6y6VT7LSNMQwbRpwFM1zzjE7mX02zGPL8vRtitoesrlkKlnhUw/3f2yz8yk7wyqaNbVvo/XZwab/mzVfsPfVfUzTaoyksMMffO2yayK+qnD4mk0dFsMddFBOqhs3b4jN3MbUlUAGHpVsbIn9SGeyrPHa1tSkDgUEGQwbaskiuqu0/uCRzObdNfLvpZwHKo6paiqaFnH60tRusLxwnloGc1pUUCQaPpo4ISNQQIWmVAICkNOldxlnESbRtShRnAPcW366rmlKqflUkCQwXStkmhSkkiXCvrotTJ0I2fRoL+5ZIAqGUyTAoJMSlVVRWhDyC6m0lVZyaRJSaHKrr2Hr1081AjuvGO0DUTZNHa97lPtTjx1CggyuGU1UkP/A7X6HABWxzFH93NLdklvm8y3z5Hi6RKfLJcCgsxKX1UvRb2V+p6+IcjOTRReX/W6F5WeU91pQorkPYmHwYBN9D2uIl3iyxtDIsNQQJBZ6ePJOGSOffRWWrYmVS15mX/o2nvLnv2H1igoy4yrutV2SXs6TdkuwzIMBQSZpS4lg+DBh9cPZXB5vZX6zKDC0p7pkgHkV1+l/95kgZk86Ub47DQhDz1SPZNqtmpniOA5pYA8dQoIMgt9zd2TnRYbaFWF0jYNXdWpMssLJg89ss5ZayduyOAPevnTedhP+tr0OftpH/uUZhQQZNLqVBE1nWjtgsuu5QiDs9ZO3PD3ZciWDMrmewrzND348PqGEdxVac1rsA1rEJx36dUb/h5GiGelA+cte/Z3qiqS8dC3KIdMuWtfn/34Q7XJrr37D02JUbXfrt1S29i970BhtU7ZsdINtvBoQ/Z5l1592Mpl6fPPjvvILs5TpxTVZxdd6Z8CgoxeWe+aorUOyrapO2V03bR0UVQ9UifAlfXsqXPNuqzHnG5rKJpNVqZHAUFmNQioa5qzdehHGBuWdExn4HkNsemeOUNdx7569qSrxODR6rKHHlnfsFZzNrPPlhTq9ETKW9ei7LNT/g1OmQKCRNP0ib1swFjZHD1NJtUrqg7J/r3O4Kmy94vWRM5rK2iiyTXLk54HKryuqn7KC37K0KdJAUFW9ibOO990UAjy6sbTwWH3vgMcYYsSRd703H0LaQyll6q6+6p+/EVP/9n3yj5XdtzsuhZBaJCO0f4i+RQQZOnq3vRNSgFl2mYmIXPftXc/xxx91KHXVXXv2cbavAwvu55zyNTrpDVdrRMy1aLjlE30V3b9616zvC6s6fPTFBTTooAgh6xKXW7Tp9CDvnGQVvapPBu4sr1/+p7nJ13HnxWeykMa0k/oy5qPKS1UQdVdwnRVS6tjoYAgS9embr/J9l3ltSWEQVvpjK0ogy3LsLucS0hTet/paqN0aSDbFtBXGuruq04JRcZHAUFylT1Fj+XprW066lZFZXvCVI2gDdtkg8FDj6z3Notpkbz6+iB2t9A6bRdN3pPhKCBING3rqZs+cXYJHOljVX2+qOdRXsNvmwwvPYo6BJ1019G6VVNNrkfTUlzRa5kGBQTJVVT0z6uTXvbNX9Zo28QQDdLpvvawqNIJ01Gk1ZlrqGgwWnZ1uPTnsvMuKWOWJhQQZPSK6uzLRhSnt286UrlNJpoNnGHls+zgr7qy6zEE6cFoITBkz7No2u4m10PdP1eTAoKUGmMPkLGko0y6oTetTttMnfUYysYeNOnCKpJm7gXdIUZo8+bNvnPnztjJWHljyYjrpqPpiOi8tQ+q9p2dRiJUZxWtLZx3rFCyyOuhlN1PXtrK1lLIS3OfbQh1jeW3s2rM7CZ331y1nUoI0thYbubY6cirsy8SnvbLMuu8uYSK9lU28rgOrUImeQpLCGZ2DfBqd//WUlNUQiUEGVLTp9eykcHZQJEtEQTZ6q/0TKTprqrphumiPv9VpZRsA3STc+2qSylMuuujhPBB4Doz+xBwmbv/sK/EiXQ1ZIbWtIopveJYVrYraGhs/vwlLy7dLm/cQtfut3X20ef1VPXQ9BQGBHf/pJldDbwJ2GlmHwYOpt6/ouvBzexc4E+BI4H/7O7v7LpPkSp15k5qU6VSNDgs3RW1aIBa6E4aBpEFedNFZ9UdLRxzxPAUOgJIdRvCD4EHgU3A8aQCQldmdiTwXuCXgDuBL5vZZ939a30dQ+ZnyO6Qeb18skEhfbzsRG7Z5SbTo5zD0pRFs6Gmg0hRxt/HOgtVXXizo7PbUJfV6SoMCMnT+xXAZ4FnuvsPej72s4Db3f2byfE+Dvw6oIAgg6iTUaUz9AcfXq8sKaQnsUsvKgMbB46FkkF2/YPsMcM+d+87wFWve9FhGXedkkL6ddFgwlgUFMatrIRwMfASd791oGM/Abgj9fpO4NkDHUtmYhlVD2Hheni0T39eMAnVPNleQeF1eq7/ot5FYfbS9HQUYR9ZRSOU28jrirpr7/5DE/pl1ylos2+VDKanrA3h5wc+dt6s8od1eTKzLcAWgLW1tYGTJHNWJ6NKL0BfZ2bTPOm1E/JkSx3Z0kN6Yfvs0pNAo8xambM0EXMcwp3AaanXTwTuym7k7tuAbbDodrqcpMnYDdEbpqhraPp42cCQzfirJrELbQ+79u7fUJVz9uknHtYYHUoEeZPWpaet6HotwjmFBXdClVffpRAZv5gB4cvAU8zsScC3gZcCvxMxPbIiYmVUIePNW00s/e9jNx3FQ4+sbyitpKuxsj2R6lDmLHVECwjuvm5mfwBcy6Lb6ZUDtleIFGqzqEv4TLanUfb9tLwn/fQCN6GEkG6DSDdQZz/fVy8eVStJEHXqCne/BrgmZhpk9RT114fDxxJUZZJ11zkI+0ln/NnlJfPaDKqqpET6pLmMRBJlJYOqLqt1Gp1Dxh9KFUWN2kVzGmUbufteCU0lA1FAkMlrOutpyNTPu/TqDT2CiuYmSr9XpKxkkF6bOV1SKNtH2UI4Y6OqpvlQQBCpoc969qoMPhtE0t1Us0EtbD9UZqzMfrUoIMhkNZ0iIWSoocfOQT+8Hj+9bdU+606Cl+4dVLZ9mbrrJS+TpqiYHwUEkQaaTnbXVFEQya6CNnTmu8xZUWU8FBBkstpkjnW7mDZZ6jK7UllQ1AOpSXrTXVCLjlNX30EkOztsn/uWOBQQZNTqZGJjq04pWrCmjmwmm522IrtdnbTUXdWtKB1pqh6aN62pLKNWJ+NpmznVWcWrbN/ZQWnZKS+KpsJIB4qiBuayabaz6zAXpS/sZ9fe/Rt6U7VZM7ruOWkltHHSmsoyaXUaLMfaqFm1YE2Zvs4h21OpbLK9Mnmzoqb/Hns6bemXAoI0MpZMtw916vTLzrNoSus8YZvd+w4cmuo6PTah7LN51T1NR1TDo9Nk9PndadqLeVFAkFFqkll3qVZqW8ee3mddYe2Esn1VZfB120vS1yZUGeWdY9OMXBn+vCkgSC1jrJ7pq7oib92DJudV1eCdfpIPsmMSsucSgkf6s2EW1PR0Ftl2iLalna6qqpZkGhQQZNSadCXNU9WfPruWMFRPIJfX7TQvUw7v79p7eC+hhx5ZP7QwTvocshPtZT+TDgZVpYXsOgdlXWebTv+hjH6eFBCkliHriuuO+C3KcMvmIKorvRRml+Uj8xxz9MYn+7LeRSGTTzcCpwenhXM+IllvMNuOkQ1wy54dtUmgUXAZHwUEmbWqQBaqZLK9cMoy0+wUGNnG4aLBb3UmqwvHzBt7AItunXVLMulptbPnXtUTqqgKaExVhtI/BQRpZIiSQdWU0lXvh0yvaa+btPQi9+kn+C5jHEJjbjqNbfaVzvjTJYi8ksyyMuqi6qomc0AVVblJPAoIshLKBpbBIhiEoFAnc6qT8WV79tSpn08HETi8CqnuqOd0usqOWzf4atzBatBIZYmuaRtCmfR8P1Wfy44ADvXyn7/kxbVGMeelLTsgDMrbDLL7KqqqyhshnDdLa1G6qhSdbxD2UzZquq509VmX/Uh9GqksK2fr9h2H9d5JvwcbM5zQINtl0rg+MrBsELllz/5Dg8iKBr8VPam3qesva5TPHqeqjUOmTQFBoqtbPVMmZKqhi+UFl127obG3aY+btnXyeZPTtQkaZQPJ0nMTde0RVXRtmrTfNDlWep+hgVztB+OhgCCTl1dNkx4RHN7LZp7h322nk26TznDstGwQgcPbMerW4Zdl2GXVW9lrk25fkdWhgCCzkO6jH+rs4dF69rZVHG2fXNMBpw95k9V1mZsoL4imSwphJbmi9o82x9S8R+OngCCjVbfapaiaZuv2HbndR7t2UW1zHlCvXj8Er65VQXklg7xBa3lBNH1tjrDxrTchwzkidgJEipRNBlckHTyqulwOJR1wbtmzn937DnTOVC+/8Llc9boXcfbpJ3LspqM4+/QT+fwlL+bMU05oVRIJ1ynsL33d0mkNYzOqjnHBZdfWrmbKfi/hekl8KiHI6OT1uukyeKmqT/3QQaNoKc2sodJVtd9sEM1e/yaN8dmMfehSl6qd+qWAIJ31fXPmlQzaLvASQ9W0EH3tu49qr6JG57rVddlxCWEiv7PWTiz8TPZ4mg5jPBQQpJEu3Q3rfibb172q8bQoYwmqRt8uS9Ouq7GOH7ZtE8BCFdNQmbyCyLAUEKS1oW7O8PnzLr2ag75x+ubQB39ZXUW7qBrVnKfuNeyzeqnLYLbwPZx36dUAG6bnrlJVklImv3wKCFJLWU+VJp+B+jd6GHFcVV1UlbEMMS3C1IJTV1XfXfiu6raXtKWgMawoAcHMXgK8Bfg54FnurgmKJqjrjKBV0tM27Nq7/7CRyOltitSdKjoY4lzqBMa2I4P7SGeTTLYondlpQLqMCld1UDyxSgi7gAuA90c6/qwNcSO1eTJb9tNcUcZSZ4bQuvXlITiFqpEHH17nvEuvnkRJoeh7KDv3ptcyBIJlfdfSrygBwd1vAzCzGIeXni3z5m9aMoBFph3WNi5Ka3b2zbIg1qSePL2Psn1WVXst46m5zj6LqoT6SKeqg+IbfRuCmW0BtgCsra1FTs24LSPzmMJN2qQvfbhGdbu1pvcdZladSsmgaS+s9L/LutDWaU+SaRgsIJjZF4BTct662N0/U3c/7r4N2AaL9RB6Sp5MVN3Mt0lf+rRQJZK3fTZjBRqNpG5SzZZ9Heupuag6qWhwWx/pnMJDx1wNFhDc/Zyh9i35YmceWX2VWLoMumo6O2gT2UVkxqrqd1Hn+tYJkCopTN/oq4ykm7EEh1i6jNwtei/GNY1VMmgazLW2wbTF6nb6G8CfAScDV5vZze4+7orYCRnLDVknA63z3lgaVKeu6fWvu79Vf+iYk1i9jD4NfDrGsVeF+nQPZxWuoTL71aQqIzlkqJu/7dO/MqXh9PnAMMbvRb+ZdhQQZkqZ6bhN5XvpO31TOe9VpYAgUaqXmgQsZR79m+sDg6pKu1FAmDndCOOyqhlWlxlV535txkQBQSbZnVKZRT/mdv3mWvJZFgUEmayq+YnGaFUzrHCeYS6qNjOqrsq1ikkBQQ6Zwg2XzSymGBTmYswZ9RjTNAUKCDI5TWcyHaMppbUP2UkE28z82vaYq3atu1BAkElpMpOpDENVOvOlgCBL1ddsmG1mMpV40k/9oYRX9b11LRkoYDWngCCtxbzR6sxkKsMYqmFcGXd8CggjMfebYYintrleq7lKtyPcsmf/4FOlzP2eGoICggDNbh4VyQX6Lxno9xSfAkJkq3Iz6KlN6izHOcTxpD4FhAnqM1NtE5CUuUuf+v496XfZngJCZMvIXMd0g4whDRLXGH8DY7pHYlJAmJAhG2bb7GvVbx5pp+i31lfJYO7Vr0NSQBiJIUsGukFE8uke2UgBYUKGrF5a1RtAlmfozLeP+yM9LcoqUkCYMTX+yhCW/Xsa8nhD93SaGgWECVLGLlO0rAeULiWDOpPvzZkCwoDG8qOKfXyZh2XXty/zeGGCxHCsVaWAMFFjCTYiTfU5fmZM7Q9zoIAwAPVckDlqm2m2/f2PPZMea7q6UECYGAUbWWXL6qm0qhQQBjD2J5sp0LUbr6Ylg66Z99h+A3N+KFNAmBgFG5krzaMVnwLCgPRjbW7OT1+rZq6Z91zPCyIFBDO7HPhV4BFgN/Bv3P2+GGmZqjn9CGW1NXkImGMmPCaxSgjXA29w93UzexfwBuD1kdIiIzLnp69VNdfvcI6/1SgBwd2vS728EfjNGOmQw83pxy3TUCdjHaIqUb/1w42hDeGVwCeK3jSzLcAWgLW1tWWlSSLTTSpjN8f2LnP3YXZs9gXglJy3Lnb3zyTbXAxsBi7wGgnZvHmz79y5s9+EdjSHHwEc/uM++/QTgemfl8xL1/tt6/Yd7N53gDNPOSH3t95lbfE690ys/MLMbnL3zVXbDVZCcPdzyt43s1cA5wMvrBMMJI65BDyRvi2rDWGZ92CsXkbnsmhEfr67/yBGGrqaW3Fx2Qugi7TRpWQAGyevO3bTUZx5ygmHSgZbt+8Y7H6eSn4Rqw3hPcAm4HozA7jR3S+KlBbJMZUfsEhsbafbDjOsFolxD8bqZfTkGMft0xDFxTFkuioZyBxV3a9DV/+k9x+CwRgfrsbQy0hGaI59rEViCsHgwYfXuWXP/sp7K8Y9qIDQUZ8lA1XPiAyr6p4a+p5L924aIwUEKaWgJNIn3zAyAAAEgklEQVSPtk/8y7wHFRBGQNUzIqtl974DbN2+Y3T3+hGxEyAiskouv/C5lT2MYhlspPIQxjhSWUSkrlgzAtQdqawSgoiIAGpDEBFZmrG3F6qEICIigEoIkzTWpwsRqWes965KCCIiAqiEMCl1RjSr9CAibamEICIigMYhTFJZyUArnolIlsYhiIhIIyohzIzaEEQkSyUEERFpRL2MZkYlAxFpSyUEEREBFBBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAATG6lsZt8B9jT4yEnAdwdKzhjo/KZrzucGOr+xOd3dT67aaFIBoSkz21lnuPZU6fyma87nBjq/qVKVkYiIAAoIIiKSmHtA2BY7AQPT+U3XnM8NdH6TNOs2BBERqW/uJQQREalJAUFERIAVCAhm9nYzu8XMbjaz68zs1Nhp6pOZXW5mX0/O8dNm9rjYaeqLmb3EzG41s4NmNpsufmZ2rpl9w8xuN7N/Hzs9fTKzK83sXjPbFTstQzCz08zsf5jZbclv8w9jp6lPsw8IwOXufra7Px34HPCm2Anq2fXAWe5+NvB/gDdETk+fdgEXADfETkhfzOxI4L3AecDTgN82s6fFTVWvPgicGzsRA1oHXuvuPwc8B/j9OX1/sw8I7n4g9fJYYFat6O5+nbuvJy9vBJ4YMz19cvfb3P0bsdPRs2cBt7v7N939EeDjwK9HTlNv3P0GYH/sdAzF3e92968k/74fuA14QtxU9WclltA0s3cAFwL/APzzyMkZ0iuBT8ROhJR6AnBH6vWdwLMjpUU6MLMzgGcAX4qbkv7MIiCY2ReAU3LeutjdP+PuFwMXm9kbgD8A3rzUBHZUdX7JNhezKM5+ZJlp66rOuc2M5fxtVqXWVWBmxwGfAl6TqYWYtFkEBHc/p+amHwWuZmIBoer8zOwVwPnAC31iA0safHdzcSdwWur1E4G7IqVFWjCzx7AIBh9x96tip6dPs29DMLOnpF7+GvD1WGkZgpmdC7we+DV3/0Hs9EilLwNPMbMnmdnRwEuBz0ZOk9RkZgZ8ALjN3a+InZ6+zX6kspl9CngqcJDF1NkXufu346aqP2Z2O7AJ+F7ypxvd/aKISeqNmf0G8GfAycB9wM3u/qK4qerOzH4FeDdwJHClu78jcpJ6Y2YfA17AYnroe4A3u/sHoiaqR2b2POB/Al9lkacAvNHdr4mXqv7MPiCIiEg9s68yEhGRehQQREQEUEAQEZGEAoKIiAAKCCIiklBAEGkpmfny783sxOT145PXp8dOm0gbCggiLbn7HcD7gHcmf3onsM3d98RLlUh7Gocg0kEyjcFNwJXAq4BnJLOYikzOLOYyEonF3X9oZluBvwJ+WcFApkxVRiLdnQfcDZwVOyEiXSggiHRgZk8HfonF6ll/ZGY/FTlJIq0pIIi0lMx8+T4Wc+LvBS4H/kPcVIm0p4Ag0t6rgL3ufn3y+j8CP2tmz4+YJpHW1MtIREQAlRBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAAKCCIikvj/3yrkc5LSNHAAAAAASUVORK5CYII=\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-4              # 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": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XeYFdX5wPHvu43eWXpZkBaahQVBwE5Ro5jY0MjPGqMGTSxRTCyoibEkxpholETsvSUoKAKKoqiwgEgRlCYsdellWZbdPb8/Zu4y9+7cPnfv3b3v53n2Yfp5uXf3nTNnzpwRYwxKKaXSS0ayA1BKKVX9NPkrpVQa0uSvlFJpSJO/UkqlIU3+SimVhjT5K6VUGtLkn4ZEZLaIXJ3kGHqKyCIR2SciN0aw/UQRecme7iQi+0UkM/GR1l4i8gsR+SjZcajk0ORfS4nIOhE5aCfJrSLyrIg0jPIYeSJiRCQrASHeBsw2xjQyxjwezY7GmPXGmIbGmPIExFUruX2XxpiXjTEjE1Re0isYKjRN/rXb2caYhsBxwEDgziTH49QZWJbsIJwSdJJTKiVp8k8DxpiNwAdA38B1IpIhIneKyI8isk1EXhCRJvbqz+x/d9tXEENEpJuIfCoie0Rku4i8HqxcETlHRJaJyG67JvgTe/nHwCnAP+3j9nDZt4tdzj4RmQG0dKyrrMWKyFgRKQjY9yYRmWJP1xGRv4jIevsK6CkRqWevO1lECkXkdhHZAjxrL79NRDaLyCYRudouq1sUx7vF/iw3i8gVjrjqichf7c96j4h87th3sIjMtT+rxSJycojPtZ2IvC0iRSKy1tlsJiKDRKRARPba8T0a4ru8XEQ+d+xrROR6EfnB/tzvF5GjRORL+3hviEiOvW0zEXnfjmGXPd3BXvcnYLjj+/2nvbyXiMwQkZ0islJELnSUfaaILLfL3Sgitwb7/yuPGGP0pxb+AOuA0+3pjli17Pvt+dnA1fb0lcAqoCvQEHgHeNFelwcYIMtx3FeBP2BVHOoCw4KU3wM4AIwAsrGaeVYBOYExBNn/S+BRoA5wIrAPeCkwLqC+va67Y9/5wFh7+jFgCtAcaAS8B/zZXncyUAY8ZJdTDxgNbAH62Md+0S6rWxTHu8/+P58JFAPN7PVP2P/v9kAmcIJdbntgh719hv2Z7QByXT6XDGABcDeQY39va4BRjs9tnD3dEBgc4ru8HPjcMW/s/1tj+/9/CJhll9EEWA5cZm/bAjjP/owaAW8C/3Ucy+/7BRoAG4Ar7O/tOGA70MdevxkYbk83A45L9t9Qbf9JegD6k6Av1kr++4HdwI/Ak0A9e13lH6b9x329Y7+ewGH7D9QtYbwATAI6hCn/LuANx3wGsBE4OTAGl3072Um0gWPZK7gkf3v+JeBue7o71smgPiBYJ6CjHMcZAqy1p08GSoG6jvWTsZO5Pd/NLqtbhMc7GPB5bQMG2///g8DRLv/f27FPuI5l032JNmD58cD6gGV3AM/a058B9wItA7Zx+y4vp2ryH+qYXwDc7pj/K/BYkO/sGGCXY97v+wUuAuYE7PM0cI89vR74FdA42X876fKjzT6127nGmKbGmM7GmOuNMQddtmmHdXLw+REr8bcOcszbsJLgPLtJ58og2/kd1xhTgVXzax9B3O2wEsmBgLiCeQW42J6+BKsGWgzkYp0EFtjNKbuBD+3lPkXGmJKAsjc45p3TkRxvhzGmzDFfjFUDb4l1pbTaJf7OwAW+Y9rHHQa0DbJtu4Btf8+R7+sqrKuuFSIyX0R+6nKMULY6pg+6zDcEEJH6IvK03YS1F+uk01SC98DqDBwfEPcvgDb2+vOwrnx+tJv7hkQZt4qS3uBSm7D+MH18te6tuCRqY8wW4JcAIjIMmCkinxljVrkct59vRkQEq/lpYwQxbQaaiUgDxwmgE1bN1M1HQEsROQbrJHCTvXw7VsLqY6z7Hm4Cj7kZ6OCY7+iYjuR4wWwHSoCjgMUB6zZg1fx/GcFxNmBdaXR3W2mM+QG4WEQygJ8Db4lIC4J/drG6Besq8XhjzBb7s1+EVTHApbwNwKfGmBFB4p4PjBGRbGA88Ab+n73ymNb81avATWLdYG0IPAC8btdei4AKrDZfAETkAt+NPWAX1h+5W5fLN4CzROQ0+w/6Fqw25LnhAjLG/AgUAPeKSI59kjk7xPZlwFvAI1ht8TPs5RXAv4G/iUgrO/72IjIqRPFvAFeIyE9EpD5W27qvnFiO59x3MvCofcM2077pWger2epsERllL69r3zzu4HKoecBesW5S17O37ysiA+14LhWRXLu83fY+5bh8l3FqhHUi3C0izYF7AtZvDSjrfaCHiIwTkWz7Z6D9OeeI9cxBE2PMYWAv7r9TykOa/NVkrJuanwFrsWqnNwDYTSd/Ar6wL9UHY3UZ/VpE9mPdHPyNMWZt4EGNMSuBS4F/YNV6z8bqeloaYVyXYLVv78RKLC+E2f4V4HTgzYBml9uxbjR/ZTdPzMSqsboyxnwAPA58Yu/3pb3qUCzHC3ArsATrhvROrBvNGcaYDcAYrOabIqxa8u9w+fs01rMNZ2O1sa/F+mz/g3VDFqwb1svs7+fvWDe+S4J8l/F4DOsG+XbgK6zmL6e/A+fbPYEeN8bsA0YCY7GuCrdw5EY7wDhgnf2ZXov1u6MSSIzRl7koFYxY3VOXAnUCTipK1Wha81cqgIj8zG6KaIZVO31PE7+qbTT5K1XVr7CaX1ZjtT1fl9xwlPKeNvsopVQa0pq/UkqloZTt59+yZUuTl5eX7DCUUqpGWbBgwXZjTG647VI2+efl5VFQUBB+Q6WUUpVEJNTT8JW02UcppdKQJn+llEpDmvyVUioNafJXSqk0pMlfKaXSkCZ/pZRKQ5r8lVIqDWnyV0qpGBljeGtBISWHa97rBzT5K6VqnfU7ismbMJUvVm1PaDkfr9jGrW8u5i/TVya0nETQ5K+UqnW+XrsDgHcWRvu2zejsLTkMwPb9h8JsmXo0+SulVASMMbz/7SbKK7wfCfnL1Tt4dMb3nh83FE3+Sqm08On3RXy4dEvM+z83dx3jX1nEc3PXVS7zakT8i//9FY/P+sGbg0UoZQd2U0opL102eR4A6x48K+p9//fNRu59bzkA2/aVVFkvIjHHlTdhasz7xkNr/kopFcacHxJ74zgZNPkrpWodr1rlV23bR96EqSzftNe9HJeCjDFUJOC+gNc0+SulUt6m3Qd5bd76qPeLozUGgI+WbwVg+Wb35F9ZjmP6pte/oevvp1X2BAq0cfdB9hw8zOOzfmDR+l1+6174cl0c0UZHk79SKuWNe+ZrJryzhN3FpQkv6/MfttNv4nT2HyoLu62vfr9k457KZf/9ZhMA/Sd+xOY9Byk5XE7+H2eQN2Eqe4oPM/TBjzn63o94dMb3/OzJuX7Hu/t/yxLSm8iNJn+lVMrbecBK+tWRFx+dsZJ9JWWs3BKktu+IwdjtPj9s2++66abdJXz+w3a277fi9z1/EMqvXqyeNxhq8ldKpZXn565j3tqdEW0rxNduVFxaRnGUQz/M/G5bXGVGSpO/UqpW2VtymNve+tZv2SpHzfyeKcu48Okvwx6ntMzwvKNPv0/h7oMUl1ZtEpr13VYe/GCF37Jxz8yLMOrqp/38lVK1yoxlW6ssG/fM11WWLSncQ5fcBjSs458Gv9u8D7AevHIz9dvNbNlTwtvXneDXq+iq592ba1J10Det+Sulai1fo83h8qo3C87+5+cc/6eZfsuWb9rLwQiS9YIfd1FcWsZLX/0YdtuNuw4eiSfe7kce0uSvlKoxTEDH+n0lhyktq4j5eAdK/RP9mY/PiXjfP7y7lG8L94TfMAb3/G9pQo7r5EnyF5HRIrJSRFaJyASX9ZeLSJGIfGP/XO1FuUqp9BCsxtxv4kdcEqR5BuDNBYWUlYc/ObxZsIHJn6+NKqZ3F0U2Yqgz9EiHmH7+y/BXFPGKu81fRDKBJ4ARQCEwX0SmGGOWB2z6ujFmfLzlKaWUU8GP/g9KBXa7fG7uurAPe/0u4Aaxl5w9hp5zuYGcLF7U/AcBq4wxa4wxpcBrwBgPjquUUkH9uOOA6/KnPl3tN79x90EOh6j9z18XWbfPWO2qhgfTYuFF8m8PbHDMF9rLAp0nIt+KyFsi0tHtQCJyjYgUiEhBUVGRB6EppWqr8/41N/xGwLNfrGN3sftQCwAXPBW+22c8Uqm27+RF8ne7oAq8tf4ekGeM6Q/MBJ53O5AxZpIxJt8Yk5+bm+tBaEqp2sr31K+KjRfJvxBw1uQ7AJucGxhjdhhjfO85+zcwwINylVJpIpJEX15hWLoxMb1vaiMvHvKaD3QXkS7ARmAscIlzAxFpa4zZbM+eA3znQblKqVomb8JUju3UlHevH8r6HcWc/9RcmtTLrlxvgLLyCrIy/euts1du49vCPdX+KsSaLO7kb4wpE5HxwHQgE5hsjFkmIvcBBcaYKcCNInIOUAbsBC6Pt1ylVO20aP1uAE585BMAtu078nL0U/4ym30lZcy8+SS/Qd6uf3khw7q1rNY4azpPhncwxkwDpgUsu9sxfQdwhxdlKaVqv++CjJ+/r8QaU+fjFf5DOBSXpuYQCqlMn/BVSqWcMf/8IuT6B6atqLJMTwDR0eSvlEo5FW7vRwzj8wifnlUWTf5KKZWGNPkrpVJO6r/+vObT5K+USppV2/bx8IcrqozWGTivvKfJXymVNOOemceTs1dTtO+Q3+ib1fQO87SmyV8plTSVA66JDtdQ3TT5K6USrqy8gn/NXs3BgO6YvtYdQRj0wKwkRJa+NPkrpRLu3UUbeejDFTw268jwC99s2M0Ou7afQm83TBv6AnelVML5XmJ+4FAZT326mvcWb8J5TzeWfv0qPpr8lVLVxhh48IOqT+f+/p3Ev7NW+dNmH6VU4oVp15n53daQ65X3NPkrparNobLwL1NX1UOTv1Iq4Xz1/nXb3d+7q6qfJn+lVML5Wn105M3UoclfKZVQu4tLOXDIGod/eZBx+lX1094+SqmEOua+GckOQbnQmr9SSqUhTf5KKZWGNPkrpVQa0uSvlFJpSJO/UkqlIU+Sv4iMFpGVIrJKRCaE2O58ETEiku9FuUoppWITd/IXkUzgCeAMoDdwsYj0dtmuEXAj8HW8ZSqlaoa9JYeTHYIKwoua/yBglTFmjTGmFHgNGOOy3f3Aw0CJB2UqpWqA/hM/SnYIKggvkn97YINjvtBeVklEjgU6GmPeD3UgEblGRApEpKCoqMiD0JRSSrnxIvm7jdVa+WYGEckA/gbcEu5AxphJxph8Y0x+bm6uB6EppZKlcFdxskNQIXiR/AuBjo75DsAmx3wjoC8wW0TWAYOBKXrTV6nabdhDnyQ7BBWCF2P7zAe6i0gXYCMwFrjEt9IYswdo6ZsXkdnArcaYAg/KVkqlmMJdxby7cGOyw1BhxJ38jTFlIjIemA5kApONMctE5D6gwBgzJd4ylFI1x9XPF7Biy75kh6HC8GRUT2PMNGBawLK7g2x7shdlKqVSkyb+mkGf8FVKeearNTuSHYKKkCZ/pZRnxk76KtkhqAhp8ldKxWXgn2Yy5okvkh2GipK+yUspFZeifYco2nco2WGoKGnNXynlidveWpzsEFQUNPkrpTzxRkFhskNQUdDkr5SKWUWFCb+RSkma/JVSMVtdtD/ZIagYafJXSsVM6/01lyZ/pVTMjGb/GkuTv1IqZvsP6Zu6aipN/kqpmF30tD7RW1Np8ldKxaxMe/vUWJr8lVIqDWnyV0rFZMuekmSHoOKgyV8pFZMftum4/TWZJn+llEpDOqqnUipi+0oOc7C0nCUb93DV8/oa7ppMk79SKmKj/vYZm7Stv1bQZh+lVESMMZr4axFN/kqpsGYs30qXO6YlOwzlIU3+Sqmwpi/bkuwQlMc8Sf4iMlpEVorIKhGZ4LL+WhFZIiLfiMjnItLbi3KVUom3css+3lqgL2qpbeJO/iKSCTwBnAH0Bi52Se6vGGP6GWOOAR4GHo23XKVU9XhrwYZkh6ASwIua/yBglTFmjTGmFHgNGOPcwBiz1zHbAB0GXKka499z1iY7BJUAXnT1bA84qwaFwPGBG4nIr4GbgRzgVLcDicg1wDUAnTp18iA0pVQ8fvemvpS9tvKi5i8uy6rU7I0xTxhjjgJuB+50O5AxZpIxJt8Yk5+bm+tBaEqpeLypbf21lhfJvxDo6JjvAGwKsf1rwLkelKuUSqBte7VPf23mRfKfD3QXkS4ikgOMBaY4NxCR7o7Zs4AfPChXKZVAr83XG721Wdxt/saYMhEZD0wHMoHJxphlInIfUGCMmQKMF5HTgcPALuCyeMtVSikVO0/G9jHGTAOmBSy72zH9Gy/KUUol3q4DpWRnZbDnoL6ftzbTgd2UCvDuokKGdG1JmyZ1kx1KtSorr+D7rfs58/E5NK6bxd6SsmSHpBJIh3dQymH/oTJuen0xl/wn/V5M/sj0lZz5+BwATfxJVjc78alZk79SDuX2C8mL9h1KSvnGGHYXlyal7EUbdielXFVVn3ZNEl6GJn+lUsh/5qzlmPtmsGFncdhtyysMpWUVUZexatt+1/fv7inWNv50oslfKackDzwy47utAGzcfTDstpc/O48ed34QdRmnP/opg/88q8rylVv1nbzpRJO/Ui7cHltPNXN+2J7sEFSCdMttmPAyNPkrFaFte0tYtmlPssPw3OINuznlL7OTHYZyyM5KfPVDu3qqlPWv2aupMIZfn9Kt2so0Idp9hj38CaVlFax78KxEBlCprLyCsgpDdmYGGQIi3ieEon2HGPPEF54fV6U+Tf4qZT304QqAmJP/jv2HqJeTSf2c6H/N3RJtLDdXI/Xb1xaxbkcxWRlWuQL87Mm5LNloXWkM796SF6+qMliuq9KyCrbtK+HJ2av55fCudGnZwHW7gX+aSdcg61Ttp80+qlpt21fCE5+swhjv76ye9tfZjPzbp5XzA/44k7Me/9zzckJZuH4Xi4N0mdx5oJS8CVN54ct1Vdb995tNfLNhNwU/7qpc5kv8ELp9/6s1O/zme9z5AcMe+oRXvl5f2ZyzYWdxld48RfsOsS1JXVrVET1aJ759340mf1UtKioMJYfLuen1b3hk+kq/xFZWXsGUxZswxlBeYfjzB9+xff+RpDRz+daQxz5UVs7eksOsLjrA91v3+61bu/1AVHH6zkm+iv+5T3xB3oSpfLJym2Ob4Ceunz85N2gzysZdVg+eNwq8HTBt7KTwD6QNf/gTTnecGH2i/XyU9+plZyalXE3+Ndy67Qd4O8SY6xOnLONNj5NNLO6Zsoxed33I3oPWk6O+h6kAJs1Zw42vLmLK4k189n0RT3+6hjvfXVq5/q8zvg957LGTvqL/xI88idMXla/R5xu7Fv+PWUcGol1ceOTE9cb8Dcx2nBhCH9s6+tKNe9m8xzoRzFy+le9dulgeKK36hO07Cwsrm8IC+Zqk8iZMDVp+sh5cU8Etu3dU0noXa5t/DXfm43MoLi3nvAEdKC2rIDNDyMw40l793Nx1AJx7bHuyM5N3rn9t/nrAvynDZ9teKym9s3Ajn35fBMCHy7a4HqfkcDkfr9jGmf3aVi5btN77J1OLS8uDJtLyCsO8tTvZuLuY297+FiCim8DOC4Yhf/6YdQ+exdUvFLhue+VzVZff/Ib1Vq0z+7alXwf/J0B73PkBE88OfHW25cChIyeSSZ+tDhunqj4N6iQvBWvNv4YrLi2vnO5x5wdcHKQJYMueEpZv2std/13qeXt7yWErUb701Y9MW7KZvAlT2bHfv5YpIXrO++LxJf5Q7nt/Ode/vJCCdTujjnPOD0VB2+MDYzkUcHPXeQP4Vy8WcOHTX3LT65G94nD/oTIe/GAFh8v9jzk/hv8DwLlPujcrTXxvuevyPvdMr5x+YJr7lYNKLaH+XryiNf9aZl6IhPJ/k+exff8hbjitG60aRT5i5c4DpRwur6B14yP7XPfSAhau38XXvz+dHQessWie/GQVHZvXB+D7rfsZ0rDOkYME/C4fOFROeYXxu0oJxrnFJvvJ130xDDw27pl5AKx+4MyIynVybr59f9Wxd874+xyGdG3huu+jH33P5C/WsnzzXr/lsTbDOJvMVPU7rVcr6uZkMvXbzckOJS5a808j2/dHl2y+27yXwl3FHHf/DI5/wH84gA+WbmGr3Vzz6crwNfbAVHvpM19z53+XWHEdCD2QWYWxmlkitS7gJuaeg4f9ero8Psv9RXIHS8tdEzu4J3yn7zbvZfIXa13XlZZbV2efBVzZXP/ywpDHDKX/xOnhN1IJUTc7kycuOS7ZYcRNk38tlDdhKq/OW0+Zo5lh+MOfVE4LwodLN4dtAjnj73MY9tCR/YpLy1i7/YDfSz4+XrGV379rJfFNe0oq29+37z/E63Y7P0CGS7/5V+dtIG/C1LA1qBVb9nHh01/y5w++Y7Z9otlbcpgnZ6+iwqUWfHLA06pH3/sRZ/1jTuX8Jyu3sf9QGRt3H+THHQd45WsrzrGTvmTUY5+5xhBLr5jvNu9l5N8+ZfU273vU6JDLsclrUT/ZIVQRzVW4l7TZp5a6450l3PHOkqDrr33JqnUGu1FZcri8yrLed1u1TeeDQb6btT6l9gnnhlcXAXDCUS3p2Lw+B12OF62nP11TOT1xyjJ2FR+mV5tGEe1buOvIQGnfFu6h7z3+Neez+rf168UTL/+bxfuDbqeq1yu/HMwJD34cdrvebRtXaabzCfUUeCz+esHRvL9kE2/M3+Dp72A4WvOvJaKpmd7wqntzw9rtB3hv8SYe+nAFn4d4qGiNo6wJIU4wAG8uKHTtyhivXXYzjluvGIBnPndvgglKm9HTQqRf8yu/jOxp6ni0b1oPgCb1s/nF8Z2rvTee1vxrqK17S/yaUv4Wpi+801drjrSf3/rmYt5aUMib1w7hgqe+9DRGsNrXg7WxJ9L977v3fAlmaS0csC3dNamX7foe4r7tG7N0o3ut3qc6etskm9b8a6CF63dx/AOzGPinmZXLpizeFNOx3rIfEEtE4q9JfvGfr5MdgvLAortGVE6LQPdWVYdOeP+G4dUZUsrS5F8DfLh0M5+ssJ4iXb5pLz9/cm6SI1IqNTVrkOM3P+Pmk1gycWT0B/K44n9sp6Zht6nulkdPkr+IjBaRlSKySkQmuKy/WUSWi8i3IjJLRDp7UW66uPalhVzx3HwA1kfwej+l1JEeZo3qZlcuy47w+Y46WbGlxtN6taqcXnz3kZPOu9cP5ZfDu/ht+9jYY2IqwytxJ38RyQSeAM4AegMXi0jgc+aLgHxjTH/gLeDheMtVSimnt687wW/eLc23ahy+W+XRHZtSN8Rga51bWL3dPr7lpCrrIn3nQtsmdRmY1zzo+stOyIvoOPHwouY/CFhljFljjCkFXgPGODcwxnxijPFVWb8COnhQbtpxGwpYKWUZ0LmZ33yLhkeagG44tZtfrTyUDnYvnGBuHtEDgK5hXrUY2CV0VJ82ldOhriwmjRtAN5d7FV7zIvm3B5zDRhbay4K5CnB967SIXCMiBSJSUFQU/qnRdHP3/5ZxqCz+/vJK1WTtwyRnH+dwJLeM7Mkzlw+MqpxHzu/P/3491G9Zj9YNQ3bJDFXxz3fU9ENdITQPuG+RKF4kf7f/heu9CxG5FMgHHnFbb4yZZIzJN8bk5+bmehBa7fOb175JdghKecrtQcPBXYM3iUw8p08iw6l0QX5Hju7YlD//vB+d7SeDzzm6XbWUXR28SP6FQEfHfAegSr9DETkd+ANwjjFGBxZXSgUVqp/9iN6t+fC3/t0127i05ceaqDsHDAFx8aBO/LR/2yBbR3YMn0njBgDw6IVHV1k3ondrANpGeGUTLy8e8poPdBeRLsBGYCxwiXMDETkWeBoYbYyJ7M0XSikVRK82jSun7zzrJ5xzzJFE37pxHbbuPcTw7rG1HjjfFRGtMce0Y8byrbx45fG4jZw+sk+boEOq/OrErlw8sBNN6me7rvda3MnfGFMmIuOB6UAmMNkYs0xE7gMKjDFTsJp5GgJv2m1d640x58RbtlKqdgrWJP7y1VWHXbh6eFe/+ab1cti69xBRjtrtiZ/2b8dP+1snol32aLVNI0zmIlJtiR88Gt7BGDMNmBaw7G7H9OlelJOOQo2xo1S6GdqtZdhtnrtyIDOWb42oW2c60yd8U9jc1du59BkddkClnwi7y7tq26Qe/zckz7NYaitN/inskn9r4lc1zzEdww9l8Mj5/eMuZ0Tv1jTICf4wVjReuup46mZb6dCrt5zm2H35B3RqFmbL5NDkn6Jmr9T74qpmeui88Ik91NOtgYZ2c3895r//L59l942O+DihDOveMuSDVbGM8tmgThbTbhzOPy45Np7QEkaTf4ravKck2SEoFZOeEb5gJ5SjHE/Pvnz14LiPF0okXUJ/dVJXxg7syBVDuwTdJrdRnSrLerdrTP2c1Bw5PzWjUp5deiqVLBkCgW/Z7NyiPj/uCD84YXUlzCUTR1IvxDg+Po3qZvOgyxXNPWf3pmfrRlxSA4cE1+Sfogp+jPyF5Uqlmi8mnEq97EyOu3+G3/JIKzUXDezIU5+urpx//OJjE9J10zniZyyuGNqFbXtr5lW6Jv8U9OHSzbyzcGOyw1AqZuHG3wnXm6eL4z3RUL3DKnj9jt5Upck/BUVyWaxUbTS4a3NO62UNc5Ah0K9D+J5DXkmHVzc6afJPQelR71CqqteuGVI5vebP7sMgpKqadp9Oe/ukoP0lZckOQamYPPCzfskOofrV0AsGTf4pSF/VqGqqcOPY/PHcvnRr1ZC2Tapn5EoVnCb/FFRe064fVVoaaQ9BHI0Te+Qy8+aTKp9+TSW+fvqpGFsipMf/soZZvmlvskNQKqyHz+9PdmYNbfNw8bcLj+Gh8/r5DRcdnZpVadPkn4LWbj+Q7BCUCqtp/ZyoHsY6uWdqv52vSf1sLhrYKer9amovIU3+KebdRYXJDkGpiBltoqyxNPmnmJteX5zsEJSKWGDq13NBzaHJXykVO4+SfZN61fcGq0SpaSc+fchLKRW1v11U9QXkEPtLWGbefBLb9tXMMXLiefFMMmljMLmgAAAU/ElEQVTNXykVNV9NPVyzz2e/OyWi4+U2qkOfdk08iCwyf7ngaKbdOLzayktFWvNPIYfLK5IdglKe6tSiftB1b107hK/W7KjGaI44f0CHpJSbSjT5p5AVm/clOwSlIuKr4fdo3ZCF63dXLo+mCSQ/rzn5UbzRS3lLm32UUhFZeu8oTjjK/5WKmQGD7Lvd9Jw0bkAiw0oZNex+rzfJX0RGi8hKEVklIhNc1p8oIgtFpExEzveizNqoaH/NvOGl0kPDOllce9JRAPTrYLXP+5J9y4ZVX2Hok51Zu+uYNfR+b/zJX0QygSeAM4DewMUi0jtgs/XA5cAr8ZZXm135XEGyQ1AqpBN75LLuwbNo1agucKS2m5WI12yphPKizX8QsMoYswZARF4DxgDLfRsYY9bZ6/SOplK10MRz+jB39XZG9ol+sDeVHF5cj7UHNjjmC+1lURORa0SkQEQKioqKPAhNKZVIQ+17AL3bNua+MX1rfRNPKDVtqAsvvim3672YPgVjzCRjTL4xJj83N7UHgVJKwW9O78Gc204J2aWztpMa+pSXF8m/EOjomO8AbPLguGllyx692atqnswMoWPz9E38NZkXyX8+0F1EuohIDjAWmOLBcdPK9GVbkh2CSmN3/zSwj4blJ21jHdtepbq4k78xpgwYD0wHvgPeMMYsE5H7ROQcABEZKCKFwAXA0yKyLN5ya5sM7S2hkqh5gxzX5Tec2q2aI6l5fM86tAjR3TUVefKErzFmGjAtYNndjun5WM1BKoi7/rs02SGoNBas6aZLywYA1EmTVxvGokm9bB4+vz/Du7dMdihR0eEdlFK49dGYfevJtGlSNwmx1DwX5ncMv1GK0dO5Usp1WIa8lg0qx+rp1qph9QakEk6Tfw3xxCXHJTsElUJ+fmxMj9IEFaxvdp2sTF68ahAvXnV8HMeuWf3f04Um/xSw4MddYbfJzBA++E16jz9enR676Jio9/n72Oj26dWmUdRl+NwVpHdONN4bPyyi7YZ3zw16Qzga2qUhtWjyTwHn/Wtu2G36tGtc414TV5Od2CP4Q4atG9cht1HVnh2j+7aJqox4noZtVPfI7brbRveMeD/fFcONp3X3G5GzcV3r5Sw6Rk/60OSfop67YmDl9MWDOtKxeX29fA6hOnPW178/vcqywKGNIxHt9zms25HeJFmOE4dEUad+9KJjWPfgWdw8ooff8p5tGvH8lYP4+9hjo4pJ1Vya/FPQJcd34uSerXjk/P4AnD/A6kngq/nnaLe7KiL5THp7+MBS4FWYuCyL9hjhRHPTdd2DZ1UZez9Qwzr+nf1O6pFb2aWzZ+vYm6RUzaBZJAVde6I1ZvoF+R1Z9+BZDOjcLMkRWU7uGfl4SyN6V9/ojr89vXtE2x0VRfIMrEu/f4N/+/jQbv6JtX5OZsTH9klUM949Z0d2P6BTi/rcProXc26r+p7dDs3qeR2WSjGa/JPMbSTAYINkdc1tQIsGOdw+uleVdesePMt1n3vP6UN2ZuhmgcAbj09dGvublxrVyeLLO051XTeka+iaaDi/PuUopt5Y9SblDad2T/j9kL7t/V8u/vD5/f0efGrXtB4ZAQN8De6a2FcUBmtpumJoFyCyk8t1Jx/l94CXNiymD03+Sfb2wo0Rb1s/J4sFd43gpB7+TxI2rZ8ddJ/LTsjjuE6hrxwCa8Sj+7Zh0rgBVZJ4JMnkwoEdadukHq0Cboi+ee0Q2ja1Hhga3LW531XEoC7uSTLwZPG7Ub3o065Jle0yMySipBXNkLvhtqyTlVn59KtPTlYGr18zmH/94ki33Aln9OL/hnSunD+6w5H44020vvff1s0+8mcceEUCcOngTkBkN3N9n1ENHahSRUGTf5Ld/va3fvMvXx19f+pRvaPrZeKz+J6RvHDlIB4+r3+VZouRfdrQton7pf+Z/YKX97tRVs+TO8/qzV8uONq1yer8AR158hfHHWlaMPAzl37rXiSgJ39R9fmIrrkNXLYM7s1rh0S87fFdW9DEPhkbA9eedBT3jenLNSd2BeCMfm0rt71kUHRPhQY7SfZq05iHz+vP4rtH8vLVgyuXt2psnYDPO64Di+8ZyTf3jIyiNM3+tZ0m/yQrr/Cv/0XSvt+68ZFH7p+9YiB//Fnfynm3vuO3jOxJC5d+2k3qZXNij1wa1Mli+X2jI465XZCTAkDdbOskUi8nk/MH+A/ndPWwrjSpl81JPXKpn5PFXy84GrB6vTxyfn+W3juKaTceeZbhgZ/1CxvLKb4riCDV6DMdyda3yW9P78GaB86sXD7z5hNDljEwL3zzzYQzqjbFuXGm1HFD8iLax6d9U+tz79/B/+pHxLriahJwBfinn/XjrxcczbGdmtGkXnaVG7xutNknfWjyTzG+5BlKo7rZlV31+rRrXNlf/Lv7RvOES013UJfmLLhrREzdEd3EkiAE6N2uMYvvGVnZR943kqkxVtfFhnWy6N3uSI+cvJYN6BRmrPhnrxgUtuw/ntu3sueULxbnKKrdWgXv2dLMkVAbBUme024czsk9W1VZ7rxy8TW9Bd47iEa0n3vDOlmcNyC28RS12af20+SfRJM/XxvzvuNP6cbiu0dWvkgbrNp2qHbdcF3/wvE1zfRoHXmvmbb2wGD1XHrD+CKtCNEWP8weKfGVX8Y+vMClgztzQX7HuKu1s245iSnjh8a07+i+bZj3h9MY2i3+kR99n9uoPlbzm5c9c3zNf83rx/9Er0ptOqpnEt33/vLwGwWRkSFVLvPDeXrcADbtPsjpj34WU5nnHtuec49tz6pt+wA4++h2XH5CHnNXbeevM7533efB8/ozondr1xu1vtplYE5+/4Zh7C4+DFi9la476aiwb4tyPjD16IVHc/Mbi4PWXmOt1bZqXJdWjWMf5dJ5oo5Fph247+rwyqF5XJjfgUZ1o/s9CGVYt5b88dy+rvdgVO2iyb+WaVIveCKon5NFt1aNmH3ryawu2l9l/dQbh1XpruimW6tGzLrlJPJaNCAzQ9i4+2DQbRvWyWLMMe6JxPfu08CKv7NpJDszI6LXBDqP4Ws6CvyfRPNEbTwv4+7Vxir/qmFdXddPPLt3le68+Z2bURBmjKe+7Rtz84geXDTQulEsIp4mft8xLx3cOfyGqsbT5F/LNK2fw7zfn8agB2YF3SavZQPyWlbt8eJWOw/mqNz4h/htbI9P09njl393D9GGD9ENhxDLy7mbN8gJ+twFwOV2P3ynQV2aB03+Azo34+lxAxARbjwtsgfalApHk38KGX+KN6/Ma9W4LuNP6Ua/DrHfXAw046bQPWJi0a1VIyZfns/gOB/+gsia81N5YDyDNbxCZobw8YptfuvqZGXQsoa9IlClPr3hmySrtlVtdrl1VOSjM4Zz66ielTcEvdA9QWO9nNqrNfVzqrcOEqoyH+5paKe7f9qbzi3qR/3cQDDPXzmIyZcPDL+hUh7Qmn+SvLOwMNkh1Fj3n9uXfu2bcJQj6bq10Qd2bY2k5v/57aeyff+hiOI4oVtLPv1d1XFxYqE9K1V10+SfJE/OXu03f3THpkmKpOYZ53JD8vguLfhyzQ7ASvo3ntadUX38B5fLsmv1oR53aN24Lq0b12VfidXbqE8770YCDeW6k4+qsqxZ/Wx2FR/m4kGdqiUGlV40+aeId647Idkh1Gj/viyfvvdMr5wPHK8eYOI5fchtVIfTfhJ+xNFGdbN5+7oh9GxTPcnfrdfO1BuH066pjq6pEsOT5C8io4G/A5nAf4wxDwasrwO8AAwAdgAXGWPWeVF2TbRhZ3GVZV49fZsM3VPg5d6RDF3QsmEd7jm7T8THHNA5saNyArx93QkU7qr6+wBUGRxPKS/FnfxFJBN4AhgBFALzRWSKMcb5BNNVwC5jTDcRGQs8BFwUb9k1UXFpGcMf/sRvWTzvck0FP/HwJSnVad4fTuPQ4YqkxjCgc7OUeV+DSi9e1PwHAauMMWsAROQ1YAzgTP5jgIn29FvAP0VETDxP0gRhjGHqks3UycokQ47UqA+XW0WVV5jKNl8DZIhwuLyCeWt30q99E+rlZFbeGCzaV8L2/aU0rZ/Nxt0HefaLdV6HC0Q3aqTyTrxP3CpVk3mR/NsDGxzzhUDgQCyV2xhjykRkD9AC2O7cSESuAa4B6NQptptcRfsPMf6VRTHtmwzv3zDM86c0vda6cfjmhynjh6b8/0MlR14Lq1eW2+B3Knm8SP5ujdWBNfpItsEYMwmYBJCfnx/TVUGz+jl88JvhFJeW+y339d/OzBC/Ln8VxnC43LDnYCmtG9clM0MQBINh466D1M/JosIYSssr+HrNThrVzeJQWQXfbNhNnawMurZswOqi/WRlZNCkXjYrt+6jeYMcBCguLado/yGO7diUcmPYe7CMutkZvP/tZsYN7sz95/Yl1b193QlhR9YE6N9Beyt5LZani1NR19yGLLprRMiXDqnq50XyLwScb6XoAGwKsk2hiGQBTYCdHpRdRXZmhmdt0L0Cenqc4lHN5Z+XeHKYaqHt0coLzVzeJ6GSy4snfOcD3UWki4jkAGOBKQHbTAEus6fPBz5ORHu/UvE4tZc2S6j0EXfN327DHw9Mx+rqOdkYs0xE7gMKjDFTgGeAF0VkFVaNf2y85SoV6KObTmRJ4Z6Y9l1x/+iI3nGrVG3hST9/Y8w0YFrAsrsd0yXABV6UpVQwPVo3okeMYxBF8gY1pWoTHdhNqRSl1yEqkTT5K6VUGtLkr5RSaUiTv1JKpSFN/koplYY0+SuVomrJA74qRWnyV0qpNKTJXyml0pAmf6WUSkOa/JVSKg1p8lcqRdWWIZ1VatLkr5RSaUiTv1JKpSFN/koplYY0+SulVBrS5K+UUmlIk79SSqUhTf5KKZWGNPkrlWLqZOmfpUo8T97hq5Tyzns3DOOz74uSHYaq5TT5K5Vi4nkRvVKR0utLpZRKQ3ElfxFpLiIzROQH+99mQbb7UER2i8j78ZSnlFLKG/HW/CcAs4wx3YFZ9rybR4BxcZallFLKI/Em/zHA8/b088C5bhsZY2YB++IsSymllEfiTf6tjTGbAex/W8UfklJKqUQL29tHRGYCbVxW/cHrYETkGuAagE6dOnl9eKWUUrawyd8Yc3qwdSKyVUTaGmM2i0hbYFs8wRhjJgGTAPLz8008x1JKKRVcvM0+U4DL7OnLgP/FeTyllFLVQIyJvYItIi2AN4BOwHrgAmPMThHJB641xlxtbzcH6AU0BHYAVxljpoc5dhHwY8zBQUtgexz7J0IqxgQaV7RSMa5UjAk0rmh4FVNnY0xuuI3iSv6pTEQKjDH5yY7DKRVjAo0rWqkYVyrGBBpXNKo7Jn3CVyml0pAmf6WUSkO1OflPSnYALlIxJtC4opWKcaViTKBxRaNaY6q1bf5KKaWCq801f6WUUkFo8ldKqTRU65K/iIwWkZUiskpEgo0yGm8Zk0Vkm4gsdSxzHd5aLI/b8XwrIsc59rnM3v4HEbnMsXyAiCyx93lcRCSCmDqKyCci8p2ILBOR36RIXHVFZJ6ILLbjutde3kVEvrbLeF1Ecuzldez5Vfb6PMex7rCXrxSRUY7lMX3nIpIpIot8Q42nSEzr7M/4GxEpsJcl9Tu092sqIm+JyAr7d2xIsuMSkZ725+T72Ssiv02BuG6yf9eXisirYv0NJP13qwpjTK35ATKB1UBXIAdYDPROQDknAscBSx3LHgYm2NMTgIfs6TOBDwABBgNf28ubA2vsf5vZ083sdfOAIfY+HwBnRBBTW+A4e7oR8D3QOwXiEqChPZ0NfG2X9wYw1l7+FHCdPX098JQ9PRZ43Z7ubX+fdYAu9vecGc93DtwMvAK8b8+nQkzrgJYBy5L6Hdr7PQ9cbU/nAE1TIa6Av/0tQOdkxgW0B9YC9Ry/U5enwu9WlVhj2SlVf+wvabpj/g7gjgSVlYd/8l8JtLWn2wIr7emngYsDtwMuBp52LH/aXtYWWOFY7rddFPH9DxiRSnEB9YGFwPFYTzJmBX5vwHRgiD2dZW8ngd+lb7tYv3OgA9Y7KE4F3rfLSGpM9rbrqJr8k/odAo2xEpqkUlwBsYwEvkh2XFjJfwPWiSTL/t0alQq/W4E/ta3Zx/fB+xTay6pDsOGtg8UUanmhy/KI2ZeOx2LVspMel1jNK99gDfw3A6vmstsYU+ZyrMry7fV7gBYxxBvOY8BtQIU93yIFYgIwwEciskCsUW4h+d9hV6AIeFasZrL/iEiDFIjLaSzwqj2dtLiMMRuBv2ANd7MZ63dlAanxu+WntiV/t/a4ZPdlDRZTtMsjK0ykIfA28FtjzN5UiMsYU26MOQartj0I+EmIYyU8LhH5KbDNGLPAuTiZMTkMNcYcB5wB/FpETgyxbXXFlYXVzPkvY8yxwAGCv7WvOuOyCrPaz88B3gy3aaLjsu8vjMFqqmkHNMD6LoMdp1o/K6falvwLgY6O+Q7Apmoqe6tYw1oj/sNbB4sp1PIOLsvDEpFsrMT/sjHmnVSJy8cYsxuYjdXe2lREfEOKO49VWb69vgmwM4Z4QxkKnCMi64DXsJp+HktyTAAYYzbZ/24D3sU6WSb7OywECo0xX9vzb2GdDJIdl88ZwEJjzFZ7PplxnQ6sNcYUGWMOA+8AJ5ACv1tVxNJWlKo/WDWUNVhnXd/NkD4JKisP/zb/R/C/yfSwPX0W/jeZ5tnLm2O1ozazf9YCze118+1tfTeZzowgHgFeAB4LWJ7suHKBpvZ0PWAO8FOsWprzBtj19vSv8b8B9oY93Qf/G2BrsG5+xfWdAydz5IZvUmPCqiU2ckzPBUYn+zu095sD9LSnJ9oxJT0ue9/XgCtS4Xce637WMqz7W4J1o/yGZP9uucYay06p/IN1R/97rHblPySojFex2vMOY52Jr8Jqp5sF/GD/6/vlEeAJO54lQL7jOFcCq+wf5y9vPrDU3uefBNxoCxLTMKzLv2+Bb+yfM1Mgrv7AIjuupcDd9vKuWD0pVtl/GHXs5XXt+VX2+q6OY/3BLnsljl4X8Xzn+Cf/pMZkl7/Y/lnm2y/Z36G93zFAgf09/hcrSaZCXPWxholv4liW7N/5e4EV9n4vYiXwlPh9d/7o8A5KKZWGalubv1JKqQho8ldKqTSkyV8ppdKQJn+llEpDmvyVUioNafJXSqk0pMlfKaXS0P8DuLBPAWlOwtQAAAAASUVORK5CYII=\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 0.4084425498386879 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = mg.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": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAELCAYAAAARNxsIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VFX6wPHvmxBq6KAgRYIiCgpI6IKwKyBSfwoKKggW0FVWWcUVFF0FdVnLiq5YUFFWUHBFBRRBLGABFBBUikhVQu8tCWnv7497EyZ9gJncycz7eZ55csu5977JTPLmnnPPOaKqGGOMMYEQ5XUAxhhjwoclFWOMMQFjScUYY0zAWFIxxhgTMJZUjDHGBIwlFWOMMQET1KQiIt1EZL2IbBSRUQWU6yciKiIt3PV6IpIkIqvc1yvBjNMYY0xglAjWiUUkGpgIdAESgGUiMltV1+YoVx64G/g+xyk2qWqzYMVnjDEm8IJ5p9IK2Kiqm1U1BZgO9Mmj3DjgKSA5iLEYY4wpAsFMKrWAbT7rCe62LCJyKVBHVT/O4/g4EVkpIotEpEMQ4zTGGBMgQav+AiSPbVljwohIFPAcMCSPcjuBuqq6X0TigY9EpLGqHsl2AZFhwDCAcuXKxV944YWBit0YYyLCihUr9qlq9UCdL5hJJQGo47NeG9jhs14euBhYKCIANYDZItJbVZcDJwBUdYWIbAIuAJb7XkBVJwGTAFq0aKHLl2fbbYwxftl22KlUqVOxTiElw4+I/B7I8wUzqSwDGohIHLAdGADckLlTVQ8D1TLXRWQhMFJVl4tIdeCAqqaLSH2gAbA5iLEaYyLYoA8HAbBwyEJvAwkDQUsqqpomIsOB+UA0MFlV14jIWGC5qs4u4PDLgbEikgakA3eo6oFgxWqMiWxjLh/jdQhhQ8Jl6Hur/jLGmFMnIitUtUWgzmc96o0xEW/zwc1sPmg17IEQzDYVY4wpFm6ZdQtgbSqBYEnFGBPxHuv0mNchhA1LKsaYiNexXkevQwgb1qZijIl46/etZ/2+9V6HERbsTsUYE/Fu//h2wNpUAsGSijEm4j15xZNehxA2LKkYYyJeuzrtvA4hbFibijEm4q3es5rVe1Z7HUZYsDsVY0zEGz53OGBtKoFgScUYE/Ge7vK01yGEDUsqxpiI17JWS69DCBvWpmKMiXirdq1i1a5VXocRFuxOxRgT8UbMGwFYm0oghE9SWb8eOnXKvu266+DOOyExEbp3z33MkCHOa98+6Ncv9/6//AX694dt22DQoNz777sPevVyrn377bn3jxkDnTvDqlUwYkTu/U8+Ce3aweLF8OCDufdPmADNmsHnn8Pjj+fe/+qr0LAhzJkDzz6be//bb0OdOjBjBrz8cu79778P1arBW285r5zmzoWyZeGll+C993LvX7jQ+frMM/Dxx9n3lSkDn37qLI8bB198kX1/1aowc6azPHo0LFmSfX/t2jB1qrM8YoTzM/R1wQUwaZKzPGwY/PZb9v3Nmjk/P4CBAyEhIfv+tm3hn/90lvv2hf37s++/4gp4+GFn+aqrICkp+/6ePWHkSGc55+cO7LNXzD57E2KPOetvdbLP3hkKavWXiHQTkfUislFERhVQrp+IqIi08Nk22j1uvYhcGcw4jTGRrdmxWJodi/U6jLAQtEm6RCQa+A3ogjNf/TLgelVdm6NceeAToCQw3J1OuBHwLtAKOAf4HLhAVdPzu55N0mWMOV3Lti8DIrPBvjhN0tUK2Kiqm1U1BZgO9Mmj3DjgKSDZZ1sfYLqqnlDVLcBG93zGGBNw9y+4n/sX3O91GGEhmG0qtYBtPusJQGvfAiJyKVBHVT8WkZE5jl2a49hawQrUGBPZXuz+otchhI1gJhXJY1tWXZuIRAHPAUNO9VifcwwDhgHUrVv3tII0xpiLz7rY6xDCRjCrvxKAOj7rtYEdPuvlgYuBhSKyFWgDzHYb6ws7FgBVnaSqLVS1RfXq1QMcvjEmUizetpjF2xZ7HUZYCOadyjKggYjEAduBAcANmTtV9TBQLXNdRBYCI92G+iTgHRH5N05DfQPghyDGaoyJYA9+4TxWHQn9VJJT0ykRJZSIDs49RdCSiqqmichwYD4QDUxW1TUiMhZYrqqzCzh2jYi8B6wF0oC7CnryyxhjzsSrPV/1OoSA2nskiWrlSyOSuyWhz/h51Ktenof6Nad86ZiAXztojxQXNXuk2BhT3KzZdoDoqCgurFWJXYcSSc9QalUpd0rnUFW+XbeL1hecxcFjJ9hxMJFRU7/nnh6X0L15XW6e+BU7DiQCUKlcSQ4dT8l2/GeP9AzoI8Xh06PeGGNO06KtiwDoWK9jkV3z971HufetJbm2zxndjZjoKI4mpZKankHV8qVJS89g2ca9PDHzR+LOKs9/bmsPwPETqTw7+2e++3UXFcrEcCQpNes87y/ZzNkVy2QlFCBXQgkGSyrGmIj3j4X/AALfpvLSvDXMWraV+Q/3yLVv2Ctf53lMr3/Oy7Y+Z3Q3Jnz8C1/8sh2A33Ye5sCxZK5/LvvwM74JBWD7geM8+E7RN0VbUjHGRLzJfSYH5byzlm0F4I99x0hOSeOCcyoB8L8lm/w+x31TlrDvSHK2bTkTSiixpGKMiXj1K9c/7WOTUtJY/ccBdh1K5ItftjPh5stylRn6slO9NqZvc9pdWIPXP//V7/P/tuPwacfmBUsqxpiI9/nmzwHoXL+zX+UTT6Sx/2gydarF8tK8NXz208mRiLcfOE7lcqUoWyr3n9fHZ/5Ip8bnBCboM9CsXlVWbd1feMHTYJN0GWMi3uNfP87jX+cxxL+PE6npTP16Aylp6Vz91Hxue3kRySlprNl2MFu5WyYu5Oqn5ud7noVrcvXjLnL/GtSG3i3PDcq57ZFiY0zE23bYGaawTsU62bbvOZxEpXIlKVkimjsnfcOm3Uf8Puezg9ty35TcT3fl9MH9XSlbqgQ3TPiCA8dOnFrgp6Br09p89lMCr/+lI3WqnRzmP9CjFFtSMcaYHLbsPkKFsiW5YYLTIH5f7yY8O/vngF9n3pjuWR0UT6Sm03u88+TXzPu7sudwEn+Z9E1W2XdGXEFKWgbRUcKgF74EoOX51RnTtzkK/N+/nLuj4VddTOKJVCZ/uT7r2P6XncfAyxtQskR0rhgCnVSsTcUYE/Hm/PoJGQp9LnIe/b3D5485EJSE8u7frsjW471UTDQvDe1AxbIliS0dQ6xPb/erW8dRtXzpXOd4/PqTM4LEli7BwMsvoFcLp1qrb5v69HjSmQHzlj9fGPD482NJxRgT8QbPcCamfaNbI175bG0hpQOjSmzuJHFejQp5lr29y0XZ1kdffSmp6RnZts28P/sEuSWio7imdRzN61ejKFlSMcZEtP1Hk2kiDwCcUUKJLV2CY8lpubY/NagNf397aR5HFG7uQ91R1VxjeHW62L8nyG7v2ui0rnsm7OkvY0zEylDljS9+pZRUppRUPqNz5bxTyHRhrUqM6decuLPKA1ChTAwP92vu1zmjgziacLDYnYoxJmJd9fhcAPbo9wCcJa0LKp6nnvF1OeoOkfLPG1sz98c/+GbdToCs4Vk6XFSTVuefxZert9OtWZ08Rw8OF5ZUjDER58CxZJZt3Ju1/rt+CPifVMYOaMEj052nTf/a/ZKs7c3rV+PSuKrcNyWZu322g9MQf9Wl4T9DrSUVY0zEyTl2VlMZfUrHn12xbL77RIR/D2l3WnGFg+JVWWeMMachJS2dK8d9whtf5D3mVkmpSEmpmGv7DR3OByD+vOr847p4ujStTYkood5Z5bmjayNm3OvfsC6RJKh3KiLSDXgeZ+bH11V1fI79dwB3AenAMWCYqq4VkXrAOiCz985SVb0jmLEaY8JX/387Y3u9t3gT/S87L9f+3bqY69ufT/0yl1OhTEn+Pcfpl9K6wdkM7tQwq1y7hjUY2bsp4PQdMbkFLamISDQwEegCJADLRGS2qvo+s/eOqr7ilu8N/Bvo5u7bpKrNghWfMSZyJJ44+ahv36c/y7X/D53Np39UZeGQmwHYfSiJad9s4Px8+o2Y/AXzTqUVsFFVNwOIyHSgD8688wCoqu9AOuWA8BgzxhgTMlZt2VdomY33f0101MnWgJs6XcBNnS4IZlhhK5htKrWAbT7rCe62bETkLhHZBDwF3O2zK05EVorIIhHpkNcFRGSYiCwXkeV79+7Nq4gxJoKpKg9M/b7QclXKVqZi6dxtKubUBTOp5PUgdq47EVWdqKrnAQ8AY9zNO4G6qnopcC/wjojkug9V1Umq2kJVW1SvXj2AoRtjwkHC/uP57mtQ00kiZ1csw4zVM5ixekZRhRXWgln9lQD4jiNdGyhoIoHpwMsAqnoCOOEur3DvZC4AbBhiY4xfdh1K5DZ3xsW8vHhbew4nplCxbEk6vdUJgP4X9y+i6MJXMJPKMqCBiMQB24EBwA2+BUSkgapucFd7ABvc7dWBA6qaLiL1gQbA5iDGaowJA2sTDlK7SjkqlC3JTzlmNpw3pju//HGAxet3UyLKqUipWLYkAHNvnFvksYaroCUVVU0TkeHAfJxHiier6hoRGQssV9XZwHAR6QykAgeBwe7hlwNjRSQN53HjO1T1QLBiNcYUf1MX/cbbXzv/o868vysf/bA1a19s6RhEhCbnVqXJuVVzHVs2Jv/OjObU2CRdxpiwcOW4T/Ld9+rtl1PPHdAxL1N/ngrAwCYDAx5XqLNJuowx5hT0iK9bYEIBeP3H14HITCqBZknFGFOs/br9IPdMXpzv/iub1cl3X6YFgxYEMqSIZknFGFOsFZRQAKpXyD3DYk4x0TGFljH+sQEljTFhrUKZkoWWeWvVW7y16q3gBxMB7E7FGFNs5RyC5ZK6Vah/dgXWbT/Ifb2ack6Vsn7NnJiZUIY0GxKEKCOLJRVjTLF0NCk12xAs59eowNM3tTmtWRUXDlkYwMgimyUVY0yxtGn34WzrE4fmOUSgKWLWpmKMKZY+W5UQsHO9tuI1XlvxWsDOF8ksqRhjiqVdhxKzls87+8zmPZmxZgYz1tiAkoFg1V/GmGJpzbaDWcv39mpyRuf6/KbPzzQc47I7FWNMsZKhmm1Ilsf6t+D8mjYXSqiwpGKMKVbWJRzMtt7ivDOfS+mlZS/x0rKXzvg8xpKKMaaYSc/IPgiuP/1QCjPntznM+W3OGZ/HWJuKMaYYSUlL5/7/Lg34eT+98dOAnzNS2Z2KMabY6PXPednWZz1wpUeRmPwENamISDcRWS8iG0VkVB777xCRX0RklYh8KyKNfPaNdo9bLyL2yTHGZDNndDdKlwxMZcvzS5/n+aXPB+RckS5oSUVEooGJwFVAI+B636ThekdVL1HVZsBTwL/dYxvhTD/cGOgGvOSezxhjAChZInB/Er7Y8gVfbPkiYOeLZMFsU2kFbFTVzQAiMh3oA6zNLKCqR3zKlwMyW+D6ANNV9QSwRUQ2uudbEsR4jTEhbMXmvVnLc0Z3C+i5Z18/O6Dni2TBTCq1gG0+6wlA65yFROQu4F6gJPBnn2N9W+MS3G3GmAj14LQfspYDeZdiAiuYbSp5DRWquTaoTlTV84AHgDGncqyIDBOR5SKyfO/evXkcYowxhXtm8TM8s/gZr8MIC8G8U0kAfOfxrA3sKKD8dODlUzlWVScBkwBatGiRK+kYY4q/9Axl656jQb3GkgSrWQ+UYCaVZUADEYkDtuM0vN/gW0BEGqjqBne1B5C5PBt4R0T+DZwDNAB+wBgTcbo/MTfb+nv3dQn4NWZeNzPg54xUQUsqqpomIsOB+UA0MFlV14jIWGC5qs4GhotIZyAVOAgMdo9dIyLv4TTqpwF3qWp6sGI1xoSmw4kp2dbrVoulYtnCpwc23glqj3pVnQvMzbHtEZ/lewo49gngieBFZ4wJdXOWbc22/uJt7YNynfHfjgdgVPtc3enMKbJhWowxIWvf0eRs66VigvPU16pdq4Jy3khkScUYE7Lqn+HkW/6a3m96kVwnEtjYX8aYkDVx3hqvQzCnyO5UjDEh6UTqyWdz/j2kbVAb6MctGgfAwx0fDto1IkWhSUVELsDpP3K2ql4sIk2A3qr6eNCjM8ZErN7jT45I3LhOlaBea/3+9UE9fyTx507lNeB+4FUAVf1ZRN4BLKkYY8LC1Gumeh1C2PCnTaWsqubseJgWjGCMMQZyz+5oig9/kso+ETkPd+wtEekH7AxqVMaYiPbfhSero56/pV3Qr/fIV4/wyFePFF7QFMqf6q+7cMbXulBEtgNbgBuDGpUxJmLtP5rM9O82Ac4jxRfWqhz0a247sq3wQsYv/iQVVdXOIlIOiFLVo+54XsYYE3A3TDg5WdZ17eoXyTXf7PNmkVwnEvhT/TUTQFWPq2rmUKHvBy8kY4xxtG5wttchmFOU752KiFyIM51vRRG5xmdXBaB0sAMzxpiypYqmK93oz0cD8M/O/yyS64Wzgt6xhkBPoBLQy2f7UWBoMIMyxkSmD5Zuzlp+5Nr4Irvu/qT9RXatcJdvUlHVWcAsEWmrqjaDjTEmqHYeTOTVBeuy1i+7sEaRXXtSr0lFdq1w58+95Up3HvnG+FR7qeotQYvKGBNxhrz4Vdby9e3P9zAScyb8SSpvA78CVwJjcR4nXlfgES4R6QY8jzNJ1+uqOj7H/nuB23A6U+4FblHV39196cAvbtE/VLW3P9f0lZqaSkJCAsnJyYUXNkWidOnS1K5dm5iYGK9DMSFscKcLivR6Iz8bCcAzXW2e+jPlT1I5X1WvFZE+qjrFHaJlfmEHiUg0MBHogjPn/DIRma2qa32KrQRaqGqiiPwFeAro7+5LUtVmp/Td5JCQkED58uWpV68eInImpzIBoKrs37+fhIQE4uLsqXRz0sFjJ7KW2zQ4q8h/X5NSk4r0euHMn6SS6n49JCIXA7uAen4c1wrYqKqbAURkOtAHZ4pgAFT1K5/yS4GBfpzXb8nJyZZQQoiIULVqVfbu3et1KCbEDHju86zlxwa0LPLrT+wxscivGa786acySUQqAw8Ds3GSwlN+HFcL8O2mmuBuy8+twKc+66VFZLmILBWR//PjenmyhBJa7P0wBWl5fnWvQzBnqNCkoqqvq+pBVV2kqvVV9SxVfcWPc+f11yPPUeJEZCDQAnjaZ3NdVW0B3ABMcMcfy3ncMDfxLLf/fvO2cOFCevbsmWv7qlWrmDt37mmd88knn8xa3rp1KxdffPFpx2eMrwevae7JdUfMG8GIeSM8uXa4KTSpiEglEblbRP4tIi9kvvw4dwJQx2e9NrAjj/N3Bh7CmaMlq2JVVXe4XzcDC4FLcx6rqpNUtYWqtqhevfj+h5OWVvSDPheUVAqLxzepGBNIRdXZ0QSPP+/gXJz2jl+AjFM49zKggTtO2HZgAM5dRxYRuRRnnpZuqrrHZ3tlIFFVT4hINeAy/KtyCznjxo1j2rRp1KlTh2rVqhEfH8/IkSPp1KkT7dq147vvvqN3797069ePW265hb1791K9enXefPNN6taty5AhQ+jZsyf9+vUDIDY2lmPHjrFw4UIeffRRqlWrxurVq4mPj2fq1KmICPPmzWPEiBFUq1aN5s1z/+eXkpLCI488QlJSEt9++y2jR49m3bp17Nixg61bt1KtWjW6du3K8uXLefHFFwHo2bMnI0eOZN68eSQlJdGsWTMaN27ME088QXp6OkOHDmXx4sXUqlWLWbNmUaZMmSL9OZvia89h7xvJJ3Sb4HUIYcOfpFJaVe891ROrapqIDMd5UiwamKyqa0RkLLBcVWfjVHfFAv9z69ozHx2+CHhVRDJw7qbG53hq7PR06pR723XXwZ13QmIidO+ee/+QIc5r3z5w/7BnWbiwwMstX76cmTNnsnLlStLS0mjevDnx8Sd7CR86dIhFixYB0KtXL2666SYGDx7M5MmTufvuu/noo48KPP/KlStZs2YN55xzDpdddhnfffcdLVq0YOjQoXz55Zecf/759O/fP9dxJUuWZOzYsdmSxqOPPsqKFSv49ttvKVOmDG+99Vae1xw/fjwvvvgiq1atApzqrw0bNvDuu+/y2muvcd111zFz5kwGDgzoMxcmjA164UuvQzAB5Fc/FREZCnwM+FZPHSjsQFWdi3On47vtEZ/lzvkctxi4xI/YQtq3335Lnz59sv5r79WrV7b9vn/wlyxZwgcffADAoEGD+Pvf/17o+Vu1akXt2rUBaNasGVu3biU2Npa4uDgaNGgAwMCBA5k0yb/ewr179z6tO4y4uDiaNXOe/o6Pj2fr1q2nfA5jujat7dm17/rkLsCeAgsEf5JKCs4dxUOcbGhXoGjGpA6kgu4sypYteH+1aoXemeSkWvDsdeXKlct3X+ZTUiVKlCAjIyPrfCkpKVllSpUqlbUcHR2d1RZyuk9Y+cbje12gwA6kOeNISvK+OsMUP3d0beTZtcvEWHVtoPjzSPG9OB0g66lqnPsqfgnFA+3bt2fOnDkkJydz7NgxPvnkk3zLtmvXjunTpwMwbdo02rdvD0C9evVYsWIFALNmzSI1NTXfcwBceOGFbNmyhU2bnEmO3n333TzLlS9fnqNHj+a5L/O6q1atIiMjg23btvHDDydnlI6JiSk0DmP84dueUq60d6MsPNP1GetNHyD+JJU1QGKwAwlHLVu2pHfv3jRt2pRrrrmGFi1aULFixTzLvvDCC7z55ps0adKEt99+m+effx6AoUOHsmjRIlq1asX3339f4N0NOMOgTJo0iR49etC+fXvOPffcPMv96U9/Yu3atTRr1owZM2bk2n/ZZZcRFxfHJZdcwsiRI7M1+A8bNowmTZpw4402Aag5M9aeEn6ksCoaEfkQZzDJr8jepnJ3cEM7NS1atNDly5dn27Zu3TouuugijyJyHDt2jNjYWBITE7n88suZNGlSnk9kRZJQeF+M96Z9vYH/LvoNgErlSjLj3i6exTJszjAgMkcrFpEVbp/AgPCnTeUj92VOw7Bhw1i7di3JyckMHjw44hOKMZkyEwpAG49neKxapqqn1w8nhSYVVZ1SFIGEq3feecfrEIwJede1yzVgRpGyGR8Dp6DphN9T1etE5BfyGF5FVZsENTJjTNh6ZPqyrOWPH7yKmGh/mndNcVDQnco97tfcA0cZY8wZ+H5D1gAaIZFQbp51MwBv9nnT40iKv3zfTVXd6S7eqaq/+76AO4smPGNMOGvX0Nu2lEx1KtShToU6hRc0hfKnob4L8ECObVflsc0YYwr14qers5bH9IsvoGTRGfunsV6HEDbyvVMRkb+47SkXisjPPq8twM9FF2Lx9vzzz3PxxRfTuHFjJkw4OWjdgQMH6NKlCw0aNKBLly4cPHgQgJkzZ9K4cWM6dOjA/v37Adi0aRMDBgwo0rgDMaR9bGxsgKIx4UJVmbP896z16CibXyfcFFSZ+Q7QC5jlfs18xauqjRboh9WrV/Paa6/xww8/8NNPP/Hxxx+zYcMGwBmY8YorrmDDhg1cccUVjB8/HoBnn32WpUuXctNNN2U9OTZmzBjGjRvn1zXT09OD880YEwBTv96Qtfzq7Zd7GEl2Az8YyMAP7M9aIBTUpnJYVbcCY4BdbltKHDBQRCoVUXzF2rp162jTpg1ly5alRIkSdOzYkQ8//BBwhlwZPHgwAIMHD84akTgqKooTJ06QmJhITEwM33zzDTVr1swaIDIvsbGxPPLII7Ru3ZolS5awYsUKOnbsSHx8PFdeeSU7dzrNY6+99hotW7akadOm9O3bl8REZ6CE3bt3c/XVV9O0aVOaNm3K4sWLAbKGtG/cuDFdu3bNGtNr06ZNdOvWjfj4eDp06MCvv/4KwJYtW2jbti0tW7bk4YcfDsJP1BR3vkml3lnlPYwku4ZVG9KwakOvwwgPqlrgC1iF0/ZyPrAJeA6YW9hxRf2Kj4/XnNauXZttveObHfXNlW+qqmpKWop2fLOjvv3T26qqejzluHZ8s6NO/2W6qqoeSjqkHd/sqDPXzlRV1b3H92rHNzvq7F9nq6rqzqM7c10vr+s3aNBA9+3bp8ePH9c2bdro8OHDVVW1YsWK2cpWqlRJVVU/++wzbd68ufbs2VMPHTqkXbt21QMHDhR4HUBnzJjhfF8pKdq2bVvds2ePqqpOnz5db775ZlVV3bdvX9YxDz30kL7wwguqqnrdddfpc889p6qqaWlpeujQId2yZYtGR0frypUrVVX12muv1bffdn5Wf/7zn/W3335TVdWlS5fqn/70J1VV7dWrl06ZMkVVVV988UUtV65cvj8XE5m6jv0462VCA85UJAH7W+xPQ32GOnOjXANMUNX/iMjKIOW4sHLRRRfxwAMP0KVLF2JjY2natCklShT8I+/SpQtdujjDVUyZMoXu3buzfv16nnnmGSpXrszzzz9P2bJlsx0THR1N3759AVi/fj2rV6/OOkd6ejo1a9YEnOq4MWPGcOjQIY4dO8aVV14JwJdffsl///vfrHNVrFiRgwcP5jmk/bFjx1i8eDHXXntt1vVPnHBG7/nuu++YOXMm4Azf/8AD9iyHOel48slBSF+49TIPIzHB5E9SSRWR64GbcNpUAPwaTlREugHP40zS9bqqjs+x/17gNiAN2Avcok41GyIyGKfqDeBxDUDP/oVDFmYtx0THZFsvG1M223rF0hWzrVcrWy3beo3YGn5d89Zbb+XWW28F4MEHH8ya/+Tss89m586d1KxZk507d3LWWWdlOy4xMZEpU6Ywf/58unbtyqxZs3jnnXeYNm0aQ4cOzVa2dOnSREdHA86dZ+PGjVmyZEmuWIYMGcJHH31E06ZNeeutt1hYyFD+eQ1pn5GRQaVKlbIm6crpdIfdN+EtPUO55unPstYbnhNaNegD3ncehJneb7rHkRR//vQ6uhloCzyhqlvc6YGnFnaQiEQDE3EeP24EXC8iOSdMWAm0UKd3/vu4UwaLSBXgH0BroBXwD3eK4WJnzx6nk9cff/zBBx98wPXXXw84E2JNmeLkySlTptCnT59sxz311FPcc889xMTEkJSUhIgQFRWV1Q6Sn4YNG7J3796spJKamsqaNWsAOHr0KDVr1iQ1NZVp06ZlHXPFFVfw8ssvA86dzZHymAbqAAAgAElEQVQjR/I9f4UKFYiLi+N///sf4CSxn376CXBGNvYdvt+YTF/+st3rEArUrEYzmtVo5nUYYaHQpKLONL4PAD+661ty3nHkoxWwUVU3q2oKMB3I9pdTVb9S1cy/kkuBzKnfrgQWqOoBVT0ILAC6+fMNhZq+ffvSqFEjevXqxcSJE6lc2cmNo0aNYsGCBTRo0IAFCxYwatSorGN27NjB8uXLsxLNfffdR5s2bZgyZQo33HBDgdcrWbIk77//Pg888ABNmzalWbNmWQ3v48aNo3Xr1nTp0oULL7ww65jnn3+er776iksuuYT4+PisJJSfadOm8cYbb9C0aVMaN27MrFmzss4zceJEWrZsyeHDh0/9h2XC1jOzf8panv9wDw8jyduo9qMY1X5U4QVNofwZ+r4X8AxQUlXjRKQZMFadueQLOq4f0E1Vb3PXBwGtVXV4PuVfxHnK7HERGQmUVtXH3X0PA0mqmu8sOqE69L3Jzd6XyJKhylWPO7OK/63nJXS7tK7HERlfgR763p/qr0dx7joOAajqKpxHiwuTV+V6nhlMRAYCLXCmLfb7WBEZJiLLRWT53r17/QjJGFOU1CehACGbUPq+15e+7/X1Ooyw4E9SSVPVnHUZBd/eOBIA38F0agM7chYSkc7AQ0BvVT1xKseq6iRVbaGqLapXr+5HSMaYouQ7cGSflvW8C6QQbWu3pW3ttl6HERb8efprtYjcAESLSAPgbmCxH8ctAxq4DfvbgQFAtgYBEbkUeBWnmmyPz675wJM+jfNdgdF+XNMYEwJS0tKZ+vUGZny3KWtbozqh+6zNyHYjvQ4hbPhzp/JXnOmET+AM3XIYGFHYQaqaBgzHSRDrgPdUdY2IjBWRzPaYp4FY4H8iskpEZrvHHgDG4SSmZThtOAdO6Ts7GcfpHGaCxN6PyDBpwbpsCQWgU+NzPIrGFKVCG+qLi7wa6rds2UL58uWpWrWq9Z8IAarK/v37OXr0KHFx/jTLmeLqynGfZFsPxSe+fPV+1/k/d/b1sz2OpOh5MUd9sVW7dm0SEhKwRvzQUbp06awOoMaEiivirvA6hLAR1kklJibG/iM2poilZ5ys/SgujxDf0+aewgsZv4R1UjHGFK1Nu45w52vfZK0Xh4RiAivfpCIi/6GAR4dV9e6gRGSMKbZ8E0pxctW0qwD49MZPPY6k+CvoTmV5AfuMMSabSQvWZlu/t1cTjyI5db0u6FV4IeOXfJNKIEYFNsZEjplLt2Rb73BRTY8iOXV3trzT6xDCRkHVXwU+W1fY2F/GmMix+9DJ0bOfu7kdjWqHbkdHE1wFVX+1BbYB7wLfk/d4XMYYw03/+SpruTgmlM7/7QzA5zd97nEkxV9BSaUG0AW4Hmd4lU+Ad1W14HHRjTER5c5JJxvnh3YunqNP92/c3+sQwkZBbSrpwDxgnoiUwkkuC0VkrKr+p6gCNMaErhOp6WzafXJSt75time/sKHxQwsvZPxSYD8VN5n0wEko9YAXgA+CH5YxJtRs3HmY175Yx9j+LZn+7UbWbj/Iqi37s5Wx4ZBMQQ31U4CLgU+Bx1R1dZFFZYwJKau27OOBqd8D0Hv8vDzL/GtQ66IMKaA6vdUJgIVDFnoaRzgo6E5lEHAcuAC42+c/EAFUVSsEOTZjTBGZt/IPypWKYfbyrdSqUo4RPZ0+Jgn7j1GrSrmshJKfUB8wsjBDmg3xOoSwEdajFBtjCnbDhM/Zf/REnvvu7n4xL8wtvIJi8p2dqFW1XKBDM0XERik2pphLSkmjTMmTv3rrdxwiLT2DxnWqAM4EVzsOJFLvrPKFnislLZ0oEdIzlFIx0Xy7bieVypXivilLAHjjzo58v2EP+44kc3vXRvR9ej6NalemcZ0qfLl6e74JBfAroQBhkVBS01MBiImO8TiS4i+odyoi0g14HogGXlfV8Tn2Xw5MAJoAA1T1fZ996cAv7uofhXW2tDsVEyiPv7+Cb9btYtaobpSOiT7j8x1PTuW+KUuYOLQ93/26mydm/kjdarGM6HkJjetUyZp7ZP7DPbINyDh+YGsujauWdZ75q7bx7zk/53udV4Z14I5JwRt7q3bVcgy/6mJGuVVhz9/SjuSUdJr5xFhcRXKbSqDvVIKWVEQkGvgNp69LAs4Mjter6lqfMvWACsBIYHaOpHJMVWP9vZ4lFZOfLbuPUKlcKSrHliIlLZ2PV/xBn5b1iI7K/qTSleM+YcBl5zHdZ8bC+Q/34HBiChkZSuXYUlnbU9LS+XHzPtpccDbJqem8Mn8Nn67cBsDbd/+Zj1f8ToeLarJozQ7+t2RzvrFN/1tnBjyXd4e7mOgoBlx2Ht+t383LwzrkmviqqH06pjtRInyzbidN61WlQpmSnsYTSFN/ngrAwCYDPY6k6BWnpNIWeFRVr3TXRwOo6j/zKPsW8LElFZOcms6Pm/bS7sIap3xseoZmJYqUtHRW/3GQZnFVuerxuQBc1+483lt8MmHUq16eCbe0o3RMNHe+9i2bffpb5CWzMfp4cirXPP0ZANe2rV9g0ihOujatzda9R/nPre35YOlmXl2wLmvf1a3juKNrIw+jM8FSnJJKP6Cbqt7mrg8CWqvq8DzKvkXupJIGrALSgPGq+lFB17OkEhp8q3JyenDa90RFCY9f34oM93MXlaNfw1/f+JbfdhzmpaHtOa9GRQC27TvGo+8tp/9l59G1aR1OpKZTskQUIkJKWjpPffQT36zbCcDgThew61Ai81clBPx7iz+vOpfUrcJ/F64nI8Sfb7mzW2Nemlfw4BdnVyrD7kNJTLrjcs6uVDbPqr7F63exdttBbiumPeX9lZjqjF1WNqasx5EUveLUUJ9XL6hT+VWsq6o7RKQ+8KWI/KKqm3wLiMgwYBhA3bo2GVBhMv+BEBEyVFm5ZR/N46ohIrz55a9M/25TtmSQ+R95THQUHz94FZt2HaF0TDRfr9vJvJV/0LtlPfq2qU9KWjq/7z3G+TVOPmV+6PgJKpUrxbHkVMqVKoGIsGLzPgBWbd3HA28X/Ijqna99y/yHe3DTf75k96EkAJ6d/TPPzj7ZplC5XCkOHs/e0Dxl4W9n9kMqwIpNe1mxKfhTUz9xQyseeueHfPdPvqsTa7cd5JnZP+W5P/M9zCupNK5Tmc5NahMTHUWXpoVP69yuYQ3aNTz1u8bipvu07kBktqkEWjCTSgJQx2e9NrDD34NVdYf7dbOILAQuBTblKDMJmATOncoZxhsWklLSWL/9EOfVqEjZUiU4dPwEN0z4IluZx69vyZh3l+V5/LOzf6JVg7OIr189q4onNT2DqYt+4+2vN2QrO2nBOib5VJEM6nhB1nL/f+c/MF9hCSXT3B//yEooecmZUAKpoJ/Rmbo0rhpP3tgqq1qu/YU1GNzpAo4kpdK4TuVCe6XXqlKOKrGlsiWVi2pXYl3CIc47+2Rin/vQVXR/4lOm3XMF1SqURlWtx3s+/tLiL16HEDaCWf1VAqeh/gpgO05D/Q15DUiZs/pLRCoDiap6QkSqAUuAPr6N/DlFcvXXd7/uYuz/VvD23X9m0Atfeh1OSHt2cNusx21v79qIVz/L+yM1/+EepGdk0P2JwmcCfPz6ltSoVJbbXl6Ute2eHpfQusFZHEtO5Zwq5Vi0ZgdPz3KSwNyHriI6KqrAc2ZWI74yrANHklJJ2H8s6xHfzDuREW9+R8/4c2l/YQ2io6O4ccIXvHJ7B6rEli40ZmMyFZs2FQAR6Y7zyHA0MFlVnxCRscByVZ0tIi2BD4HKQDKwS1Ubi0g74FUgA4gCJqjqGwVdK9KSitdPAhW1QDSI39urCVc2q5Ot3cd3+WhSKjO+25it/WDz7iNULV+aIS9+RZNzq3L5RTV4yk0Ot3e5iNnLf+et4X/KKp+Slo4qlMqjfeKdbzbQ4aKa1Knm9/MnWdIzMhjz7jJu6NCAS+pWOeXjTcEOJx8GoGLpih5HUvSKVVIpSpGSVJJS0vi/f80vsuvFli7BseS0XNtbnV+dHzaeevtCziewAGbc25kNOw/z3Mc/s//oCT564Mpc3+P8h3tkdfTr8WTuu4e8HgxQVfYcTuKn3/fTtWmdXPsBXp6/hpS0DO7pcYlf8aemZ9DzyU+Jr1+NJ28svmNdmeysn4ollVwiJan4c4fycL/mjHv/xzz3/W9kF5as301aegY94s8FYNALX7LncBJzRndDRJi5ZDNvfrUecP5YJ6WksftQEre/+jUAT9/UhgtrVaJkiWh+2LCHh6efbHsoVSKK9+/vSq9/5h50cN6Y7tnq9N/9diMNz6lE8/p5d54b8eZ3rEs4xJzR3ShZ4uR//oP/8yW7DiXxzxtbM3ra97wyrANxZ9tQdOb0fbDOGXz9mouu8TiSomdJJR+RkFS6PzGX9DyeZb2vdxPmrdzGmm0HAScRbN1zlPEfruTyRjW5tt15REcJJ1LTsw0PcqoWrdnBibT0XP/1Z6hmNTpn3jGkZ2SQlq7ElIjiWFIqx0+kUbNyYB7XTE5J41hyGtUqWNuBMWfKkko+IiGp+N6lNDynEut3HAJg2j1XUKV8KYa+tIi/9riYZvWK/7AZxhSlfYnO4+7Vykbe705x6qdiAuiQz+OzpWOieeHWyzhwLJl5K7dl/cf+xl2dPIrOmOKt33v9gMhsUwk0SyrFhG+/j1mjugFQJbY0N3Ro4FVIxoSN+9re53UIYcOSSjFwNCk1a/mx/gG7SzXGuHo17OV1CGGj4B5YJiT0e+azrOU2F5ztYSTGhKddx3ax69gur8MIC3anEuJ8H6Rofxoj9xpjCjfg/QGAtakEgiWVEOc7sODD18Z7GIkx4WtU+1FehxA2LKmEuNT0DACub3++x5EYE766nd/N6xDChrWphLiffz8AwI2X21NexgTLtsPb2HZ4m9dhhAW7UwlhhxNTspZjoi3/GxMsgz4cBFibSiBYUglhv7k95uuexqi2xhj/jbl8jNchhA1LKiHswDGnF/0/rrMGemOCqXP9zl6HEDasTiWE/bHvGDHRUdSsXM7rUIwJa5sPbmbzwTObr8c47E4lhP2wYQ81KpUhOsqmgDUmmG6ZdQtgbSqBENQ7FRHpJiLrRWSjiOR6EFxELheRH0UkTUT65dg3WEQ2uK/BwYwzVGVkaJ4zCBpjAuuxTo/xWKfHvA4jLATtTkVEooGJQBcgAVgmIrNzzDP/BzAEGJnj2CrAP4AWgAIr3GMPBiveUJOhSsKB47Q6v7rXoRgT9jrW6+h1CGEjmHcqrYCNqrpZVVOA6UAf3wKqulVVf8aZi97XlcACVT3gJpIFQET1TvrF7Z+ycdcRjyMxJvyt37ee9fvWex1GWAhmm0otwLc3UQLg76TeeR1bK2chERkGDAOoW7fu6UUZov7+9lIAhvypoceRGBP+bv/4dsDaVAIhmEklr9Zlf6eZ9OtYVZ0ETAJn5kf/Qwt959eowMZdR+jatLbXoRgT9p684kmvQwgbwUwqCYDvZOa1gR2ncGynHMcuDEhUxURyajqVypVExJ78MibY2tVp53UIYSOYbSrLgAYiEiciJYEBwGw/j50PdBWRyiJSGejqbosYCfuPc+h4SuEFjTFnbPWe1azes9rrMMJC0O5UVDVNRIbjJINoYLKqrhGRscByVZ0tIi2BD4HKQC8ReUxVG6vqAREZh5OYAMaq6oFgxRqq/nzxOV6HYExEGD53OGBtKoEQ1M6PqjoXmJtj2yM+y8twqrbyOnYyMDmY8YWqY8nO9MFh1UhkTAh7usvTXocQNqxHfQhav90ZSDK2dIzHkRgTGVrWaul1CGHDxv4KQZkDSV7dKs7jSIyJDKt2rWLVrlVehxEW7E4lBO06lIgA1SuW9joUYyLCiHkjAGtTCQRLKiHo858TUKBkCRv3y5iiMKHbBK9DCBuWVELQviPJXodgTERpVqOZ1yGEDWtTCUFpGUqbBmd5HYYxEWPZ9mUs276s8IKmUHanEmKSUtIAWL0t4rrlGOOZ+xfcD1ibSiBYUgkxuw8lATC080UeR2JM5Hix+4tehxA2LKmEmN2HEwE4t3p5jyMxJnJcfNbFXocQNqxNJcRk3qmcXamMx5EYEzkWb1vM4m2LvQ4jLNidSojZfTiJkiWiqFyulNehGBMxHvziQcDaVALBkkqI2X0okbMqlrEh740pQq/2fNXrEMKGJZUQs/tQEmdXKut1GMZElIbVbIbVQLE2lRCz+3ASZ1e09hRjitKirYtYtHWR12GEBbtTCSHJKWkcTkyxpGJMEfvHwn8A1qYSCEFNKiLSDXgeZ5Ku11V1fI79pYD/AvHAfqC/qm4VkXrAOmC9W3Spqt4RzFhDwfYDxwGoWdmqv4wpSpP7ROTUTUERtKQiItHARKALzpzzy0Rktqqu9Sl2K3BQVc8XkQHAv4D+7r5NqhpRA/L8se8YYH1UjClq9SvX9zqEsBHMNpVWwEZV3ayqKcB0oE+OMn2AKe7y+8AVEsGPPe11B5I8y6q/jClSn2/+nM83f+51GGEhmNVftYBtPusJQOv8yrhz2h8Gqrr74kRkJXAEGKOq3wQx1pDwxhe/AlC2lDV1GVOUHv/6cQA61+/scSTFXzD/euV1x5Fz2vX8yuwE6qrqfhGJBz4SkcaqeiTbwSLDgGEAdevWDUDIxphI9PbVb3sdQtgIZvVXAlDHZ702sCO/MiJSAqgIHFDVE6q6H0BVVwCbgAtyXkBVJ6lqC1VtUb169SB8C0VHVSkVE03XprW9DsWYiFOnYh3qVKxTeEFTqGAmlWVAAxGJE5GSwABgdo4ys4HB7nI/4EtVVRGp7jb0IyL1gQbA5iDG6rnDiSmcSE2n/tkVvA7FmIgzb+M85m2c53UYYSFo1V9uG8lwYD7OI8WTVXWNiIwFlqvqbOAN4G0R2QgcwEk8AJcDY0UkDUgH7lDVsJ5gJPNx4lpVynkciTGRZ/y3Tm+Hbud38ziS4i+oLcKqOheYm2PbIz7LycC1eRw3E5gZzNhCTcJ+N6lUtaRiTFGb3m+61yGEDXvMKEQk7D9OiSihhg15b0yRqxFbw+sQwoaN/RUitu8/Rs3KZYmOsrfEmKI2Z/0c5qyf43UYYcHuVELEtv3HqV011uswjIlIzy55FoBeDXt5HEnxZ0klBKRnKDsPJtK6wVleh2JMRHr/uve9DiFsWFIJAXsPJ5GanmGN9MZ4pFrZal6HEDasAj8EbNvvDCRp1V/GeOODdR/wwboPvA4jLNidSggY8+4yAM6xIe+N8cQL378AwDUXXeNxJMWfJZUQUiW2lNchGBORZg2Y5XUIYcOSisfSM06OsRnBo/4b46mKpSt6HULYsDYVj+06mAjAvb2aeByJMZFrxuoZzFg9w+swwoLdqXhsyx5nNP+4s2y2R2O88vLylwHof3H/QkqawlhS8dic5b8DULeaPflljFfm3ji38ELGL5ZUPLbDrf4qXdLeCmO8UjbGnrwMFGtT8VB6RgbHklLpEW+zVhrjpak/T2Xqz1O9DiMs2L/HHvpx8z4SU9Joem5Vr0MxJqK9/uPrAAxsMtDjSIq/oN6piEg3EVkvIhtFZFQe+0uJyAx3//ciUs9n32h3+3oRuTKYcXrl/SXOZJbN6xfvqZCNKe4WDFrAgkELvA4jLAQtqbjTAU8ErgIaAdeLSKMcxW4FDqrq+cBzwL/cYxvhzALZGOgGvJQ5vXC4OJKYwqqt+6leoTTly8R4HY4xES0mOoaYaPs9DIRgVn+1Ajaq6mYAEZkO9AHW+pTpAzzqLr8PvChOD8A+wHRVPQFscacbbgUsCWK8eVJVMjTzq7OckZG5rGRkKKqQoUp6hmaVz9yXoc7+lLR0Ek+kcfxEGjsOHOf1L34F4G89rX+KMV57a9VbAAxpNsTTOMJBMJNKLWCbz3oC0Dq/Mu6c9oeBqu72pTmOrVXQxTbuOkKf8fNOO9goEaKiICMDUtMz3GTgJIhguerSOsSfZ1VfxnjNkkrgBDOp5DXmSM4/0fmV8edYRGQYMMxdPTF79FWrTylCb1QD9gF8BvzN21jykxVjiCsOcRaHGKF4xBn0GOXmgAyVVNx+lucG8sTBTCoJQB2f9drAjnzKJIhICaAicMDPY1HVScAkABFZrqotAhZ9kBSHOItDjFA84iwOMULxiLM4xAjFI85gxhjMp7+WAQ1EJE5ESuI0vM/OUWY2MNhd7gd8qarqbh/gPh0WBzQAfghirMYYYwIgaHcqbhvJcGA+EA1MVtU1IjIWWK6qs4E3gLfdhvgDOIkHt9x7OI36acBdqpoerFiNMcYERlA7P6rqXGBujm2P+CwnA9fmc+wTwBOncLlJpxOjB4pDnMUhRigecRaHGKF4xFkcYoTiEWfQYhSntskYY4w5czb2lzHGmIAJi6RS2HAwQbjeZBHZIyKrfbZVEZEFIrLB/VrZ3S4i8oIb288i0tznmMFu+Q0iMthne7yI/OIe84KcxpSQIlJHRL4SkXUiskZE7gnROEuLyA8i8pMb52Pu9jh36J4N7lA+Jd3tpzy0T6A+HyISLSIrReTjEI5xq/uerBKR5e62UHvPK4nI+yLyq/v5bBuCMTZ0f4aZryMiMiIE4/yb+3uzWkTeFef3ydvPpbqd/IrrC+chgE1AfaAk8BPQKMjXvBxoDqz22fYUMMpdHgX8y13uDnyK0/emDfC9u70KsNn9Wtldruzu+wFo6x7zKXDVacRYE2juLpcHfsMZLifU4hQg1l2OAb53r/8eMMDd/grwF3f5TuAVd3kAMMNdbuS+96WAOPczER3IzwdwL/AO8LG7HooxbgWq5dgWau/5FOA2d7kkUCnUYswRbzSwC6c/R8jEidMhfAtQxufzOMTrz2XQ/vAW1ct9U+b7rI8GRhfBdeuRPamsB2q6yzWB9e7yq8D1OcsB1wOv+mx/1d1WE/jVZ3u2cmcQ7yygSyjHCZQFfsQZeWEfUCLne4zzNGFbd7mEW05yvu+Z5QL1+cDpK/UF8GfgY/eaIRWje+xWcieVkHnPgQo4fwglVGPMI+auwHehFicnRySp4n7OPgau9PpzGQ7VX3kNB1PgkC5Bcraq7gRwv57lbs8vvoK2J+Sx/bS5t7mX4twFhFyc4lQrrQL2AAtw/js6pKppeZw729A+gO/QPqcS/6maAPwdyHDXq4ZgjOCMPPGZiKwQZ8QJCK33vD6wF3hTnKrE10WkXIjFmNMA4F13OWTiVNXtwDPAH8BOnM/ZCjz+XIZDUvFrSBcPnepQNAH9fkQkFpgJjFDVIwUVPcV4AhanqqarajOcu4FWwEUFnLvI4xSRnsAeVV3huzmUYvRxmao2xxkd/C4RubyAsl7EWQKn6vhlVb0UOI5TjRRKMZ68uNMe0Rv4X2FFTzGeQHwuK+MMvhsHnAOUw3nf8ztvkcQYDknFryFdisBuEakJ4H7d427PL76CttfOY/spE5EYnIQyTVU/CNU4M6nqIWAhTp10JXGG7sl57qx4xL+hfQLx+bgM6C0iW4HpOFVgE0IsRgBUdYf7dQ/wIU6SDqX3PAFIUNXv3fX3cZJMKMXo6yrgR1Xd7a6HUpydgS2quldVU4EPgHZ4/bk8k7rGUHjh/OezGSdbZzYmNS6C69Yje5vK02RvwHvKXe5B9ga8H9ztVXDqliu7ry1AFXffMrdsZgNe99OIT4D/AhNybA+1OKsDldzlMsA3QE+c/wx9GxvvdJfvIntj43vucmOyNzZuxmloDOjnA+jEyYb6kIoR5z/V8j7Li3HmIwq19/wboKG7/KgbX0jF6BPrdODmUPz9wWl7XIPTFik4D0D81evPZVD/8BbVC+fJi99w6uIfKoLrvYtTh5mKk81vxamb/ALY4H7N/OAIzmRlm4BfgBY+57kF2Oi+fD+4LYDV7jEvkqNR088Y2+Pcqv4MrHJf3UMwzibASjfO1cAj7vb6OE/HbHR/SUq520u76xvd/fV9zvWQG8t6fJ6kCeTng+xJJaRidOP5yX2tyTxPCL7nzYDl7nv+Ec4f25CK0T1PWWA/UNFnW0jFCTwG/Oqe522cxODp59J61BtjjAmYcGhTMcYYEyIsqRhjjAkYSyrGGGMCxpKKMcaYgLGkYowxJmAsqZiIJCILRSTo84iLyN3uSLzTcmxvISIvuMudRKRdAK9ZT0RuyOtaxgRbUGd+NCYciUgJPTm2UmHuxHnuf4vvRlVdjtNXA5y+L8dwOisGIoZ6wA04IyrnvJYxQWV3KiZkuf9xrxOR19w5Iz4TkTLuvqw7DRGp5g6hgogMEZGPRGSOiGwRkeEicq87eOFSEanic4mBIrLYnYuilXt8OXHmy1nmHtPH57z/E5E5wGd5xHqve57VIjLC3fYKTke02SLytxzlO4nIx+5gn3cAfxNn3o4OIlJdRGa6MSwTkcvcYx4VkUki8hnwX/fn842I/Oi+Mu92xgMd3PP9LfNa7jmquD+fn92fRxOfc092f66bReRun5/HJ+LMd7NaRPqf2btqwt6Z9C62l72C+cL5jzsNaOauvwcMdJcX4vZaBqoBW93lITg9hsvjDAFzGLjD3fcczsCamce/5i5fjjvkDvCkzzUq4fQmLueeNwG3B3WOOONxelGXA2JxerNf6u7bSo6h6N3tnTjZM/9RYKTPvneA9u5yXWCdT7kVnJw/oyxQ2l1uACzPee48rvUf4B/u8p+BVT7nXozTI7saTk/yGKBv5s/JLVcx5/diL3v5vqz6y4S6Laq6yl1egZNoCvOVqh4FjorIYWCOu/0XnGFhMr0LoKpfi0gFEamEM3dGbxEZ6ZYpjfOHHWCBqh7I43rtgQ9V9TiAiHwAdMAZfuZ0dAYaycmJACuISHl3ebaqJrnLMcCLItIMSAcu8OPc7XESBar6pYhUFZGK7r5PVPUEcEJE9gBn4/zMnhGRf+Ekpm9O83syEcKSigl1J3yW03EGnQTnDiaz+rZ0Acdk+KxnkP0zn3OMoszhvvuq6i9zEgUAAAFnSURBVHrfHSLSGmeY9ryc8jSwhYjCmUwpyXejm2R8Y/gbsBto6h6T7Me5CxrOPOfPuoSq/iYi8ThjQP1TRD5T1bF+fRcmIlmbiimutuJUOwH0O81z9AcQkfbAYVU9jDPr3V/F/QsuIpf6cZ6vgf8TkbLiTDh1Nc5IvP46ilNdl+kzYHjminsnkpeKwE5VzQAG4Ywsm9f5csZ6o3veTsA+LWCeHRE5B0hU1ak4E0I1z6+sMWBJxRRfzwB/EZHFOG0Ap+Oge/wrOCNNA4zDqVb6WURWu+sFUtUfgbdwRn79HnhdVU+l6msOcHVmQz1wN9DCbUxfi9OQn5eXgMEishSn6ivzLuZnIM1tXP9bjmMezTw3ToP+4EJiuwT4QZyZOR8CHj+F78tEIBul2BhjTMDYnYoxxpiAsaRijDEmYCypGGOMCRhLKsYYYwLGkooxxpiAsaRijDEmYCypGGOMCRhLKsYYYwLm/wEosyCN3K0mkQAAAABJRU5ErkJggg==\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
}