{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "measured-strengthening",
   "metadata": {},
   "source": [
    "# MCMC in EasyVVUQ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "chief-morocco",
   "metadata": {},
   "source": [
    "EasyVVUQ provides support for MCMC sampling with multiple chains in parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mechanical-deposit",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import easyvvuq as uq\n",
    "import numpy as np\n",
    "import chaospy as cp\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "leading-aging",
   "metadata": {},
   "source": [
    "We define a Rosenbrock function in 2 dimensions for testing purposes. This will be a stand-in for our probability density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "allied-earthquake",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rosenbrock(directory):\n",
    "    json_input = os.path.join(directory, 'input.json')\n",
    "    if not os.path.isfile(json_input):\n",
    "        sys.exit(json_input + \" does not exist.\")\n",
    "    with open(json_input, \"r\") as fd:\n",
    "        inputs = json.load(fd)\n",
    "    x1 = float(inputs['x1'])\n",
    "    x2 = float(inputs['x2'])\n",
    "    output_filename = os.path.join(directory, inputs['outfile'])\n",
    "    y = (1.0 - x1) ** 2 + 100.0 * (x2 - x1 ** 2) ** 2\n",
    "    with open(output_filename, 'w') as fd:\n",
    "        json.dump({'value': -y}, fd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "royal-camel",
   "metadata": {},
   "source": [
    "Next we define a helper function to create a campaign, sample the search space and return the corresponding DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "robust-intensity",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mcmc(tmp_path='.'):\n",
    "    campaign = uq.Campaign(name=\"mcmc\", work_dir=tmp_path)\n",
    "    params = {\n",
    "        \"x1\": {\"type\": \"float\", \"default\": 0.0},\n",
    "        \"x2\": {\"type\": \"float\", \"default\": 0.0},\n",
    "        \"out_file\": {\"type\": \"string\", \"default\": \"output.json\"},\n",
    "        \"chain_id\": {\"type\": \"integer\", \"default\": 0}\n",
    "    }\n",
    "    encoder = uq.encoders.GenericEncoder(\n",
    "        template_fname=os.path.abspath(\"rosenbrock.template\"), delimiter=\"$\", target_filename=\"input.json\")\n",
    "    decoder = uq.decoders.JSONDecoder(\"output.json\", [\"value\"])\n",
    "    campaign.add_app(name=\"mcmc\", params=params, encoder=encoder, decoder=decoder)\n",
    "    vary_init = {\n",
    "        \"x1\": [-1.0, 0.0, 1.0, 0.5, 0.1],\n",
    "        \"x2\": [1.0, 0.0, 0.5, 1.0, 0.2]\n",
    "    }\n",
    "    def q(x, b=1):\n",
    "        return cp.J(cp.Normal(x['x1'], b), cp.Normal(x['x2'], b))\n",
    "    sampler = uq.sampling.MCMCSampler(vary_init, q, 'value', n_chains=5)\n",
    "    campaign.set_sampler(sampler)\n",
    "    action = uq.actions.ExecutePython(rosenbrock)\n",
    "    iterator = campaign.iterate(action, mark_invalid=True)\n",
    "    for _ in range(1000):\n",
    "        next(iterator).start()\n",
    "    df = campaign.get_collation_result()\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "amino-charger",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = mcmc()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "included-announcement",
   "metadata": {},
   "source": [
    "Finally we plot the the five different chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "legislative-brief",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x124a943a0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABAoElEQVR4nO3daXRc553f+e9z7629UFVYCzvBBSTFfZWoXaa10JbUXmR3293udKY7cTInOadzJjOZybvJvMqrTOacSc6Jx92xO922u9st2dbSomVtpChSEinu+waA2HfUvtx7n3lREEiKpMQFYBWA/+ccHBZQl1X/Wyz+6sFzn0VprRFCCFG5jHIXIIQQ4otJUAshRIWToBZCiAonQS2EEBVOgloIISqcNRcPWldXpzs6OubioYUQYkE6dOjQqNa6/mb3zUlQd3R0cPDgwbl4aCGEWJCUUt23uk+6PoQQosJJUAshRIWToBZCiAonQS2EEBVOgloIISqcBLUQQlS42xqep5TqApKAA9ha621zWZQQQoir7mQc9Ve01qNzVokQi5jraoa7EpiWgT/swRe08PhMlFLlLk1UgDmZ8CKEuDNjfSl6To5d9zPTY+ALevCHPPjDHvwhC3/Igy/kweM1y1SpKIfbDWoN/FYppYH/prX+0RzWJMSiM9KTxB/2sPLBRnKpIrlMkVyqSD5jk57MMz6Qhms2+bC85nRol8L7sy9fyMLySIgvNLcb1I9prfuUUg3AW0qpM1rrPdceoJT6IfBDgPb29lkuU4iFK5sqkBrP0fZAzUzgfp7ruOSz9kx451JFcukiybEcY72p6461fOYN4f3ZbdOS8QPz0W0Ftda6b/rPYaXUK8CDwJ7PHfMj4EcA27Ztk/29hLhNIz1JlKGobQ3f8hjDNAiEvQTC3hvucxyXfNomly6ST5cCPJcuMjWSZfRK8rpjPf7PQvtqN8pnYW6aEuKV6kuDWikVAgytdXL69rPA/zXnlQmxCLiuZrQ3RSwexOu/u0tGpmkQjHgJRm4S4rY7E9yfhXkuXWRiKIOdd64eqBRev1nqCw96Zi5olv70YBhyUbOcbuedEQdemb76bAE/01q/OadVCbFITA6WArO+rWpOHt+0DEJRH6Go74b77KJzXXjn0qV+8fGBNHbh+hD3BazrwnsmzAMWSkJ8zn1pUGutLwEb70MtQiw6I1eSeAMW0frAfX9uy2NixUxCsRtDvFhwrnajpIrkpgM9NZHDKbozxylDlcL7mm6Uz7pVvAFLhhfOEhmeJ0SZ5DOlfuSWzljFtUo9XhOP1yRc7b/u51prinnnuguan/WNJ0ZzuM7VEDdMNTO88POjUzx+GSN+JySohSiT0Sul0Rp1bbe+iFhplFJ4/RZev0VVzU1CPOdc15XyWdfK1EgG17k6xsAwDfwh65pWuEcm+nwBCWohykC7mpErSaL1AXzBG4fjzUdKKbwBC2/AIlJ3fVeOdjWFnD3ThfJZKzybLDA5lEG7V0P8uok+oek+8UU+0UeCWogymBrJUsjatK+pLXcp90WpL7s0guTz/fHa1TNjxG9nos/MBc1FNNFHglqIMhi5ksTymcQag+UupeyUoW490cfV5K8J78/CPDmWY6zvcyG+gCf6SFALcZ8VcjaTQxniSyMyPvlLGIa6w4k+9hdM9LGuCfH5NdFHglqI+2ysN4V2NfXtczN2erGY1Yk+n7ugWWkTfSSohbiPtC5dRAzX+G/aShSz404n+uTTdkVP9JGgFuI+So7nyKWKLN0UK3cpi9ZtT/RJf7YIVvkn+khQC3EfjfakMD0GNU2hcpcibuK2J/pcc4Hz2ok+ltdky3NLZr0uCWoh7hO76DA+kKKutWrej0JYbG53oo99Tat7NklQC3EfTI1kGbw4ieto6pfIRcSF5NqJPnNFglqIOaa15uyBAQCCt7jAJcQXkd+/hJhjQ5cTM7dlSJ64GxLUQsyxazetrW2Ri4jizklQCzGHrh3SVddataDXoxBzR4JaiDnUfU1ruq59/ixnKiqLBLUQc+izNSe8gRuHdQlxuySohZgj2VRh5na8IyKL4Yu7JkEtxBy5+OnIzO3aVun2EHdPglqIOaC1JjOVByDWGMTrlykL4u5JUAsxB6aGszO3G5ZEyliJWAgkqIWYA+cPDs3cjn5u/0Ah7pQEtRCzzLHdmc1am1dW39d1i8XCJEEtxCwb7r66DVR9m1xEFPdOglqIWXblVGmSizdg4QveuGGrEHdKglqIWZTPFGdut6+pLWMlYiG57aBWSplKqcNKqdfmsiAh5rO+c5Mzt2ONwfIVIhaUO2lR/zlweq4KEWK+01rPTBmvaQ5XzA7WYv67raBWSrUCzwM/nttyhJi/kuO5mdstq2LlK0QsOLfbov7PwL8DbrkhmFLqh0qpg0qpgyMjI7c6TIgFq+fk+MztQNhbxkrEQvOlQa2UegEY1lof+qLjtNY/0lpv01pvq6+vn7UChZgPHNudmTLeJhcRxSy7nRb1o8DvKaW6gF8AO5VSfz2nVQkxz4wPpGduN3TIdltidn1pUGut/73WulVr3QF8D3hHa/2DOa9MiHnk8pFSd5/pMTBNGfUqZpe8o4S4R9eOne7cHi9jJWKhuqO1F7XW7wHvzUklQsxTIz1Xp4xHagMUnAIewyMbBYhZIy1qIe6B1pr+85MARBuCpAopfnz8x3w8+HF5CxMLigS1EPfg2rHTyzbV88qFVwAIe2QxJjF7JKiFuAeDF6dmbhseSBZK3SBr69aWqySxAElQC3GXHNtlcigDQOvqat7rfQ+AFbEVZaxKLEQS1ELcpYnBq2OnG5dFOTt+FoCd7TvLVZJYoCSohbhLl4+Oztw+O1kK6aAVxDJkI1sxuySohbgL+UxxZrstgHd63gXgW53fKldJYgGToBbiLlw7dtqqdfhsyHTUFy1TRWIhk6AW4g5dO3Ya4EDhfQB2dewqU0VioZOgFuIOXTt2uugW0ZECAEujS8tVkljgJKiFuEP912y31Rs8jzJgU/0mmTIu5owEtRB3wLFdEqNZAFztMhrsA2B70/ZyliUWOAlqIe7AtetODxUGIejQEGzAY3jKWJVY6CSohbgDn607rYHe8HmUgmc7ni1vUWLBk6AW4jZdu+70ZG4CaksXFSPeSLlKEouEBLUQt2ngmgWYLjhnUB7Ni8teLGNFYrGQoBbiNmitGe5KAJCxM9BQuqDYWtVazrLEIiFBLcRtuHbs9KXBM9R0nWdH3TYZkifuC1k9Rojb8NkCTAW3SPjiCcK+QdZGV5W5KrFYSItaiC/h2C75dOlCYn/3KYpOhJivGl+VrOsh7g8JaiG+xMiV0gJMjutQvHwZgNX/8n8pZ0likZGgFuJL9JwYA2Dk4kkAfOEENQ3t5SxJLDIS1EJ8gc/GTruOQ/LKJQA2/4vvl7MksQhJUAvxBbpPllrTE6ePAaBrbZa1rCxnSWIRkqAW4ha01kwOZtDFIuPD3QCs/tYzMiRP3HcS1ELcwtRIaVLL5NFPAcjXBdm8ekM5SxKLlAS1ELdw7qNBnHSa0cQAAGuffhTTNMtclViMvjSolVJ+pdTHSqmjSqmTSqn/cD8KE6KcHNtFA8ljRwDINMZormopa01i8bqdFnUe2Km13ghsAnYppXbMaVVClFnf2Qmc8XFGsqVlTZcs28TAxcnyFiUWrS8Nal2Smv7WM/2l57QqIcps4OIkqbOnAUi31hLzx2jokOVMRXncVh+1UspUSh0BhoG3tNYf3eSYHyqlDiqlDo6MjMxymULcP7lUEXtwcKY13bFkPQDta2vLWZZYxG4rqLXWjtZ6E9AKPKiUWneTY36ktd6mtd5WX18/y2UKcf8cfesyqa6LAKSW1BHzxQjFfBiGDMsT5XFHoz601pPAu8CuOalGiDLTWlPo6WE0O4prmSxpWQPA8i0NZa5MLGa3M+qjXikVm74dAJ4BzsxxXUKURe/hKyQHewFIt9VS7a8BwB+SzWtF+dzOetRNwE+VUialYP87rfVrc1uWEOVx6e0TTOTGsQNeWuOdKKBze7zcZYlF7kuDWmt9DNh8H2oRoqzSF7pIjpUmtxSqQ9R6qwGIxYPlLEsImZkohC4UKPT2cuAn+5nKlzawjY/a2F3drHmsRdb2EGUnW3GJRUW7Ls74OMWhIeyhYeyhQeyxcbSrSRW8M8eFPCGMSBXhal8ZqxWiRIJaLFhaa9xkEntoiOLQMPbQEPbIMLpoA6D8PjzxOIGWFi7tO0mqWJrXFfSGMb1+try0sZzlCzFDglosGG4uNx3K063l4SHcTGkFPGWZmHV1+NeswYrHsRoaMGMxlFKk9u3jwoTCtUwM2yFkBfE0xgnEAmU+IyFKJKjFvKRtG3t09LpgdiYnZ+43a6rxLlkyHcpxrLpa1E1WvnOmphj86AAJXwBP2sVr+vB6/Kx5/oY5XUKUTUUFtZvLobxelCHXOMVVWmuciYnrW8tjo+C4ABihEFa8Af8Dq7EaG7Hq6zF8t9e3nNq3j3MTaRx/BE8qR9hTeqxYs6zrISpHRQX15N//Em/HEsKPP17uUkQZOak09vAQ9uBgqW95eBhdKACgvF6shgaCmzaVWsvxOGY4fFfPU+zrY/DkIQbDcXzpKUzlwe8N0v5Ip4z0EBWlooLajEXJnTxJcPt2DL+/3OWI+8AtFLCHpy/0TV/0c1PTizUaCquuHt+qlXg+C+VYbFZ+49Jak9i7l/PZPPmaML7BNCFvGE99PfGVslaNqCwVFdTeJUsodHWTO3WK4JYt5S5HzDLtONijY6XW8nQ3hjMxCbq0aq4ZjeJpbsYTbyi1luvqUJ65mbqdP3OGwe5TDNa04U8VUCiCniCRznZ8gYr6byFEhQV1Rwe8v4fssWMENm2Svup5TGuNOzU13XUxND00bgRtOwAYwQBWQxxfZ2eptdzQgBG4P6MsdKFA4sMPuGDkUcFGPF1dBD1BvPUNtKyT6eKi8lRUUJuRCGZtDc7YOIVLl/CtWFHuksRtcjOZ64bFFYeG0Lk8AMpjYdXX41+/Yaa1bFRVla0fOPPpYfqHLzLUvIxwzkA7DiFfCE9zE9WNMl1cVJ6KCmoodX9kx8bJHj0qQV2hdLGIPTJCcXBoprXsJJKlO5XCqq3Bt3w5VkMcT7wBs7a2Yn47clIppg5+xPmIS9hYhj3cQ9D04YvWEGutxjAro04hrlVxQe3r6CD76WGK/QMUh4bxxGUd4HK6Ycr18BD22Bi40/3KkSqsePxqa7m+HuX1fsmjlk/6ww+5krzC+OqVNGT9ZManCPlr8bS30bq6ptzlCXFTFRfUVlMTyudD5/Nkjx7B8+yz5S5p0Zp67XWKvb3oYhEA5fPhiTcQ3Lp1prVshEJlrvL2FYeGmDh5hEstFtWZFUyMXiFsePBbPsxQiFCscj9gxOJWcUGtDANvezv58+fJX7iA88ijmOH5EwYLiRmpwlzzwA1TrucjrTWpvXvpLg6RaVpL1TD4ugaJBurxrVxJrDE4b89NLHwVF9QA3qUd5M+fB8cld+IEoR0PlbukRSn8xBPlLmHWFC5cYKzrLENVLjVvOwy2D7I04OA1vZixGB3r6spdohC3VJFXTrzt7TDdusmdPIG27TJXJOYzbdsk9+3jsjlOIB3AxiF6fpCGUAwABXhl7LSoYBUZ1EYggKexNJ7VzWTJnztX5orEfJY9epTh4cuMNPgoDrWTDBl0xPJkCiGCmzexYWdbuUsU4gtVZFDD9OQXwAj4yR47hp6evSbEnXDTaZKffMzFSI7qKS8Ft0BoMkc8FKPomqzY0SIb14qKV7lBvWQJAGYshj0ySrGvv8wVifko/dHHDEz1MtYSJnM5TsHQdNbYJPIBQuseoFr2QxTzQMUGtVlXhxEOo3y+Uqv66JFylyTmGXt0lOSJI5xv0jQOeMg7eYLaQ3PYSyLvo351k0xwEfNCxb5LlVJ4O5ZQ7OvHt2oVhctdOFNT5S5LzBNaa1IffEBvfphEaw3Jc1VorVlZq0kWfJjxRuraqspdphC3pWKDGsC/ahW6WCztzKEge+x4uUsS80ThchdTl89xrsNLwwWTgl3aFCAe0kzm/ITam2TjWjFvVHRQe5qb8TQ3kTt3Dt+yZeROncKdXkBeiFvRjkN63z66jQkKrQ2kzxqY+STLPSPYrkHGCFO/skEmuIh5o6KDGiC4dStuMoUKBNCFAvnTp8tdkqhwuePHGR/u5mJnmKqjLnZ2koD2UONLM5X3YdXHqWu5u11hhCiHig9qz5IlWHW1FPv6sOINZI/KUD1xa24uR/rjT7gUTOPWNFC44GA4mlBxmGKwgfFsgOjyRvxhGZIn5o8vDWqlVJtS6l2l1Cml1Eml1J/fj8KueX4CW7bijE9gRqI4U1MUurruZwliHsl8/DGjU/10P1CD71ARN5vG51oE1BSJYhgnGKN1fWO5yxTijtxOi9oG/q3Weg2wA/hXSqk1c1vW9XydKzCjEZzJSYxwmOzRo/fz6cU8YU9MkD52lLP1RbRRje4q4CnahHSKlK+WsUyI2s44sQYZOy3mly8Naq31gNb60+nbSeA00DLXhV1LGQaBzVuwR0awaqopXuktrYksxDXSH+xjMD/KwKo6rCMOOp3Bpy18up8h6im4mmVPrSp3mULcsTvqo1ZKdQCbgY/mpJov4H9gNUYwiJvNoSyT7NFj97sEUcEKV66QvnSeM20KJ1mF1ZPGZzv4VZEelSXtKNo7YwSq5SKimH9uO6iVUmHgH4B/o7VO3OT+HyqlDiqlDo6MjMxmjaXHtywCmzdhj4xgVteQP3sGN5ud9ecR8492XdIffEAfE4wtqcdzWqEyOfzaIGL0kjdXgcqx5GHZ2k3MT7cV1EopD6WQ/hut9cs3O0Zr/SOt9Tat9bb6+vrZrHGGf9260u4vjo22HXInT87J84j5JXfqNMnBXk4t96KHQvi6JgkUHTymQY9dIOeLEfFp/MuWlrtUIe7K7Yz6UMBfAKe11v9p7ku6NcPrJbB+Hc7EJEZVmOzxE2jHKWdJoszcQoHMRwfoDqRJ1dRhXTYwcgW8GMSNK0z6V1BlTvL0sw+jPDIkT8xPt9OifhT4Y2CnUurI9NfX57iuWwps3IgyDXBc3FSK/MWL5SpFVIDswYNMTQ5zpjMIvWH8F4cJ2Q4+y6C/6CMXCtEWnqR63QPlLlWIu/al21porT+gtAlGRTCCQfxr1pA9cQLlscgePYp/5cpylyXKwJmaInP4COfrihTMOowBjVlw8SqDBqufT3QN4fAUD8TX4GltLXe5Qty1ip+ZeDOBzZuB0kUke3CI4uBgmSsS5ZDev5/xwgSXl0fgSpjAhWHCjo3P52MoHyMTNWnK9hBdsxFlzMu3uhDAPA1qMxLBt3IlOC6ADNVbhIr9/WTPneN0i8ZJRmG0iLeo8ZoGUWuCPq2pc/tobFtB6MHt5S5XiHsyL4MaILhly8zt/IXzOKlUGasR95PWmtTeDxgmQX9bNQwECF0cI+QU8PmDDKQ9GGqSSKjA+u/9zyivt9wlC3FP5m1QW7W1eJdOD7dyNbnjslb1YpE/e5bsYB/HlxrowSr0aBa/rfFZJobOMq6zGPoSlzYtwVddW+5yhbhn8zaoAYLbts7czp08iS4Wy1iNuB90oUB6/wH6/FkmqmvRIz5C3QkCTo4sPq5kMli6i4srm3ju0T8td7lCzIp5HdSexkY8LaVlR9xsjty5c2WuSMy1zOEjZKZGOb7Mgp4w9uUxwnZpLL3r5HBJYOo+sps6WV0nIz3EwjCvgxoguPVqX3X26FFZq3oBc1Ipsoc/pavGZrQYZexyltqJDF43R9BjklJpPO453n76Mb695quyg4tYMOZ9UHva27Gmp6w7Y+MUe3vLXJGYK+n9+xlNTfBqqMjUMZe27gEiTp6w1yJBAZxBMmGbluY61tevLne5QsyaeR/USqnrW9VHZK3qhSjVN8C5fZ/yCgmMXs2KSwNUafCZJq5ySdtZvPoi+x97mK0tKwh7ZZU8sXDM+6AG8C5fjhmLAVDo6sKemChvQWLWOK7mSM8Ev/vpr7iQGoXoGGvOO3iVF58y8VoGQ3YW5V5itLqWztZqtsQ3l7tsIWbVgghqZRgEtlz9zylD9eY/rTWXRlL89YFuPn7vEA3jl7Dil1h6wo92fWizHp9lMuHm0W4CU/dz5qG1/MvN/4S2qrZyly/ErFoQQQ3gX7UKIxQCSkP1ikNDZa5I3K3hZI6XP+3j10f6UYUsz/a8SZV7FmPAg63bcFWUamecgrLJahvDOUfaW8//+bX/jWp/dbnLF2LWLZigVpZFYNMmALTtkHjtdZxksrxFiTuSytv89uQgP/uoh+Fknmda8ny752/wTfTQh0mBDRS1SdBxSFFgyi2g7H5s08H46nepD4fKfQpCzIkvXT1vPvGvW0vm0EEMfwA3kyHx2mtEX3oJQ6YQV7Si43Koe4JD3RM4rmZzaxUPqxNYXQcYO3mFgaoGMhkfhbwPbZvknFG00ngwMNwLOP5mfvDN75T7NISYMwumRQ2fbSywAWdyktCOh7DHx0nu3o123XKXJm5Ca82p/gQ//bCL/RfHWFIb5E/W+Xgy8Rre/oMkLrjk42vpL2TI2Y2k8zbYozhoLMPAUefIeKLUbHoCn7Wg2hxCXGdBBTVAYOMGlMfCHhkh/ORTFLq6Se/dW+6yxOdcGc/ws4972H1ykKDX4rtbmnih6iLR0z8HO0cu/CBFt5a+zCC5YhOpnMZyc3gxMQ2F49Oo4jCGGeGFP/xBuU9HiDm14JohRiCAf+1asseOUfPQQwQ2byZ7+DBmLEZg48Zyl7foTaQLfHBhlAvDKar8FrvWNbK6Ko8680tIDkJ8LXbNVpJ//wqpYoqhdJpkrhmcHMFIEGMqA0qR8B3FU4ixfM2DWCHpmxYL24ILaoDApk1kjx0je/gwoccew5maJLX3A4xIBN9S2eC0HHJFhwOXxjjWO4VpKB5ZXsuW9iie/kNwaA+YHlj7LWyjjomf/4K8k+fU6GmmskvJF4tE41G8A1OkDYdkfZ7gcBHDU82Tf/pn5T41Iebcggxqs6oK/6pV5E6dIrhtG5FnnmHylV+R3P1bzJe+PTPlXMw9x9Uc7Z3ko0vj5G2Hdc1RHl5eS8hNwvFfwOQVqOuElbuwE1kmX36Z0ewol6cuk8yFSeYDeOMh/BM5JlUeUykc5xJFM8jOBx/FDATKfYpCzLkF10f9mcCWLWjbIXvsGMrrJfL88yifj6nXXsdJpctd3oKntebCcIq/2t/F+2dHiEd8/NFDS3j6gQZCo8fgk7+A1BCsfh7WvYSdyDL6yi85P3SKy1OXyeQdxnNNTMRNfI5NNpvCUIrhWi/BcYe2aD3L/uCPy32aQtwXCzaorZoafMuWkj12HLdQwAyHiL7wPDqfJ/Haa+hCodwlLljDiRy/PNTLq0f7MQ3FNze38K3NLdR78nDs7+Dcboi0wPZ/Bk0bKI6M0PWLn3D8yiEm85M4jslwrp7RpioMo0BgKEkOB8fV+D2XUZaH9Q9sw/D5yn2qQtwXCzaoAQJbt6LzeXInTgJg1dcT2fUc9ugoid++JcP2ZlkyV+TNE4P87OMextMFdq5u4AcPLWFpbRA1fAo++TFM9UDns7Dxe+CPkh3s5/BP/2/ODZ7AMiya/MsYmMoz1rIKLxlaJ83SynhAZukQvqEpfCpGzbe/V+azFeL+WZB91J/xxON42lrJHjlCYMN6lGXh7egg/MTjpN7fQ3rfh4Qff6zcZc57BdvlYPc4n3ZP4GrYtqSGbR3V+D0mFNKlFvTIWYi2wOoXIFgDQH/XCU799X+lmE3TGGokbNZzqO8yPWsfpDqZROVczGSeIqDI43VToExitSvwBaVvWiweCzqoAYJbtzL1q1+TO3OWwLq1AAQ2lCbFZI8cKQ3bW7+uzFXOT66rOTWQYP/FMVJ5m1WNVTy6vI5o0FM6YPQ8nH0D7DwsewraHgLDwHEdDp74LaOv/BK/Y7CqZjWu4+fUwBiHNzbSOgGZQpG2UYccpd1bfGsH6bWzrDDq2fAv/gzTkE0BxOKx4IPa09qK1dBA9tNP8a95AGWUentKw/amSO15HzNShXfJkjJXOr/0jGXYc36EkWSepqif5zc00RybbuUWc3DhdzB4HMINsPH7pT+BsewY73/6Mt639hO3qmmrbieddzk7mOBCkyKS9JPNZaiesMm7pT0wg5FPOd2wHp2/yJql7bRX2eU6bSHKYkH3UcP0xgLbtuJMTZG/cOHqzw2Dqueew6qpIfHmbuyxsTJWOX+Mpwv8+kgf//BpL3nb5evrm/iD7W1XQ3qiCw7+BQydgCWPwNZ/CuEGXO1yZPgIv97/l4TePsjK8FI6okvJ5F3ODibJUqAnojCTNp60jT9dCuOAuxfPA61c8oXw1f0eq+ur4NO/grGLZXsNhLjfFmyLWhcK5M6WNrtVlglAcvdvMfx+lGWBZaEsi/CTTzL58itMvvwy1d//PkYoJHvt3US2cHXCimUqHu+sY1NbDMuc/qx3inDpPeg9WOqD3vzHpT5pIFFI8Hb324x3n+OBg8N0RFfjMTyk8jZnBpPYjssny5LUj8YxHE1stIAD+N2P8K1p5/2aTjwKdm34Gr64F078Eo7/femCZHVHuV4SIe6bLw1qpdRfAi8Aw1rredOZWxweJvXeezf8fOrXv7np8TqXZ/y//wQA5fGUwt2yUKaF8lhgmiir9HNlmjNBr6zP3WdZYFrX3DZLj2ea131AXPcYRuX+YmM7bmnCyuVxirZmfWuEHctqCXqveetM9cGZ1yAzDq3bSv3RpgetNWfGz7Cvfx/ekQSPHC9SG+pAoUjnbY6M5jFdGK/S+Fwvhquo688CYOk+Bh7wYjXYJI0cG2oeZXNbLRgKNv0Aej+GqGwQIBaH22lR/wT4f4G/mttSZpenpQXf8mXkL16i6rlnMWMxJv/270Apoi++gHYcdNEGx0Y7Dvlz5yj2DwDgX7sGbdvgOGjbRttO6Tjbxs0U0HaxdF+x9HdL39/DUD9D3STores/LD73wXHzD4trPwQ8N/+wME3weL70t4bPJqzsPT/KVLbI0roQj3fWURu+ZuyyY0P3B9BzAHxVsOn7My3cTDHD+73vc3nqMu2ZIJtPW3itGADpgs2hSZdsQzORngtcqOsllFhBZDyPmt5Evr9ziELch+FfzhLPDv5w08MYn11AtLzQIaN1xOLxpUGttd6jlOq4D7XMKqUU4Z07KQ4PkzlwgNgf/AHhxx8jtfcDlNeLt6lp5tjusTSp2naWDl4ivW9fqUvk8cfv6Pm01mCXwrz0IVC8Pujt4nSoO1eDfvp4Zv6OjXbsaz4ESve5ufz0B4c983il4+7hw8E0Sh8OpnnDbwxTBZczI1nGcg71IT9PtVXTkAqjjveQsSz8Gzdi5CfgzKuQGoGmDbDiabBKIX5p6hLvX3mfglPgEaOTlmMXwS391pAp2Bya1AxtepiW/W/TU9VNKLEKf8omkC6N8Jhs6sKusXj8QJopj4v1rM243UWNuwLLWLC9dULc0qy965VSPwR+CNDe3j5bD3tPDL+fyHPPMfnyy6TefY+qnV8hc/AgmUOfEn3h+ZnjjvdNcX4oxabWZrasWUPm4CHMaBT/mjW3/VxKqVJL1eOZi1O5Ke261384XBv6n92+5kPg2qC/LvSnPzgymTzn+icZGk/hMzQbawLEg3nUSB+5fnvmw8FflYShg2D5Yf13Smt1AAWnwAd9H3Bm/Ax1gTqeD2zHeOuD0nMDiVyRwwlF346nqT95EOwxJquqCEw6RMdLk1rywQyWb4QdR4qMG3ES9iVSdQ309FyiLlBHXaDuvr2+QlSKWQtqrfWPgB8BbNu2Tc/W494rT1MToYceIr3/AJ62VvwbNpD56GPssTGs2loAvr6uiQ/8oxzqnmA4spSnmqZIvvsuRiSCt7W1zGdwa8owwOtF3eMONnnb4WDXBJ92T0AMti6pZmtHNb7pi7AzMuPok7+G/gMQfwA6nwNvEIC+VB/v9LxDqpBia3wrG4uNpN74x5mQHknmOZWz6Nuxk2ghTaTvIueifQQn1hMbzZfOB5vq9BGq+hySLhxb3szKZ75KgydBqpii2if7IYrFqXKvYs2iwNateNpaSe/di7elBeXxkDl0aOZ+w1A8sbKe5zc0MZpx+E31alK+EIk3/hF7YqKMlc8t19Uc753iJ/u6+PjyOJ3xMH/yaAePrKi7PqS1Lo3mOPgXqPwkau03YM03wRvEdm0+7PuQ31z4DYYy+Fbnt9jstMyEtEZzeTTFqbxF30M7iUZDBI98xHjxBEZuPbHRq2uueN1jVHkd2q0qrLYn2PjUn/O9tc+SttMsiSzBNMwbT0KIRWBRdPgppYg88wwTv/gFyfffx7d6FbmTJ3EeeggzGp05bmW8ipqQl9eO9vNGw0aeOL8P9ZtXqf7972IssOU0u0bT7D0/wmiqQEsswDc21dMY9d94YG4KzrwOE91QuxxWfa104RAYzY7yu+7fMZ4bZ23tWh5pfgTdP1ha9Mp2sF2X0wMJxn0R+h78CsrvJ3XuMAy8Rd6zlcjENSGtz9Lgz7HB18B4VtOz/ff5o9VxhjPD5OwcHZGO+/TKCFF5vrRFrZT6ObAfWKWU6lVKzcuV2o1QiKpnnsEZG8dNpkApMocP33BcXdjH9x5sp7W9gb1tWzh3eYiJ114v9fkuAKOpPK8c7uWVw30UHc0LG5r47rbWG0Naaxg4VlpIKdEPq3bB+u+CrwpXuxwaOsQvz/2SvJPn+WXP82Tbk+i+gZmQTuVtDnVPMOaP0vfQTlyvD7dYJHbsTVy9GbN49alMPUbEGmKFN0pRGbz3wv/K1hVNxIJeLicuYyiDtogMxROL1+2M+vj+/SjkfvC2txPcuoXMoU9Rfh/506cJbd+O8bmtnPwek9/b2MxHET/H8xkyBw+w1v9b4i98bd5OhknnbQ5cGuN43xRey+CJlfVsbI1enbByrXwKzr1ZWqsj1lZaMzpQ6h+eyk/xds/bDKYHWR5bzhOtTxCwAhS6u5l6/XW04zCUyNM1miYXraX/wadwPaU+9Lpzh8mnm0ApDH11xEqV1UWzJ0ykOsJrX/n3BAwv25aUnq9rqouWcAs+U5Y0FYvXouj6uFbwoYco9PVhDw4BkD12jNDDD99wnFKKHctqiUce5sNcmsJ7B1nlD7HsmSfvd8n3xHZcDl+Z5OPL49iOZmNbjB1Lawl4b9HfO3ymFNJOEVZ8FVq3g1KlHcPHTvFh/4copXh6ydN0xjpRSlHo6mLq9dexbYfLo2nGUgVysVr6t18NaTM1Tv7EeUCj9PXXmvPqQaINGc79wZ8z3p3mpS1xLNNgIjfBZH6S9XXr5/hVEqKyLbqgVqZJ5LnnmPjF36LzebLHjhPYsuWWi9AvrQsR+8NdfPCTDCffeI9xM8DWr2yv+Ja11pqzQ0n2XRgjkS2yrD7E45311IRuMUKkmIXzv4WhU1DVCA+8CKHSULh0Mc27V96lJ9FDa1UrO9t2EvaGAWZCOp0rcn4oSa7okq2uo3/bU+jpoYq54iiNu3+J1iE0Crg+qBvrDWr/5J/zj+czrG6sor22NJKkK9EFQEe0Y9ZfHyHmk0UX1ABmJELVzq+Q+Mc3S2uCnDhBcOvWWx5fHfLx7D/7Dh/9t7+m//XdjGgPO5/ccOPwtQryac8ke86NUF/l4ztbW2mrCd764LGLpeVICxlY+ji0PwzTIywuTFzg/d73cbTD4y2Ps65u3cyHVP7yZRKvv8FwIkvXaBpXc11I11d5Odv9ASv37KPo1KHIA9d/IIabczzxRIBfjQSwzBxPrLy6n2XXVBd1gTqqvFWz/voIMZ8syqAG8K1YgX/dWnInTpI9cpTAhg1fOFnF5/Pw6D//Hif+v/9B12938/eOxdcfXXXrFmqZrW2O4PcYPNAYuTr1+vPsPFx8B/qPlFrP674DkdKMzbyTZ2/vXs5NnKMh2MBX279Ktf/qOOb8pctMvPYaXaNpRpKlcdDZ6nr6tz1Jc0OEoEpz8M3/yMohk0KhDksPonCwVcvMY+iOMV54qIkr1Q/SM5Zh5+oGQr7SWzJTzDCYHmRr/NYfoEIsFos2qAHCjz1GcWAAZ2yc3JkzBNZ/cV+oGQiw7p/8PuG/+jmn9r3N32qDZza3s6Kh8lp8fo/J2uborQ+Y7CkNu8tNQftD0PEEmKW3w5XkFd7peYeMnWF743a2NGy5bgxz/tIlhn71KueHkmQKpWnf2Zp6+rc+yUMr6hk48gb9773OMquBQrYWSw/j1WcoqmUzj1Fc2c9XNtSRXfY8r3dHaI55Wd9ytd6eZA8aLd0eQrBIJrzcivJ4iOzahfJYZA8fvq09FM1YjPaXfo8ttRbLTuzn1cO97LswiutWzGTML+bYcOFtOPKz0veb/giW7wTTougW2du7l1cvvorH8PBS50tsb9x+fUhfvMilv32FE31TN4T0zkiaC3/zvzP63uv4rFrcbD2WHsOrT+ESxlalpQWGOy/R8oBJy4pv8EpfDJ9l8PyG5uta/l1TXYQ8IeoDV7tChFisFnWLGkq7lYefeILku+/iplKYkcgXHu84LlM6SvUzT7PprbcIDZziY7WOoUSOr61ruvVoikqQGCgtR5oehebNpYC2Sl03Q+kh3u55e2aUxY7mHXiM67uC0mfPcfJnLzOUyM/8LFvTwMTyNWy++Crneo6TztoErGq8uXbyegCvPoFLjLyxEUNB/7JBdEuOR+Nf5dcjLeTtIt/d1krYd/WtaLs2V5JXWFm9suIv2gpxPyz6oAbwr1mDd8UKjNtYMyObKHLp8DCxeIym7Q+y6uAnRJ063p1o5mcf9/DihiYaIjeZ4VdOrgPdH5a+vCHY8PulWYaA4zp8OvwpB4cOErSCvLj8RdqqbpxcMnr8NMf+5h9I552Zn9m+ALZyqD38IwaDPs5FQ6y2ozjpGDk9jE8fx6GBgrEay1B4VBJl9fJY1Ro+Mh5jeCrPixubaai6/vXqS/VRdIvS7SHENAnqabcT0gChmJdoQ5DJoQxufTuNK6ZovHCclx6r4x9TJn/7yRW++kCcNc1f3DK/b9KjcPpVSA5CfC10PgOe0nT4idwEb/e8zXBmmJXVK3m89fGbTiw5/9ERzv/db3Cmu3eUKq0Tkk13YaoUE+va6B5Ls7W7hkxGUWAMv3sUW7VRVB0ETQufqXDtU0SGl+BseIHT4/DEynqW14dveL6uqS48hoeWcMsN9wmxGElQ3yGlFEvW1XLi/SyJkSy6eg0NjUl8B/bw+8+/yO5hze6TgwwlSkPNyrZbtutC7ydweQ+YHlj3bahfBZTGWB8fPc7+/v1YhsVzHc+xPLb8hodwXM3+tz9m4o03Z34W8lnknSz9ubO4a1vIrt9IzRv9bE36mSgUUNj43GMUVCeOimOZioCRxdH9JPWDhD0dHO2vZsNDUba0x254Tq01XYku2iJtsva0ENPkf8Jd8Ic8NHfG6D0zQXKigK7fSn16D8W3dvPNl77D/mE/h7onGEnm+fqGpuv6X++L7ERpRMfkldJa0St3ga/Uck0VUrxz5R16k720R9r5SttXCHlCNzxEIlfk3Tc+xNj7HgAKaIx56U32ciaQxX5qMzGfSfXfHsO200y6y1A4+NwTFIz1uISxTEWTlWXMCWHrpaTqGnA7NtG6rpanVjXctP95JDtCuphmaWTpHL5AQswvEtR3qXFZlNHeNLlUgVTSxW1+hPpL75F8/TUe+853iEf8vHVqkJ9/1MPzG5qu7tI9l7SG/sOlsdFKldboaFw/MwX83MQ59vbtxdUuT7Y+yZraNTcNy8ujaT7YvZ/YwX0A+C2DqnCSMxOnObB2CbXLlrD2SDeTF4fBzmPoVYCB1z1FztgGgNc0qDZ9pFwPtrbpranFs+NJotVBvrm55Za/aXRNdaFQtEcqY/MJISqBBPVdMkyDjvW1nNk/gD/sJZspMtS8g/ruPSTffJOVL7xAzfZ2XjvWzy8P9fLkyno2tEbnbhRDLgFn/xHGL5X2LVz9dfCXxiVn7Sx7evdwcfIijaFGvtr+VaK+G8dYu67mwKUxTu79hNoj75PWaSxvGuWxGZoo8u72lTxxvo941xiTCQ+RYoYcqwALj75I3tgIgGUo/IbJpJvHweBysxe9bT0+j5cf7Fhy84WgphXdIq1VrQSshbWsrBD3QoL6HkTqAtS1VjHWn6J1VTW9Z2Go4UHiXR9h7NlD3VNP8f0H29l9cpB3zgwzmMixc3UDni8IqjumNQydLK3ToR3ofBZatpRa1EB3opv3rrxH1s7yUNNDbG7YjKGuf37HdeiaHOBXx0+QP7Sf+gtnGVcu1SEvAY/BaMElbbt84+gA7VWddI0NYGSOkzM3AgFMt5+CWnW1JCBlFykoL1eWTpJeu502Twud8TB+zxcPX3y05dHS/pNCiBkS1PeobU0Nk8MZJoezdG6Lc+EQ9GfW03zsOGYsRnDzZn5vYzMHLo1z4NIYo6k8L2xoJhqYhb0VC2k4txtGzkK0BVa/AMEaAIpOkQ/7P+Tk2Elq/DU8v+z5mf0GbddmKDNEf6qfgfQAZ0eucOHKCBv3nqHasQj6fIT8GtOAiUyRXN6mM7ocD1VcHDhNIXcJx9qEq0NAGltdv12Z42pcY5wrywwSq7fzL7c/R0dN5LbHmMvYaSGuJ0F9jzw+k9bV1XQdG6W+PUzn9kbOfQK9mSx6zwHMaBTfsmU8vLyWeMTHmycH+fnHPXx9XdPMKnF3ZfR8aSElOw/LnoK2h8AotZQH04O83fM2iXyCTfWb2BzfzGhmlI8mP6I/1c9QZgh3ej3oVMILR9N89XwCn7mMQFUSmykAJjNFikUfDd7lpDJQyH2Ca49SNJZj69h0ITcOr3OsXtSa1XzjxR+wrUOG2Alxr9Rc/Jq5bds2ffDgwVl/3Eqlteb0vgFy6SIbvtJKJlHg3IF+7LMnWeIfpO4738QTbwBgIl3gtWP9jKULPLqijm1Lqu+sBVnMwYXfweBxCDeUliMNlx7bcR0+GfqET4c+BSAejAMwnB1Ga41SivpAPc2hZqq9DZz/oAv748NY2TTKSmJ4x3G0TcF2GU8XqDIbCZm1OM4UCesY/qk8RVVHzrMBw7nxfeNVI4TbPTz8wp/SsHGTtIyFuANKqUNa6203vU+Cenakp/Kc+qCfurYqlm6oIzme4+wHPRRPHqOjNkXD91/CrCot3lSwXd46NcS5oSSd8TDPrInf3pKp45dLreh8srQUacdjYJhk7SwnRk/wyeAn1x1uKIN4ME5TuInmUDONoUa8ppfek+c5/PJuGB+j6OZI2P1g5ig6Gq01lvIRs1oJaNDZ00yocXw2FJSfnLHjuucwDYWhDHyh8zz3/X9KZOvTEtBC3IUvCmrp+pgloaiP+NIIgxenqGsNU1XjZ/Xj7Zy2bS6dPI7xqzeo/4NvYXi9eC2Dr69vpDHq44PzY4ynr/DChuZbL5nqFOHiu9B3CII1pNe9RL8JA/376E32MpmfvO7w7Y3baQ430xBswGN4cFzNeLrAuVO9nP6H1/FPjeNqh5QzQtoZQ6HxYBDymmi7inrXTzDTxbgzREYVMYGiMskZ22eeQ6nSEDxMD9t3pFn5yA9h2fza/UaI+UJa1LPIsV2Ov9eL5TFY83gLhqFIT+Y5+eZpCmdOs3J9FfXf+jrKuDrq4sp4hjeOD2C7mufWNrKi4XN9vlO9JE6+zECyl/5YE/3hWqaKqRueO2AF+O7K7+I3g4ylCwwn8gwncwwn80wMj9P44e/wZNNorcm7SRLOAFUBhc8y8JgKC4Ue19QW0uR1hgHS2KrUjz1V7aWoHyKYKH2uGwZ4DINgvJ6dO3JEm1ph/e/P9JELIe6cdH3cR+MDaS4cHKJtTS1Ny0tjldNTeU786gjFSxdZ/WgLdc9e3/JM5Iq8fmyAwakcD3ZU80CrxWCyl/7Lb9M/dJiUUlC3El84TlOoieZwM4lCgpOjp8kWXVaEtuHXbYyk8oylCjNrcoTyKVbt3w2uSyJnY+sCPY0O8dwwtY6Bq12y2STeXJ66rIuFwaAJCXccgIvrIuSCFjVdK4iOlmYvKgVB06Kuo5EH12eoamku9ZN77+HCqBBCuj7up+rGILF4kL5zE9Q0h/AFLEJRH+u/tZljP89xel8/ayNHqNmxCShdiCzqBGuWTNJ7/gx/caKb6JkEy1Q/YadAU90aNi9/ltpQG64d5srkJG+eeZ/LU5ex3Bra/Du4OBki4E3TUOVjS1uM+sQQwT1v4/MYDAU8XB5N0dtgUTQnWTaZwqN9TKSGcfNTNDgeIirElLeWbkYxvQlGG0KYjqapK4OZfgBHXZ1iXqU8tK5Zyepdm6latqy0jogQYk5Ji3oO5DNFjr/XR7jGR+fWOKan1CWQSeQ4+t/fxx4fZ/13tzPYbPNh/4fk7BwAISuIOTyK1X2Bek816x75Psdycfons0xkCkwW++jJfYJh2Gys286W+CbikQANER9BO0/+1Ckm9h0ia1ukCyZXJh3GMw5FfQRljBM1/ajsFAVnHD8GDUTJeWJ06wm8ThplFNEGWEWNi0nGegK0iass/G6RagV1Gx5g5bdfoLb5xmF5Qoi7J10fZTDSk+TysVH8IYsVW+MEI6ULhdmJNId//DZOJkdgYxUnqwdJhXJY5FmaGGalrQlE1/GbzFomixagcbRNyjhOIDhOfTjC15Y+Q7W/mlwiT+L0ZZInL5LqHSdvWzgu2LbNZCZD3h3B1eP43W6iukjGKGCZBkEslG0wTu66mm2PwVQ0Qt7agNZBYiN5FFClXZZFNPaSVSx9fhcNK+L3/wUVYoGToC6TxFiWi5+O4BQdlqyro64tjFKKzPAkx/7HHrKTWQAKVpYruotCnY2uDUHEh21muJIaI1vIYbgOlm1Qo4NUqxAq68VIeIh0ZaefycUkBTqBXZwiHcvjzyTwFRxMQ+Ex1Mz6GgnHQ1GbaAy0UigzzWRtPdlAA9gtM1PPa8dsqtIujWaKYLQaWlppXx2h6dlvluGVFGLhk6Auo2Le4eLhYRIjWepaq1iyvhbTMtBak+kbZvLtV5k6fJyE66W7qLEJYvst8lUeHJ9BRmnyjiZke6nKW3iLmjojS1BnUGSxnSx5O0/O0UzEggwEx2juHyaQz1EXqmKsYJMzguTworSLqYsYriYTqmWssR6rGEO5BjPvAgUBt4nW4Sx1ZgpvdQj/us0s8x0g9siLULeinC+nEAuWXEwsI4/PZNWDjfSfn6Tv/CTpqTwrtjYQqPISbGnAs+tZfB2t6HdOstz2kXGDOFi4Rcjm8qTtLFm7QLboYGOA8jNsRfDoAoYqkm1sJF8ToOboL/ANX6TTdahVEQpWA32FNAoFxTRB7aCJUbTaGWqqw/X40LZFzvSgTYNVtoNtRMkVLOKj40TNLJ6GGqofe4rl29vwFlsgeuMWXUKIuXdbLWql1C7g/wFM4Mda6//4RcdLi/rmpkayXDw8jGtrqmr9ZKYKFPM2AMpQeJ0MnD+O30nht2z8lo2hIJFPc+hKH4liDltbKIJYqgYvIUynQMFN4qIJaptGI4dtjDDh9FPUUWxqsVUdI/FaXMMi5/eilYHSLoa2MbRNB2EiKGzHxds3TthnUdtZQ/y5r9PQ2YQq1y41Qiwi99SiVkqZwH8BngF6gU+UUr/RWp+a3TIXvmh9gHVPtNB1bJR8xibaECAU9RGKeQlEvJimgbZXMf7Tn+Jm7Jm/ZykfdYE49QHNeHYKt9CLXTgMOo/GR0hXEVV1OKqRCRoYCDZAcRXevINGMdQWRikHQ+fwOzksQ+G1DKosP08ub8QbiDEwHAZPAKetSOeDzdStkx1WhKgUt9P18SBwQWt9CUAp9QvgG4AE9V3w+i1WPth4y/uVZVH7Z3+Gtm2yx46TOXiQoFJse2obgbWrcI/8nMH+IgP9XlLJJMGcQU24gaLycz7pI+Vm8SdLvyUVfQZTDT6CPoXHEyLqj5DGhrYMj21+iHXNW2fW5YjnbLqOjZbW2F5246YCQojyuZ2gbgGuXPN9L/DQ5w9SSv0Q+CFAe7tso3SvlGUR3LIZ/wOryRw8RPb4MQqXLxFoidO+62m0ruGNSzbe4W5Wn97N5Q1hckcu4RlKYZpFzLYY+Z3reTL2FaLZOL1Dg1xQJwg0G3xt+XeIh64fYvdlHyBCiPKZtYuJWusfAT+CUh/1bD3uYmcEAoQff4zA+nWkDxwgc/4CuZHjRLauZ1lbkoPH3sAd7eeQWkd283I2mFGeCcRZ8dhzhIJRbLc0qeaccYKWcAvPLHmGoEemewsxn9xOUPcB117ub53+mbiPzFiMyK5d5Df0cfT1nzL2yz3kI37WosjWN9ERfYRvrdvM6njtzN9JFVLs7trNUGaITQ2b2NG044ZtuIQQle92gvoToFMptZRSQH8P+MM5rUrcUrEuyskdjTStaGDV2UmieZOqLesJ7/zqdetA9yZ7eav7LWzX5rmO51geW17GqoUQ9+JLg1prbSul/jWwm9LwvL/UWp+c88rETYW9Yf6n9X8KgH7CIX/hIp6mxpmQ1lpzdOQo+wf2E/VG+eaKb1Ltry5nyUKIe3RbfdRa6zeAN+a4FnGHlGniX7Vy5vuCU+DdK+9ycfIiy2LL2Nm2E695i80IhBDzhsxMXCAmchO82fUmk7lJHml+hI31G2VLLCEWCAnqBeDS5CXeufIOpjJ5cfmLtFa1lrskIcQskqCex1ztcmDgAEeGj9AQbGBXxy7CXlknWoiFRoJ6nsoUM7zV/RZ9qT7W1q7l0ZZHsQz55xRiIZL/2fPQUHqIN7veJGfn2Nm+k9U1q8tdkhBiDklQzyNaa06Nn2Jv715CnhDf7vw29cH6cpclhJhjEtTzhO3a7Ondw5nxM7RVtfH0kqcJWIFylyWEuA8kqOeBRCHB7q7djGRG2BbfxrbGbTIVXIhFRIK6wl1JXOG33b9Fo/na0q+xNCrrRAux2EhQVyitNYeGDvHJ4CfU+GvYtXQXUZ+sEy3EYiRBXYHyTp53et7h8tRlOqs7ear1KTymp9xlCSHKRIK6wiQKCV69+CqJQoLHWh5jfd16mQouxCInQV1hAlaAmC/GzradNIWbyl2OEKICSFBXGI/h4fllz5e7DCFEBZExXkIIUeEkqIUQosJJUAshRIWToBZCiAonQS2EEBVOgloIISqcBLUQQlQ4CWohhKhwSms9+w+q1AjQPesPPLfqgNFyFzELFsJ5LIRzgIVxHgvhHGB+nMcSrfVNdwKZk6Cej5RSB7XW28pdx71aCOexEM4BFsZ5LIRzgPl/HtL1IYQQFU6CWgghKpwE9VU/KncBs2QhnMdCOAdYGOexEM4B5vl5SB+1EEJUOGlRCyFEhZOgFkKICidBfQ2l1HeVUieVUq5Sal4N5VFK7VJKnVVKXVBK/R/lruduKKX+Uik1rJQ6Ue5a7pZSqk0p9a5S6tT0e+nPy13T3VBK+ZVSHyuljk6fx38od013SyllKqUOK6VeK3ctd0uC+nongG8De8pdyJ1QSpnAfwG+BqwBvq+UWlPequ7KT4Bd5S7iHtnAv9VarwF2AP9qnv5b5IGdWuuNwCZgl1JqR3lLumt/DpwudxH3QoL6Glrr01rrs+Wu4y48CFzQWl/SWheAXwDfKHNNd0xrvQcYL3cd90JrPaC1/nT6dpJSQLSUt6o7p0tS0996pr/m3cgDpVQr8Dzw43LXci8kqBeGFuDKNd/3Mg/DYaFRSnUAm4GPylzKXZnuMjgCDANvaa3n43n8Z+DfAW6Z67gniy6olVK/U0qduMnXvGuBisqllAoD/wD8G611otz13A2ttaO13gS0Ag8qpdaVuaQ7opR6ARjWWh8qdy33atHtQq61frrcNcyBPqDtmu9bp38mykAp5aEU0n+jtX653PXcK631pFLqXUrXD+bThd5Hgd9TSn0d8AMRpdRfa61/UOa67tiia1EvUJ8AnUqppUopL/A94DdlrmlRUkop4C+A01rr/1Tueu6WUqpeKRWbvh0AngHOlLWoO6S1/vda61atdQel/xPvzMeQBgnq6yilvqWU6gUeBl5XSu0ud023Q2ttA/8a2E3p4tXfaa1PlreqO6eU+jmwH1illOpVSv1ZuWu6C48CfwzsVEodmf76ermLugtNwLtKqWOUGgJvaa3n7fC2+U6mkAshRIWTFrUQQlQ4CWohhKhwEtRCCFHhJKiFEKLCSVALIUSFk6AWQogKJ0EthBAV7v8HG0iK/b4q2B4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(df[(df['chain_id'] == 0).values]['x1'], df[(df['chain_id'] == 0).values]['x2'], alpha=0.5)\n",
    "plt.plot(df[(df['chain_id'] == 1).values]['x1'], df[(df['chain_id'] == 1).values]['x2'], alpha=0.5)\n",
    "plt.plot(df[(df['chain_id'] == 2).values]['x1'], df[(df['chain_id'] == 2).values]['x2'], alpha=0.5)\n",
    "plt.plot(df[(df['chain_id'] == 3).values]['x1'], df[(df['chain_id'] == 3).values]['x2'], alpha=0.5)\n",
    "plt.plot(df[(df['chain_id'] == 4).values]['x1'], df[(df['chain_id'] == 4).values]['x2'], alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mineral-barrier",
   "metadata": {},
   "source": [
    "Finally let us plot a histogram of this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "affecting-collapse",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD7CAYAAABDld6xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAL5klEQVR4nO3d3Y8ddR3H8c+n26VLoYYLQZASy4UhaUigcVNUuLGKqdVANDGBRK8weyMGExOfboz+AcYbYmiQaOIDGtTE4APWUNIQFSkPmkIhVoKhSKwECEUN0O7Hiz2VFg7s2c7MznfnvF/Jht3uOZPvsN13fp09vx0nEQCgrnV9DwAAeGuEGgCKI9QAUByhBoDiCDUAFEeoAaC49ZM8yPaTko5KOi7pWJL5LocCALxmolCPfCDJs51NAgAYayWhntgZ3pA5ndXFoYE1z+uaX3HM4mILk6CSo3r+2STnjvvcpKGOpN/ajqRbkux+qwfP6Sxd4Q+ucExgOsycvanxMY4fPdrCJKjkd7nj72/2uUlDfVWSp22fJ2mP7ceS7Dv5AbYXJC1I0pw2nvawAIBTTfRvsCRPj/57RNLPJW0f85jdSeaTzM9qQ7tTAsAUWzbUts+yvenE+5I+LOlA14MBAJZMcunjHZJ+bvvE43+Y5DedTgUA+L9lQ53kCUmXrcIsAIAx2JkIAMURagAojlADQHGEGgCKI9QAUFwnv+sDAKbJzKbmvxZAL775p1hRA0BxhBoAiiPUAFAcoQaA4gg1ABRHqAGgOEINAMURagAojlADQHGEGgCKYws5sMq4g/jwdP01ZUUNAMURagAojlADQHGEGgCKI9QAUByhBoDiCDUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQ3MShtj1j+yHbd3Y5EADgVCtZUd8k6WBXgwAAxpso1LY3S/qopFu7HQcA8HqTrqi/JemLkha7GwUAMM6yobb9MUlHkjywzOMWbO+3vf9VvdzagAAw7SZZUV8p6RrbT0q6XdIO299//YOS7E4yn2R+VhtaHhMApteyoU7ylSSbk2yRdJ2ku5N8qvPJAACSuLktsCLr33lB42Mc+8czLUyCabKiUCe5R9I9nUwCABiLnYkAUByhBoDiCDUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQHFvIgRVg+zf6wIoaAIoj1ABQHKEGgOIINQAUR6gBoDhCDQDFEWoAKI5QA0BxhBoAiiPUAFAcW8gxNWY2bWp8DJ9/XuNjHPvr3xofA9OFFTUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQHKEGgOIINQAUR6gBoLhlt5DbnpO0T9KG0ePvSPK1rgcD2sb2b6xVk/yuj5cl7Ujyku1ZSffa/nWSP3Y8GwBAE4Q6SSS9NPpwdvSWLocCALxmomvUtmdsPyzpiKQ9Se4b85gF2/tt739VL7c8JgBMr4lCneR4ksslbZa03falYx6zO8l8kvlZbWh5TACYXit61UeSFyTtlbSzk2kAAG+wbKhtn2v7nNH7Z0q6WtJjHc8FABiZ5FUfF0j6nu0ZLYX9J0nu7HYsAMAJk7zq4y+Stq3CLACAMdiZCADFEWoAKK6Tu5B73TrNnH36d3w+fvRoi9MAS159Rwt3If9rC4MAK8SKGgCKI9QAUByhBoDiCDUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQXCdbyLO4yDbwQmY2Nds6PZSv5cy/X2l8jMUW5gBWihU1ABRHqAGgOEINAMURagAojlADQHGEGgCKI9QAUByhBoDiCDUAFEeoAaA47kI+BYby/zNXXd7o+f7z39oZBFhlrKgBoDhCDQDFEWoAKI5QA0BxhBoAils21LYvsr3X9qO2H7F902oMBgBYMsnL845J+kKSB21vkvSA7T1JHu14NgCAJlhRJ3kmyYOj949KOijpwq4HAwAsWdE1attbJG2TdF8n0wAA3mDinYm2z5b0U0mfT/LimM8vSFqQpDltHMxuONTxr20bGz3/vHv5O4m1aaIVte1ZLUX6B0l+Nu4xSXYnmU8yP6sNbc4IAFNtkld9WNJ3JB1M8s3uRwIAnGySFfWVkj4taYfth0dvuzqeCwAwsuw16iT3SvIqzAIAGIOdiQBQHKEGgOIINQAUR6gBoDhCDQDFEWoAKK6Tm9s2NbPp9G+MewJb2Idn/X/7ngDoBytqACiOUANAcYQaAIoj1ABQHKEGgOIINQAUR6gBoDhCDQDFEWoAKI5QA0BxJbeQs/17eI587v2Nj3H+vS80ev5i4wmAfrCiBoDiCDUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQHKEGgOIINQAUV3IL+fp3XtD4GMf+8UwLkwxDG3d1b+o/zb+k8qGnmh8EWINYUQNAccuG2vZtto/YPrAaAwEATjXJivq7knZ2PAcA4E0sG+ok+yQ9twqzAADG4Bo1ABTX2qs+bC9IWpCkOW1s67AAMPVaW1En2Z1kPsn8rDa0dVgAmHpc+gCA4iZ5ed6PJP1B0iW2D9u+ofuxAAAnLHuNOsn1qzEIAGC8klvI29j+3XQbepUt6G1s/256V/cnf3xZ4xnO2N/4ENydHlOLa9QAUByhBoDiCDUAFEeoAaA4Qg0AxRFqACiOUANAcYQaAIoj1ABQHKEGgOJKbiFvQ5Ut4E21sW163batjZ5/xxW3NJ7hS5+5uvExjjc+ArA2saIGgOIINQAUR6gBoDhCDQDFEWoAKI5QA0BxhBoAiiPUAFAcoQaA4gg1ABQ32C3kTTXddi1Jz77nnOaDtOD+b3y70fN3XdLC9m/uIA6cNlbUAFAcoQaA4gg1ABRHqAGgOEINAMURagAojlADQHEThdr2TtuP2z5k+8tdDwUAeM2yobY9I+lmSR+RtFXS9bab7wYBAExkkhX1dkmHkjyR5BVJt0u6ttuxAAAnTLKF/EJJT5308WFJV3QzTh0+9NTyD1rGef98vvExDt14ceNj7LrkqkbPZ/s30K/WfteH7QVJC5I0p41tHRYApt4klz6elnTRSR9vHv3ZKZLsTjKfZH5WG9qaDwCm3iShvl/Su21fbPsMSddJ+kW3YwEATlj20keSY7ZvlHSXpBlJtyV5pPPJAACSJrxGneRXkn7V8SwAgDHYmQgAxRFqACiOUANAcYQaAIrj5rZvoo3deDMtzLHlq79vfIzjLcwBoD+sqAGgOEINAMURagAojlADQHGEGgCKI9QAUByhBoDiCDUAFEeoAaA4Qg0AxTlJ+we1/yXp760feGXeLunZnmfowlDPSxruuQ31vKThnlsf5/WuJOeO+0Qnoa7A9v4k833P0bahnpc03HMb6nlJwz23aufFpQ8AKI5QA0BxQw717r4H6MhQz0sa7rkN9byk4Z5bqfMa7DVqABiKIa+oAWAQBh1q25+0/YjtRdtlfoJ7umzvtP247UO2v9z3PG2xfZvtI7YP9D1Lm2xfZHuv7UdHfw9v6numNties/0n238endfX+56pTbZnbD9k+86+Zzlh0KGWdEDSJyTt63uQpmzPSLpZ0kckbZV0ve2t/U7Vmu9K2tn3EB04JukLSbZKeq+kzw7ka/aypB1JLpN0uaSdtt/b70ituknSwb6HONmgQ53kYJLH+56jJdslHUryRJJXJN0u6dqeZ2pFkn2Snut7jrYleSbJg6P3j2rpm//CfqdqLkteGn04O3obxA+7bG+W9FFJt/Y9y8kGHeqBuVDSUyd9fFgD+KafFra3SNom6b6eR2nF6PLAw5KOSNqTZBDnJelbkr4oabHnOU6x5kNt+3e2D4x5G8RqE2uf7bMl/VTS55O82Pc8bUhyPMnlkjZL2m770p5Hasz2xyQdSfJA37O83vq+B2gqyYf6nmGVPC3popM+3jz6MxRme1ZLkf5Bkp/1PU/bkrxge6+Wfsaw1n8YfKWka2zvkjQn6W22v5/kUz3PtfZX1FPkfknvtn2x7TMkXSfpFz3PhLdg25K+I+lgkm/2PU9bbJ9r+5zR+2dKulrSY70O1YIkX0myOckWLX1/3V0h0tLAQ23747YPS3qfpF/avqvvmU5XkmOSbpR0l5Z+KPWTJI/0O1U7bP9I0h8kXWL7sO0b+p6pJVdK+rSkHbYfHr3t6nuoFlwgaa/tv2hpAbEnSZmXsg0ROxMBoLhBr6gBYAgINQAUR6gBoDhCDQDFEWoAKI5QA0BxhBoAiiPUAFDc/wDaTZegbGllawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist = plt.hist2d(df['x1'].T.values[0], df['x2'].T.values[0], bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "finished-prior",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "decimal-ordinary",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}