Chain_Log_Files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<form action=\"index.ipynb\">\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: July 19, 2019\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": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.9.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from pymcmcstat.MCMC import MCMC\n",
    "from datetime import datetime\n",
    "import pymcmcstat\n",
    "print(pymcmcstat.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a simple model and sum-of-squares function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "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.\n",
    "\n",
    "To generate a set of results and save them to a specific directory, the following code can be executed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "datestr = datetime.now().strftime('%Y%m%d_%H%M%S')\n",
    "savedir = 'resources' + os.sep + 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)\n",
    "\n",
    "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`.  In this case, we are going to reference a saved solution from the `resources` directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the folder `resources/20190517_073038_serial_chain` matches the pattern that was specified for `savedir`, and we display its contents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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 resources/20190517_073038_serial_chain/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, there are log files saved in both binary (h5) and text (txt) format.  Note, if you run this simulation on your machine, the results folder (`savedir`) will be different because of the date/time stamp."
   ]
  },
  {
   "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": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymcmcstat.chain import ChainProcessing as CP\n",
    "from pymcmcstat.chain.ChainStatistics import chainstats\n",
    "from pymcmcstat import mcmcplot as mcp\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the plotting class and define the directory in which to find the log files.  If you want to look at results you generated, then chage the value of `savedir` accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define directory where log files are saved\n",
    "savedir = 'resources' + os.sep + '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": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary: 0.05322208404541016 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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text: 0.4920275926589966 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": 24,
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\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",
      "==============================\n",
      "Acceptance rate information\n",
      "Chain provided:\n",
      "Net    : 20.32% -> 5079/25000\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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAEKCAYAAAAmUiEiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dfZQlRXn/P8+9ihhgVjkq6M6K4C4EASVhkRc5vKzyQ5FETXbBI0cXf3uyC0aM6AaEjUg0os7KYn7BXSG7gRhfsuOJB876RuQwGg2EwZeoqCirrM6sihpkBlRAd+v3R3Vl6vZUdVf37b597536nlNn5nZXd9frt556nqeqRClFRERERF1oNZ2AiIiI4UYkmYiIiFoRSSYiIqJWRJKJiIioFZFkIiIiakUkmYiIiFoRSSYiIqJWRJKJiIioFU9oOgFNQEQEeBbwcNNpiYgYYBwA/ETlePQuSJJBE8x004mIiBgCjAK7syIsVJJ5GGBqaoqRkZGm0xJRIXbvhh/8AJ77XFi8uOnUDC9mZ2dZsmQJBMwGFirJADAyMhJJZoiwbRusXQt790KrBTfcAGvWNJ2qCFmICyRFZASYmZmZiSQzJJiehkMO0QRj0G7Drl0wOtpYsoYWs7OzLFq0CGCRUmo2K260LkUMBe67r5NgAPbsgZ07m0lPxBwiyUQMBZYt01MkG+02LF3aTHoi5hBJJmIoMDqqdTDttv7dbsP118epUj8g6mSiTmaoMD2tp0hLl0aCqRNFdDIL2roUMXwYHY3k0m+I06WIWjE9DRMT+m/EwkQkmYjasG2bNiuvWKH/btvWdIoimkDUyUSdTC2Ifit+TE9rk/uyZYNbFtFPpgZ0I/YP25QhJD/Rb8WNBSndKaUWXABGADUzM6NCsHWrUq2WUqD/bt0a9FjXz/YjQvMzNTUXz4R2W19fqBimMpmZmVGAAkZUXn/Li1B3AC4H7kYvtPo5cDNwRM4zfwF8CfhVEm4DXljgm8Ek003DGKZGpVTx/Gzdqu+beINOsN3i9ts7y86EiYmmU1YcRUimH6ZLpwEfBE4EzgSeCPy7iOyX8czpwMeBM4CTgKnkmcrX3XYj9g/blKFoftasgTvvhE2b9N9+XqzYiyntgvVKzmOhXgfg6WiGPLXAM21gFnhdYPwoyZRAGUlmEKaKvUznsEh3AzVdmpcgWJok/ugCzxwA/BY4JzB+YZ1M2YYxLI3KIDQ/g0KwTaRzakpPkZooi6kpPW3r9tsDSzJoa9engC8XfG4z8ANgX8/9JyXEYsLiIiSjlL9hhFRak42qDoTkp5f6h246zjDpSfJQpcQ2yCSzBdgFjBZ45m3Ag8DzM+JclRRIRyhCMi7083SgqhGrm+/3QkLotg7qTGfTdZBOS5X5HEiSAa5DK3APLfDMeuAhYHlOvK4lmTT6eTpQRceronPUPVWsqg7qSGe/DUBVS2wDRTKAJASzG1hW4LlLgRngxBLfLKSTcaFfxexuO17VnaPOqWKVdVBlOvtROlrQkkyiT3kIbco+2ApPtuJ8GHiP9fsy4DHgz1PP7B/4za5JZnKyukqrUqwO6Xi+7/WzdOZClXUQipC6qmsA6nYAqFJiGzSSmacrScIFVpwvADdZv3d5nrkq8JtdkczWrUqJzG/cZSqtiAdtCBHlEUXW9/pVOnPBzke3dVDmm732dq7qnVVJbANFMk2EIiST7txTU/MJBvSoWhSuhtNqKbV9e2cjKDqCbdw4F7/dVmpsTOchb+Tvd0nG1IUrH61WuTpwvd+V3157O6fT0m8DQCSZikjG7twiSq1fr9Tmze7KHh/Pr5g0fA3HJhNf456cdHeIdJpXruz87Uu7eVe/+vWk81V1h8sj8jKd3CU1ZE1V7TpIp6XfBoBIMhWQjKtSs0IZksn7RrutpRofCaU7RNE0mw7ratBFROq8qdzUlM5HWkKrqpzyOtzkpFLXXOOXdEI6cLedfGpKD1KueksTqGsqvnFj5/VulPJV6AAjyVRAMllShqujlu1A9tTGR16hHaxImn0Nuui0I08CSOuvRIp3jpB8rV6t46Y70OrVnfFWrpzfwUKllLJS3tat/npzTf2yBhX7d7ftrRuiiiSTl+kAkhkbK95Ri1RaemQz73I1wvXr5xq3r0EayaOIJOObdoTmI2t0N9KL6xtFxfxQSSbdga64wh8/TwL0deIiUyBzz1fOoNSmTeXryTVVy0qLq02XnXJFksnLdA7JFOmsLmLIqzSXdcq8y+4kZ501Xye0Y8f8Z0XmvplW+q5ePUdQISNkaD6mpvQUxPVsmjyzOohLse7TNeURbfp6VudO59ElcRmFedY00M6reca+//a3Z3/fJcmIZNeZCRs3+ttVWmLMIrsyuqxIMnmZziGZ0GlHllSRRpZlpNtgRt303N40+Kkp97QrPfqHNj6X+dhFlFlhcnL+VGv16mxRfnJSj/w7drgljzJlZ5TeWfUSMg1Md/60DsX1TkNiGzfOn4aFtBNb2pqczLZ6+tp02WlXJJm8TJeQZFotpS68sLMxjI25dRp5Fp8qCcbuLHmKSZ9OoYhTWxnlso+MikoaNgHZVjNTF2XSFVoveab+ouHiizvzMzamSd1YDX0K/3RYt06TVFY5+6xTMF8aCkUkmbxMB+hkfB3SnpO7RFB76mLiVy25uILPtJ729N28Wakrr5yv3A1VahZVLncb8nRNRmKbmtIdzu64daSlrjJIS5UhJBwaWi2lbrxRqXPPnU9sZRFJJi/TgX4yhlB8PikhVom8Bimi1DnndN8p8iQZn5XHNi9PTnYqNW39SJ3Tvaxgyj6vDNMd8oILynXSEB1VnjK3bHANWCG6maLhrLP8bToUkWTyMl3A49flkGfPg7OsK3mdMlR/kRfMN7OkL5+i2WdezvLdOOSQahp7SCcyI24viC2LOOxpRV0k4wpXX13Pe7s1Y/eEZICnAucBb0nCq4Gnln1fL0MoyfgIIq2ktDv71q3ZCs2qO0R6apOezt1+u55GhL63iO9G0ZDWaRlJanw8u9OaqUQdo3q6PH1h3bq5Mu71lLHuUMaMXTvJAGuAe4FN6D1d1gPXAt8F1lRBBHWGUJ1MkdHKOLH5fEdcpueQcPXVbjNp2l8j7VWb5+SXFUJ8N8oEMyUbH5/vvJjnl2TyOjHRXd584UUvyo+zbt0ceZf5flWSax2hqBm7FyTzPWA/x/X9ge9XTQpVhyr9ZNKd02cV2LKlHMn43MztDpplXSga6pRk0vlJi+obNvjTlB5pe60bstM/NhburOl6NsSPqJehXyWZe4GDHdefCXyvzDt7Garyk8nrROnrvvhlvmNGnjKN3YT16+evh3H5btSpf7AbuK/czznHrYxOO+itW6elpFWr6kuvCUWmoHYwEu/kpM5XXluqI+22LqzsIthekMw5iTTzb8D/S8Ink2tBJwY0GaqQZFat8usIQsXiVktPo+xlAyENy1Yud9u5jU5kfHy+W77tiwJKHX549Q0eOr1/s5TkLmc9n5v/5s1KXXKJUu9+dz1SQzcEkGeevuSSORJNx7viinwdVkgYHy+/Mb5SvSGZxUk4Cb073Z8n/7ebJpDA9He9dsl4nl50kft+3ihlgpEetm/XFR8yDTBKyCIS14oVfmWxUuFTxBUrsvOS/r1lS/Z70wsyfVMmX2d1+fu43PTr9J+pMoTs7zM5WV6SSn/Drn/fKnEXaiMZ4EXA/cCeJDwAjIV8qJ9C1auwfZ2naIM2FWtPA7LiFTGl+tbXVJVfUGr58s7fZ589f0rmCsY1IM8V3/es7e+TJQkZEs9aGGqv9cpKS9b0tyyRpck/aylA0TKyvaNt1wYzNXblJ2vJQZ0kcw/wGeCP0YewvQ74KrATWNw0eRTIR+X7ybiCPQ0q0tDsacy6ddkrmdNbGRT9hhGNd+zoLq9NhtCOvWWLrluXlGpLVGaNlM+LGrTE5frmhg3Za7uywpVXzrcWdtMGWy3dhkxd29OjUEL37ZNUJ8n8Fjg8dU2ATwCfaJo8qiQZpeYrFs89N7yCjVi7eXPx0Xl8PMwC0c3c/MIL619P1W/BXlltk4Q9uocQxIYN89uH/Q2zWr7sWipb91Rkqvfyl2crdM2gUmRKWgXJSNLpgiAiXwPepJT6cur6kcCkUuqA4Jc1CBEZAWZmZmYYGRnxxpuehjvugAcfhAMPhG9+E979bv97Wy19IH2rBWeeCZ///PwD6tNx56dN/82rlnYb/uEf4A1vyI4X4UarBe99Lxx/vD7wfnQU7r4bTjzRX2cAF18Mr3oV7L8/3H8/7NwJf/M38+ur1YLXvhY+8hHYs8df3yEQgbe/Xaf1Fa/wv2d8XKdDBJ7zHHjkkbl03n47/OM/FkuDCPz4x7ps0pidnWXRokUAi5RSs5kvKigBvAn4OrAkdf1EYHfTEkqBfAQ54xUdia68UksIISt6jVOavaivyKK4U0/N1i30WjqpQ6FaZMQtE+y9kkOnD0UdNHfsmFv/VsS1wfe+LO9xey+ZKurf9nLuRpLJvOnonHuT8CjwMfQBa5cD3wbOb5o8qiKZOldOG2tSXQ5ZrZZ+d52d04R3vENPB8fHs/UXxoelaMPftKl+oqmbjG0rTdokXWZdliHG006rv36zLEx1TpcOAo4FXpD8PRZYlnzsu8C3gG8C31RKfS74xT1G3nRpYgJWrKjz+7oaq8Yll2hx/9BD88X+KmCmAFn5abXgllvg+9+Ht7613Df+7M/gk5+sPz91od2GXbv0tGN6Gu68E/7nf/QU/OST4dZbYd06Pa0Kwfg4nHdePW0oDTvtNmqbLnmkgn2B44G/QB83+2XgoaallbolGbNRUl3rfMqMyHUvyuxmROzm+aJrv8qYnusO9tog35EnIW0p6wSLXqS9jCSTeXNYQx7JpBtCupEas6BSza2hSafJtb6p7KLMfgyhZG78ZnybZp9/fj3ln7WXcp6Dncjc4tE8gsza5a5MaLfne3b70h5JpkAIIRml/COG8Wco6w9RNlx8sXtjI1/8KtzPzTeWL+8cfUPeW9aE6yLMUDK3dyZML5PI0guJKPWa18wNIJOTSl17bbauyTxn/GvsVeLpbVqNP1KW02OI06Kt3+l2bZlRTPvKtSqdTObNYQ2hJJM3YvRaSjDOVOPjYe763YjVy5d3WivS64SyNrUyviJlpLy0VGAU2emjYbKCcbozdVhEyW535JA2sHq1e72Pa1sKkXIK8HS9pkmtG2k6S0JsxON3WEIoyShVXFrJMivn3c8ijrSXbt7aFeN4VpWkZUjLt7H32NjcQkvf5uWhnrn2eVOuqUWeU6RNFGWmFq4N2EOWEaRH/rqslC4Jo4xk02opddRR2XF8+8xEksnLdAGSUSpcIsjqRGkXb9vj1t5Rz7eEIL3jXmgjqnpHOZ/ZNet87vRufekphStknZwZ4gcScqpmVlnanStPV+L6rlL17qCXtcixKn+frH1mIsnkZbogyWTtkWuL9XkVZxqu63wkn0LQJqYyo+K5585tEF7VjnJ50lreCl67XF16o26nenZ5++rO6COy1oWZuio6vTH1PDmZHc9V3+lpp33An+87dnnm1a+RmEPysX69v+4iyeRluiDJKOV2pLL3qM06KdDE95FFu51/pEk3o6Ld6c1oV+d+uemOmgUXCRhlaqhy+cYb3WVqvu9zTpyYcJ9TbZ+rnaf38imp885OMvGMTsWW9NLTztBzsfJWbRuJOKQtZeljlIokkxvKkIxSnUpXo58IncIYPYGvgvNG1G7n9+lGYysNJyb8x99edlnYu32dOF1+RQ67D931z5CFb3Ggj9izSD9rqmUUuNu3z7ck2ebgLJLMkhJcSB8/7JIUs/afSe8HnVWeIbvlRZLJy3RJklEq26qSF4qYYl0Vne5IRU9RzGvYvlXFWem89tqwkdblgKZUdidXKv/AszQRu3Z7y/p+Fsm50tZqda7gNjqqiYnwJRB5UkJem8s6lC2PbO2tS9PpOuMMf/mlEUkmL9NdSDLd6jTSo67vfb4l9kaaMpJUumNlSVYhU5gi5tC07qKMJJH3bDrPrnOjQ1CW5NJpy1J6+wjZpeQPRV76fM+kycLnZbxli16Dlt5hMA+RZGoimW6tBa5Rt8g51Er5O4uNqSm/mTfv6IvQ3dh8ZOAaCUNO2iwihRjJIVQaqILkzPd8efH5m9i6uCJpLlJ23ea/DCLJ5GW6R5KMa9tDl14i9BzqIo2lbMPKm9en5/dly62IYrjbDtItyYWkxycBlj3QPu97Vee/KCLJ1EQySs2f6lxxhV9MNhYCW9TPctzKa+RFG0soeVX1XJF32q72Wei3kdxY51zTn3TbCDnQ3jXopNFtffh0S0WnSDYiydRIMkq5jw5JdwLTiLvVkdgo01jKiullnwt5ZxbZup6pgiCqIE7fSQhFBwvX+0LKoZv6cCn1Q/2ZXIgkUyPJ+Bq9a8/YkOlVUZG16sbSa4SQRnp0r0qy6qajZtVlGdKrQ0+Sh6L6vywUIZlW5mYzEfNw333zN0/aswde8hL40Y/0hle7dsGaNe64Ntptvb9sEaxZozc9alk1t3ev3vRoerrYu6rC9LTOd8j3feW3c6f+f9s2OOQQvWnYIYfo32vW6DK1y7YMRkfh9NPnb8AUkv6surTTH1oWeeVQBx55pPffBKIkUxTdKl/tZ8qOyHUo8sqiiMivVHnHuKbTHyLJFJ3+9DqvVX4zTpdqJBmlionvLoWnT2QPUQKaeL1uoFWmw1d+PvL0+Qz1Ov2uqWrW1DivLOpQsOehqm9GkqmZZJQqNr8PiVtUImiigabRjUTlKhOftFCXzqlM+m3/Jjv9ZcvCfl/IAFMFqlDqR5LpAclUiW58Wqq2ABVBHRKVzxpX1BU/BFWbtsu+q+gAUwdCpWiDSDIDRjL9pGMpijokKt+K9DqmTVWmv8y7+mHqW4bkIskMGMn0Q0PrBlVLVL4tEurUzVSV/qLvanqAKdv2BsqELSKXi8jdIvKwiPxcRG4WkSMCnlslIveKyKMi8i0RObsX6a0Do6Nwww3apA367/XXu48H7Uf4TMNlcfLJc8f1GojASSdV8/40qkx/0XctW9bpjgDlXBvKohem9MZJBjgN+CD6qNszgScC/y4i+/keEJGTgY8D24A/Am4GbhaRo+tPbj2oyhekiM9Kv2J0VJ/bbDpfq6V/DwrpFkHTA0xPSC5P1Ol1AJ6OFsNOzYizHfhU6tp/AR8K/EZfTZeqQj8oEKtE04rtXqLJvJbRJdV2TG0vICJLgfuAY5RS93ji/BjYpJT6gHXtb4FXKqVe4Ij/JOBJ1qUDgGnfMbWDiOlp7SFri76+I0YjItKYntZTpKVLw9pLkWNq+2G69L8QkRbwAeA/fQST4GDggdS1B5LrLlwOzFhhgCcTbjThph5RDfphilu1Xs1GX5EMWjdzNPDqit/7HmCRFYZqbJ+ehl/8olkFYkQ5uNZqDRv6hmRE5DrgHOAMpVQep/8MOCh17aDk+jwopR5TSs2aADzcdYL7BKaRnneeNkAaohk0C9VCxPQ0rF07J4E2vdC1LjROMqJxHfAqYIVS6v6Ax+4EXpy6dmZyfcEg3UiNem18vDsLVURvsFCmuE9oOgHoKdJrgFcAD4uI0avMKKV+CyAiHwZ2K6UuT+79PfBFEXkr8Gn09Go5sLanKW8Yrka6dy88/elRghkEGPNxWlk/bFPcxiUZ4CK0nuQLwE+tcJ4V59nAM80PpdQdaGJaC3wDWIm2LGUpi4cOTTtyRXSHpn1keoW+M2H3AiIyAswMgwl72zY9j9+zZ66RxmnSYKGo+bgfUMSEvaBJZmpqauBJBmD3bvjhD+Gww2Dx4qZTE7EQMDs7y5IlSyCSjBsispgh9JWJiGgAo0qp3VkRFirJCPAs8k3ZB6DJaDQgbkQsr6IY9PI6APiJyiGRfrAu9RxJoWSyL4DMLQV+OE8kjIjlVRRDUF5Bae4H61JERMQQI5JMRERErYgkk43HgL9N/kbkI5ZXMSyI8lqQit+IiIjeIUoyERERtSKSTERERK2IJBMREVErIsl4ICJ/KSK7ktMQ7hKRFzadpqohIqeKyA4R+YmIKBF5Zeq+iMg7ReSnIvJbEblNRJal4hwoIh8VkVkReUhEtonI/qk4zxeRLyVlOSUilzrS0venT4ScrCEi+4rIB0Xkf0TkERH5NxE5KBXn2SLyaRH5TfKejSLyhFSc00XkayLymIjsFJELHOkZjDaatwnwQgzoFeCPAa8HngfcAPwKeEbTaas4ny8D/g69l49Cr2S3718GPITehuP5wC3AD4F9rTifBf4bOAE4Bb0/88es+yPozcQ+AhyF3pbjN8BaK87JwO+BvwaOBN4FPA4c3XQZpcrjc8AFST5egN5m5EfAflacLcCPgRXAceg9jv7Tut8GvgV8Hjg2qYNfAFdbcQ4Ffg1ck5THG5PyOWsQ22jjCejHANwFXGf9bqE9hN/WdNpqzHMHyQCC3nJjvXVtEfAo8Ork95HJc8utOC8F9gLPSn5fBDwI7GPFeS9wr/W7q9MnGiyzjpM1kvJ5HFhpxfnDJM6Jye+XAXuAg6w4F6L3nt4n+f0+4J7Ut/4V+NwgttE4XUpBRPZBj0C3mWtKqb3J75qOF+tLHIremN0uhxl04zblcBLwkFLqK9Zzt6FJ5gQrzn8opR634twKHCEiT7Xi3EYnbqX/y3tR8vfB5O9x6HPD7DK7Fy3Z2GX2LaWUvRH+rWiJ7ygrjrc8Bq2NRpKZj6ehRdoipyEMI0xes8rhYODn9k2l1O/Rnc6O43oHAXH6trw9J2scDDyulHooFT1dZmXLY0REnsyAtdEFuUAyIqICmJM1Tmk6If2OKMnMxy9J5syp697TEIYUJq9Z5fAz4Bn2zcRKcmAqjusdBMTpy/LOOFnjZ8A+IvKU1CPpMitbHrNK73s9UG00kkwKie7gq1inISSi8YtZWKch3I9usHY5jKB1LaYc7gSeIiLHWc+tQLeru6w4p4rIE604ZwLfU0r9yorT96dPBJys8VXgd3SW2RHoPartMjtGRGxyPhO9bcJ3rDje8hi4Ntq05rkfA9o8+CiwGm1BuR5tHjyo6bRVnM/90WbUY9EWkEuS/5+d3L8syfefAscAN+M2YX8NeCHwIuD7dJqwF6HJ6sNoxeZ5aPNs2oT9O+CtaGvMVfSnCXsz2qR/Glr3YcKTrThb0GbtM9DK2TuAO6z7xoR9K9oMfhZar+UyYY8l5fEG3CbsgWijjSegXwPaN+FHaF+Eu4ATmk5TDXk8PSGXdLgpuS/AOxOSeBRtvTg89Y4DgY+hd3abAf4J2D8V5/nAl5J3TAOXOdKyCvheUt73AGc3XT6ONLrKSgEXWHH2RetrHkyI4pPAwan3HAJ8Bu0v9Avg/cATHHXz9aQ8fmB/Y9DaaFyFHRERUSuiTiYiIqJWRJKJiIioFZFkIiIiakUkmYiIiFoRSSYiIqJWRJKJiIioFZFkIiIiakUkmYiIiFoRSSYiIqJWRJKJiIioFZFkIiIiakUkmYiIiFoRSSYiIqJWRJKJiIioFZFkIiIiasWC3EhcRAR4FnqjpYiIiHI4APiJytmUakGSDJpgpnNjRURE5GEUfaicFwuVZB4GmJqaYmRkpOm0RCxw7N4NP/gBPPe5sHhx06kJw+zsLEuWLIGA2UDjJCMiF6GPMn1OcunbwDuVUp/NeGYV+rzk56DPXr5MKfWZot8eGRmJJBPRKLZtg7VrYe9eaLXghhtgzZqmU1Ut+kHxOw28Db2z+3LgduAWETnKFVlETgY+DmwD/gi9g/7NInJ0b5Lbf5iehokJ/TdicDA9PUcwoP+uWzd89dg4ySildiilPqOUuk8p9X2l1AbgEeBEzyN/hT54fKNS6rtKqbejj+R4Y6/S3E/Ytg0OOQRWrNB/t21rOkURobjvvjmCMdizB3bubCY9daFxkrEhIm0ReTWwH/5Dqgb1cPbK4RoJ166F8fHhGw2HEcuW6SmSjXYbli5tJj11oS9IRkSOEZFH0OfHfAh4lVLqO57ohQ9nF5EniciICWjT28DDNRLu3QvnnRelmkHA6KjWwbTb+ne7Dddfr68PE/qCZNCHeh2LPgJ1C/DPIvK8Ct9/OfrgMROGYpx3jYQGwzq/HzasWQO7dmmd2q5dw6f0hT4hGaXU40qpnUqpryqlLge+gda9uFDmcPb3oI9LNWEoxor0SJjGMM7vhxGjo3D66cMnwRj0Bck40AKe5LlX+HB2pdRjSqlZExgiT18zEo6PD9/8PlrNhgONk4yIvEdEThWR5yS6mfegzwH+aHL/w8k1g78HXioibxWRPxSRq9Cm7+t6nfZ+wegorFo1XPP7hWY1a5pQa/1+04dxo/1ddqGVvj9HW47OtO5/geQAeOtaV4ezAyOAmpmZUf2EqSmlbr9d/+3mHRMT3b2jaUxNKdVqKQVzod0ejDyVqcOtW+fy22rp371Eme/PzMwoQAEjKq+/5UUYxtCPJNN0Q+sn3H57J8GYMDHRdMqyUaYOmybUst8vQjKNT5ciBt/zs2pRe//93df326+a96dRRfrL1mHTDnm9+H4kmZKosmM13dC6QR26k0cecV//9a+7f3caZdOfrv+yddiEQ56d9p58P0/UGcZAl9Olqqc2ZUXWKnQ43aAuUX9ysjdTiLLpd9V/N2WxdauOa56pc6rsSnuZ70edTI0kU6Yx2WTgI4aiFd0POpxudCdZ5eAq3zryVyb9WQTYDVn0QmGf1XaLfj+STI0kU7Rh2p1GRAcfMYRWdNPKwm7T4SNI1/taLd2x+yH9W7fO1Z+v/ot21l5Ko1Uq1CPJ1EgyRRqmK24VxFCksdTdiIuO3lnl14RVKTT9WXVZth57LY1WOThFkqmRZJQKb5i+TtNtBwrVWfSqERcZvbOIJE+cr4ssQ9LvS3fZcm1KGq1K/xNJpmaSmZpSavt2pcbH5xqFubZ9e+e1qiWZUJ1Fv0yp0shLl6sT9IP+qeqpXBNSmyHqycnu9T+RZGokGVeD37ixc64uMtcRNm7s1MmY/8uMIkUaereNOC05VOWNfPvtukyyRlNbsugnsiwytcorq17nq2qijiRTE8nkSSZ2EFFqw4ZOghkb686KUBW/9P8AABjCSURBVFQX041J1W6Qq1d35mP9er8Oyte50u8cGwsrhypH/F4s2yjSmXtluq6D0CLJ1EQyITqWrJBVsbYo6+sIRRuLLUWFNmKXvidEF5HVubpp5FV1kG5G8lByKqNTKmJRTE/HQ+Frt+Pj5Uk3kkxNJFNEkvEF1wjs0rP4OkLo6Jc2nY+N5ecvy0SbRZp5RNCNNJJOk5GCinSOEKIK8dtJm9vT8devd+dz/frupirpMrCn4778pqe66fzbU/fQ9mEjkkxNJKNUZydvtcp1SqU6JZciptG00tnV2Mt4zJYl0ImJfBIp28nzOkdoh81Lny3x5fnttNvu+FNT7rbQanUnifne63uHjxRD2u3GjWFpUiqSTG6owrpkRFy78vIIxq5wexTJ68QGWboSnxI6T3IwnXv79vw8+Bp6CIlkSWD2SG2P0j7JII+4XcTsS9/YmP+dWWbr9O+3v90dd/ny8LpwIWuKnn5HXj2YdptV16HTsUgyeZmueKuHqSktWbg694YNnXPuIhJDuoF0M1Vz6QZc07TQkJ62hEzjXPoH10gtoqWxUClxYiJb57J16/xnVq/2SwnmnT5JKrSMfHFbLb8Umta9FJFkfIR09dVKXXPNnBUyK9+u8nMhkkxepguSjD21cfnCZEkCmzZ1NihfPNcIZHfibpTO9kI4+ztFOkz6fbblzJ4yFLWc+crjyivd19NpbrWU2rEjW+Hq66S+b4vM5SFUUnW9f906//tdUqhL97J6tTvPPhLPGzRWr9Zxr7giP/1Z9RhJJi/TBUjGN9qnG4q9LskV154OpOO123MOUpOTWipat26+ybeM1GGkgiqU1iatrrSEWM5c0xjfNGPLFn+5h+YjS18EfunTJhmlwi1uJlx77ZxEUqSeQkjf1GfR9moHY8EMKT8f4qZVFSG9EZENpeCf/3nunuYu9xElptrS/xtcdhk885n6+he/qM9Nuv76zg2QLr8c3vc+/8kEPoyN6XePj7vzEYpWS79j1y5Yvjx875T0fi0bN+q9TN7/fv37Xe9yf+u44+Atb5lfnqb8RLLTa/ZEWbbMHVcEfvlLXb9pKNWZl0ceCS+7dhtWrtR7K598cn467W+m24Uvnm9fnelpOOwwuC5nt+ubbvJvDGbjK1/JjxOEPBYaxkCgJNOtX0yREDKKmanIpk1h79ywobhZ2oRVq/w6Fp8eJa00LCM9tdudEmKZ4NLJZJWBazpi58OlHLafzdJFdaP38pWPS9HtW+2fVUarV4d/K404XcrLdCDJVDXFqCrs2KHTFaoU3bKlHMGsXDn3nU2b9HfTfhdZU8MrrihnsXrHO4pPTXzvMek0OrTJSV0evjT7/FDy2oCZvrgU2vYeQuvXl9PtpDu9z0LpIv28chRR6tJLs+P4pkyRZPIyXYBkfI0jrZPpRWi1NAGEfDOLXI48Mv8755/vbrjr14eTR6iewQ5r1xbPTzpMTrqVqCEmcbsMtm9XavPm/Ljj453EkuV341No29895RT3vS1b5t5TZbvLKluf70wkmbxMB5BMWsw991ylLr54/qZTZfxL6ghpJ7WylqMqGmY6XtER3EVuoc+uXJntGFdHuYhkk7+9c16e6TjL6nP11bpt1jGN96XLN2VqlGSAE5omkW5JJnSkaLX0SNf0lMpMM4wncK+ILzTf4+Na7F61KvzdNpnnEYMtmeVJLHWTb1YZuMqr3S4mnW7dWsyHqIrgmjI1TTI/bppEuiWZoiNFEbNqXcHulGXN3WW/e/bZ/vvGj2X79uJp2rLFb+IepJDlk+Mz1WeVdy/bWmOSDDDuCZ8AHmmaRLolmTJz3lYrX9Gap1Mpo8PwNQx7z5ZeNMTzz6/+vb0YrU1dGP+fc86pPg/r1/sV2j6HvX4ot6yFmL0gmQeBlwOnpcLpwANNk0i3JKNUOdPj+Hi+ctFYXy680O2yXpUyeWKiUwFZd/BNB/o52E6QWd7B3XZU+286uBZQlg3r1s3lJWTKfPbZ+fltzIQNfBI41XPv802TSBUko5QuYJsMQkyCRRufq+G73OSLhFZLqRtvDPcgrcIvZXIy3++in4LLp6XoNPlP/qQaUrrwwjmpM0+azbpnewJXZYEylrM0ekEyi4HFTZNF2RBKMmlJ4JWv7I0Ibxb82Uvzzzqr0+nLtabFhGOOyf/GO97RuUdxt/4pVU316gxmSmRLLnn7rviCSwIR0eRephzXrevcvmN83K2vybPWGR8l12r8PKnKFXpOMsCLgPuBPUl4AHhfyIf6KYSQTJaXZ90dwTQ0Yy0ylqPNmzuJoawVyafMy5NE6lQ6imgXgTLvD5Hc1q3r9EhOd0KzpahN7lnvO/fcYtdD82FLV1k72vnWXYWUc6jeKe35bKNOkrkH+Azwx8BS4HXAV4GdgyTZhCh+8yowtMPlWV/SnX9srHMnNddoZDt3FW1oxvrkynNIfuy0VRXsqcvkZLl3hC61MMr3rLLZvFmpl7wku0596eyWiPO29zD3u3FTCGkzeds91EkyvwUOT12TxKr0iabJo0A+KjFhb9miPTjzGlWWr4chFrPyOqQB2A0x1Iv19a/vJC5DNEU2rTL6ote+NuyboZ3e1iX4JMgsy5xvu4e6QhbJhDybV8dmiuLyWh4bc++PUzScddb8dJmy9G0Ub6NOkvkacIrj+pHAw02TR4F8dG3Ctkes9P+hFX3ttbqxrl9fXCIxDlKhjd31/lWrii2qO/zwYmksmpcsycx23a9iHVC3IVRycoW8+jYWyLpI01jVNm+ec+QsuhdQnSTzJuDrwJLU9ROB3U2TR4F8dG3Cdk1jtmwppkBdtarcvDpkk+4yocwSgCLvLpsXl25gx46wlcZ1dtKy+qM6fGNc5f2a12hyTntan3yye21VEdRJMnuT8CjwMeBS4HLg28D5VZFA3SHUulR01zZ7h7OqGrhLwjDnN+VtRF4mGCtH1R201dLvtXfUS5uRfVMls0jPXoDYlCXLTnNoOspIuVWGvO+WOWKmTpI5CDgrIZePAd8Bfgc8DnwD+Ehy76V1EUQVoYifjEvxFnLCQLfrh+y58caN2Q1o9epOi8hrXlOeJOzFfFVLNeee69YL+coaNCkpVf2+LGWC2fGubHrKWoN6kbeih+X1dO0SsC9wPPAXwHXAl4GHun1vnSGUZJTyb5Cd1ch8m1CHBEMuO3bozZ+vuSasIRrdjktP5AunnZbd4Iy/Rp3TjrypUjdlWWWHTlt9etX5Q50qq6qHUMStHvIyXUCSyTqgPO98o5DRzp4bG8mljOesy3HLrKfySQiuzpu29ph82JLS8uXVNfw8BXbofrTdemNn7ciX1lv0cnuPvP1nygbflDUUkWTyMl1Q8ZulHMs7CiRramVMrzaBlTWN+hqjy3vYnqa4pkWu/BrrQ5U6oKokGTP1MukrurwjawHjm988n3TrIBmfB3FdUozJd1EJxiCSTF6mS5iwjYXDd1ZO3iHsPv2GXdlTU0pddFG5RpN1LEheGoucOOkjg9e9rtjWDIaQTXlu2OCO4zqaxOV3ZHtKh7oEmGlmVr7M9/L2OC5KKBs2dA5Oxgkz9HRS22embFqyPHrzEEkmL9MlnfGWL3dLN74jP2yYTu5zGut21EpLLFkHrKXT6huZfSdOZklmeef5gNYFpTe/doX01o9mz2HflphFvZFFwtctuU7EzFpZ/fKXZ3/bLA9Jr5Y3g06oc6QhVt/+xWYanlVmZRBJJi/TFSwrMBXoOuQsC934tYT4mmRJLK4poK+zZCkDsySzdjt/B7xQPxPfEb0uQi67ZYL5RohEkCZdXx4MOU5OKnXJJf62YyQ5lxQZ6sFs0pS1zimrTZdR+ioVSSY3hOhkimw6XaTSilhKTj9dWxeuvVY3Opena6jizvVdnwUq5J1TU36vV6MbufZa/wia5zGbt4bHdhxst8vVV9ZZ3llxbUkwS9+lVPagkrVjXkgbCV3nlNemi5qvlYokkxtCSKYbs2lepYX6WNg+K2lx2mfx8uUlxBRugm95f0gZpRt2WV+jkNXIW7ZosvKdkJlnkjbfyPI0tuP6jAE+6dH4S2XVdYiLgC8f6alO3nQ5pL5CEUkmL9OBJuwyDlehlRa6BsfVCEOkJdehX1Wm3y6jPD1Qlq9RWvnp66xZUpg99Ut/J+0OsHKl+4wkHxHae88UqYeQw9Zcjo9Z+q70NddxJUWMEGXN10oNGMkkyxLuBh4Gfg7cDByR88wFSQbt8GiBbwY74xVZkFem0rK2yMwSp33SUtGTBF0jdlGELK7LGu1DJLI865Lx73G9z0zd8s6Qzup8WSZ2V15d1kl7B7z0N7JcBNrtfOtbERRdDOnCoJHM5xLSOAp4AfBp4EfAfhnPXADMAAdb4aAC3wwmGQO7Ebg6rr2ZVJF35k0ZioygZaZ46RHbvCfPWtYETB1k6THKjsz2+0PLtqiZ37w7r4O7CK8IyfUCA0Uy8xIET08S79xDWM2RTOmlC2VIxoZLR1KmQ/oajr1GxnwvRMQNsVylFaZZZzd322nrQh45lx3d84g1tB6q0H2kyahKfUoVGHSSWZok/uiMOBcAv08kningFuCoAt/oimSUmptGddMhi0opIVOSEGVn1tSlnxpyFvLWjpV9V149lpnedaP7qPudZTGwJAO0gE8BX86JdxJ6689j0Uex7EimT6Oe+E9KiMWExVWQTBUdsuqG4zNTpqUjF+oQyeucehXxVM5KX13EWoXuo+p3VlUfg0wyW4BdPrLIeO6J6H2G3+W5f5VDUdwVyVTZIatsjN10mqo7XC+mXt2SdL/pOupElfUxkCSD3iZiCji05POfAD7uude3kkwd6KbjVSVZ9bJ8uiHpfq7HKlF1PgeKZNAbkV8H7AaWlXxHG7gX2BQYv2udjFLFFIG9ttZ02/Gyng3JzyBJCP2k66gLVdfHoJHMZuChRLdim6SfbMX5MPAe6/eVwP8BDkMfz/Jx9EkKzwv8ZikTtqtjhTg/9bu1pghC8zNoEkId+pOsbzUx6CxkSWaeriQJF1hxvgDcZP2+NrEsPQb8DO1b80cFvlmIZMoSxaB1tDyU8XoddgmhKJocdKqsjyIkI0p3ugUFERkBZmZmZhgZGcmMOz0NhxwCe/fOXWu3YdcuGB3N/s7EBKxY4b5++ulFU908yuTn7rvhy1+GU06B44+vNXl9j27aUpVp2LkTli7t7puzs7MsWrQIYJFSajYr7hPKf2Zh4L77OhsFwJ49uqLyKmnZMmi15jeqpUurT2cvUDQ/27bB2rU6fqsFN9wAa9b0Jq39iG7aUlUYHe3dtwxavf3c4MF0LBuhRDE6qjtWuz333PXX976Sq0KR/ExPzxEM6L/r1unrCxXdtKVBRiSZHHRLFGvWaHF4YkL/HfSRPDQ/WaP2QsWwDTqhiDqZHJ2MQVVz2YWCftA/9CuGoS1FnUwNaGIuO8gwo/a6dVqCWSijdggWWluKkkygJBNRDsMwakfMR5RkIvoGC23UjpiPqPiNiIioFZFkIoYK09Pa8rWQTeX9hkgyEUODbdu0RWvFCv1327amUxQBUfEbFb9Dgmgy7y2KKH6jJBMxFIjOf/2LSDIRQ4GF6rI/CIgkEzEUWKgu+4OABa2TmZqaijqZIcPu3fDDH8Jhh8HixU2nZngxOzvLkiVLIEAns1BJZjEQjZwREd1jVCm1OyvCQiUZAZ6FPho3CwegyWg0IG5ELK+iGPTyOgD4icohkQW5rCAplEz2BdBcBMDDeSJhRCyvohiC8gpKc1T8RkRE1IpIMhEREbUikkw2HgP+NvkbkY9YXsWwIMprQSp+IyIieocoyURERNSKSDIRERG1IpJMRERErYgk44GI/KWI7BKRR0XkLhF5YdNpqhoicqqI7BCRn4iIEpFXpu6LiLxTRH4qIr8VkdtEZFkqzoEi8lERmRWRh0Rkm4jsn4rzfBH5UlKWUyJyqSMtq0Tk3iTOt0Tk7HpyXR4icrmI3C0iD4vIz0XkZhE5IhVnXxH5oIj8j4g8IiL/JiIHpeI8W0Q+LSK/Sd6zUUSekIpzuoh8TUQeE5GdInKBIz2D0UbzzrFdiAE4D63xfz3wPOAG4FfAM5pOW8X5fBnwd8Cr0OcavzJ1/zLgIeAVwPOBW4AfAvtacT4L/DdwAnAKcB/wMev+CPq88o8ARwGvBn4DrLXinAz8Hvhr4EjgXcDjwNFNl1GqPD4HXJDk4wXoM9h/BOxnxdkC/BhYARwH3An8p3W/DXwL+DxwbFIHvwCutuIcCvwauCYpjzcm5XPWILbRxhPQjwG4C7jO+t1Cewi/rem01ZjnDpIBBPgpsN66tgh4FHh18vvI5LnlVpyXAnuBZyW/LwIeBPax4rwXuNf6vR34VCo9/wV8qOlyySmzpyf5P9Uqn8eBlVacP0zinJj8fhmwBzjIinMhMGPKCHgfcE/qW/8KfG4Q22icLqUgIvugR6DbzDWl1N7k90lNpasBHAocTGc5zKAbtymHk4CHlFJfsZ67DU0yJ1hx/kMp9bgV51bgCBF5qhXnNjpxK/1f3ouSvw8mf48Dnkhnmd2LlmzsMvuWUuoB6z23oiW+o6w43vIYtDYaSWY+noYWaR9IXX8A3ekWCkxes8rhYODn9k2l1O/Rnc6O43oHAXH6trxFpAV8AD0Vuie5fDDwuFLqoVT0dJmVLY8REXkyA9ZGF+QCyYiICvBB4Gi0HioiA1GSmY9fksyZU9cPQiswFwpMXrPK4WfAM+ybiZXkwFQc1zsIiNOX5S0i1wHnAGcopex9iX4G7CMiT0k9ki6zsuUxq5T6LQPWRiPJpJDoDr4KvNhcS0TjF6MtBQsF96MbrF0OI2hdiymHO4GniMhx1nMr0O3qLivOqSLyRCvOmcD3lFK/suK8mE6cSZ+Vd2LSvw5tjVuhlLo/FeWrwO/oLLMjgGfTWWbHiIhNzmeit034jhXHWx4D10ab1jz3Y0CbBx8FVqMtKNejzYMHNZ22ivO5P9qMeizaAnJJ8v+zk/uXJfn+U+AY4GbcJuyvAS8EXgR8n04T9iI0WX0Yrdg8D22eTZuwfwe8FW2NuYr+NGFvRpv0T0PrPkx4shVnC9qsfQZaOXsHcId135iwb0Wbwc9C67VcJuyxpDzegNuEPRBttPEE9GtA+yb8CO2LcBdwQtNpqiGPpyfkkg43JfcFeGdCEo+irReHp95xIPAx9M5uM8A/Afun4jwf+FLyjmngMkdaVgHfS8r7HuDspsvHkUZXWSngAivOvmh9zYMJUXwSODj1nkOAz6D9hX4BvB94gqNuvp6Uxw/sbwxaG42rsCMiImpF1MlERETUikgyERERtSKSTERERK2IJBMREVErIslERETUikgyERERtSKSTERERK2IJBMREVErIslEDDSS7Sff3HQ6IvyIJBMRDBG5SURuTv7/goh8oIffvkBE0vu0AByP3noyok8R95OJaBQiso/q3DWvEJRSv6gyPRHVI0oyEYUhIjehVyL/VXLKgRKR5yT3jhaRzyY79T8gIv8iIk+znv2CiFwnIh8QkV+iVyMjIm9JTin4dXKiwWZz6oGInA7cCCyyvndVcq9jupScBHBL8v1ZERm3TwsQkatE5L9F5LXJszMi8q8ickC9pbZwEUkmogz+Cr1vyT8Cz0zCVLJZ0+3o1cPL0ZuKHwSMp55fjd7K4UXoTbRB7wv8JvR2EKvR+9KMJffuAN6M3nPFfO/96UQle6rcgl4Zfhp6D5bD0BuV23gu8Er0xlPnJHHfVqgEIoIRp0sRhaGUmhGRx4HfKKX+dyc2EXkj8HWl1BXWtf+LJqDDlVLfTy7fp5S6NPVOW7+zS0T+BvgQ8Aal1OMiMqOjqayd316M3vfmUKXUVPL91wHfFpHjlVJ3J/Fa6K0THk7i/Evy7IaiZRGRj0gyEVXiBcAZIvKI495z0Rtagd7VrQMi8hLgcvQmTSPotrmviPyBUuo3gd8/EpgyBAOglPpOojA+EjAks8sQTIKfktpGNKI6RJKJqBL7AzvQO+ql8VPr/1/bNxJ9zqfQu8ptQG/4dAqwDdgHvblTlfhd6rciqg5qQySZiLJ4HL2VpI2vAX+OlhR+X+Bdx6E7+VuVPj8IETk34HtpfBdYIiJLrOnS84CnMLd/bkSPEdk7oix2ASeIyHNE5GmJ0vWDaKXrx0XkeBF5roicJSI3ikgWQexEH4p2sYgcJiKvZU4hbH9vfxF5cfK9P3C85zb0/rkfFZE/Ts6G/jDwRdV5AF1EDxFJJqIs3o8+luM76H1qn62U+gnaYtQG/h3d4T+A3nx7r+9FSqlvAG9BT7PuAc5H62fsOHegFcHbk+9dmnoNSu8l+wr0htr/gSadH6I33Y5oCHGP34iIiFoRJZmIiIhaEUkmIiKiVkSSiYiIqBWRZCIiImpFJJmIiIhaEUkmIiKiVkSSiYiIqBWRZCIiImpFJJmIiIhaEUkmIiKiVkSSiYiIqBWRZCIiImrF/wdhslO9t4kCbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 300x300 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEKCAYAAAC2QVjtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaQklEQVR4nO3dfbBcdX3H8fd3V001JLEVBWUTqSUyKiqWXB7iQ9IgjdK0lE4DGaRD6g0Pqe2gkoBh1FpFHkKMpJLRBG8n4ihDbDUdfAIfAmpIIPFhIio2WB6yUYmC5iaOoN777R+/3d69e/d5z+45Z8/nNXNmuWfPOfnlkHzz+31/T+buiIjEKRd3AUREFIhEJHYKRCISOwUiEYmdApGIxE6BSERip0AkIrFTIBKR2D0j7gLEwcwMeBFwOO6yiAy4GcBPvcnI6UwGIkIQKsZdCJGMKAAHGl2Q1UB0GGD//v3MnDkz7rKIDKTR0VFmz54NLbQ8shqIAJg5c6YCkUgCKFktIrFTIBKR2CkQiUjsFIiaKBZh+/bwKSK9oUDUwMgIvPjFsGhR+BwZibtEIoPJsrhCo5nNBA4dOnSobq9ZsRiCz/j4xLl8Hh55BAqFvhRTJNVGR0eZNWsWwCx3H210rWpEdezbNzkIAYyNwUMPxVMekUGmQFTH3LmQq3o7+TyccEI85REZZApEdRQKsHlzCD4QPjdtUrNMpBdiD0RmttLM9prZaOnYaWZvbnLPUjN70MyeMrPvm9nZvSjb4sXw6U/D1q0hNzQ83ItfRURiD0SEyafvAk4B5gFfB/7bzF5R62Izmw/cBowArwG2AdvM7KQoC1XuMTv/fFi2DO68M8qni0ilRPaamdmTwGp3n9Jhbma3A9PdfUnFuV3A99z9shaf37DXTD1mIt1Lba+ZmeXNbBkwHdhZ57IzgK9WnbuzdL7ec6eZ2czyQVgjpS71mIn0VyICkZm90syOAE8DHwPOdfcf1rn8WODxqnOPl87XswY4VHE0HCetHjOR/kpEIAJ+DJwMnAZ8FPiEmb08wudfB8yqOBo2sDrpMdNUEJHOJWI9Inf/HVBu+HzbzIaAy4FLa1z+c+CYqnPHlM7Xe/7ThNoWAGGl2MaGh0Ov2UMPhZpQoyA0MgKXXBKac7lcCGLqYRNpXVJqRNVywLQ63+0Ezqw6dxb1c0odKxRg4cLmNaFyEILweemlqhmJtCP2GpGZXQd8CXiMkES+AFgILC59fytwwN3XlG7ZANxjZlcAXwCWEbr9L+lvyYNGiW31sIm0JvZABLwAuBV4ISGRvBdY7O5fKX0/B/j/v+rufq+ZXQBcA1wL7AP+1t0f6GupS8qJ7equfiW2RVqXyHFEvdbK7Pt2jIyE5tjY2ERiWzkiybp2xhEpEEW0eH6x2FpiWyQr2glESWiaDYRCQQFIpFNJ7TXLFI1BkqxTIIqZlqMVUY4o1g0WNblWBllqJ71mQWUzTJNrRQIFoj6qboZ97WtQPdtEY5Aki9Q061PTrFYzrJrGIMkgUfd9AtVqhlXK5WDnThgammi2zZ2rXJFkg5pmfVJrjaNK4+PwrW/BunXqRZPsUdOsj71mlcuFtEq9aJJW6jVLqOFh2LVraoK6EfWiSRYoEPXZ0BDccsvE6o/NgpJ60SQLFIhiMDwcmlvbt8N99zXOHY2NwW239a1oIrFQIIpJefXHoaHJ62PXCkpXXhmS2CKDSsnqGKd4VCovI3LwYNjUsVouB48+GgKYuvclDZSsTqFyDWn+/Np5o/HxEKgqR2fPnq31sWUwKBAlTKEAN9ww9Xw+D9OnT+3+37wZ5szReCNJNwWiBFq9Gm68cSJflM/DO94Be/bUHoPkrpqRpJtyRAnJEdVSLMKGDfChD4Vg08zWrbB0ae/LJdIK5YgGyPr1rQUhCEluNdEkjRSIEqzZRNlqaqJJWikQJVizibK1aEqIpJECUYIVCpMHO7ZCU0IkjRSIEq5yOsj994fPM8+sf/2mTfUHOWq3EEkqLYyWAtV7pk2fDqeeOvW6O+6AJUtqP6NyCZJcLtS0tBKkJIVqRCk0NAQXXTT53EUX1Q9CxeLkgZDj40pqS7KoRpRSW7bA294GO3bAa18bglM9jXYL0Vw1SQIFohQbGqodgKonxZZ736r3T1NSW5JCTbMBU71l0bp1E/PXKqeMNEpqi/SbpngkeIpHu+ptWXTGGWGHEAgz+2+4IcxnE+klTfHIqHojsctBCMLo6zVrlKiWZFGOaIDUygXVUjn6et8+OOooOHJEC61JfFQjGiD11jKq5bOfDesYLVoUxiRpHzWJk3JEA5QjKnv3u+GDH+zs3solaUW6oRxRxl1zDaxd2/6EWQjNug0boi+TSCOqEQ1gjaisWAyJ6vPOa+++fD7c9/DD8MQT8LznhbW0VUuSdrRTI1KyeoAVCnD00e3fNzY2dS6bWdgYUvPTpBfUNBtwnaxpVIt7mK+mbn/pBQWiAdfJmkb1lLc0EomaAlEGDA9PHtTYqVxO89OkNxSIMuLhh7t/xvXXK2EtvRF7IDKzNWa228wOm9lBM9tmZie2cN/bzezHZvZbM9tvZh82sz/qR5mzasaMuEsggyoJvWYLgI3AbkJ5rgXuMrOXu/tvat1gZhcA1wNvBe4FXgpsARx4Zx/KnDrlray7Ga2xciX88Idw/PHw+tc3XgNJpB2JG0dkZs8HDgIL3P0bda65GXiZu59Zce5DwGnu/roWfo1MjCOqVrlcLEwEplbmp9Vy0UVhgTaRWtI+snpW6fPJBtfcC5xiZqcCmNlLgLOBL/a4bKk2PBymb2zfDvv3w2OPhd1hO/236BOfgN27oy2jZFMSmmb/z8xywE3ADnd/oN517v5pMzsa+JaZGeH38TF3v7bOc6cB0ypOZTbbUb0Q/9FHd9dc27FDTTTpXtJqRBuBk4BljS4ys4XA1cA/AX8O/B3wV2b2njq3rAEOVRwallfS7YDHvXu1RZF0LzE5olLe5xzgDe7esLPZzL4J7HL31RXnLgQ2A0e5+3jV9bVqRMWs5YjqGRkJu3qMjXX+DG1RJNVSlSOy4GbgXGBRsyBU8hygOr1a/mtk1Re7+9PuPlo+gMNdFXrAlDdxXL++82doiyLpRuyBiNAcuxC4ADhsZseWjmeXLzCzW83suop77gBWmtkyM/tTMzsL+ABwh7t38e96dhUKsHRpd820sTG49lolsKV9SQhEKwk9ZXcDP6s4zq+4Zg7wwoqfrwE+VPr8ITAC3Alc2vviDq7qeWk2pW7Z3Ec/GmbuL18eadFkwCUmR9RPWR1H1KpiMUxuPXIE/vqvO3/O/ferRy3LtB6RdKWyi3/BArjnns6es3mzFuWX1qhGpBpRQ8UizJ7d/XM+/nH1qGVNqnrNJNkKhRBEul3PaMUK9ahJfQpE0lS5e3/r1s4S2GWf/GRkRZIBo0AkLSl3799yS+fPuPpq7ZsmtSkQSVuGh+Hyyzu/f8UKeO97J481KhY1TSTrFIikbW95S3f3f+ADYazReeeFGtKLX6ydZrNOvWbqNevI8uVhGZCyo4+GX/6y++fm8yEfpe7+9FOvmfTcli1hwOKHPxw+v/vdaJ47NqadQrJIAxqlY0NDEyOnt2+P7rnTp0f3LEkH1YgkElFt5AihaSbZokAkkShPmI0iGJ1/vpLWWaNAJJEZHoZdu7ob9Ahh6drKtY127w5rJWl5kcGlHJFEamgoDHpcsaK754yNwYYN8ItfTO6d084hg0nd9+q+74nPf767JUQa0fIi6aDue4ldL3u+duzo3bMlHgpE0hNz5/bu2Sec0LtnSzwUiKQnCgVYu7Y3z969O6wEoLlpg0M5IuWIemrdOli9uvl1nTALiXEtuJZM7eSIIg9EZnaau98X6UMjpkDUX8Ui7NwZpm5cfXW0z9bctOSKe83qzxB23RABJtYyAnjBC+Dii7vb5rpSeW6aAlG6dRSIzGxrva+AP+m8ODLohodh8eIQPKZPh9NO6z4o7dkDCxdGUjyJSUdNMzN7EvgH4Ej1V8Dt7n5MBGXrGTXNkqFYhDlzug9Eap4lUz+aZncDh939G9VfmNneDp8pGbNvXzRNtLEx+MxnQvNPwSidOu2+/xfgJ7W+cPezOi+OZEmUM/bf+c5Qu1q3LprnSX+19cfAzF5rZg8DjwGPmdnjZnZDqakj0pYoZ+xDqF2tXg033hjN86R/2v0jsAn4ETAEnAisBt4IfMfMjou4bJIB5Rn7UbrqKg14TJt2A9GfAW939++4+0PufiswD/gucFPkpZNMOFLd5dEl97CmkRbjT492A9GPgBdUnvDQ7fZe4E1RFUqyJcpcUaXxcbjkEq1jlAbt/u/fAnzEzKp3Q58FNOyeE6mnnCsqb2vd7cJqlcbHw1gl1YySra1xRGY2XvrP3wGfBb4H5IELgWvd/VORl7AHNI4omYrFMNCxPLt+505YtiwEk26ZwWOPqXu/n3o5juiFwMnAq0ufy4G5gANXmtmbgb3AXnf/cpvPlowrFCYHiqVL4dFHo5k06x4CW3mqiSRLW4HI3R8H7iwdAJjZHwGvZCJA/Q1wNfDc6IopWXXKKdE96z/+Q4Eoqbqe9OruTwG7S4dIpMqJ7CiaZ1/+Mpx+Ovznf6qJljRaGE0SrTqRnc/Dxz/eeU3pvvtg9myNwE4aLYymZHUqVCayC4Xw8/ve111v2Nq1vVu0TeJfj0gkcpWJ7JGRMD6o2+balVfCUUfBvHlhUOXcuWqyxUU1ItWIUqVYDCOmo8gZVcvlQjNQS89GQ9sJycDat692EIpiEOT4+OQdZqV/FIgkVWpNB8nnQxJ61arunz82FsYbSX8pEEmq1OpF27Qp7Px6+eXR/Brnn68pIf2mHJFyRKlU3YsGsH07LFoUzfO1/Gz3lCOSgVcohAXzKwNFlLP41UTrr9gDkZmtMbPdZnbYzA6a2TYzO7GF+55rZhvN7Gdm9rSZ/Y+Znd2PMksyRb3io5po/RN7IAIWABuB04GzgGcCd5nZ9Ho3mNmzgK8AxwN/T1gt8mLgQK8LK8lWXvEximDkrl60fol9QKO7T1pQzcyWAweBU4Apu4SUvJWwf9p8d/996dwjPSqipMzQUKgZXXppaGJ1Qxs49kcSakTVZpU+n2xwzd8AO4GNpQX8HzCzq80s3/viSRoMD4dk8/r13T0nl5tYH6lYDAlx1ZCil6hAZGY5wtrXO9z9gQaXvoTQJMsDZwMfAK4A3l3nudPMbGb5AGZEW3JJovJW19000049NTxnZCSM6F60SGth90KiAhEhV3QSsKzJdTlC8+0Sd/+2u98OfBC4rM71a4BDFYf+TcuIegnsefNau3/XLtiyBVasmBjRrRHY0UvMOCIzuxk4B3iDuz/c5Np7gN+7+xsrzr0Z+CIwzd1/V3X9NGBaxakZQFHjiLKjWJzojj/++DDJde3asEZRp1at0h5qjaRqHJEFNwPnAouaBaGSHcAJpaZc2UuBn1UHIQB3f9rdR8sHcDiSwktqlJtpo6NhcbRFi+Cuu7p75vr1qhVFJfZARGiOXQhcABw2s2NLx7PLF5jZrWZ2XcU9HyX0mm0ws5ea2V8Rlqfd2M+CS7oUi5OXDxkf726y7Ph46FGT7iUhEK0k9JTdDfys4ji/4po5hIX7AXD3/cBiwo6ze4F/BzYA1/elxJJKtWbuu8O553b2vHx+okdNupOEcURN/01y94U1zu0kDIIUaUmt9a/zeTjrLPjc59p/3ooVcO+9MH++xhl1Kwk1IpG+qDdzv1ObNoVpILNnqzu/WwpEkinlgY7bt4fP4WF43vO6f+6KFUpcd0OBSDKneub+/PnRrPCo2fqdUyCSzCsU4JZbJppsuRycd177z1EPWucSM6Cxn7QwmtRSudgatL9IvxlccUVYKVLJ65QNaBRJisomW3ViuxXuYePG2bM14rpdCkQidVQmtu+/v7080pVXwmWXKYHdKgUikQbKtaShoZBHamcm/6ZNMGeOuvZboRyRckTShvLk2SeeCD+vXNn8nlwObrstewMf28kRKRApEEmHikXYsCHkhVpRuZNssRimnAzyNtdKVov0WHmhtHXrQu7oVa9qfk95HaMbb9Qia9VUI1KNSNpULLbftV/JLPSwlQ3qHmqqEYn0UK1Z/BACTLOetVxuchCCiQX6s0yBSKRN9TZydIf3vKf+ffk83HDD1Hu1nIgCkUjbCoUQUKrl87BkydRAk8vB1q2h+bVq1dQ1tMfH4c47e1rkxFMgEunAqlVhzetyQCkvKVLeU61yqZHNm8MyteUc0KteNbl5po0claxWslq6Ujk/rTLZXH2+3F2/Zw9cddXUPBGEEdwLF/at6D3XTrI69hUaRdKsPC+t0fmRkclrZdeS9TyRmmYiPVS9YH89Wc8TKRCJ9FC9rv5qtfJEWdriWoFIpIfqdfXXUjmeKGtbXCsQifRQrQX7b7wxdOdXD34s54lq7b826L1qSlaL9NjwMCxePLV3bXQ0BJixsYnu/0IhNMeqm3Pl2tKgTQMpUyAS6YNavWuVAWr6dDhyJNR69uyZev+g96opEInEqFAIvWXlplituWgA118/uLUhUCASiVWtfFAt8+b1r0xxULJaJEatdO8PerMMFIhEYlWrez+XmzqHbZCbZaBAJBKrWt37mzfDo49O3hZ70GnSqya9SgLUmzybZpr0KpIy9SbPNjMoi/CraSaSUoM0DURNMzXNJIWKxbB5Y+Vf31wu5JaSUjPS4vkiA+7ee6cOfBwfD7vRpnFOmgKRSAp9/eu1z7///elspqlppqaZpEytZlm1JOyVpqaZyADbt69xEIL07ZWmQCSSMq0stpbPhxn9aVnhUYFIJGVqjca+6KLJP194IZx+enq69pUjUo5IUqrWlkXltY1OP33yZNo4ckYaWS2SAdWjscs/p3GFRzXNRAZMrRxS0pcSUSASGTC1ckhJX0ok9kBkZmvMbLeZHTazg2a2zcxObOP+ZWbmZratl+UUSZPh4ZATSstSIknIES0ANgK7CeW5FrjLzF7u7r9pdKOZHQ+sA77Z4zKKpE6nM/rjEHsgcvc3Vf5sZsuBg8ApwDfq3WdmeeBTwL8Crwee27tSikgvxd40q2FW6fPJJte9Fzjo7k1HSJjZNDObWT6AGd0WUkSik6hAZGY54CZgh7s/0OC61wHDwMUtPnoNcKjiSMFYU5HsSFQgIuSKTgKW1bvAzGYAnwQudvdftvjc6wg1rfKRkpazSDbEniMqM7ObgSXAG9y9UY3lz4DjgTtsYvPwXOkZfwBOdPefVN7g7k8DT1f8WtEVXES6FnsgshAVPgKcCyx094eb3PIg8Mqqc9cQ8j6XA/sjL6SI9FTsgYjQHLsAOAc4bGbHls4fcvffApjZrcABd1/j7k8Bk/JHZvZrgEZ5JRFJriQEopWlz7urzv8jsKX033OAJvthikhaxR6I3L1pwsbdFzb5fnlU5RGR/ktar5mIZJACkYjEToFIRGKnQCQisVMgEpHYKRCJSOwUiESkLcVi9NsUKRCJSMtGRsL2RFFvU6TthLSdkEhLisUQfFrdpkhbTotI5Pbtq79NUbcUiESkJb3cpkiBSERa0sttipQjUo5IpC3VW13Xoy2nRaRnerFNkZpmIhI7BSIRiZ0CkYjELtM5otHRhvkzEelCO3+/stprdhzaZFGkXwrufqDRBVkNRAa8CDgc0SNnEAJbIcJnDiK9p9YM0nuaAfzUmwSaTDbNSi+lYYRuR8WGjYebjZfIMr2n1gzYe2qp/EpWi0jsFIhEJHYKRNF4Gvi30qfUp/fUmsy9p0wmq0UkWVQjEpHYKRCJSOwUiEQkdgpETZjZGjPbbWaHzeygmW0zsxNbuG+pmT1oZk+Z2ffN7Ox+lDcunbwnM3uFmf2XmT1iZm5mb+9XeePS4Xu62My+aWa/Kh1fNbNT+1XmflAgam4BsBE4HTgLeCZwl5lNr3eDmc0HbgNGgNcA24BtZnZS74sbm7bfE/Ac4H+BdwE/73kJk6GT97SQ8OfpL4AzgP2le47rbVH7R71mbTKz5wMHgQXu/o0619wOTHf3JRXndgHfc/fL+lPSeLXynqqufwS4yd1v6nXZkqTd91S6Jw/8Cvhnd7+1l+XrF9WI2jer9Plkg2vOAL5ade7O0vmsaOU9SWfv6TmEmtTAvFsFojaYWQ64Cdjh7g80uPRY4PGqc4+Xzg+8Nt5TpnXxnm4AfsrUf+xSK5OTXruwETgJeF3cBUk4vafWtP2ezOxdwDJgobs/1auC9ZsCUYvM7GZgCfAGd2+2ltHPgWOqzh1DBhKybb6nzOrkPZnZKkJi/43uvreX5es3Nc2asOBm4Fxgkbs/3MJtO4Ezq86dVTo/kDp8T5nT6XsysyuB9wBvcvc9vSxjHFQjam4jcAFwDnDYzMp5nkPu/lsAM7sVOODua0rfbQDuMbMrgC8QqtLzgEv6WvL+avs9mdmzgJeXrnsWcJyZnQwccfcINjJOpE7e01XA+0v3PVJxzxF3P9LX0veKu+tocABe51hecc3dwJaq+5YCPybMoH4AODvu30vS3hNwfJ177o7795Ow9/RInXveF/fvJ6pD44hEJHbKEYlI7BSIRCR2CkQiEjsFIhGJnQKRiMROgUhEYqdAJCKxUyASkdgpEIlI7BSIRCR2CkSSOGZ2bGkx/cvN7LulDQh+YGZa32hAKRBJEp1c+nwr8PbSz48BnyqtaigDRv9TJYleDfweOMfd73H3B4F3A3MIS4UsMbMfm9k+M1sRa0klElqPSJLoZOCz7v5IxbnRiv9eT9ha5xDwbTP7nLs/0cfyScRUI5IkOhn4XtW5M4BfArOBH7j7AQ+Lgn0J+Ms+l08ipkAkiWJmzwbmAvmKczlCrugTwIuAAxW3HAAGZqPBrFIgkqR5JWH1wQvN7AwzexlwO/Bc4JpYSyY9o0AkSXMy8CBwLfBfwB5C7WiBu/+asJ9XZQ3ouNI5STEtFSuJYmYbgT929wvqfP8M4EeE/eAPAd8G5itZnW7qNZOkORm4o96X7v6H0u4o2wk1+rUKQumnGpEkhpkZoZazzN2/GHd5pH8UiEQkdkpWi0jsFIhEJHYKRCISOwUiEYmdApGIxE6BSERip0AkIrFTIBKR2CkQiUjsFIhEJHYKRCISu/8DjkkLWZ/gAJoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "settings=dict(fig=dict(figsize=(3, 3)))\n",
    "mcp.plot_chain_panel(chain[burnin:, :], settings=settings)\n",
    "mcp.plot_pairwise_correlation_panel(chain[burnin:, :], settings=settings);"
   ]
  },
  {
   "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": 27,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------\n",
      "Display log file: resources/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: resources/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
}