{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<form action=\"https://github.com/prmiles/pymcmcstat_examples\">\n",
    "    <input type=\"submit\" value=\"Return to Index\" style=\"background-color: green; color: white; width: 150px; height: 35px; float: right\"/>\n",
    "</form>\n",
    "\n",
    "# Save Chain to Log File\n",
    "Author(s): Paul Miles | Date Created: June 18, 2018\n",
    "\n",
    "Many models are time consuming to evaluate.  As MCMC simulations required many model evaluations, it can be useful to periodically save the chain elements to a file.  This can be useful for a variety of reasons:\n",
    "\n",
    "- Chain visualization while simulation continues to run.\n",
    "- Chain is saved in the event that simulation ends prematurely.  \n",
    "\n",
    "This is important when working on remote systems where you may have limited computation time.  This tutorial demonstrates the following:\n",
    "\n",
    "- How to specify a log file directory\n",
    "- Format to save log files in (binary or text)\n",
    "- How to read in log files for analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run Simulation & Export to Log Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import required paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pymcmcstat.MCMC import MCMC\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a simple model and sum-of-squares function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define test model function\n",
    "def test_modelfun(xdata, theta):\n",
    "    m = theta[0]\n",
    "    b = theta[1]\n",
    "    nrow, ncol = xdata.shape\n",
    "    y = np.zeros([nrow,1])\n",
    "    y[:,0] = m*xdata.reshape(nrow,) + b\n",
    "    return y\n",
    "\n",
    "def test_ssfun(theta, data):\n",
    "    xdata = data.xdata[0]\n",
    "    ydata = data.ydata[0]\n",
    "    # eval model\n",
    "    ymodel = test_modelfun(xdata, theta)\n",
    "    # calc sos\n",
    "    ss = sum((ymodel[:, 0] - ydata[:, 0])**2)\n",
    "    return ss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize MCMC object:\n",
    "- Add data\n",
    "- Define model settings\n",
    "- Define model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize MCMC object\n",
    "mcset = MCMC()\n",
    "# Add data\n",
    "nds = 100\n",
    "x = np.linspace(2, 3, num=nds)\n",
    "y = 2.*x + 3. + 0.1*np.random.standard_normal(x.shape)\n",
    "mcset.data.add_data_set(x, y)\n",
    "# update model settings\n",
    "mcset.model_settings.define_model_settings(sos_function=test_ssfun)\n",
    "\n",
    "mcset.parameters.add_model_parameter(\n",
    "    name='m',\n",
    "    theta0=2.,\n",
    "    minimum=-10,\n",
    "    maximum=np.inf,\n",
    "    sample=True)\n",
    "mcset.parameters.add_model_parameter(\n",
    "    name='b',\n",
    "    theta0=-5.,\n",
    "    minimum=-10,\n",
    "    maximum=100,\n",
    "    sample=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define log file directory and turn on flags in simulations options\n",
    "The following keyword arguments of the simulation options allow you to setup the log files.\n",
    "\n",
    "- `savedir`: Directory in which to store log files.  If not specified, but log files turned on, then saves to directory with naming convention 'YYYYMMDD_hhmmss_chain_log'.\n",
    "- `save_to_bin`: Save log files in binary format.  Uses `h5py` package for binary read/write.\n",
    "- `save_to_txt`: Save log files in text format.  Uses `numpy` package for text read/write.\n",
    "\n",
    "By default the feature is set to `False`.  You can save to either format or to both.  Regardless of what format is used to save the chain, a text log file will be included which appends a date/time stamp with corresponding chain indices.  This will be explained in more detail later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "datestr = datetime.now().strftime('%Y%m%d_%H%M%S')\n",
    "savedir = str('{}_{}'.format(datestr, 'serial_chain'))\n",
    "mcset.simulation_options.define_simulation_options(\n",
    "    nsimu=int(5e4), updatesigma=1, method='dram',\n",
    "    savedir=savedir, savesize=1000, save_to_json=True,\n",
    "    verbosity=0, waitbar=False, save_to_txt=True, save_to_bin=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcset.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process Simulation\n",
    "At this point, the simulation is either running, or has completed running.  You will observe a folder in the working directory that matches the input argument for `savedir`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m20190517_073038_serial_chain\u001b[m\u001b[m/ \u001b[31mChain_Log_Files.ipynb\u001b[m\u001b[m*\r\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the folder `20190517_073038_serial_chain` matches what was specified for `savedir` and display its contents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "savedir=20190517_073038_serial_chain\n"
     ]
    }
   ],
   "source": [
    "print('savedir={}'.format(savedir))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20190517_073038_mcmc_simulation.json  s2chainfile.h5\r\n",
      "binlogfile.txt                        s2chainfile.txt\r\n",
      "chainfile.h5                          sschainfile.h5\r\n",
      "chainfile.txt                         sschainfile.txt\r\n",
      "covchainfile.h5                       txtlogfile.txt\r\n",
      "covchainfile.txt\r\n"
     ]
    }
   ],
   "source": [
    "ls 20190517_073038_serial_chain/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, there are log files saved in both binary (h5) and text (txt) format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing the log files\n",
    "We start by importing several modules from the [pymcmcstat](https://prmiles.wordpress.ncsu.edu/codes/python-packages/pymcmcstat/) package.  We note that this operation should be done from a separate script file, and possibly from a separate computer.  For example, if running a long simulation on a remote server, you can periodically copy the log files from the remote server and analyze the chains on your local machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymcmcstat.chain import ChainProcessing as CP\n",
    "from pymcmcstat.chain.ChainStatistics import chainstats\n",
    "from pymcmcstat.plotting.MCMCPlotting import Plot\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the plotting class and define the directory in which to find the log files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize classes\n",
    "MCP = Plot()\n",
    "# define directory where log files are saved\n",
    "savedir = '20190517_073038_serial_chain'\n",
    "# For testing purposes we can repeatedly read in the data to see how binary versus text is processed.\n",
    "ns = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in binary data files and print amount of time it takes to process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary: 0.06698329448699951 sec\n",
      "\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "for ii in range(ns):\n",
    "    results = CP.read_in_savedir_files(savedir, extension='h5')\n",
    "end = time.time()\n",
    "binary_time = end - start\n",
    "print('Binary: {} sec\\n'.format(binary_time/ns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text: 0.5604548215866089 sec\n",
      "\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "for ii in range(ns):\n",
    "    results = CP.read_in_savedir_files(savedir, extension='txt')   \n",
    "end = time.time()\n",
    "text_time = end - start\n",
    "print('Text: {} sec\\n'.format(text_time/ns))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is clearly seen that the binary files are more quickly processed.  In either case, the results extracted from the log files are identical, and we can proceed with the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "We extract the following from the results dictionary:\n",
    "- `chain`: Sampling chain for model parameters\n",
    "- `s2chain`: Observation error chain\n",
    "- `sschain`: Sum-of-squares error corresponding to each row of `chain`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain = results['chain']\n",
    "s2chain = results['s2chain']\n",
    "sschain = results['sschain']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the burn-in period for the chain as half the simulation run time.  Display statistics for burned-in portion of chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "---------------------\n",
      "name      :       mean        std     MC_err        tau     geweke\n",
      "$p_{0}$   :     2.1135     0.0391     0.0017    38.1465     0.9989\n",
      "$p_{1}$   :     2.7027     0.0984     0.0044    37.5644     0.9975\n",
      "---------------------\n"
     ]
    }
   ],
   "source": [
    "# define burnin\n",
    "nsimu = chain.shape[0]\n",
    "burnin = int(nsimu/2)\n",
    "# display chain statistics\n",
    "stats = chainstats(chain[burnin:,:], returnstats=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Chain\n",
    "- Chain panel\n",
    "- pairwise correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAEJCAYAAACgxlOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXuUX0WV7z/714g4QkdZalA6IJiICD5mQAF1AYbBCIOjjomwhCGZGyXRmXEEWxByBUYHHBMBrxfTgsnIZRRMs3ThCj6YYaV9DUj7wiXjoGE02o2IT7pBBSSp+0edM119+jyqzuN3fr/u/V2rVvfvnDrn1PNbu3bt2iXGGBQKhaIpdNpOgEKhmN9QklEoFI1CSUahUDQKJRmFQtEolGQUCkWjUJJRKBSNQklGoVA0CiUZhULRKJRkFApFo1CSUSgUjUJJRqFQNIq92k5AGxARAZ4FPNR2WhSKPsN+wM9MwKbHBUkyWIKZbDsRCkWfYgi4zzfyQiWZhwAmJiYYHBxsOy2Kirj+enj728EYEIEPfxjOPrvtVM0/TE9Ps2TJEgicAchCdPUgIoPA1NTUlJJMn2NyEg4+GPbsmbk2MAC7dsHQUGvJmpeYnp5m0aJFAIuMMdO+z6niV9HX2LlzNsEA7N4N997bTnoUc6Eko+hrLFsGnUQrHhiApUvbSY9iLpRkFH2NoSG49lpLLGD/XnONTpV6CaqTUZ3MvMDkpJ0iLV2qBNMUyupkFurqkmKeYWhIyaVXodMlRWOYnISxMftXsXChJKNoBFu32qXl5cvt361b206Roi2oTkZ1MrVDbVfmYnLSLrcvW9a/ZaB2Mg2jiug/n6YNPnlR25XZWPBSnTFmwQVgEDBTU1PGB1u2GNPpGAP275YtXo9VfrbX4JuXiYmZeHEYGLDXFxrmU1lMTU0ZwACDJqS/hUSeLyGEZKo0kvnUwELzsmWLvR/H62dyrYIdO2aXWRzGxtpOWTjKkkzr0yURuVBEviEiD4nIL0TkZhE5rOCZt4jIV0Xkt1G4TURe2kT6qoj+82naEJqXFSvghhtgdNTqYtaubTyJwejGNFYtkntDJ3MC8BHgWOBk4AnAv4nIk3OeORG4EXglcBwwET1zYN2Jq9JI5lMDC8lLrIM4/XQ44wy49dbupDEE3dKTqEUyvTddAp6OFcmOD3hmAJgGzvaMH6yTSRP9JyasOJw3/ZlP0wafvPTDFLGNNE5M2ClSt8vBp436Yt7oZIClUUaODHhmP+APwGkZ958YEUscDgwhGWPmNpIQhW63G1idDSvt3Xl56aYOomw+m0xjk2UfiroXHeYFyWCnb7cAXwt8bjPw38A+GfcvjQpnVgghGRe9PFq3vZrVrbKpks+m0th22btoIo/zhWRGgF3AUMAz7wZ+A7wwJ05lScZFr64Y1NGw6hiJm54i1pHPutPYJLmWqZMm2mjfkwxwNVaBe0jAM8PAg8DRgd8K0skkMT5eX4OqU7yu2rDqHImbnCLW1YFC0lhUT00NPGXrRCWZ2R1eIoK5D1gW8Nz5wBRwbIlvliYZt9LdyivTIUOM23yIyKdhZb2rl6eASdRJ8j7wqacmyq/qO+uW1vqZZDZH0sgJwAFOeJIT53rg/c7vC4BHgTckntnX85ulSCat0sGYDRuCXpP5rk7HmG3bZjei0JFs06aZ+MmGlfeuXp0CJlEnyfsgpKNX7dTJAaCOOqlTouxnkjEZYY0T50vAdc7vXRnPXOr5zSCL323bZkJapXc69a1wuARQxso2ji9izMaNMw03bfR3Sa2XJZmiPIyP1/P+tLyGdvSsTp32Dfda2gDQa3XStyTTRvAlmS1bbGfNIoOq+oA0qchtTFmkNjqa3mCT7xOZTTpFpFZmJM7roC45lNU7JYnTt+zdwaHIjilPUqy6rWTHjtnSpVvWbr6SeRsYsM+516voyerQ/SnJ1EwyExP+BJM2zfGB26mzyCSPONxGlycZheRhfNxfvM7roGnTmjqWmrOkyKRU4NZdLNX5EHMagZQh37T8u2kuylccLy2voUgjuTJQkgnJtAfJ+HZadxTyrUB3hMsisrixuw2800kf8bKmO2XJxjcPWR00Ly0h4r5PHYjM7kRFA4MPMWdJRiErUVXrIlRqy5JSNm6sVgculGRCMu1BMmmV49MYiiowb4RzO8Lw8Mw0I5YuNm/Ob3jJEctXEivTCLM66OioMVdckf/+uLO6HSNLZ1FmxPfNX9b7163Lz787FUtOBScmivOfJsmI5Eu1cdi0aW57iutZZPYAkSeNl1HoK8mEZLqAZKo27jwdgc97kw1DxJgTTsiOPz6ervTNm44VNeq8RpiVF3cql1dmSTJcvTp/2pUnyZWVGGK91oYN2elMk+iy9HSdjjErVxanR8QOIJs2zZ2CpSm2074TE9r4eHpailanyk67lGRCMl1AMlmVs3797IaxcWN6o0+b09cxlcnrMHlTl7Gx9EY9MZH/bBqKFJY+oegZ9/vJkdrtyHEdlClbn+lVmo1RWekw+d1Ox6Z9bGxGGspS9CfDunW2PrPur18/k960sklKQ75QkgnJdAlJJtlp4//TpI4yIn+VUDSNitOxebMxF188d8nXV7HZjby4ac/7nrtEv27d7M5bdzpiVFWuJ0O8glSVtNPKZts2Y7ZvN+aNb5xLbGWhJBOSaQ+dTFJMHx6eO7r7KA6L7GHOPLN6wyqSRtJWW2JJxtUtuIpNV0cS/+870tYRzjyz+HtpnXLNmno6alw/bp376ulCCSH520c3ExpWrKhmSmCMkkxQ8LWTmZiw5OKONC7ZZJm3u5WZJWL76C98QtoqVNLnTdb304jHmHqmREUdyyXwEBLpdnBJt1tS3OWXN5PvqsvYSjIhmfa0k8lS1CaVlW6HT1NiDg/X21jiBpic2iSncjt22KmE73tjgqy7M7mE6uqDYn1EXodqamT3DcPDtmzrnir51G+TocwydtdJBngqcDpwXhTOAJ5aNyE0EYpIpoyittOxc+AsyaZMxz33XGNGRtKNskZHZ4+ymzcb85732G+58/zQcOWV9Tfos86yhD0ykm60WDQNGR2dUV53S5pww/h4dyWZboXQZeyukgywFrgHuBLrbmEYuAr4L2BtkwRRR8gjmSqN6eKL06+Pjpabz8e6oLwGkrfKUOZ7mzc3I8kkv+Nrz5HUizQhafmEiy4qr5NZuXL2tLsXQs9LMsAPgCenXN8X+GHdpFB3yCOZKmJxlg6h0zHmnHPqbyBVFJHDw3N1Mu5ScdMdwm3keWXu6opcS2lXp7NunSXyVauaTXOerVJe6HRmpKGzzipuQ02k3Z12lt213m2SuQc4IOX6M4EfdJs0SqS/siRzzDF+ClWfhhtXvu9zmzZVs9lwl+NHR9OnZCFh7dpyz8XSWFFeBgas0ZyrjL7oorlm/hMT9azW5YWy5RTr8bLSdu656fuuwJLoyEg9q5BJnV3y/14imdMiaebTwIej8JnoWqoz714KTehkqjTazZttA/CdCoyNhUlcy5dnK4uNqa7UTCPWIrIVmbHX2bKl3Hdjso0td4u+V+ZeMrzmNfW3AVeqy7LRGh8PU+LnfSNr9bBo1anbJHNgFI7DOo56Q/T/QNsE4pn+wtWlbtqEuBXsLkVnddS4c/kSYXK7QdyY49Fr+/bu5tUNriFdaPCd1l1ySbbB4llnzbWJWrMm+11Z76my3O46PcvbClC2bGLL6G3bjLnssmK3H1kSTVdIBng58GNgdxQeADaGfrTt4LuEXVWaKbNxb3zcNob164t3aK9eHZ6mNCvTKnnshzAyUmyd7VMmK1Zkt4tVq8qbKsQDTMj+trx3xSuP7paSkHoeHW2XZO4GPg/8GfZ8pLOBbwH3Age2TR51kkzc8FxlaIjiL54GZa045T3n23HKEsRCIJZkfrM2JRrjN011CSlN4V51ZS7pYMydwuQ9t3z5bLuspGV6kR1SWmibZP4APDdxTYCbgJvaJo+AfHhtK0gb+Xw6b3KlJu09VZStvWAJ22SoonvIC+6mxKwtF3lpii25s6ZMcd3UvY+qaHqTjCdi07tuXXg7Se69a4Nkvg28IuX64cBDbZNHXSRTx27brMYQz483bzbmtNOa+U63CChePq5bST462lwe3JW1qlOTvI4aK/LrWIHshklBHGIL5zpJJnGEeiGuA/6viCxJXF+EPYt6XmDnTlvkdWLPHli3Dt79bjj/fHjb2+CWW+r9TqcD551Xf9rTcPnlcOON8KIX2bxlQcSGEPzqV3DRRdXSl4Xdu+Gyy+Dgg+H00/PTXhbGwPg4vOQl8LGPzc5/p2OvfexjMDDg/75PfQpOOKH+tCaxeHH97xQT0CJFJK6Sx7BL1ndhD7s/C7jcGPPJ2lPYAERkEJiamppicHBwzv3JSTjooLmdVcQ2kt27bQN5//thr71sx24TIrbjPOc5cMghcOyxzXQeF51O8TfizlWG9Dod+Ku/gk9/ujukWTcGBmDXLhgasu3pjjvg17+G/feHl71s5vpNNxW3n4EB+/wxxzRfFm66k5ienmbRokUAi4wx/kJF4DRjMbACe7DaDcD3gT9iSee7wCeie69ue0pUZbpkTLZOJjmnb8vM3Q2uI6eszZv9GAYG8vUfybq56KLe0lW5e4OyjqvZtq04zZs2dXeDZtaeptZ2YQP7AC8B3oI9CfJrwINtE0lVkokJJG3fjWtE1q3OXOQ3Ntk5P/7xZtLRbfLy1c/ECssyWy1iRbxrVexbzlmE7rr8yBqIfAmxyIFXaP19/OPZ78rb06SuHkIy7UkyWaNHPBJ1s8ONjYXZxdS1mzo234+dWvlKbmUU0GmEvn373H1WWWFkxNZbmd3arq3K2Fj6jno3rFs3+2C84eHZy+N1SpOxY3MfQ82qbSPPNaeSTEimA5xW1UkkRasERcZ3sUuHs88u/lboqJzX+ZLuGfLIzrXVKDpXKi3/eY7P163LtztyN1OWmcYmR3Gf5W3XFN/1kVP3AJQkwTTS8X3Pq1+dfT/P/YOSTEimA46pzWtoISOJuxM3y11m2i7i2GgsdGrW6cw2PquzkSfTIWLN1V0fNzHSHJnnhTz7odjgrWgzZbxlokw+3U7ma8iWJKem9CdpU5lYivKdUhbFyTvyV0kmJNMBJGNMtvLRPb3AZ3oQN+AkccVGYmnOqao4THrjG2emOXU4fBKxjrGy7hdtsHNH+jRdSyw1+aQjr6yLdBh510ON9JJ1a4yfXqis7UtS0vBpGxs22PIu834XSjIhmQ4kmTTjvLhBjo/bOW7RKkgcP0sSyOq8oTuu8zp+UnfQRAg5Jzpt9c5Xeuh00n3hut/PMv/P0rfE+ohQUvdR8rrfT+4tSkqwRxyRrUz2dWYfl2e8IdanDRXVnZJMSKYDScaYuTt148PTfBtjrC8IIYy6jqBN7qxNShRppFfFmbXPUap5Jz34rhCNjeUf55L3jTRJpej4Wtdc3/2ma0JQNJXLkvTGx4256qqZ6crEhJWUXUvxtGezBq14U6gbr6itFTmyUpIJyXQJkjGm+hk5If5+k5We7EyrV889aC7vvXnm4mkK2k7HmFe9Kj+Nl1ySrltKjoYuGefpdtxni8o66YMl7ZzqvG+E3ut05jrOGh7OPoEyLWRtPMyqkzQXHUX1l+ZcPib3LL89acf9pEFJJiTTJUimrpWm5OibFq66Kr3SY6Wxu3Sa3OxXxv7BmHTy63SK31d0MFxeZ/Z5Npa4kkvEvu4jkythq1fb60VnZiXTlqUz8x1oQnzqFhFw1jNJos0i95GRGSf1IX5+lWRCMl2CZOpYMUiOvlWOiM1Ssk5MWIVvXgcKyV+aDifvKBbf97quN7OeTeY5lh7q6Kw+HdlNW0j9J6WeUJ+6PocGVsl7WSjJhGS6Jkkmb8oUm/cnR8OkL9UqR8RmNZqyI2Feh3QN8kIaatbUI2+ptEo+XIRKK0UrY76SnqtADi2vOvLtk/cyUJJpmGSMST+69rLL0isztj6NV5/ckS0phfg0xtBGE9KBqjzjgyydT9H7q3aUUGml6F2uWwvXfim5IOCTriLH3VXrogq5Z0FJpgskY8zMMnCeUjJezamiI0n7bmijKTOSlh19i5B1pG9Z6coXdRBn2kqUSyYhZeYz5S3z3qxvlSH3LCjJdJFk0qZNyRHOR1EcKrrW3Wi6CR+pJG2Er4Mk8lafiiSKNHueMmQXv6tuPUkRypB7FpRkukQyWZ3FPdMmL17Viq6z0dQBn44ax8tLd94I34R05StR5NWjq7z2KYMm9CRFqPObSjJdIhnf0ShPkskzrCpqrG001CyEiP5x/DSppNsjfKgSPU+SCZ3+dHuAqPObSjJdIhlj/EX4NFuLrBHZt7G20VDrTEeaVJJFnFde2Uy+yijRk3vNfAwKs95VtOJYN+pS6CvJdJFkjAlblSiKF9pYm1oFCkGdElWe1NeEzqnsEv/o6Oyd5mXLIG4TrlVz07q1tHboO82LoSTTZZKpE2Uaa1OrQL6oW6LKs4JuQlKrS6FctgzalkhDp7rGKMn0Ncm03eDKom6JKm8HdhM6pzqIumwZtKlbK9veypLMXihax9AQXHutPTIlPgnhmmvSPcb3EtauhRUr4N57YenS6ul9+GHb5JPodOz768bQUPU0ly2DZcvmnvgwMNBMPpPYuXPuSRO7d9s8NNLmQhipiQBcCHwDeAj4BXAzcJjHc6uAe4BHgO8BpwZ8s6ckmRhVR9bQOXavIUs3k+d3tp/Rlm6t25JML5DMF4E1wBHAi4DPAT8BnpzzzMuAx4F3YU+vfB/2WJYjPb/ZkyRTBWXm2L2IMmb6/Yy2dGtlCK4syQQd7tYNiMjTsRLNCcaYr2TE2YYlodOca18H7jLGrPf4Ru7hbv2GyUl7ImJS9M46pKvXMTlZ3xRMkY3Qci57uFsv6mQWRX9/kxPnOODKxLVbgdelRRaRJwJPdC7tVzp1PYiuz7EbRh26EkUxulXOoWdhNwoR6QAfAv7DGHN3TtQDgAcS1x6IrqfhQmDKCZMVk9ozmJyEX/7SKhFddEuJqCiHyUkYG7N/5zt6imSAjwBHAmfU/N73YyWkOMyLcXLr1pmD442ZIZp+WZ1aqIjrbfly+3fr1rZT1Cx6ZrokIlcDpwHHG2OK+P3n2HO5XSyOrs+BMeZR4FHnWxVS2huYnIRzzpmZJhljD7gfHYXjjlOC6VUk623PHmu6sGLF/K2z1iUZsbgaeD2w3BjzY4/H7gBOSlw7Obq+IJCmh9mzB57+9PnbWOcD8vRn8xW9IMl8BHgT8FrgIRGJ9SpTxpg/AIjI9cB9xpgLo3v/B/iyiLwTu+R9BnA0cE5XU94i2jTmUpTHQqy31iUZ4K1YPcmXgPudcLoT5yDgmfEPY8ztWGI6B/gusBJ4XYGyeF4hthIeGLC/VQ/TH1iI9dZzdjLdQGwnMzEx0fd2MvfdBz/6ERx6KBx4YNupUfiiH+ttenqaJUuWQKCdzEIlmQOZR8vYCkWXMWSMuc838kIlGQGehd0vlYX9sEQ0VBBPoWUVgn4vq/2An5kA4ugFxW/XERVQLhM7y9wPhYiGCxFaVv6YB2UVnOZeUPwqFIp5DCUZhULRKJRksvEo8I84lsKKTGhZ+WPBldWCVPwqFIruQSUZhULRKJRkFApFo1CSUSgUjUJJRqFQNAolmQyIyN+KyC4ReURE7hSRl7adpjohIseLyHYR+ZmIGBF5XeK+iMh7ReR+EfmDiNwmIssScfYXkU+KyLSIPCgiW0Vk30ScF4rIV6NynBCR81PSskpE7onifE9ETm0m1+EQkQtF5Bsi8pCI/EJEbhaRwxJx9hGRj4jIr0XkYRH5tIgsTsQ5SEQ+JyK/j96zSUT2SsQ5UUS+LSKPisi9IrImJT391y5DvI4vlIDdAf4o8DfA84Frgd8Cz2g7bTXm8RTgn7B+fAx2F7t7/wLgQawLjhcCnwV+BOzjxPkCcBdwDPAKYCdwg3N/EOtI7BPY0yjOAH4PnOPEqXTyRBfKqfA0DWAE+CmwHDgK69foP5z7A9hje/4deHFU9r8ELnfiHAL8DrgiKoe/i8plRb+3y9YT0IsBuBO42vndwW5DeHfbaWsov7NIBhCsu41h59oi7BlXZ0S/D4+eO9qJ82pgD/Cs6PdbsQ7h93bi/DNwj/N7G3BLIj1fBz7adrlklNXTo3wf75TLY8BKJ87zojjHRr9PAXYDi50467H+pveOfn8AuDvxrU8BX+z3dqnTpQREZG/saHRbfM0Ysyf6fVxb6eoyDsE6ZXfLYArbyOMyOA540BjzTee527Akc4wT5yvGmMecOLcCh4nIU504tzEbt9K7ZZ08TeMo4AnMLqt7sJKNW1bfM8a4zu9vxUp6RzhxMsuhn9ulksxcPA0r3oachjDfEOczrwwOwJ6P9T8wxjyO7XxunLR34BGn58o64zSNA4DHjDEPJqIny6psOQyKyJPo43a5IHdhKxQlEZ+m8Yq2E9JPUElmLn5FNH9OXM88DWEeIs5nXhn8HHiGezNaLdk/ESftHXjE6amydk7TeKWZfZrGz4G9ReQpiUeSZVW2HKaN9XXdt+1SSSaBSH/wLZzTECIx+SQWzmkIP8Y2XLcMBrG6lrgM7gCeIiJHOc8tx7apO504x4vIE5w4JwM/MMb81onTsydPeJym8S3gj8wuq8OwfqndsnqBiLikfDLWN8v3nTiZ5dDX7bJtzXMvBuxS4SPAauwqyjXYpcLFbaetxjzui11OfTF2JeTc6P+DovsXRHn+S+AFwM2kL2F/G3gp8HLgh8xewl6EJavrsQrO07HLtMkl7D8C78SuylxKby1hb8Yu5Z+A1X3E4UlOnBHssvYrscrZ24HbnfvxEvat2GXwFVh9VtoS9saoHN5G+hJ237XL1hPQqwFrp/ATrF3CncAxbaep5vydGJFLMlwX3RfgvRFJPIJdxXhu4h37Azdg3UhOAf8C7JuI80Lgq9E7JoELUtKyCvhBVNZ3A6e2XT5O2tLKyABrnDj7YPU1v4mI4jPAAYn3HAx8Hmsn9Evgg8BeKXXynagc/tv9Rj+3S3X1oFAoGoXqZBQKRaNQklEoFI1CSUahUDQKJRmFQtEolGQUCkWjUJJRKBSNQklGoVA0CiUZhULRKJRkFApFo1CSUSgUjUJJRqFQNAolGYVC0SiUZBQKRaNQklEoFI1CSUahUDQKJRmFQtEolGQUCkWjUJJRKBSNYkGeuyQiAjwL65tWoVD4Yz/gZybAb++CJBkswUwWxlIoFGkYwp7B7YWFSjIPAUxMTDA4ONh2WhQLDNdfD29/OxgDIvDhD8PZZ7edqmJMT0+zZMkSCJwBtH5agYi8FXgr8Ozo0n8C7zXGfCHnmVXA+6JndmKP2fh8wDcHgampqam+J5nJSdi5E5Ytg6GhtlOjKMLkJBx8MOzZM3NtYAB27er9+puenmbRokUAi4wx077P9YLidxJ4N/ZQrKOBHcBnReSItMgi8jLgRmAr8KfYQ8duFpEju5Pc3sHWrbbBLl9u/27d2naKFEXYuXM2wQDs3g333ttOerqB1iWZNIjIb4B3GWPmdBsR2QY82RhzmnPt68Bdxpj1nu/ve0mmn0fEhYx+rrd+lmT+ByIyICJnAE8m+3zf47CnGbq4Nbqe9d4nishgHLAa8r5G1oh40022ISt6E0NDcO21lljA/r3mmt4nmCroCZIRkReIyMPYozc/CrzeGPP9jOgHAA8krj0QXc/ChdhjVOPQ991w2TLopNTeeefp1KnXsXatlVzGxuzftWvbTlGz6AmSwZ6D/GLgGOzh5f9PRJ5f4/vfjz38PQ59P24kR0QXe/bAunUq0fQyhobgxBPntwQToydIxhjzmDHmXmPMt4wxFwLfBf4hI/rPgcWJa4uj61nvf9QYMx0H5okRXjwiXnnl3Hv9rkycnLQjvRJl/6MnSCYFHeCJGffuAE5KXDuZbB3OvMbQEKxaNXfqNDAAS5e2k6aqWEirZm2Sade+bYxpNWCnMsdjbV5eEP3eA5wc3b8eeL8T/2XAH4F3As8DLgUeA44M+OYgYKampkwvYWLCmB077N9QbNlizMCAMWD/btlSf/q6gYkJYzodm484DAyUK5NuoWy9bdkyk9dOp7t1VubbU1NTBjDAoAnp4yGRmwhYe5ddWKXvL7ArRyc7978EXJd4ZhVWj/MocDdwauA3e45k6mhwExPGjI11t0NWIcY0DA/PJpg4jI3V8/66Ubbe2iTTst/uW5JpI/QayfTj6G1M/SPxxIQxInMJpqmyqEqQVeptx472yLTst8uSTK/qZHoedc5n+9EKdHISzjlnJt11rGjt3GmbexLnnlv/KkwZvU+yzqvUW5oJQpN6NDft3f5261JFG4GKkkzICO6OllkjZ5kRse5pSiiaGInHx+eWQ6dTfx7LlHdanVeVQLulR0tLe5lv63SpSyQT0rDcyhWZmQqkEVNIpbepMIxRtoNlkaObJ/d9TeQtlCDz8lqmsyYHnib1aHlpD/22kkyXSMa3gaZVblGH9Kn0kM7dtLQT2sGyyDEtT52OlWyaQChBFimjQzprtweIOiVOJZkukYxvA82q3KoV7dtoutWYfTtYXrm1oQT1Jcg6ldFtKPjr/KYqfruI884r3uCWtbcoRllF2777FivtmlDKZsHXPD5PSVqkiGzCaMx3/1Cdyug2FPw9sSEzhJHmS6CkJJPUsQwPz4wI4+PGXHHFbBHfHS1FZp4tq2vw1Vu0uTyahaIRNUuyaFv/lDWVS5vqFk1N25Jkduyw7bKq7kenSw2TTFYDGR83ZsWK2ddf9arZ5HPllfZvFSVfiN6iSmNO6yxFHci3gw0P509RkuVTV6esqpsqmlqFEGE3LbPrJmglmYZJxkfH4gYRY1avDq/krA4RKp1s2hQmOcUkkLbUmSW9GePXkIvekYVt28LynIYqHc1nFais+YGvHmvbNhvq0v+Mj5cnXCWZhkmmaLXIJxStArnEkOwQVZbON27Mz9uWLenKzU4nPc8hdiJlpZGsNIVIMr7pK1pST66EufHzyL+KBJXMv4i/PVZeuqpINkoyDZOMMbNF3bKE4y57JomlqEMlpZONG9OnNiGduix5Dgz4SRpFElghETPfAAAXSElEQVTW9CyL3HylIJ9vhyypDwykDwIbN6aXzYYN1SSoEIL1NRRMK88QAlSSCcl0xSXssbF061SfjhkbcPk+m9YhRIxZuTK9Yb3nPfnvSebliivC8uCG0dFqkow7WrsjdZZdSpYxY5bEUGSIFrqknqYTS6vHv/iL7Lr3Qd7U3MceK81QMKu9vf71/tKWkkxIpmvaIOlWYt5oEVe87wjj0yHSOmDa6JfVwEOILitf27bZ0b1IkZmm7EwbrUUseWflI6sjZUkMeVOOPAO7tPL2SVNR3NHR7LaU1P/4SjJZhPSOd9j6iVeVisrVR9pSkgnJdCDJxA1g+/a5y9Tx6tHmzemVd/HFtnEVzZXTKj2eGoQqnZMNM6lP2L69GsEkt0hs3OhnqezGyZpqXXyxfzouvjhMUgFLilmdWCR9ST0k5JVrXFZJySFJlKtXz01fFgkUDUAusa5bV9xW8upQSSYk0wEkkzXir16dvTcpraJjwshrFAMDxlx0kTFvfOPsRrdxY3lSiEfPqpJLnJaRkXCdT5rOJWtal/b+0DA2lj/lySI4l2SMCZ8Sx3oyH6knrtdt24rjxxJeaBtN1k9WvpNllwUlmZBMe5JM0SiRJlIXxU+OUiKWWEZH7UiTJSa705KQ6db4uG1cdayMbdkStpSeHKE3bpy9TJ5WPuPj+XF80pk35QBLZFn3fJTWaeGqq2YIKk0ZXDXkbd7csSNbko5DbKuVFydJskkoyYRk2pNkqkxTQoLPyBdPNcbGihtU/M40sTu0wyanQnm6FBehq1YDA7PtirLC2Wdn30tOKbI6e6djFedp191O5ksWaZLcpk3Vyj7r/a5kmJSki94TD3J5bUZJpqbQlCTTZNi+fSZNWVMNN1x+eblGHtvVjI7OKA7d6U6ePmN4eCa+j2juEkdVPREYc8klM2mMjdg2bMhO75lnZiuGfXQdMQFk2dAkDRzLBDdNRdPzvCm7G+eyy7Lv63SpphBCMq4pvBtWrqxHz+Eb4lHI53tVRtC0zueSiC95+DR43/f4xh0fT19RetObwsp527b8jgizFfpZ1tIxfCTPvPz7rjDWETZtyu4PSjIhmfYgmSSBrFtnzN///VxbjZARu6mQXO2pS0zP+55vvDIrNKHfAUv6eXqYpsphxYr85eYsy2VXf1eU5jxFdt1lm6fEV5IJyXQByWSNGmnLips3d3falAwnnzy7sZ5zTntpSQsjI3aUD+38PmR5+OGz66bX8p5mrAj++qc4bNrkbz9UR8iaMvUMyQDHtE0iVUmmzGbINonGDVlWqG2G0M5x7rnGnHZa++muEvK2XYQu03ezfeVtNeglp1U3NfDOrmLZMhDxj2+M/XvBBfnxjj4635FVHdizxzrVCkl/04jLxxdXXQW33NJMWgDWr5/txOmcc+r/ximnwCGHpNf3XXfNdV6VB2OK49dV3294Q/0OrcSEtgBAREazbgGnGGP2rZSqhiEig8DU1NQUg4ODqXE2bYLzzw97b6fj13hOOAG+8pW5nU/EvmP37rDvuhgYgDvugGOP9UuLSDgJuIg7UUinaRMDA9YTHliPdEuXwv33w0tf2r00xIRQpdxjnHkmvPnNNh+33w6nn149bT/9aTrRTE9Ps2jRIoBFxp4p74W9Sqblz4G/Bh5OphF75Gzf413vsgV+wQW2Aw0MWJeLH/xg9jO+He3LX06/bgzceKP9/4wzynXcN78ZvvlNv5HvAx+A/fef7aozBJ2OlZoWL7aEXEenaRJJ15Px3507/d9xySXw+ONw+eXl85s2uKRdd+9n3bvhBntulDEzklMVwjfGDlKrVpV/R8pLS+k0PgMcn3Hv39vWuVTVyRgzsyzp6hPOPLPZ+XCWO8p4Xh7HqWpk51qnGlNuR7mbpm6u6JTVM7jbOnxdS6SFNH2KSLVl6g0bbBpGR7Pf77NKl7WKFZqerI2cXVX8AgcCB7ZNFmVDEclkLTuWCb7vcf3DxAZtsYl98h1Vls6zFHt5lqBN5l3EmgaUVWwWGR2uWzfjWS5plvCmN82UhY+Frkh2uWdtCQmtk6xFh9FRG+quI9/2YUyXSAZ4OfBjYHcUHgA+EPrRtkMeydRp9BTvSyraYj8yYhupr3uBeE9SWekjzbduU0Z+ec+41rJllmjzVnDS0pH1/pUr/Qgmy0lVXI8bNpQvu1h6SHt/3PGLdlFXDUVuWrtFMncDnwf+DFgKnA18C7i3nySbPJKp2+ip07Ejdda9LGOuojA2ZhuEbwdLa7Tx1MGno3Y6xqxdW2/ZuGJ5Xl7Wr58racQWyFllW3eIN2/6bDVIC0XL1kW77U84oXoezjprbhrT9qe1TTJ/AJ6buCbYZeub2iaPOkimiiST9VzdOgufncZFYcWK2enttl4luekvr8xdV50jI8a85jXt2AJdeWW1vOZJOnXVQZ67Eddo86KLwk/O6BbJfBt4Rcr1w4GH2iaPgHwU6mTKmMNffnn9jT8p5ruOj+re0tAU0STLMimW50mProd9353NTeSj7BQ1lhR8rZ7rSvtRR+W/t8zRMt0imbcD3wGWJK4fC9zXNnkE5MNrdSnNLFwke8dwyCZG3waafN+pp86dNtTVkUTqcRqVfKfrVDve5e0iS88Ran7fVHBJMc/xe12hW/kNPeyvWySzJwqPADcA5wMXAv8JnFknETQZfEgmhs/JhskGWdZtQbzUun27Fc2vuMKPRKpsQkxreGUluZBOWzRVWr++vHK77uB6F+zGtDJeCGjyG92UZEKN3J8JnAJcHP1eA7wXWAacLyKfEJHzReTVge/tSUxOwqGHWuOk5JnJa9fOGM652L3bnlftnj+chdgIq9OB4WH4yU/gec+D177WGrm98522SRThfe8L264gAitXzr0enz8dnxM9OtrMNgj3/Oe086HBluHDD3fHkrjTyd4SMjAAxx1n28Jb3uJXH3Xguc+t/51xXXb9POwapIJ9gJcAbwGuBr4GPFj1vU0GAl09hDhxTo7SecrC7duLj2X1CVnOty+5JH26F3uxS54jlHamURNSTbKMslwl5EkyIna0j73xb99ePi1ZUzK33rvlJTH+bhPvTba3UPTMLux+CEUkU3Rej2sxWnS2cR5xuErc8fFyZyDFpJGVXp/05TlciuOMjeU74A7pQPGRKHE5vupV6XGTU7eBgfTzpuI0Js8kL0pH7Mkvy62Hqzuqw1eNq8SPFcKuEV9djr6Kyr4slGRCMl3S1cP69ekNPD4WJcujfJM6jthzWh6ZxCSR5kozzcl43nG6VUfZDRv8fNO6aYjLN03XlTzZ0TfEdVW0spU8uqRK3k87bUbyctMc4nWw0yleZcsjrTyL3iIoyYRkuqTTqrQK9z2OtGjqVCW4diRZ4nDa2T4+tilp76lCmL7+boaH56Y7q2OVIb44f0XOwpPlkGeZHE/xRkby859mdBeyYBDXcdZ3YiLLIq24bEOhJBOSaQ+dTNljLeo6d9odjdasyX7OZ5WgzMkBRQe1ZRHmyIjdgOmSbxki8vVrW+Ybvidz5p09XWT7k3VCZV6affLipsnnnPE816ChUJIJybQHyVRR9OXZH/hKAhs2zBWr0xpL0Rx7YiJM1+Pzzvi9RdOsNH1HmiTj7jJ2v59VB64/49CD73zen1YOaScRxFPQpPI+1rEVbSPIS2eWtORKIT51kHccbyiUZEIy7WmMV0YJ5ytZ+BBIVifNWyXIOpenKHQ6s4/T9UGRUjkrTtbZ2Mk8pXUi1zw+np6mKYd98pf2/rR4PiuNafGypqS+h/WlXU+eJlBm4UElmS4EH5IxJv14Dd9R0gdFJJA1JcnTl7j6ixA3E2VXHfL0QHlxfJ6L8+R2xjRF5vj43Pe5BJ6Xv7o6aVa82F1H2jdcaSjt2bS9TlnfzitLn8HAB31LMliL4W8ADwG/AG4GDit4Zk2UWTc8EvBNb4vfeHvB6Gi2CBxPbUJGBx99QN7SdOj7fKWXNKdObaNoCT3PjsmnXvLi+R7L66MfCSWCkCOBi+BbFnnoZ5L5YkQaRwAvAj4H/AR4cs4za4Ap4AAnLA74pjfJJOE2hvgYjjIV56sP8B2FfHRIWboPN28+04K2kEekZaYAPoRaVZIJHXiKjDPLTnXqQN+SzJwEwdOjjKS69zQzJFPaqrgKyRjjZ8Dm8w4ffUAc12dK4uMgKutddTfopiSivClmyAgfQqi+RF/XtKTpd5bFfCKZpVFGjsyJswZ4PJJ4JoDPAkfkxH9iRCxxOLAqydTRIetuQFkrCUmfvmmoUzRvWiIKmUamoUz91TH1Kosq76yT7OcFyQAd4BbgawXxjsN65XsxcAKwPZo+DWXEvzRFh1OaZHptruy+q2znq4s4uyXiVyHoOuuvl1E32c8XkhkBdmWRRc5zT8C6AH1fxv2elGSaQJXOV4dk1c0OXJage7n+6kITeex7ksHu4J4ADin5/E3AjZ5xK+lkjAnrkN1esakqXuc9W5SXfunAvaTraAJNkH3fkgzWR/DVwH3AspLvGADuAa70jB9MMmmdy6cz9/qKTQhCjNL6oQM3oT/J+1a3BxqVZGY6/GbgwUi34i5JP8mJcz3wfuf3xcCrgEOxJyfciHVy/nzPbwaRTFmi6JdR3QeheYnti+IzjxYy2hpo6ib7fiaZOQrZKKxx4nwJuM75fVW0svQo8HOsbc2fBnwzyBivLFHMJwVjaF7mkwRXBW0PNHVKa2VJpuxZ2LXBGCMecU5M/D4XOLepNLlIcw8Zu48scl+4bNncs4ljF5f9hpC8TE7OPl97zx5Ytw5WrOiiy8ceQZX2UweGhtov8wY8uM4vxJ3LhS9RDA3N9vXbdd+qNSIkL3kda6GhSvuZL1CSKUBVooidcicdkfcjfPOiHWsG82mgKQsxVkexoCAig8DU1NQUg4ODXs9MTtqReOnShdVAymLrVjtF2r17pmP1M8FWxXxoP9PT0yxatAhgkTFm2vc5JRlPklGEYz50LMUMypJM64pfxfxFLygdFe1DdTIKhaJRKMko5gUmJ61CenKy7ZQoklCSUfQ9tm6Fgw+G5cvt361b206RwoUqflXx29eYnLTEkjQS3LVL9UF1o6ziVyUZRV9DDf96H0oyir6GGv71PpRkFH0NtajtfSxonczExITqZOYJ7rsPfvQjOPRQOPDAtlMzPzE9Pc2SJUtALX6LISIHArrYqVCUw5Ax5j7fyAuVZAR4FvZAuSzshyWioYJ4Ci2rEPR7We0H/MwEEMeC3FYQFVAuE1seAuChENFwIULLyh/zoKyC06yKX4VC0SiUZBQKRaNQksnGo8A/Rn8V+dCy8seCK6sFqfhVKBTdg0oyCoWiUSjJKBSKRqEko1AoGoWSjEKhaBRKMhkQkb8VkV0i8oiI3CkiL207TXVCRI4Xke0i8jMRMSLyusR9EZH3isj9IvIHEblNRJYl4uwvIp8UkWkReVBEtorIvok4LxSRr0blOCEi56ekZZWI3BPF+Z6InNpMrsMhIheKyDdE5CER+YWI3CwihyXi7CMiHxGRX4vIwyLyaRFZnIhzkIh8TkR+H71nk4jslYhzooh8W0QeFZF7RWRNSnr6r12GHDe5UAJwOnaJ8W+A5wPXAr8FntF22mrM4ynAPwGvxx49+rrE/QuwZ5S/Fngh8FngR8A+TpwvAHcBxwCvAHYCNzj3B7HHCH8COAI4A/g9cI4T52XA48C7gMOB9wGPAUe2XUZR+r4IrInS/yLskcg/AZ7sxBkBfgosB44C7gD+w7k/AHwP+HfgxVHZ/xK43IlzCPA74IqoHP4uKpcV/d4uW09ALwbgTuBq53cHuw3h3W2nraH8ziIZQID7gWHn2iLgEeCM6Pfh0XNHO3FeDewBnhX9fivwG2BvJ84/A/c4v7cBtyTS83Xgo22XS0ZZPT3K9/FOuTwGrHTiPC+Kc2z0+xRgN7DYibMemIrLBvgAcHfiW58Cvtjv7VKnSwmIyN7Y0ei2+JoxZk/0+7i20tVlHAIcwOwymMI28rgMjgMeNMZ803nuNizJHOPE+Yox5jEnzq3AYSLyVCfObczGrfRuWS+K/v4m+nsU8ARml9U9WMnGLavvGWMecN5zK1bSO8KJk1kO/dwulWTm4mlY8faBxPUHsB1vISDOZ14ZHAD8wr1pjHkc2/ncOGnvwCNOz5W1iHSAD2GnQndHlw8AHjPGPJiIniyrsuUwKCJPoo/b5YLcha1QlMRHgCOx+ieFJ1SSmYtfEc2fE9cXY5WYCwFxPvPK4OfAM9yb0WrJ/ok4ae/AI05PlbWIXA2cBrzSGOM6PPs5sLeIPCXxSLKsypbDtDHmD/Rxu1SSSSDSH3wLOCm+FonJJ2FXDRYCfoxtuG4ZDGJ1LXEZ3AE8RUSOcp5bjm1TdzpxjheRJzhxTgZ+YIz5rRPnJGbjZHqkrKOl/Kuxq3DLjTE/TkT5FvBHZpfVYcBBzC6rF4iIS8onY32zfN+Jk1kOfd0u29Y892LALhU+AqzGrqJcg10qXNx22mrM477Y5dQXY1dCzo3+Pyi6f0GU578EXgDcTPoS9reBlwIvB37I7CXsRViyuh6r4Dwdu0ybXML+I/BO7KrMpfTWEvZm7FL+CVjdRxye5MQZwS5rvxKrnL0duN25Hy9h34pdBl+B1WelLWFvjMrhbaQvYfddu2w9Ab0asHYKP8HaJdwJHNN2mmrO34kRuSTDddF9Ad4bkcQj2FWM5ybesT9wA9aN5BTwL8C+iTgvBL4avWMSuCAlLauAH0RlfTdwatvl46QtrYwMsMaJsw9WX/ObiCg+AxyQeM/BwOexdkK/BD4I7JVSJ9+JyuG/3W/0c7tUVw8KhaJRqE5GoVA0CiUZhULRKJRkFApFo1CSUSgUjUJJRqFQNAolGYVC0SiUZBQKRaNQklEoFI1CSUbRl4hcUL6j7XQoiqEkoyiEiFwnIjdH/39JRD7UxW+vEZGkrxaAl2DdTyp6HOpPRtEKRGRvM9tjXhCMMb+sMz2K5qCSjMIbInIddjfyP0QnHBgReXZ070gR+ULkrf8BEflXEXma8+yXRORqEfmQiPwKuyMZETkvOqHgd9FpBpvjEw9E5ETg48Ai53uXRvdmTZei0wA+G31/WkRG3RMDRORSEblLRP46enZKRD4lIvs1W2oKJRlFCP4B67vkY8AzozAROWzagd1BfDTWofhiYDTx/GqsG4eXYx1pg/UJ/HasK4jVWJ80G6N7twPvwPpdib/3wWSiIr8qn8XuCj8B64flUKyTchfPAV6HdT51WhT33UEloAiGTpcU3jDGTInIY8DvjTH/441NRP4O+I4x5iLn2v/CEtBzjTE/jC7vNMacn3inq9/ZJSL/G/go8DZjzGMiMmWjmTzvbydhfd4cYoyZiL5/NvCfIvISY8w3ongdrPuEh6I4/xo9uyG0LBT+UJJR1IEXAa8UkYdT7j0H68wKrGe3WRCRPwcuxDpqGsS2yX1E5E+MMb/3/P7hwERMMADGmO9HCuPDgZhkdsUEE+F+Ei5EFfVDSUZRB/YFtmO96SVxv/P/79wbkT7nFqxnuQ1Yp0+vALYCe2MdPNWJPyZ+G1Rl0DiUZBSheAzrTtLFt4E3YCWFxwPedRS2k7/T2DOEEJE3enwvif8ClojIEme69HzgKcz40FW0BGVxRSh2AceIyLNF5GmR0vUjWKXrjSLyEhF5joisEJGPi0geQdyLPRjt70XkUBH5a2YUwu739hWRk6Lv/UnKe27D+tD9pIj8WXQ+9PXAl83sw+cULUBJRhGKD2KP5vg+1lftQcaYn2FXjAaAf8N2+A9hHXDvyXqRMea7wHnYadbdwJlY/Ywb53asInhb9L3zE6/BWB+yr8U61f4KlnR+hHW8rWgZ6uNXoVA0CpVkFApFo1CSUSgUjUJJRqFQNAolGYVC0SiUZBQKRaNQklEoFI1CSUahUDQKJRmFQtEolGQUCkWjUJJRKBSNQklGoVA0iv8PV0QdnCJMIzgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 300x300 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEJCAYAAAAw1SpDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGypJREFUeJzt3X+UVdV99/H3d0ZFBQRZix/RATVCiFErrRmigwqSEi2OUTQKi+VTwDGithYXaBISta2NsbYRSB4oiXEUeFafSEzQCBg0dY3VDFpHU0UpRiFOZIxhSDTMTCtEZ3b/2HfCZZw7c++ds+85997Pa62zLvfcc/bZxx9f9vme/cOcc4iIxKki7gqIiCgQiUjsFIhEJHYKRCISOwUiEYmdApGIxE6BSERip0AkIrFTIBKR2CkQiUjsFIhEJHaHxV2BOJiZAccB7XHXRaREtbscBrKWZSDCB6GWuCshUsKGAW3ZHlyugagdYPfu3RxzzDFx10WkZLS1tTF27NiczyvXQATAMccco0AkkgBKVotI7BSIRCR2CkQiEjsFoj7U1cGYMf5TRMKxcpyz2syOAfbt27cvY7L68MPhww8Pfj/sMPjgg8LUT6RYtbW1MWzYMIBhzrmsX9+rRdSLurpDgxD472oZiYShQNSLzZtz2y8iA6NA1IuLLsptv4gMjHJEyhGJREY5ooh98AFcfTWMHu0/FYREwinrIR79qa+PuwYi5UEtIhGJnQKRiMROgUhEYqdANABNTbBsmf8UkfwpWZ2n+fNh7dqD3y+4AO67D6qqYquSSNFSiygPTU2HBiGAxx+HsWMPvmlraYGGBv8pIn1TiygPzzyT+beFC+Hdd+ErX4GuLqiogHvv1Tg1kb4oEOXh3HMz/9bZCV/+MnR3WO/qgmuvhaFDoaZGj24ivdGjWR6qq2HevN5/q6g4GIS6dXXB7NlwwgnqJCnSGwWiPK1ZA88/D1deCWZ+X2Ul3H23D0a96eryj27KG4kcSoNeI1jFo6UFdu6E8eP9o1d9vQ84nZ29H3/qqfDqqwO+rEjiaNBrjKqqYNq0g/mfujpobobVq3s/fvt234pas6ZAFRRJOAWiQKqqYOLEvo9ZsMC3okTKnQJRQBMmZM4Xddu1Sy0jEQWigKqqfB+i/vz4x+HrIpJksQciM7vezLaZWVtqe9bM/qKfc64ws9fMbL+ZvWJmMwtV31zV1cHu3X0fc/rpGrMm5S32t2ZmdjHQCbwBGDAPuAX4U+fc9l6OrwGeBpYCm4C5wJeBP3POZfUuKuq3Zv0ZOhQ6OrI7dt48PapJ8Srat2bOuY3Oucecc2845153zn0N6ADOynDKImCLc+6fnXM7nHO3AT8H/rpQdc7F7bdnH4TAj2Grq4PTTvPnipSDRA3xMLNK4ApgMPBshsPOBpb12Pc4cGkf5Q4CBqXtGjqAauZkw4bcz7n/fv+5fTssXw7t7dHWSSRpYm8RAZjZ6WbWARwAvgPMcs79V4bDxwB7euzbk9qfyVJgX9pWsL7Nl102sPM7OtQyktKXiEAE/AKYBHwGWA2sNbNPRVj+XcCwtK1gQ0/vuAOGDBlYGQ8/HE1dRJIqEYHIOfcH59xO59yLzrmlwMv4XFBvfgOM7rFvdGp/pvIPOOfaujegoA877e1w220+71Nb23/fop5OOilMvUSSIhGBqBcVHJrTSfcs8Nke+2aQOaeUCHfcAa+8Ahs3wnPP5RaMfvlLvdqX0hZ7IDKzu8zsPDM7MZUruguYBvxr6vd1qX3dvgVcaGZLzOyTZvZ3wKeBlYWue76qq31Hx2yD0fbtMHmyH582bJjvcyRSSmIPRMAoYB0+T/QkUA1c4Jz7aer3ccDHug92zm3F9x26Fv8I9wXg0mz7ECVFXR2cfHLu57W1wZIlMGpU9HUSiUvsHRrjUOgOjZkcdRTs35//+ffcA4sXR1cfkYEq2g6N5WzGjIGdv2SJf7Wvzo9S7NQiirFFBL1PLZuvI4+EO+/0c2pXV0dTpkgu1CIqUl1dcMop0ZS1f79vJU2e7NddEykWCkQJ0HONtKjK1Ct/KRYKRAnQ16ogA9HYGH2ZIiEoECVE96ogy5fD8OHRlLlhg6YUkeKgZHXMyepMli3zj1fbtg28rJNP9quMiISmZHWJWbwYXn45mkc2zYstSadAlHBr1sC//MvAy1mwQMFIkkuBqAi880405Wj5IkkqBaIi8OST0ZWlxzRJIgWiInD55dGWt2ABfOITcNNNcMMNsGlTtOWL5EpvzRL61qynUaNg796D382iGxoCUFOjfkcycHprVuJaW/1o+ylT/Gf3BPtR2bpVLSOJj1pERdIi6s348T7nE5W/+itYWTTTy0kSqUVUhnbu9K2jwYOjKe/jH9eKsxIPtYiKuEVUXw/XXutH8EedM9KKs5IPtYjKTEvLwSAE0QYh0Oh9KSwFoiL1xhsHg1Aoy5eHLV+kmwJRkZowIff10XL1/e9DVcGWopRypkBUpKqq/JJElZX+e2Ul3Hdf9Nd5+23fx0gkJCWrizhZDT5XtHOnf5VfVeX7GW3dGv11nn9e82BL/5SsLlNVVTBt2sFHqKiHg3RbsAAmTtRqIRKGWkRF3iLqTc/hIFEbMgTa28OVL8VLLSL5o/ThILW10Zff0aGWkURLLaISbBH1dPbZ8Nxz0ZY5ejRs3Ki8kRxKLSLJ6Gtfi77MPXu0fppER4GoDNTWhnsFv3atH5/W0gINDf5TJFcKRGWisREWLQpT9pIlMHYsTJ8OJ5zgx8CJ5EI5ojLIEaWLeuqQ3lRWQnOzemWXI+WIJCs7d8IDD8Cll8Kxx4a5Rmen1lGT3CgQlaH58+Hhh+Hxx8Nd4+67lS+S7CkQlbHq6mgWcOzNli0+bxTijZ2UHgWiMhfVAo6ZfOMbcOWV4cqX0qBAJFx8sZ/hMZSHHtIka9K32AORmS01syYzazezVjN7xMwmZnHeTWb2CzN738x2m9lyMzuyEHUuNVVV/hV8SBdcELZ8KW6xByJgKrAKOAuYARwOPGFmGaeEN7O5wD8Cfw+cAtQBs4FvBK9tiVq0KOxEa++95x8Dm5o0Qb98VOL6EZnZSKAVmOqcezrDMSuBU5xzn03bdw/wGefcOb0cPwgYlLZrKNBSjv2I+pI+GX8haIL+0lNK/YiGpT7f7eOYrcCZZjYZwMw+DswEHstw/FJgX9qmF8u9qKuDX/0KTj21MNfTBP3SLVEtIjOrAB4FhvfWsulx7N8A3wQMOAz4jnPu+gzHqkWUo02b4Ic/hHXrol8hJN2MGfDEE+HKl8IqlRbRKuA0YE5fB5nZNOCrwA3AnwGXAReZ2W29He+cO+Cca+veAE3r1Y/aWv/Y9L3vHZwXO4Sf/lQdHyVBLaJU3ucS4Dzn3Jv9HPsM8Jxz7pa0fVcB9wJDnHN9ZjnKeaxZPlpa4NOf9lN/hHDssfBuXw/iUjSKtkVk3kpgFjC9vyCUcjTQM9h0dhcZZf3Ev96f02cbdWDeew/OPFMto3IWeyDCP45dBcwF2s1sTGo7qvsAM1tnZnelnbMRuN7M5pjZSWY2A/gHYKNzrhOJ3IoVMGhQ/8fl6+c/90NCbr013DUkuZIQiK7Hvyl7CngnbZuddsw44GNp378O3JP6/C+gHngcWBi+uuVr/37fCzukO++EK64Iew1JnsTkiApJOaKBCbV2Wjqto1acijZHJMWnsRGuvjrsNe69N2z5kiwKRJKX664LW/5998H55/v+TFL6FIgkL9XVMHNm2Gs89ZTPSekRrfQpEEneNm+Giy4Kf50XXoBzztFwkFKmQCQDsmmTTywvXx62hdTYqHXUSpnemumtWaSamuDGG+E//iPcNe68E/7yL7VKSBLl+9ZMgUiBKIiQMz52Gz/ez4mtVlJy6PW9JEro1/vglyxasMAHJClukQciM/tM1GVK8amvDzvjY7pduzTBWrEL8Z/KQwHKlCLU2QkjRxbmWg8+WJjrSBiH5XOSmf0g00/AiPyrI6WmtdW/WQs9Ru3FF/3o/TfegAkTlMguNnklq83sXeD/AB09fwLWO+dGR1C3YJSsLrxCBKNuFRV+iEhdXWGuJwflm6zOq0WEHynf3tvk9ma2Lc8ypYTV1kJNTfjBsuAn/1+40C9hpJZRccg3R3QjsKu3H5xzM/KvjpSyxkbYuBHGjQt/rc5O/1ZNikNOgcjMppjZm8BbwFtmtsfM7k496oj0q7bWrxQSup+Rmc9PadbH4pBri+i7wA6gGpgI3AL8OfBzMzs+4rpJCevq8jmjULM+OgezZ/vWV319mGtIdHJKVpvZ+8AZzrnX0/YZ8AMA51xRzK2nZHWyzJ4NP8j0HjYiu3crX1QIhepZvQMYlb7D+Uh2O3BhjmWJAHDzzeGv8Rl1s020XAPRGuD/mtnYHvuHAVlHP5F01dV++emQfv1rOPpoTSWSVLkGohXAGcAbZvb/zexLZrYUP3n9lyKvnZSNNWv8dCLnnhvuGu+/r6lEkirXHNFoYBI+GE1KbRMAh39sewXYBmxzzm2JvLYRUY4ouZqafLAITZPzhxHbNCBmdiRwOocGqNOcc8MHVHBACkTJ9id/Aq+8EvYalZXw4Ydhr1GOCt2z+o+cc/uBptQmMmDbtvl8zvvvh7tGZ6fv06TJ+ZNB8xFJIv3P/8Btt4W9xubNMHEi3H572OtI/zRDox7NEq2lxS9FHdqRR4ZtgZULzdAoJamqyq9xVlkZ9jr798OnPhX2GpKZApEkXl0dNDdDQwOMCDjb1Y4dfg5sKTw9munRrOgMHQodPWfCitD69X7KEg0JyZ0ezaRstLfDqFH9H5ev2bN9Xqq6Wj2xC0WBSIrSnj2waFHYa7zwgnpiF4oCkRStFSvg9NPDX2ftWrWMQlMgkqK2bRuceWb46zQ2hr9GOVMgkqL3wgswfXrYa0yZErb8cqdAJCXhySfDvnqfPBkuuUQLOYYSeyAys6Vm1mRm7WbWamaPmNnELM4bbmarzOwdMztgZq+b2cxC1FmS6etf9zMxzpoVpvxHH9US16HEHoiAqcAq4CxgBnA48ISZDc50gpkdAfwUOBH4An7+7C8Cb4eurCRbVRVs2BB2nNquXbBsWbjyy1HiOjSa2UigFZja27ppqWOuw0/c/0nn3Ad5XEMdGsvAqFGwd2+48u+7T4s49lRKHRqHpT7f7eOYzwPPAqtSSxq9amZfNbNeRySZ2SAzO6Z7A4ZGXGdJoNbWsL2jr7lGyxVFJVGByMwq8NPRNjrnXu3j0I/jH8kqgZnAPwBLgFszHL8U2Je26T+fMrFhQ9jytYhjNBL1aGZmq4G/AM5xzmUMFmb2OnAkcJJzrjO1bzFwi3PuY70cPwhIX0FrKNCiR7PyMH++75QYwsiRvuUlXtE/mpnZSqAWOL+vIJTyDvB6dxBK2QGMSSWyD+GcO+Cca+vegPbIKi6J1z0x/5VXRl/23r0wZ0705Zab2AOReSuBWcB059ybWZzWCIxPPcp1+wTwjnPuDyHqKcWtutqPqg/xNm39+rCDcMtB7IEI/+r+KmAu0G5mY1LbUd0HmNk6M7sr7ZzVwAjgW2b2CTO7CPhqqiyRjJ58Mky5e/fqDdpAJCEQXY9/U/YU/pGre5uddsw44I+5H+fcbuACoBq/fNG3gW8B/1iQGkvRuvzycGXffz+cdJLvY3TDDZqYPxeJSlYXivoRlbfQ/YvS1dSU14DZok9WixRKayvcc48fyHrPPX4LZetWtYyyMeB1zUSK0eLFfus2YoQfRxbCD3/o11CTzNQiEsH3NTr55DBlr10Lt2bqaiuAApHIH+3cCVdfHabsO+/0rS4Nlu2dktVKVksP48f7EfahlHJvbCWrRSKycyc88ACccUaY8vfuVcuoJwUikV7Mnw8vvQQXXxym/NCDcYuNApFIH372szDlllPfomwoEIn04XOfC1e2mf9cs0bzYStZrWS19KOyErq6wpXdmTaHxMknF/ccR0pWiwTS2emXoR42DMaNi77sdLt2lWfLSIFIJAsPPgi//z386lcwZEjYa/34x2HLTyIFIpEctbfDokXhyj/3XGhoKK/5sBWIRPKwYoVfxaM74RylW27xK9eecALU10dffhIpEInkqa4O3noLVq+OdobG7sR4VxcsXFgeLSMFIpEBqKqC666DPXtg40aYGfFaw52dxf0WLVsKRCIRqa2FefOiLbOyEpqbS3/GR81HJJJghx12cJ6k1atLd8ZHtYhEIlRTE20C+8CBQ7+X6oyPCkQiEaqqgu99zz9Sgf+M+s3a2rV+9H5TU7TlxklDPDTEQwJoafFJ5vHjfXD6/Od9Mjtq8+Ylqyd2vkM8FIgUiKRAvvlN30coas8/7xeQTAKNNRNJuJtvhhtvjL7czZujL7PQFIhECujb34YvfCHaMseMiba8OCgQiRTYQw9Fu5ZaKSxVpEAkEoPFi+HUUwdeznHH+WR4sVMgEonJli0Df7VfUSL/B5fIbYgUn559jioq4Pzz4ZOfzL6MlhY477xkvcLPh17f6/W9xKxnnyPwgSXXJbCTMM2s+hHlQIFIkm74cNi3L/fzHnjAL4UUF/UjEikRy5blF4SgeKeZVSASSZgf/Sj/c084wY9Bu+YauPzy4hkgq0AkkjCXX57/uVu2wOTJforZDRv8SrVHHQW33x5d/UJQjkg5IkmgUaNg795oyxwyxE/8H5JyRCIlpLXV976eMgWuvDKa/kIdHcltGcUeiMxsqZk1mVm7mbWa2SNmNjGH8+eYmTOzR0LWU6TQFi+Gn/0M1q+H556LpsyHH46mnKjFHoiAqcAq4CxgBnA48ISZDe7vRDM7Efgm8EzA+onErrraL180ULNmDbyMEBKXIzKzkUArMNU593Qfx1UCTwP3A+cCw51zl2Y4dhAwKG3XUKBFOSIpNi0tfgHG5ub8zj/jDHjppUirdIhSyhENS32+289xtwOtzrlslqBbCuxL28pgpSgpRVVV8Ld/2/tvQ4f2f/7LL/t8U9Kmmk1Ui8jMKoBH8a2bc/o47hzgQWCSc+63ZrYGtYikjIwfD7t2Hfw+YgT87nf+z0ccAR98kF05UU81WyotolXAacCcTAeY2VDg/wFfdM79NptCnXMHnHNt3RsQ+CWmSFg7d/rhHJde6j+7gxD4vkPZWrs2GS2jxLSIzGwlcAlwnnPuzT6OmwT8J9CZtrs7oHYBE51zuz5y4qFlqB+RlKyWFhg7NvvjBw/2r/ajULQtIvNWArOA6X0FoZTXgNOBSWnbo0BD6s+7A1ZXJPGqqnJ7w/bf/+1zRnGKPRDhH8euAuYC7WY2JrUd1X2Ama0zs7sAnHP7nXOvpm/A74H21Pc/xHIXIglSVwe7d0NDA1xxRf/Hb9gQvk59ScKS09enPp/qsX8BsCb153H4xy4RyVJVld+mTfODX3/yEx+Yduz46LGXXVbw6h0iMTmiQlKOSMpJfT1cey10ZfirfORIP6QkCkWbIxKRcFpa+g5C1dUwdy6cdBLcdFNh65ZOLSK1iKSENTTA9OnZHz9oEOzfn//11CISkY+YMCG3kfsHDmTXQztqCkQiJayqCu699+BKIdno6PDLHDU0+M6ODQ3+ES8kPZrp0UzKQPdKId/5jp9WJFcVFT6g1dX1fZxW8ciBApGUs/p6WLgQOjv7PzZdZaUf9d/XyrLKEYlIVurqfEBpaPCdHocMye68zs5w66YloUOjiBRYd2dH8PNYZ7P0dWWlH/UfglpEIoJzcOKJ/s8jRnw0uV1ZCd/9bt+PZQOhFpGIAPBm2nDz7uT24MF+UGxzMzz/PIweDbW10V9byWolq0X6NGUKbN168HtNDTQ29n6sktUiErlNmw4NQuC/R72CrAKRiGT02GO979+yJdrrKBCJSEYzZ/a+/8ILo72OApGIZFRb63NC6Wpqok9Y662ZiPSpsdHnhLZs8S0hvTWLiN6aiYSht2YiUrQUiEQkdgpEIhK7sk5Wt7Vl/QgrIlnI9/+pck1WHw8EnnNOpKxpYrT+mJkBxwHtcdclzVB8cKwiWfUqhHK+dyjN+293OQSXsnw0S/0DejvueqSzgxPCtOfyN0kpKOd7B90/KFktIgmgQCQisVMgSo4DwN+nPstNOd876P7LM1ktIsmiFpGIxE6BSERip0AkIrFTIBKR2CkQFYCZLTWzJjNrN7NWM3vEzCZmcd4VZvaame03s1fMLMPEncmVz72b2alm9iMzazYzZ2Y3Faq+Ucvz/r9oZs+Y2Xup7d/MbHKh6hwHBaLCmAqsAs4CZgCHA0+Y2eBMJ5hZDfB9oB74U+AR4BEzOy18dSOV870DRwO/BL4C/CZ4DcPK5/6n4f/dnw+cDexOnXN82KrGR6/vY2BmI4FWYKpz7ukMx6wHBjvnatP2PQe85Jy7rjA1jV42997j+GZghXNuRei6FUKu9586pxJ4D/hr59y6kPWLi1pE8RiW+ny3j2POBv6tx77HU/uLWTb3Xsryuf+j8S2pkv1npkBUYGZWAawAGp1zr/Zx6BhgT499e1L7i1IO916SBnD/dwO/5qN/MZWMshx9H7NVwGnAOXFXJAblfO+Qx/2b2VeAOcA059z+UBWLmwJRAZnZSqAWOM8519/EbL8BRvfYN5oiTd7meO8lJ5/7N7Ob8Qn7P3fObQtZv7jp0awAzFsJzAKmO+fezOK0Z4HP9tg3I7W/aOR57yUj3/s3sy8BtwEXOudeCFnHJFCLqDBWAXOBS4B2M+vO8+xzzr0PYGbrgLedc0tTv30L+HczWwJsxjfPPw1cW9CaD1zO925mRwCfSh13BHC8mU0COpxzOwta+4HL5/6/DNyROq857ZwO51xHQWtfKM45bYE3wGXY5qcd8xSwpsd5VwC/wE8P8SowM+57KcS9AydmOOepuO+nQPffnOGcv4v7fkJt6kckIrFTjkhEYqdAJCKxUyASkdgpEIlI7BSIRCR2CkQiEjsFIhGJnQKRiMROgUhEYqdAJCKxUyCSxDCzManJ8heZ2X+mFg3YbmblOn9R2VAgkiSZlPq8Grgp9f0t4F9TsxtKidK/XEmSM4APgEucc//unHsNuBUYh58KpNbMfmFmb5jZNbHWVCKl+YgkSSYBG5xzzWn72tL+vAy/xM4+4EUze9g597sC1k8CUYtIkmQS8FKPfWcDvwXGAtudc287PznYT4DPFbh+EogCkSSCmR0FTAAq0/ZV4HNFa4HjgLfTTnkbKNkFB8uNApEkxen4WQivMrOzzewUYD0wHPh6rDWT4BSIJCkmAa8B3wB+BLyAbx1Ndc79Hr+uV3oL6PjUPikBmipWEsHMVgHHOufmZvj9MGAHfl34fcCLQI2S1aVBb80kKSYBGzP96Jz7MLWiSQO+Jf9PCkKlQy0iiZ2ZGb6VM8c591jc9ZHCUyASkdgpWS0isVMgEpHYKRCJSOwUiEQkdgpEIhI7BSIRiZ0CkYjEToFIRGKnQCQisVMgEpHY/S8hPsb2EMrsiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "MCP.plot_chain_panel(chain[burnin:, :], figsizeinches=(3, 3))\n",
    "MCP.plot_pairwise_correlation_panel(chain[burnin:, :], figsizeinches=(3, 3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Print log files\n",
    "Log files display a date/time stamp associated with when chain segments were appended to the correponding log file.\n",
    "\n",
    "| Date    | Time | Start | End |\n",
    "|--------------|---------|---------|---------|\n",
    "|2019-05-17| 07:30:40|\t0|\t999|\n",
    "|2019-05-17| 07:30:40|\t1000|\t1999|\n",
    "|2019-05-17| 07:30:41|\t2000|\t2999|\n",
    "|2019-05-17| 07:30:41|\t3000|\t3999|\n",
    "|2019-05-17| 07:30:41|\t4000|\t4999|\n",
    "|2019-05-17| 07:30:41|\t5000|\t5999|\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------\n",
      "Display log file: 20190517_073038_serial_chain/binlogfile.txt\n",
      "2019-05-17 07:30:40\t0\t999\n",
      "2019-05-17 07:30:40\t1000\t1999\n",
      "2019-05-17 07:30:41\t2000\t2999\n",
      "2019-05-17 07:30:41\t3000\t3999\n",
      "2019-05-17 07:30:41\t4000\t4999\n",
      "2019-05-17 07:30:41\t5000\t5999\n",
      "2019-05-17 07:30:42\t6000\t6999\n",
      "2019-05-17 07:30:42\t7000\t7999\n",
      "2019-05-17 07:30:42\t8000\t8999\n",
      "2019-05-17 07:30:43\t9000\t9999\n",
      "2019-05-17 07:30:43\t10000\t10999\n",
      "2019-05-17 07:30:43\t11000\t11999\n",
      "2019-05-17 07:30:43\t12000\t12999\n",
      "2019-05-17 07:30:44\t13000\t13999\n",
      "2019-05-17 07:30:44\t14000\t14999\n",
      "2019-05-17 07:30:44\t15000\t15999\n",
      "2019-05-17 07:30:45\t16000\t16999\n",
      "2019-05-17 07:30:45\t17000\t17999\n",
      "2019-05-17 07:30:45\t18000\t18999\n",
      "2019-05-17 07:30:45\t19000\t19999\n",
      "2019-05-17 07:30:46\t20000\t20999\n",
      "2019-05-17 07:30:46\t21000\t21999\n",
      "2019-05-17 07:30:46\t22000\t22999\n",
      "2019-05-17 07:30:47\t23000\t23999\n",
      "2019-05-17 07:30:47\t24000\t24999\n",
      "2019-05-17 07:30:47\t25000\t25999\n",
      "2019-05-17 07:30:48\t26000\t26999\n",
      "2019-05-17 07:30:48\t27000\t27999\n",
      "2019-05-17 07:30:48\t28000\t28999\n",
      "2019-05-17 07:30:48\t29000\t29999\n",
      "2019-05-17 07:30:49\t30000\t30999\n",
      "2019-05-17 07:30:49\t31000\t31999\n",
      "2019-05-17 07:30:49\t32000\t32999\n",
      "2019-05-17 07:30:50\t33000\t33999\n",
      "2019-05-17 07:30:50\t34000\t34999\n",
      "2019-05-17 07:30:50\t35000\t35999\n",
      "2019-05-17 07:30:50\t36000\t36999\n",
      "2019-05-17 07:30:51\t37000\t37999\n",
      "2019-05-17 07:30:51\t38000\t38999\n",
      "2019-05-17 07:30:51\t39000\t39999\n",
      "2019-05-17 07:30:52\t40000\t40999\n",
      "2019-05-17 07:30:52\t41000\t41999\n",
      "2019-05-17 07:30:52\t42000\t42999\n",
      "2019-05-17 07:30:52\t43000\t43999\n",
      "2019-05-17 07:30:53\t44000\t44999\n",
      "2019-05-17 07:30:53\t45000\t45999\n",
      "2019-05-17 07:30:53\t46000\t46999\n",
      "2019-05-17 07:30:54\t47000\t47999\n",
      "2019-05-17 07:30:54\t48000\t48999\n",
      "2019-05-17 07:30:54\t49000\t49999\n",
      "\n",
      "--------------------------\n",
      "\n",
      "--------------------------\n",
      "Display log file: 20190517_073038_serial_chain/txtlogfile.txt\n",
      "2019-05-17 07:30:40\t0\t999\n",
      "2019-05-17 07:30:40\t1000\t1999\n",
      "2019-05-17 07:30:41\t2000\t2999\n",
      "2019-05-17 07:30:41\t3000\t3999\n",
      "2019-05-17 07:30:41\t4000\t4999\n",
      "2019-05-17 07:30:41\t5000\t5999\n",
      "2019-05-17 07:30:42\t6000\t6999\n",
      "2019-05-17 07:30:42\t7000\t7999\n",
      "2019-05-17 07:30:42\t8000\t8999\n",
      "2019-05-17 07:30:43\t9000\t9999\n",
      "2019-05-17 07:30:43\t10000\t10999\n",
      "2019-05-17 07:30:43\t11000\t11999\n",
      "2019-05-17 07:30:43\t12000\t12999\n",
      "2019-05-17 07:30:44\t13000\t13999\n",
      "2019-05-17 07:30:44\t14000\t14999\n",
      "2019-05-17 07:30:44\t15000\t15999\n",
      "2019-05-17 07:30:45\t16000\t16999\n",
      "2019-05-17 07:30:45\t17000\t17999\n",
      "2019-05-17 07:30:45\t18000\t18999\n",
      "2019-05-17 07:30:45\t19000\t19999\n",
      "2019-05-17 07:30:46\t20000\t20999\n",
      "2019-05-17 07:30:46\t21000\t21999\n",
      "2019-05-17 07:30:46\t22000\t22999\n",
      "2019-05-17 07:30:47\t23000\t23999\n",
      "2019-05-17 07:30:47\t24000\t24999\n",
      "2019-05-17 07:30:47\t25000\t25999\n",
      "2019-05-17 07:30:48\t26000\t26999\n",
      "2019-05-17 07:30:48\t27000\t27999\n",
      "2019-05-17 07:30:48\t28000\t28999\n",
      "2019-05-17 07:30:48\t29000\t29999\n",
      "2019-05-17 07:30:49\t30000\t30999\n",
      "2019-05-17 07:30:49\t31000\t31999\n",
      "2019-05-17 07:30:49\t32000\t32999\n",
      "2019-05-17 07:30:50\t33000\t33999\n",
      "2019-05-17 07:30:50\t34000\t34999\n",
      "2019-05-17 07:30:50\t35000\t35999\n",
      "2019-05-17 07:30:50\t36000\t36999\n",
      "2019-05-17 07:30:51\t37000\t37999\n",
      "2019-05-17 07:30:51\t38000\t38999\n",
      "2019-05-17 07:30:51\t39000\t39999\n",
      "2019-05-17 07:30:52\t40000\t40999\n",
      "2019-05-17 07:30:52\t41000\t41999\n",
      "2019-05-17 07:30:52\t42000\t42999\n",
      "2019-05-17 07:30:52\t43000\t43999\n",
      "2019-05-17 07:30:53\t44000\t44999\n",
      "2019-05-17 07:30:53\t45000\t45999\n",
      "2019-05-17 07:30:53\t46000\t46999\n",
      "2019-05-17 07:30:54\t47000\t47999\n",
      "2019-05-17 07:30:54\t48000\t48999\n",
      "2019-05-17 07:30:54\t49000\t49999\n",
      "\n",
      "--------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "CP.print_log_files(savedir)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}