rl_ltl_pac.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code Supplement of _Reinforcement Learning with General LTL Objectives Is Impossible_\n",
    "\n",
    "In this notebook, you will be able to: \n",
    "* [Reproduce the results in Section 5 and Appendix F of the paper](#paper-repro)\n",
    "* [Visualize constructed counterexample MDPs for different LTL formulas](#mdp-construction)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Dict, Any, Callable, List, Optional, Union, Tuple\n",
    "\n",
    "import math\n",
    "import re\n",
    "import subprocess\n",
    "import tempfile\n",
    "import pathlib\n",
    "import dataclasses\n",
    "import textwrap\n",
    "import concurrent.futures\n",
    "import os\n",
    "import functools\n",
    "import itertools\n",
    "import collections\n",
    "import hashlib\n",
    "import json\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "\n",
    "from scipy import interpolate\n",
    "\n",
    "from tqdm.notebook import tqdm_notebook as tqdm\n",
    "import inspect\n",
    "\n",
    "from IPython.display import display, Markdown\n",
    "\n",
    "import spot\n",
    "from spot.jupyter import display_inline\n",
    "spot.setup(show_default='.bans')\n",
    "\n",
    "from num2tex import num2tex\n",
    "\n",
    "import buddy\n",
    "import graphviz\n",
    "\n",
    "plt.style.use('seaborn-paper')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproduction of Empirical Results (Section 5 and Appendix F) <a class=\"anchor\" id=\"paper-repro\"></a>\n",
    "\n",
    "We first declare the main routine for estimating LTL-PAC probability under a particular setting.\n",
    "\n",
    "The routine simply calls into (a modified version of) Mungojerrie to run multiple reinforcement-learning algorithm sessions in parallel obtain a Monte Carlo estimate of the LTL-PAC probability.\n",
    "\n",
    "For documentation on Mungojerrie, please see https://plv.colorado.edu/wwwmungojerrie/docs/index.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def canonicalize_args(f):\n",
    "    \"\"\"Wrapper for functools.lru_cache() to canonicalize default                                                          \n",
    "    and keyword arguments so cache hits are maximized.\"\"\"\n",
    "\n",
    "    @functools.wraps(f)\n",
    "    def wrapper(*args, **kwargs):\n",
    "        sig = inspect.getfullargspec(f.__wrapped__)\n",
    "\n",
    "        # build newargs by filling in defaults, args, kwargs                                                            \n",
    "        newargs = [None] * len(sig.args)\n",
    "        newargs[-len(sig.defaults):] = sig.defaults\n",
    "        newargs[:len(args)] = args\n",
    "        for name, value in kwargs.items():\n",
    "            newargs[sig.args.index(name)] = value\n",
    "\n",
    "        return f(*newargs)\n",
    "\n",
    "    return wrapper\n",
    "\n",
    "class NumpyEncoder(json.JSONEncoder):\n",
    "    \"\"\" Custom encoder for numpy data types \"\"\"\n",
    "    def default(self, obj):\n",
    "        if isinstance(obj, (np.int_, np.intc, np.intp, np.int8,\n",
    "                            np.int16, np.int32, np.int64, np.uint8,\n",
    "                            np.uint16, np.uint32, np.uint64)):\n",
    "\n",
    "            return int(obj)\n",
    "\n",
    "        elif isinstance(obj, (np.float_, np.float16, np.float32, np.float64)):\n",
    "            return float(obj.asttype(np.float32))\n",
    "\n",
    "        elif isinstance(obj, (np.complex_, np.complex64, np.complex128)):\n",
    "            return {'real': obj.real, 'imag': obj.imag}\n",
    "\n",
    "        elif isinstance(obj, (np.ndarray,)):\n",
    "            return obj.tolist()\n",
    "\n",
    "        elif isinstance(obj, (np.bool_)):\n",
    "            return bool(obj)\n",
    "\n",
    "        elif isinstance(obj, (np.void)): \n",
    "            return None\n",
    "\n",
    "        return json.JSONEncoder.default(self, obj)\n",
    "\n",
    "\n",
    "def persistent_memo(func):\n",
    "    cache_file = pathlib.Path(f\"{func.__name__}.cache.json\")\n",
    "    try:\n",
    "        with cache_file.open(\"r\") as f:\n",
    "            memo = json.load(f)\n",
    "    except FileNotFoundError:\n",
    "        memo = {}\n",
    "    @functools.wraps(func)\n",
    "    def wrapped(*args, **kwargs):\n",
    "        key = json.dumps((args, kwargs), sort_keys=True, cls=NumpyEncoder)\n",
    "        key_digest = hashlib.md5(key.encode(\"utf-8\")).hexdigest()\n",
    "        if key_digest not in memo:\n",
    "            result = func(*args, **kwargs)\n",
    "            memo[key_digest] = {\"args\": args, \"kwargs\": kwargs, \"result\": result}\n",
    "            with cache_file.open(\"w+\") as f:\n",
    "                json.dump(memo, f, sort_keys=True, cls=NumpyEncoder)\n",
    "        return memo[key_digest][\"result\"]\n",
    "    return wrapped\n",
    "\n",
    "\n",
    "mungojerrie_path = \"mungojerrie/build/mungojerrie\"\n",
    "\n",
    "pattern = re.compile(rf\"PAC Probability for tol (\\d+(\\.\\d*)?) is: (?P<satprob>\\d+(\\.\\d*)?)±(?P<std>\\d+(\\.\\d*)?)\")\n",
    "\n",
    "@canonicalize_args\n",
    "@persistent_memo\n",
    "def estimate_pac_probability(epsilon: float, \n",
    "                             model: str, \n",
    "                             min_num_samples: int = 100, \n",
    "                             max_est_std: float = 1e-3,\n",
    "                             reward_type: str = \"multi-discount\", \n",
    "                             discount: float = 0.99999,\n",
    "                             gammaB: float = 0.99, \n",
    "                             learn=\"Q\",\n",
    "                             lr_schedule: Union[Tuple[str, float], Tuple[str, float, float]] = (\"kkt\", 10),\n",
    "                             explore_schedule: Tuple[float, float] = (0.1, 0.1),\n",
    "                             ep_number: int = 100, \n",
    "                             ep_length: int = 10,\n",
    "                             tolerance=0.,\n",
    "                             sarsa_lambda: Optional[float] = None,\n",
    "                             model_params: Dict[str, Any] = None):\n",
    "    \"\"\"Estimate LTL-PAC probability for a specific setting, calls into Mungojerrie.\n",
    "    \n",
    "    Args:\n",
    "        epsilon: the $\\epsilon$ in LTL-PAC definition.\n",
    "        model: the environment MDP to solve for. \n",
    "               there should be two files with paths `models/{env}.prism` and `models/{env}.ltl`, specifying the\n",
    "               environment MDP in PRISM format and the LTL formula.\n",
    "        min_num_samples: minimum repeats to run for Monte Carlo estimation.\n",
    "        max_est_std: the standard deviation below which the Monte Carlo estimation should terminate.\n",
    "        reward_type: the reward scheme to use for Mungojerrie.\n",
    "        discount, gammaB, learn, ep_number, ep_length, tolerance, sarsa_lambda: see documentation of mungojerrie.\n",
    "        lr_schedule, explore_schedule: a wrapper interface for mungojerrie's learning_rate and exploration schedules. \n",
    "        model_params: the parameters to be formatted into the PRISM model (e.g., the parameter $p$ in Figure 2 of the paper).\n",
    "    \n",
    "    Returns:\n",
    "        the estimated LTL-PAC probability and its estimated standard deviation.\n",
    "    \"\"\"\n",
    "    if model_params is None:\n",
    "        model_params = {}\n",
    "    \n",
    "    with tempfile.NamedTemporaryFile(mode='w+') as tmp:\n",
    "        model_str = pathlib.Path(f\"{model}.prism\").read_text()\n",
    "        model_str = model_str.format(**model_params)\n",
    "        tmp.write(model_str)\n",
    "        tmp.flush()\n",
    "        try:\n",
    "            cmd = [\n",
    "                    mungojerrie_path, tmp.name, \n",
    "                    \"--epsilon\", \"1e-12\", # necessary for accurate model checking \n",
    "                    \"--ltl-file\", f\"{model}.ltl\",\n",
    "                    \"--learn\", learn, \n",
    "                    \"--reward-type\", reward_type, \n",
    "                    \"--discount\", str(discount),\n",
    "                    \"--gammaB\", str(gammaB),\n",
    "                    \"--tolerance\", str(tolerance), \n",
    "                    # Learning rates\n",
    "                    \"--alpha\", str(lr_schedule[1]) if lr_schedule[0] == \"linear\" else \"0.1\",\n",
    "                    \"--linear-lr-decay\", str(lr_schedule[2]) if lr_schedule[0] == \"linear\" else \"-1\",\n",
    "                    \"--kkt-lr-decay\", str(lr_schedule[1]) if lr_schedule[0] == \"kkt\" else \"-1\",\n",
    "                    # Exploration\n",
    "                    \"--explore\", str(explore_schedule[0]),\n",
    "                    \"--linear-explore-decay\", str(explore_schedule[1]),\n",
    "                    \"--ep-number\", str(ep_number),\n",
    "                    \"--ep-length\", str(ep_length),\n",
    "                    # PAC settings\n",
    "                    \"--est-pac\",\n",
    "                    \"--est-pac-probability-min-samples\", str(min_num_samples),\n",
    "                    \"--est-pac-max-std\", str(max_est_std),\n",
    "                    \"--est-pac-epsilon\", str(epsilon),\n",
    "                    \"--lambda\", str(sarsa_lambda) if sarsa_lambda else \"0\"\n",
    "                    # \"--seed\", str(2),\n",
    "            ]\n",
    "#             print(cmd)\n",
    "            output = subprocess.check_output(\n",
    "                cmd,\n",
    "                stderr=subprocess.STDOUT\n",
    "            )\n",
    "        except subprocess.CalledProcessError as exc:\n",
    "            print(\"Status : FAIL\", exc.returncode, exc.output)\n",
    "            raise exc\n",
    "    output = output.decode(\"utf-8\")\n",
    "    match = pattern.search(output)\n",
    "    pac_prob = float(match.group(\"satprob\"))\n",
    "    pac_prob_std = float(match.group(\"std\"))\n",
    "    return pac_prob, pac_prob_std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we define a function that scans over a grid of settings (see Section 5 of the paper) to estimate all those setting's LTL-PAC probabilities.\n",
    "\n",
    "We also define a function for visualization of the obtained results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_ps = np.geomspace(0.1, 1e-3, num=5)\n",
    "default_ep_numbers = np.geomspace(1, 10000, num=21, dtype=int)[1:]\n",
    "default_ep_length = 10\n",
    "\n",
    "def collect_pac_probs(env: str, \n",
    "                      ps: List[float] = default_ps, \n",
    "                      ep_numbers: List[int] = default_ep_numbers,\n",
    "                      ep_length: int = default_ep_length, \n",
    "                      rl_algo: str = \"Q\", \n",
    "                      reward_type: str = \"multi-discount\", \n",
    "                      skip_after_stable: int = float(\"inf\")):\n",
    "    \"\"\"Estimate LTL-PAC probability for varying settings.\n",
    "    \n",
    "    Args:\n",
    "        env: see :py:func:`estimate_pac_probability`.\n",
    "        ps: the parameteter $p$ to the environment MDP.\n",
    "        ep_numbers: an array of $N$, the number of epsidoes, to run an algorithm for.\n",
    "        ep_length: each episode's length until an environment reset. \n",
    "            This is a hyper-parameter of reinforcement-learning algorithms in Mungojerrie.\n",
    "            `ep_length * ep_number` is the total number of samples available to an algorithm.\n",
    "        rl_algo: the learning algorithm to use for Mungojerrie. \n",
    "            Must be One of `Q` (Q-learning), `DQ`(Double Q-learning), `SL` (SARSA($\\lambda$)).\n",
    "        reward_type: the reward scheme to use for Mungojerrie.\n",
    "        skip_after_stable: skip the rest of the estimation for larger $N$, when the estimated LTL-PAC probability \n",
    "            stablizes at 1. Skipping vastly reduces computation time, but assumes an algorithm cannot do worse with \n",
    "            larger number of samples. By default this is set to `inf` which means no skipping.\n",
    "    \n",
    "    Returns:\n",
    "        a dict of estimated LTL-PAC probabilities, and a dict of their estimated standard deviations.\n",
    "    \"\"\"\n",
    "    lr_schedule = {\n",
    "        \"Q\": (\"kkt\", 10),\n",
    "        \"DQ\": (\"kkt\", 30),\n",
    "        \"SL\": (\"kkt\", 10),\n",
    "    }[rl_algo]\n",
    "\n",
    "    explore_schedule = {\n",
    "        \"Q\": (1.0, 0.1),\n",
    "        \"DQ\": (1.0, 0.1),\n",
    "        \"SL\": (1.0, 0.001),\n",
    "    }[rl_algo]\n",
    "\n",
    "    probs, errs = collections.defaultdict(list), collections.defaultdict(list)\n",
    "    for p in tqdm(ps, desc=f\"<{reward_type} | {rl_algo} > Iteration over p\"):\n",
    "        err_zero_cnt = 0\n",
    "        for ep_number in tqdm(ep_numbers, leave=False, desc=\"Iteration over N\"):\n",
    "            if err_zero_cnt < skip_after_stable:\n",
    "                prob, err = estimate_pac_probability(\n",
    "                    epsilon=0.1, \n",
    "                    model=f\"models/{env}\",\n",
    "                    model_params=dict(P=p),\n",
    "                    reward_type=reward_type,\n",
    "                    max_est_std=1e-2, \n",
    "                    discount=0.99999,\n",
    "                    gammaB=0.99,\n",
    "                    learn=rl_algo,\n",
    "                    lr_schedule=lr_schedule,\n",
    "                    explore_schedule=explore_schedule,\n",
    "                    ep_number=ep_number, \n",
    "                    ep_length=ep_length,\n",
    "                    sarsa_lambda=0.)\n",
    "#                 print(p, ep_number, prob, err)\n",
    "            if prob == 1.0 and err == 0:\n",
    "                err_zero_cnt += 1\n",
    "            else:\n",
    "                err_zero_cnt = 0\n",
    "            probs[p].append(prob)\n",
    "            errs[p].append(err)\n",
    "    \n",
    "    return probs, errs\n",
    "\n",
    "\n",
    "def interpolate_and_intersect_pac_delta(probs, \n",
    "                                        errs, \n",
    "                                        delta:float = 0.1, \n",
    "                                        ep_numbers: List[int] = default_ep_numbers, \n",
    "                                        ep_length: int = default_ep_length):\n",
    "    \"\"\"Interpolates the LTL-PAC curves (the left in Figure 3) to obtain the intercepts (the right in Figure 3).\n",
    "    \n",
    "    Args:\n",
    "        probs, errs: the estimated LTL-PAC probabilities and their standard deviations, \n",
    "                     by :py:func:`collect_pac_probs`.\n",
    "        delta: the desired LTL-PAC probability to reach. \n",
    "               For example, `delta=0.1` means we would like to intercept \n",
    "               at 0.9 LTL-PAC probability.\n",
    "        ep_numbers, ep_length: see :py:func:`collect_pac_probs`.\n",
    "    \n",
    "    Returns:\n",
    "        the intercept curves and their errors.\n",
    "    \"\"\"\n",
    "    line, lower, upper = [], [], []\n",
    "    for p in probs:\n",
    "        prob = np.array(probs[p])\n",
    "        err = np.array(errs[p])\n",
    "        try:\n",
    "            line.append(interpolate.UnivariateSpline(ep_numbers,  prob - (1 - delta), s=0).roots()[0])\n",
    "            lower.append(interpolate.UnivariateSpline(ep_numbers, prob - err * 2 - (1 - delta), s=0).roots()[0])\n",
    "            upper.append(interpolate.UnivariateSpline(ep_numbers, prob + err * 2 - (1 - delta), s=0).roots()[0])\n",
    "        except IndexError:\n",
    "            # The intercept does not exist, so we set them to +inf\n",
    "            line.append(np.inf)\n",
    "            lower.append(np.inf)\n",
    "            upper.append(np.inf)\n",
    "    return  (np.array(x) * ep_length for x in (lower, line, upper))\n",
    "\n",
    "def frexp10(x):\n",
    "    exp = int(math.floor(math.log10(abs(x))))\n",
    "    a, b = x / 10**exp, exp\n",
    "    if a == 0:\n",
    "        return \"$0$\"\n",
    "    if a == 1:\n",
    "        return f\"$10^{b}$\"\n",
    "    return f\"${a:.1f} \\\\times 10^{b}$\"\n",
    "\n",
    "def plot_pac_probs(env: str, \n",
    "                   probs, \n",
    "                   errs, \n",
    "                   ep_numbers: List[int] = default_ep_numbers,\n",
    "                   ep_length: int = default_ep_length, \n",
    "                   rl_algo: str = \"Q\", \n",
    "                   reward_type: str = \"multi-discount\", \n",
    "                   title: Optional[str]=None, \n",
    "                   delta=0.1, \n",
    "                   max_intercept_num_samples:Optional[int]=None, \n",
    "                   plot_for_paper: bool = True):\n",
    "    \"\"\"Plots the figure as in the paper's Figure 3.\"\"\"\n",
    "    if title is None:\n",
    "        with pathlib.Path(f\"models/{env}.prism\").open(\"r\") as f:\n",
    "            lines = f.readlines()\n",
    "            env_title = env\n",
    "            for line in lines:\n",
    "                m = re.match(r\"\\/\\/\\s*\\@title\\s*:\\s*(?P<env_title>[^\\n]*)\\n\", line)\n",
    "                if m:\n",
    "                    env_title = m.group(\"env_title\")\n",
    "                    break\n",
    "    \n",
    "    sample_numbers = ep_numbers * ep_length\n",
    "    if plot_for_paper:\n",
    "        plt.figure(figsize=(3.5,2.1))\n",
    "    ps = sorted(probs.keys(), reverse=True)\n",
    "    lines = []\n",
    "    linestyles = ['solid', (0, (3, 1, 1, 1, 1, 1)), 'dashdot', 'dashed', 'dotted']\n",
    "    linestyles_iter = itertools.cycle(linestyles)\n",
    "    for p in ps:\n",
    "        prob, err = map(np.array, (probs[p], errs[p]))\n",
    "        err = np.maximum(err, 1e-3) * 2 # we visualize 2 * std in the plots\n",
    "        filler = plt.fill_between(sample_numbers, prob-err, prob+err, alpha=0.4)\n",
    "        eb = plt.errorbar(sample_numbers, prob, yerr=err, linewidth=1, linestyle=next(linestyles_iter))\n",
    "        lines.append((filler, eb))\n",
    "\n",
    "    plt.gca().set_xscale(\"log\")\n",
    "    major_ticks = np.geomspace(sample_numbers[0], sample_numbers[-1], num=5, dtype=int)\n",
    "    minor_ticks = np.linspace(major_ticks[:-1], major_ticks[1:], num=6)[1:-1].T.ravel()\n",
    "    plt.gca().set_xticks(major_ticks)\n",
    "    plt.gca().set_xticklabels([frexp10(x) for x in major_ticks])\n",
    "    plt.gca().set_xticks(minor_ticks, minor=True)\n",
    "    plt.gca().set_xticklabels([], minor=True)\n",
    "    plt.gca().tick_params(which=\"major\", length=4)\n",
    "    plt.gca().tick_params(which=\"minor\", length=3)\n",
    "    plt.gca().set_xlabel(\"#samples\")\n",
    "    plt.gca().set_ylabel(\"LTL-PAC Prob.\")\n",
    "\n",
    "    plt.axhline(1 - delta, color=\"#000\", alpha=0.4, linestyle=\"--\", linewidth=0.6)\n",
    "    \n",
    "    pstex = [f\"{p:.0e}\".replace(\"-0\",\"-\") for p in ps]\n",
    "    plt.legend(lines, pstex, title=\"$p$\", handlelength=1.2, loc=\"best\")\n",
    "    rl_algo_to_name = {\"Q\": \"Q-Learning\", \"DQ\": \"Double Q-Learning\", \"SL\": \"Sarsa($0$)\"}\n",
    "    reward_type_to_name = {}\n",
    "    if title is None:\n",
    "        plt.title(f\"{env_title} | {reward_type_to_name.get(reward_type, reward_type)} | {rl_algo_to_name[rl_algo]}\")\n",
    "    elif title != \"\":\n",
    "        plt.title(title)\n",
    "    plt.gcf().tight_layout()\n",
    "    plt.gcf().savefig(f\"plots/{reward_type}-{rl_algo}-{env}\", dpi=600)\n",
    "    plt.show()\n",
    "    \n",
    "    # Plot intercept\n",
    "    lower, line, upper = interpolate_and_intersect_pac_delta(probs, errs, \n",
    "                                                             delta=delta, \n",
    "                                                             ep_numbers=ep_numbers, \n",
    "                                                             ep_length=ep_length)\n",
    "    if plot_for_paper:\n",
    "        plt.figure(figsize=(2.217,2.1))\n",
    "    plt.errorbar(range(len(ps)), line, yerr=np.array([line-lower, upper-line]), color=\"#000\", alpha=0.4)\n",
    "    plt.gca().set_xlabel(\"$p$\")\n",
    "    plt.gca().set_xticks(range(len(ps)))\n",
    "    plt.gca().set_xticklabels(pstex)\n",
    "    plt.gca().set_ylabel(\"#samples\")\n",
    "    if max_intercept_num_samples is not None:\n",
    "        plt.ylim([0, max_intercept_num_samples])\n",
    "    plt.gcf().tight_layout()\n",
    "    if title is None:\n",
    "        plt.title(f\"{env_title} | {reward_type_to_name.get(reward_type, reward_type)} | {rl_algo_to_name[rl_algo]}\")\n",
    "    elif title != \"\":\n",
    "        plt.title(title)\n",
    "    if plot_for_paper:\n",
    "        plt.gcf().savefig(f\"plots/{reward_type}-{rl_algo}-{env}-intercept-{str(delta).replace('.','_')}\", dpi=600)\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproduction of Results in Appendix F\n",
    "\n",
    "Running the cell below reproduces results in Appendix F. \n",
    "\n",
    "It takes approximately 2 hours to finish all the experiments. \n",
    "For quick reproduction, you may want to run for less `reward_types` and `rl_algos`.\n",
    "\n",
    "Note that the code caches previously ran experiments, so it's ok to interrupt and re-execute any time.\n",
    "If you are using Docker, the cache is recreated for every new container launch (meaning that by default, caches are cleared if you exit Docker).\n",
    "If you are in your own environment, the cache is stored at `estimate_pac_probability.cache.json`. You may delete or rename this file to fully re-execute a reproduction.\n",
    "\n",
    "Also note that due to multi-threading's non-determinism, the results might not look exactly the same as we report in the paper. However, the difference should be minimal since we always repeat enough simulations until the standard deviation of the estimate is within $0.01$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eed1a782b34148a9974e317c454e5274",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "<zeta-reach | Q > Iteration over p:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAACQCAYAAADKtM0CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzd0lEQVR4nO2dd3xUZbrHv+9Meu8JJRAg9A6hgwRERVmxLOoq1rVdC2vZe8Wyu7ZdxfWurousyrWsorgrqyjggiK9hhJIARJKEtJ7nWT6vPePmWAISZjMnBCSnO/nM59MznnnOc9M8sxbzvv8HiGlREVFpfui6WwHVFRUOhY1yFVUujlqkKuodHPUIFdR6eaoQa6i0s1Rg1xFpZvj0dkOuEtERISMi4vrbDdUVDqdw4cPl0spI5sf75AgF0L0BjYAI4AAKaWl2bnPAR/gD1LKn4QQgcBqIAz4QEr5mbPXiouL49ChQ4r6r6LSFRFCnG3peEcN1yuBK4H9LZx7FvgdcLXjJ8CDwJfAFcADQgivDvJLRaXH0SFBLqU0SCmrWjk9BtgnpdQBdY5efBrwk5TSCqQAQzvCLxWVnkhnzMm18ue9tDVAKBAC1DY7puIieZUNPLY6mdyKBkL9PVttV1VvbvP8SEsGp3We/NHrHyzzepxnTe+e+xkuqwH4H+/fsbThbd7we6rV88+a3uVV02J+7/VFi+dv4a98Lh/gTtv7/DVsEH8rzec/ITUc8IMiEUmUrKRUhBElKyk3DuS9ilL+E1LDEKPkn5ZnMEa+RakIo7e5gkKPcG5KLmVmqo1gK9TICFbdWINn2d08sPcTqrQCD2zsHQf7xmm4/18aQvUWsIRDWDUvXTOTv+3aT0Fdf2qCRhNe9C1bZt+EzprDzbuTyY5/lqkZy9g0ZST9POOI2LOTothrCag+THF4PgV9fkNs/p+JrBuPh0ckvsX72Tv3ZmxlJhJPrudM/3uIy17O1unTafD1Yf6WHZT0u4/Igp2UxdWSGfEwt6T9iVO+M9BICzbDAY5Muh1LtSdv/P23Lv0/iI7cuy6E2A7MazYn3y6lTHQ8XwfcCXwEPCalLBVC/A34UEqZ2orNh4CHGn/v16/fxLNnW5yK9EiKa/TcuGIvd03rT2G1nrunxbXa9rN9Oa2e9y3YS+8fHmSH3zUMmP8Y5sBYPOvyzv0UNjMApuCBZH39EgN/+VKr5z3r8kje8m8mXLno3PkGQzVJlYdZW7mbGmMpT41+heSdK7jqmpeItdqoMJRQb2nAHNgXj/piLP4xeNQX8/mBL3gh8VkqDCX4aLzxCBtJcel+rAZvop57mfIlD5O+7UPGz7udEM8gPILjKfEsJCZyMsaMXVR6euOhLydl8ypG3fsUskaDh66MyPCxGPS56AeNwvNoDgc3fcmQxBloIiIwe2iQfiH4Fp/EJvvT189IodbAsbUHmHHrDZQVnkIGheChrUMXNISgukys1QbCPMPZ+e1uRj5wLdEh8dRn76PSEoi/RykmvyBsWg88LBoy1mzmisX3orPkYw4bTXT1KfJK8xEAkdHYvPzIXLWHR1e+1ObfXghxWEqZcMHxTgjyv2Gff6cC30spE4UQTwNFwFfANuAqKaXRmWskJCRIdeHNTlmdkRtX7OHmCX347dVuzHiq8+H9GXDT+zD0WkV8k1JSYaggry6PLWe3sPnsZkw2E4uHL+aekffgqWl9RHExTHl5nL3zLoJ+sYDo//kft/w0GS1899YRbvztBLQeAqvV6pa9jkCr1aLRXDjTbi3IO2p13RPYCIwFfhBCPA/cKaVcAvwZ+AzwBV50vORD7KvrS4CVzga4ys9U6Iwsen8vc4ZFuhfgZ/fBma3w6x8gapgivlUaKnlpz0vMip3Fjzk/crb2LAsHLeS+Uffh7+nvlm1rdTWFzz1H0HXXuR3gALpKAzf9dgJGswGNVYOHx+V3l9loNGKz2fD3d+6z65B3IKU0A/OaHU5ynMsH5jZrXwv8oiN86QlU6kzc+sE+RvUO5tUbRrluqPwUnPoRgnopFuAWm4V9RfvYlr+NiTETuSbuGibHTKZfUD+3bTckH0Gflkb00qX4jh7ttr28jEqO/JDLL5aMQWPV4Ovr67bNjsDT0xO9Xo/NZmuxR2+OuuOti1OuM3LHh/vpHeLL8tvHI4RwzVDZKVh9K/RJgMkPXby9kxwuOUyDuQGAEeEjWDhooSIBXrdtO/mPPYa1vFyRAAfQagVX3D4Eq9V6WfbgTdFqtU5PJS7vd6LSJiW1eh749DC+Xlr+766JaDQuBvjpLXDoI5j/Bgy5WjH/inRFLDuwjAnREwCYFDPJbZvSYqFmwwZKl71B4DVXE/n0U27bBMjYX4RJb2HM4FDMZrMiNjuS9nyZq0HeRSmobmDJ6iOYrTZW3T8VHy8X/pRSwsmNsP4pGDgb4pvPsFzHZDWx4ugK6s313D70dsJ9wt22adPr0e3YiT4llYB584j53e9cH7k0QUqJrtJIr8HBbtu6HFGDvAuSVabj+bVpVDWYWHX/FEL8XNwgeGYLHPwI+kyA698BJ+Z3zrIxeyObz25m6eSlxIfGEx8a77bNum3bKX//fSIefIDAK69EeLq+It+U7KPl9BocQp/BIYrYu9zo8kFuMpnYuXMnYWFhxMTEcPz4cby8vBg7diwHDx4EYNq0aRw4cACr1cqYMWPIycmhtraWQYMG0dDQQFFRETExMQQEBHD69GkCAgIYNGgQKSkpaDQapk6dyt69ewFISEggPT0dg8HA8OHDKS0tpaKign797PPM3NxcwsPDiYqK4sSJE/j4+DBq1Khz++unT5/O/v37sdlsjB07ljNnzqDT6YiPj0en01FcXEyvXr3w8/PjzJkzBAUFERcXR2pqKlqtll7xo7jl7f+Alz8r75lMfmYqWSYTI0aMoLi4mMrKSvr374/VaiU/P5+IiAgiIiLIyMjA19eXESNGcPjwYYSulBkn/sD+gPlY4uYwTm/mVOpx6uvrGTx4MDU1NZSWltK7d298fHzIysoiODiY2NhY0tPT8fDwICEhgf377TuXJ0+ezJEjRzCbzQTEBvDX/X9lOMMZ7TGaM2fOUFBQQGRkJKGhoZw8eRI/Pz+GDRtGcnIyADNmzGDPnj0ATJgwgYyMDBoaGhgyZAhVVVWUFhQQ8dpr1JaUkKvREH76NH369OHYsWN4enoyfvx4Dhw4AMDUqVM5dOgQFouFUaNGkZeXR01NDQMHDsRgMFBYWEhUVBTBwcFknMgga4uVOYtHsnNn6rn3otfr8fDwwM/PD4PBgM1mw8fHB4vFgsViwcvL69z/n4eHBx4eHhgMBoQQ+Pn5UV9fD4C/vz+33norUkry8vIoLi7mL3/5CzfccAM2mw2z2YynpycajQaj0YhGo8Hb2xu9Xg9AQEAAOp0O4DxftFotWVlZVFZWnvvfaxUpZZd+TJw4UfYk/vurI7L/0g1yx8lS9wzl7JHyxSApK7KUccxBvalePrH1Cblw7UKZX5uvmN3KL7+UJ2fOkseHDlPMppRS5mdWyMJTVecdM5lM0mQyKXaNYcOGyWeffVZKKeWuXbvkpEmT3LbZko/AIdlCjHT5nrwnUVprYH1qEQCz4iNcNyQlrH3Y/jxsgAKe/cyhkkMEeQXx+PjH6RPYRxGb0majft9+gm+6SbEhOkBlUT0/fXKCX/1+imI2m6PX6ykvL+fFF+1bQkaMGEFVVWtpHR2DGuRdiC+SzjI0OpCU/Br3FpwqTsOMJ6GuWDHfANLL09mau5WpvaYyJ3aOYnZNZ8/iP2M6QddeizYwUBGbUkqqiuuZedtgvP06LgzS09MZPHgwPj4+ACQnJzN27NgOu15LqPfJuxD/SSvm5gl9eOLKwe4Z2vIylGXA3BeUcQyoNdVyoOgApQ2lTO41GQ+NcoFT8trrWKuqFQtwsC+2ZR0po99w91f92yIlJYXc3FwMBgP19fW8+OKLPPWUMrf9nEXtybsIafnVFFbr+eWEWAJ83PizWUyQtQPuWquYbzZp49tT3+Lt4c2j4x4lwteNqUQz9OnpNBw8SPTSZxSzaTFbqS5tIKxPAJ7eWsXstkRKSgqLFy8mMTGR2tpann/+eWbMmNGh12yO2pN3EVYfyGX6oAj3Ahzg7B4IiYW+F+QxuMzB4oN8kPoBtcZaRoSPUMwuQMX/fYjf5Ml4x7t/C66R7JRyCjIrGT1bmTWDtkhJSeHXv/41+/fv5/jx49x5550dfs3mtDvIhRB/6ghHVFrHarWx+XgJt06Kdd9YWSYseMt9Ow7MVjPP73qe4WHDWTxiMRqhXL9hqarCVl9P+P2/VsymlJKsI2XkHq/Cy90vTCc4c+YMgwe7Ob1yE1f+In9T3AuVNtl8ogSNEMwdeoFGX/uwmCD1X+CrnCbHd2e+o1RfypMTniTIK0gxuwCG9GME33QjfpMnK2azvtpI7PAwxexdjIKCAqeSSDoSp64uhIgXQlwrhIiXUpY4+Zq3hRC7hBDvNDv+VyHEdsejynHsXiFEpuPYn9v/Nro3aw7lM294FFqtm/8spzbDqF9CpHLqWmtP2ef2oyLdyH5rAWttLYVLl6LxD1Bk62oj+7/LoqGuZ2UyX3S8IoR4BZgOHASeEkLsk1K+eJHXTAD8pZSzhBDvCSEmSSkPAkgpn3S0GQ801bN5U0r5oYvvo9tSZzCzL6uCNf81zX1jBYfBT7leLKcmh8yqTO4afpdiNhupWr0aTWAgATOVW6SyWW2E9/Gn38hwbJaeU83Xma7hKinlPCnlc8A12FVWL8Y04CfH85+AqS20uQn4psnvTwohdgohrnTCfo9hzaF8+ob6MrK3m8kTulKoyoGE+xTxCyCnNofbht7GM5OVW/kGsJlMVP/7a0JuWYRQMOWz4GQVGo2G8N4BTL5+oGJ2L3ecCfI0IUSM43kUcNyJ14RwcWHG+cAmx/Nvsau4/hL4XyFEx97X6EJ8d7SABaN7u28o+TOozgEv95RYGrHarBToCrht6G2K2GuK8dQpfEaOIPRXv1LUbl5GFcaGyz+NVGlaDXIhRJoQIhVIBLKFECeBHGCmE3argcZVmCDH701tDwYKpJQNAFLKaimlTUpZBpwEotvw6yEhxKHGR1lZmRPudE1yK+rJLKlj8RT3RRY4/h2MvNl9Ow7SytNYdXwVvf0V+AJqgrTZaDh4iNA7FqMNCFDMrsVsJSTSj7HzFPgsuxitBrmUcrSUcoyUcoiU0rfJT2dWbfZhL64Adhmo5kUWbgLO7cYQQgQ5fvoCg4FWI1dKuVJKmdD4iIx0c8X5Mubz/WcZFxtCRKC3e4aMdRAzBsYvVsYxoN5cz59m/AkPrbK3oer37qV8+XK8Biq7pz5tRwH5mZV4+yrj79ubT7arfWFhIRMmTDiXyaZ0+7ZwZuEtEnvVk6HYe9llUsrStl4jpUwWQhiEELuAFCnlASHEcmkXcgS7ntsNTV7ylBBiPvYvnWXSrhHX49l4rJglcxW4x7r/ffttM4VunVUZqnj3yLu8M/edizduJxWffEJAYiKeCn95+wV5MXRKzMUbtsDqpNwLjr2z5RTRQT7nHbujjRFXWFgYW7Zs4aabbnLqmu1t3xbOfK19CnwCvAXMwK60Ov9iL5JSPtHs9yVNnl/R7NzLwMtO+NJj2H+mnJoGMzeMVWA47BtsF4ZQiK9Pfo3ZZibKL0oxmwD6jEz0h5Pp/+VqRe3WluupLNAxaYGyo4P24OPjcy5JpZFXXnmFrVu3otFo+Pjjj2lauLOl9q7izMKbn5RyjZSyQEr5FXYpZZUO5suDeVwxJBJvTzfXIItS4eDHMGC2Mo4B32d/z3UDrlPMXiPm/DzC778f3+HDFbVbdLqaugoDHl6Xz3puWloaBQUFbN++nRUrVvD666932LWc6ckrhBBPYZ9XT8NezFClAzGarWzPLOO9xQr0vmY9TF8CCs2d08vTKdAVsGjoIkXsNWIuL6d85f/R989vKGoXQOOhYfJC12+ZPb82zanjz69NI2fZAqdsnjhxgu3bt5OYmAhAr169WLVqFR999BHz58/n2Wefddnf5jjzl18C3ArcBZzCXtZIpQPZkFpEoI8H090Rhmhk7cPwi7fdt+Pgnxn/ZErMFMW3sFatWoUxNRXPi0kZtZPirBqSN53lludcT8h57aYLJZ+fX5t2wfG25uTNGTp0KFdffTXLly8HOCcDddddym8scma4/qmU8q9SykellG9LKesV90LlPL5Ozufqka4tEp1H9m4wVEPcLPdtYVdg3ZG/g0VDlO3FAWo32rdMCIX3eVutkjFXxqJxd0uwm5jNZubNm0dKSgrXXHMNBoOBmJgYEhMTmTNnDp988kmb7ZOSkly+tjM9+VEhxPXYh+tWACmlOmTvICp0RpJzq/jjjQrsBc/YAAPnKjZUz6vLY37cfGb1UeZLo5GGtDSsFRWK2gT7Ntbso6WMnete9l5LPfTza9Pa1XN7enry008/nXdsypQpvPBCy8IdLbV3FWf++gmORyOSZmWOVJTjXwfzCPLxZGCkmxtBrFYI6gPTHlPGMSCpKIkFAxcon1VlMhH9/POYCwsVNVtZ3EBpTh3+oW7uM2gBt9V5LiFt/rWEEJOB94DfSCnnOB5qgHcg/0kvolSJLKnCZDjymT3QFUBn0vFlxpeE+SqbpimlpHbTD/hNSiByyeOK2jbUmZh2c3yHpHo+ddUQxW12FG1ta30HeAkYD7wnhHjyEvnUYymrM3CqRKeMsaIU6D1BsYIJZ2vP8uTEJ+kXqOzCWMP+JGrWrcOjt7LbY00mC9s+z8A/SDl1165KW/8BE6WU1zmyz+YAN14al3ou3yQXMLxXkPtDQasF0r+G1H8q4peUkteSXsPHQ5nNGU2pXvsNATNnoFG4wGB9pZFRs/sSFOmnqN2uSFtBLoUQoUKIMCAQ0Db5XaUD+OFYMVePiHZ/KHh6M8Qrl7GbVJREdm02E6KU2zUHYLPZqN+1m2AFtm4253RyKdFxyt7m66q0FeQW4OsmDwv2/O9/XwK/ehw1ehPHCmu5eYICc+i974KhBmYrs6Fizck1zOozC18PZTc7WkpKCLvnHvymT1fWrslK3rFKAjpgwe0c29q3Qy09PZ3p06cza9Ys7rvvPuwFT1onKSnpXHt3JZxbHSNJKZVTx1e5KN8eKWBQZAAxwW4GUl2pXQHmhncVqY5Sb65nT+EeVly5wm1bzanb9AM+I4YrvjBWU6Zn3FX9CIpQ6Evp0CcXHtuxDAKb7WVoQ5Bj6NCh5+rp3XfffRw6dIhJk1ov5dy/f3+2bt2Kj48PixcvJi0tjdEu1mFXJZkvEzallzBvhAIJH2UZMPZXipU/+u70d0T6RjIxeqIi9hqxGY2Ur1iBJiREUbsA21dnYrXYFLfrDp5Nyjt5e3sTGxvLK6+8QmJiInPnziUnJ+e89jExMecSVDw8PNBqXd93rwb5ZUCDycLRvGoWTVBAcvngShhxo/t2HHyY9iHz4y6adNhuajf9gDYkBL8xYxS1a2wwM2BMBLHDlVOkVYp169YxatQoSktLKSsrcypBJTU1lfLyckaMcF3Pvq1baKFNZJ8aj8UIIUKcMdyGWutLQogUhzLr045jgUKI9UKIPUKIu114H12a9SmF9A31pV+4myvB+irwC4deytTaqjXVUqYv49ahtypi7zzb69cTMFf5LRdnkssozqrBx9/Fmu0tkbkRNjxpf2RutD8ajzc+Njx50Xn6woULSU9Pp0+fPhw/fvxcgsojjzxCbW0tq1atIjExkWXLlgFQWVnJ448/zkcffeSW+23dt1gJvAg0rYoXjv3e+S1tGW1LrdXBb6WUTffsPQh8CfwL2CaE+KeU0tSO99Gl2ZhWzJxhCgzVd78D3sHgp0wvti13GwCRfsoKOEizGb+pUwi65hpl7dokBSeryE4pV9QuQ6+1P5py6ofzjw29ts05udFoxNvbvhAYFBREcXFxmwkqFouFO++8kzfffJOYGPfyGNoarkdLKc8TbZRSHqMN/bUmXEyt9Q0hxE9CiHFN20sprUAKdhWaHoHJbOVgTiW3TOzrniEp7QtBQ5wR03XGnGRn/k5FbDVHt3s3wssbr1gFpidNqCnX03fY5XmHd9OmTcyePZvZs2dTUlLCkiVL2kxQWbNmDQcPHmTp0qUkJiayb98+1y/eUtFyx/L+jlaO72ztNU3avADMdzyfB/yhybkwx8/BwC7H882Aj+P5H4Er2rD9EHCo8dGvX79WyrR3DdYdzZdXvLHVfUOnt0m5YpqUVov7tqSURboiufbUWrk8ebki9pqSvfhOWfL224rb3fDuUZm0PksmrTvjlh2TySRNJlPbjV4Mcusa7tKSj8Ah2ULMtNWT7xRC/EkI4QUghPAWQrwGOPP1Xk0raq3SkcEmpTzlTPvmyG4m5Lg+pYhEd8sfASR/CtEjQaOM+snbh98mtzaXx8cru5/cUlWFIS2NkEVtzvjajVFvoe+wUIZMir40muoK7UG4FLQV5C8D9dhTTbOAZKAO+zz9YrSq1tpEmTWCn9cE9gFXOvTWxwGZzr+FrovVamN/VgU3T3BzqG4xQ85umHivIn7pLXp25u9kZh9n1LfbR83XX+MdPwivvspWFD2ZVAxCEBJ9ibaxznnu0lxHAdqSZLZIKV+TUo6QUg6UUo4E3gcuWmJSSpkMNKq12qRDrdVx+k0hxB5gPXYVWIAPgcXALuBjKWWPKFa1LbMUf28PxsaGuGeoJBWu/D3EKVNS6Pus7wnzCWNCtLLbWAE8YnoR+ZvfKGrTZpOUF+jw8VOTUVrCGUlmX2AhdgkoP2CbM4ZlK2qtUsqHW2hbi12muUexLqWQKwYrMFRPWgkByimnrj+znnn95ylmrxFzcTH6Y+lEPfHExRu3g+qSemIGBBE/SVn12O5CW/fJrxNCrAI2AJFAiJTyWimlWnVUAaSU7D5dwU3u7lVvqIDQOJj124s2dYaS+hKOVRzrkPJHVV99hSE1FY1CUsON7P3mDDVlevervnZT2vpUvgbKgV9KKd8FesQQ+lKx90w5GgFTBrh5yyd7F+TuA98QRfz6KvMrhoUNo3eAsvndAHU/biboWmWlnI0NZvqPCmfkLOX9bYu/H/17u9q3tyLKJUlQAWKwFyD8QgjRAEQLIbx7yny5o/n2SCEzB0e4X3u7OBWufEkRnwD6BPThsXHKSUY1YjhzBnNBAcE3LFTU7rHdhSAhMKzjygGsObnmgmPvpbx3wSahW4a0fsegvRVRLkmCipSyRkr5sZRyAfA49sopm4UQyqjL9XB2nSpnobvVUcpOwv6/Q1icIj6dqDjB2bqzJES7Ll/cGra6OqKWPoM2MFA5mzZJfbWRoAjlxSyUxsfHh9DQ83cidnqCihBinBBioxDiGLAC+E7ayxs95PLVVABIPluJwWwlcYibC0WHP4bYKeCnzC6vtLI06kx1eGqVX6Wu++kn/CdPVtRmRb6O8D4BDBzX9fZKOFtBRYkElbaG6yuAp4Gj2O91vwPcIKXMcvlqKgB8c6SQaYPC0WjcHKrbJCTcr4hPVpsVL60XD41R/ju8/vBhqj79jMglSy7euB0c/E82wZF+Ha6pvj1ve4tbfN869BZPJzzN9rztQNvD9eY4U0GlMUHlq6++csv/toLcKKVsVHT/XgihzPKtCjsyS3nuOjfrfeXsgYZyGOZcWZ6LsfnsZr448QU3xN9w8cbtpOqL1UizGY23ckothnoz/UaEXZJePDE2kcTYxPOOvbLvFZ5OePrc+fZysQoqlypBZYgQ4m+Ox/Kmv7t1xR7OiaIaqhrMXDXcmTyfNtCVKrqN9dvT3zIuapz7C4HNMBUWUr9jh6I2AdK25aGvM+MX1IESTwrS3goqSiaoCNmK1pQQotUymFJK5f9qLpKQkCAPHTrU2W44zasbjpFd3sDH97Yu/XNRrGb46Gq4eSVEuC/yX2OsYd6aefzzF/9kUMggt+01RbdrN8aTmVjrG4j6jTLDdZvVRtK6bGIGBjFgrPI9udlsBs5Xc2nO6E9Hk3ZPy4UQLwUt+SiEOCylvGDVtC2Nt/MCWQhxt5TyM+Xc7Jlszyzj0Tnx7hmpOgsDrlAkwAG+OfUN/YP6Kx7g5pISqtesIfyxR/Edqlz2cEl2LUERPvQfFa6YzfbyyNhHOu3a7aU9Kxb3dpQTPYXsMh2F1QauG9XLPUM//QEC3BzuN2Fj9kbmD1Be4qnio48w5Z7FZ5CyXx5pOwqoKmno1CKGj457tNOu3V7a8ymd6TAveghfJxcwsX8Ivl5uzKMbKuHMVhikjHRSVnUW2TXZLBqsbKVSS1UVtevWE3b33QgFCyfodSb6DAlh4vz+itns7jgd5FLKBwGEEB90nDvdm59OlHCtu7340S/AKxCihiniU52pjuenPE+IT4gi9hqpXb8ebVgYwb9QNu/oyI+51FUa8A1QUMOtm+PKeMepjPw2hBxfFELsczyudBy7VwiR6RB37JYJMKdL6zhVqnNvl5uxzv6zvlQRn2w2G2tPr2VspDLCj+fsNjQgfH3p9cc/IryUC0ar1Yant4beg0MUs9kTcCXI2y79wPlCjoCXEKLpUvJnUsppwLWcL0DxppQyUUr5jAs+XdbYbJIPd2VjtUkCfd3YTZaxEZKUG0jl6fIo1BUSG6is1lrl6i9pSErCZ6Tru7RaIvdYBcVZtZeFjlvZ8nfb1b69FVTa274tWp0sCSHWc2FAC8CZr/2WhBwPAkgpsx3Hjc3sP+mQY35ZSrnFiWt0GfacKWdDapF7RspOQvYOWPAXWK2MRHJ+XT4Pj31Y0W2sVoOB6tWrCbntNkU3v5j0Fk4fKiX3WKX7OwXbSdW/LtxxVr5iBR5R529LDr2t9b9LeyuotLd9W7S1IuKOuFcIPy/U1QAjW2jzEtDYLX2LPQEmHPhRCJEg7cqtFyCEeIgm++f79VO2lK7SVNWbeHNTJhP7h7LjZJlrRqwWyNoO9WX2veoK6ItVGip5YtsTrF6w2m1bTalZswZpsRC6+A5F7RacrKL34BBOHihR1O6lorUKKlu3bkWj0fDxxx8TFxfXZntXaWu4/qGU8mxLDyfsVtOGMKMQ4iYgXEq5GkBKWS2ltEkpy4CTtCH73JWEHG02yfs7TlNYo2fp/KGulyTO2ADHvoHEpfa8cTf1xSw2C5+kf4LRamRIqJsVVJsgTSb0KamE3LIIbUCAYnarSxvY9dUpNB4aJi2IU8zupaa9FVSatg8Pd31PQFs9uTv7JfcBDwNfYU9u+UfjCSHEGOAxYEGTY0FSylqH1NRgwMUu7/Jif1YFqw/ksWROPCN6BzOid3D7jdRXwPbXIXoU9BqniF9bcrfw9amvFbHVFENmJn6TJxF0nXLCEFarjZLsWsbM6cugcZF4+Spbx9wZ6rZto3779nO/hz9sVzArfvHnJSX/xMQ2h+tgr6CycOFClixZcl4FFWg5QaVp+w0bNjidi96ctj6xMUKIdc2OCUBKKdvM/JdSJgshGoUcUxqFHB06b29i76l/EELUSClvAJ4SQszHPrJYJqU0u/RuLiNqGsy8+UMmY/uGsHiqG/d0j35hL0M893eK7FPPqMhg+ZHlTIiawI585XYnS4uFot//geCF1yvai2cdKSNzfxFTbxzUKQEOEDhnDoFzLizyG/Pyy07baG8FlebtfX1dF8Vo61M7Bbi82bgNIccLauNIKV/GLgHdLZBSsnLXGc5WNvCP+ybh7+3iP2fhEXt9swVvKVKltFxfztuH30YrtDwx8QlGhCu3+l27+SfMBQUEKXhf3KAzU19jJGpAEBGxyolNdAabNm3irbfeAmDw4MG8+uqrvP766yQmJiKE4Pbbb+ehhx5qtf3VV7teGaetBJVtsoUa5UKID1pSXO0sLscElQPZlTzw6UEenj2Qx+a4OA83G2D321BwCG79DLz83fLJYDHwY86PbMvbxvWDrmduP+WKDUqbjZxf3Y7P0KH0evUVxewe/D4bIWDYtF4EhF469RdnElRODBvO8IwTl8qlC2hPgkpbC29vtXL8EpSn6LrUGsz8eVMGI3sHc890N3rftDVQcQrmvex2gEspWXd6HRlVGSwYuIDZfVtNMHQJY3Y2HpGRhD/0oGI2q4rrqS5pwNPb45IGuLNEPKa8Dl5H0ZbG2/rWTnWQL10eKSUf7crmdJmOZ+YPJcDVYXpdCehK7EP06JbuPraPtPI08nR5NJgbmNZ7GlqFctABrLp6Sl5fRsgtixQrYGgxWzmbVkHfYaGMuMQqrM4SuUTZ8lEdSUdthumRHMmt5tO9Odw3I47x/VwsH2yzwfbX7OmkN68EN0Uc8uvy+ebUN4wIH8EVfa/A39O9UUFTLJWV1KxdS8Pu3cS88IJidgtOVpGyNY/5D4/C051kHhWg4zbD9DjqDGbe2JTBkOhAfj3TxWG6zQbHv4Ojq2HB225XRakz1bH8yHKSS5O5If4GYvzdkxFqiim/gJpvvqZy1ecAeA+IU8SurtpAXYWBGYviieofdPEXqFyUtobrLW6EcXIzTI/j0705ZBTX8cz8oQT6uLBN1Gqxp5Bm7wCrCca4V/XTYrOw+sRqduTv4NGxjzI+arxb9ppiyDxJ/d491G3fQdACZTTmAKRNkvxDLllHy+k9OFRxKaqeilpXRgFS86r5aHc2d0/rz8T+LgzTrWY49aN900ugIxXVw/V93yariU3Zm1iZtpJbh9zKgoHKBKKUkobkI1R89BHW2jpC77id8AfuV2wRqqKwntAYP0Zd0Qe/oMs7lfTA+vaJFre3gsolSVBRcY56o4U3fshgQIQ/D8wa2P7ex2KEY2uhINku6dQ3wa296eX6cnbk7eDj9I8xWU3cN+o+vLTuB4y0WNDt2UvlJ59gPH0an+HDCJwzB4/wcEUWoeoqDRxYn0XMoGD6j+w8WaeWOLar4IJjB7/PwT/k/C/ikbNar2vX3goqlypBRcUJPt9/lrT8GlbePZHg9qaRmhrgxHqozIaqHJjxBPSfDvFXXvSlLZFZmUl6eTopZSlE+UWRW5dLqI+LC4BNsBkM6HbtouKDldh0OqKeepLAa+ajDXBvEc9mtVFdqqeysB5DvQldlYHYYQPQena/AaaPj8+5iiiNXA4JKioX4VhhDSt3ZnHHlP5MGdDO3sdYB0c+hx9fAHMDzH7GXl/chXmo0Wpkd8FuVmesJqUshdGRo3ligjLlga01NdRs+B5rZRV+U6cS+dunCbr+ercC3GyyUpJdy8mDJRSdriZjfxHFWbWU5eqI7Ne1d7Y5izMVVJRKUFGD3EX0JitvbMygb6gvj8we1L5hur7aHuDbX7P33KNvsQ/TXaBcX866M+vIrc0l2i+afoH9mBQ9iXFR49xSFLXp9eiPHiX/N7+h7H//F3NBAQEzZxJ45ZVoXFR70deZyM+o5NSBEkpyajj0nxxsVhu940OIGXj5rqRnp5aTnVoOQNK6LBpqTQAc2pgDwOnDpeSdqGyXzaYVVB555BFqa2tZtWoViYmJLFu2DLAntKSnp9OnTx82bNjgsv/qcN1Fvkg6y5Hcat6/ayLBfu0YptdXwOF/gNYTRv0SJj3osl7biYoT/JjzI/8+9W/m9pvLjN4zmNp7KkFe9oBxRVHUqtNhPHECY3Y29bv30JB0gOCbbiLw6qvxHdX+jTkmvQVdlZGD/8mm79BQSnPrMDVY6DsslMEJUXj6eBAzMIDQGD8aakzttn8pGDAm4tzzKQt/3vCZcG0cAPET23+r82IVVC5VgopKC1isNo7kVvPW5pMsntKPaQPbMYyqK7Gvoufth/B4mPoohLdfrthgMZBUlMR3Z75jT8EeFgxcwHUDrmNc1Dg8NK79Sa3V1RhOnEB//DiWsnK8BwzAZrRXqQ5/6EG8Bzh3799qtqGrNlJfbURXbcBssCKl5PShUqLjgggK80GGSXwDvegzNJTgCF+EQ+ll8vXdd8e02Wzm2muvPVdB5bXXXjtXQaXTElTcRQjxNpAAJDfNSBNC9AY+B3yAP0gpfxJCBAKrgTDgg/YUcbhUCSo2m+RMmY6k7ArWpRRxILuSI7+/ilD/iwxdjTqoKwZdMez8C1SchllP2xfXQn5Wtakz1ZFbm8vZ2rPndNB9tD54ab3w8bD/9NZ6oxVakkuS+erkV+hMOm4behvz+s+jX5BrCjmWsjIMGRnoU1Kp37MHw4kT+E2dal85j4yg4eBBopcubfX10iZpqDPZg7rKSEOtkYZaE/7B3pTm1FJbYSB+YhTbv8jkil8NISTGj/DeAZf1LTJnElRW/NdWHntfuSSf9qJIBRV3aCrkKIR4TwgxSUp50HH6WeB3QCqwAbsG3IPAl8C/gG1CiH9KKS+bsduxghq+PJjHwexKcirqiY+y50u3GOCNQV2TD3VFUJkFVdn2vejZO+Ga12DIfAjqhdFqJLcml5y6HHJqckgqTqJQV0heXR75unz8PP3w8/BjVt9ZWCwmaiuKiNMH8vfcfzCk7zgemPgAM/vMJNCrfYtV1poaLGVllK74O77DhqENCaF+3z60ISFEPvkEPsOG4T10KB6hoQTOvfAf2Wq1oas0UlehR1dl5HRyKaExfgghCAr3IWNfMWOvjCWstz9hfQLw9LFvTR02vVe32abalRRqOmq43qqQIzAGeEJKKYUQdY5efBrwmJTSKoRIAYYCnVdoCmgwWVhzOJ9vDueTUVzHgAh/xsYGc/e0/vQJ9eXeTw6ClPaeOWc3lBwDnxCozgGNpz3QfUNB6wHVeeBvn9eVDZnHrpJ9lGeXE+4bzifHPmFe/3kEeAYQ4BlAYmwiq46vIiEmAb2uBm1eMVEphzCUFGE2VxIwKIE3/mHEf6Yn4UcOY+pdSP2YsXj2ikEbHo62hbmbubgY45kzWEpLERoN+vRjmAsK0G3diuHIEcIffpiQRYvwHjQQ7/h4NH5+F34etSYqC+uxmK1kHSkjMNwHT28tqdvyqaswEBzhC0IiNIKE6+LwC/IiMMyHgFBvvP08qS5u6DYBDp0/tZBSOr3Y21FBHkLrQo5a+fMcoQYIdbSvbXZMEf726H8TXVxLUM1uSqLuZM+QfRi9LMw5NoeheWspDBuPPsrEuL6p7D9zGxEl3yKFN6VRc/h6xkauOTaR/yooZ0jtUXb3fYjBU05RU2nllexINl31DFUhV5HXt4iU0aeZu/MBEoq3UO3Vi6yo8Uy5IZak7/MJLqrEx3Ccwl5/5uPv7uCqrAmMPh1NvO5bHgi+gpCJQfgNj8H/XxrGVnzLaNsEZKSGiEmx5O+NJqh0P6GeZhqif4XNUoo+IIGGqn54Zx0mKXA2g0uPIsxGThypY9TZf1PYdxbe4ycwcMYg9nyYxKBTX2MLDCVv4LVMSIzidHEIXgYr4UDqzOeZEtsbo284OUcqSegFxzaeIjjaj75DQ9m0Mp2J18VRmFlNZVE9gxN+XmTy8NIwYkYvR/HBYAIcQe0f7HVBCaPODor2oNVqMRqNbQ7XOxur1XpuYe5idFSQV9O6kGNTFdbGc43tDS20P4/2qrUOnz+Thuw0pPeNBJoDmDf4FqRG4Blch3/4EkLqjYR4G9HG/5Lwo7kIr4VobJKB3uG8Pm0pxn51hFVU0it4HvEVkuHzHkFfXYvPkDP4hywlsN6TuD7hzB0Qhs43i75xrxJusOBvMjJi/jWYArZh1A3Cw8tGROoxdt29n/r8cnSnzhIScgPhlVZCRgzANyqMSOshoqNHE63xh8BQQvsG039ONf5+d6CVJnrn6YiMEGw6+2dm3HMrwZH3ElBuJTTGH+HlSViJgYjI+4nIq0QT4E/YkL5M79Wf8L73YzWY6F1ST3iEFt/R9WhsFszxWqbfNomo/kFYTFaCo/0JjvRjxMw+ePpo8Q/2Zs6dw4joG0D8+CikTWKzSYZOicFmk9isEiklZqOVwZOUq83W2Wg0Gmw2G3q9Hq1We1ntoZdSYrVasdlsaDTO3QHvkIU3x5z8YSnlw0KIvwP/kFIecJz7G/b5dyrwvZQyUQjxNFCEXfhxG3CVlNLozLUuR2UYle6BzWbDam1RGbxT0Wq1LQb4JV14u4iQ45+xa6z78nMFlQ+xr64vAVY6G+AqKh2JRqNxure8nOmw++RtCDnmA3ObnasFlK2Mp6KiAqjbWlVUuj0dthnmUiGEKAMuJmQRAZQrcLnGe0sNLrZr6XhLx5Tyty1fnG3T/Nzl4G9b7S71Z3w5+dtfSnlBSaEuH+TOIIQ41NKChAt2JgJIKQ+70q6l460cU8RfZ31uq03zc5eDv221u9SfcVfwVx2uq6h0c9Se/DJE9bfj6Wo+qz35xVnZ2Q60E9Xfjqer+eyyvz2iJ1dR6cn0lJ5cRaXHoga5iko3Rw1yFZVujhrkKirdHDXIVVS6OT0uyIUQY4QQ3wohXJNIvcQIIeYJIZ4TQrzT2b44ixBiihDiv4UQf+xsX5xFCLFICPFeZ/vhDEKIJ4UQLwghbnamfbcJciFEbyFEY4qrR5PjbwshdjUGiZQyFfi2s/xspB3+/iSlfB3odO2kdvichL3MdXBn+erwyyl/hRCjgXrsqkSdhrP+AlXYY9enRUPN6DZBDlQCVwL7Gw80FZQEvIQQrhWT6hic9lcI8RR2hdvOxmmfpZRvAvmd4uXPOOvvXOy6ggkONeHOwil/pZSfSilfBcYJJ2Rruo3uupTSABiavecLBCWFEKXA1cAQIcTLnSVQ0Q5/JwCTACmESJKduHupHT7HAiOATq1c6Ky/UsrGHj1GSll4ab38mXZ8vjHYBVFNzvw/dJsgb4UQmglKOuqr39FpHrVNCBf6+wrwQad5dHFCuNDn5cA3neZR24TQisiolNL1crIdRwgtf77rnTXQnYbrLVFN64KSlyPVdC1/oev5XE0P87e7B/k+7HMcgHk0metcpnQ1f6Hr+dzj/O02QS6E8BRC/ASMBX4QQkyRUiZjn+PsAmyNirGXA13NX+h6Pqv+OuyqWWgqKt2bbtOTq6iotIwa5Coq3Rw1yFVUujlqkKuodHPUIFdR6eaoQa6i0s1Rg7yHIux8J4TwEkJ8dYmu+ZIQYtGluJbKz6hB3nMZDJzEnuiQ3sm+qHQg3T1BRaUFhBBvArMAC/Y0y1ohhAEYgj3lUgLPAZnYS0p7OV56q5SyVAhxDNgCTHH8tDrs5EopbxdCJAL/A5iAAcAfpJTrmvnwFHAz9v/BP0kpNwghXgauAszAcinlvzvmE+hhSCnVRw98AJ8CYcC/sNeK9wQOAxrHeY3jmKfj94eBZx3Ps7GnkgrgFHC14/gW7EGdCBzBHsAhQLrj/EvAIsdr/+045gMccjxPBrwar9/Zn1F3eag9eQ9DCDEbeBkYBfTHnmq5EXsQLwM+FULogVew96h/F0JEYFd52eUw0yClPO6wV4A9oAHysH9xgD1wLUC1EEInhPBv4sZIYLwQYrvj9wAhRCD23n+lEEIC/wscU/TN91DUOXkPQ0q5A7gd+Ap4CPhYSpkInAa+k1LeBewAHgUWAzuklLOxl+lpVDNonvDQ9PfGNhOFEFohRDAQIKWsb9LmBLBfSpnouPYYKWUdsFdKeS/wIfCCEu9XRZ2T91QmAYewz6kbs5oCgW8dvagX8Dj2ufbnQoirgQLsPbuzFANrgX40C1gpZboQ4oAQYqfjGmeBe4G1Qghv7EP4F117ayrNUbPQVBTHsfC2SEr5eCe7ooI6XFdR6faoPbmKSjdH7clVVLo5apCrqHRz1CBXUenmqEGuotLNUYNcRaWbowa5iko35/8B/Nz0zAYpbokAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x151.2 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJoAAACRCAYAAAA7OKHUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAP5UlEQVR4nO2dfWxV533HP1/Ar+CYYQcI9mwDWSJiSLbgvCBCIC+qmqhpM6St2ppt6qZmL1mbSVO3ttMarV2VRpO2ru2ahWqtujTRtmpJukbJQlBD3khDSLSFEJRxjzHgBBswtsH4BWz/9sc511wu9vW5L+dc38vzkax77uPnPOdn8eV5nnPO9/k9MjMcjqiZV+wAHJcGTmiOWHBCc8SCE5ojFpzQHLHghOaIhUiEJmmFpHckjUpaIGmlpFclvSLpSUnzg3ofSNoZ/FwTlN0u6Q1JL0lqDsrWSnpN0uuSro0iZke0KIrnaJKqgRrgaeBOoA4wMxuQ9A3gF2b2M0mvmdktaee+BHwSuAb4XTN7QNLTwBeASeB7ZvapMHE0NjZaW1tbwf4ux3nefvvtE2Z2edj6C6IIwsxGgVFJye/9Kb8eByaC4yWSXgH2Aw/i97AjZnYaeFPSN5P1zOwIgKT6sHG0tbWxZ8+evP4Wx/RIOpRN/VjnaJJW4Pdw24OiW8zsVuAQcD/wS8CplFPmB5+pcWaMWdL9kvZI2nP8+PHCBO7Im9iEJqkK+BHwOTMbBzCzk8GvnwbWAv3AZSmnTaZ9ph9fhJltM7MOM+u4/PLQPbsjYiIZOmdgG/786n0ASZX4c8QxYCPgmdmwpBpJi/DnaO8H554MbgwmgcEYY3YUiEiEJqkCeB64DnhB0teArUCrpAeBfwR2Ac9LGsLvye4LTv8G8CIwCvxeUPYQ8G+AgAeiiNkxM8PDw+zfv581a9ZQW1ubUxtR3Qycw5+LpVI3TdXrpzl3B7Ajrexd4Jb0uo542LVrF7t376a1tTVnobkHto6MjI2NsXfvXsbGxvJqxwnNkZGDBw9iZjQ2NubVjhOaIyOe57Fs2TKqqqryaifOu05HiTE8PMxHH33E9df7U+lc52fghObIQGdnJwDt7e1UV1fn1ZYbOh0z4nkezc3NeYsMnNAcMzA0NERvby+rV68uSHtOaI5p8TyP+fPnUyj3ixOaY1oSiQQtLS1UVlYWpD0nNMdFDAwM0NfXV7BhE2Jy2AZlXwxcsk8E70KR9BlJuyQ9K+myoMw5bIuM53lUVFTQ0tJSsDaj6tFOAncAvwCQdDlwW+CmfRe4NxDbHwG3Ao8Dfxic+9fAx4AvAV8Oyr4O/Bbwm8GxI0I8z6O1tZUFCwr39CsSoZnZaJqr9kZgZ3C8A7gZuArYG3jTdgA3S6olcNia2Zv4ViEIHLZm9iEQ2mHryJ4TJ04wMDDAlVdeWdB245qjLea8c3YQ30k7XVneDltHfnieR1VVFc3NzQVtN65/tAHOO2cvC75PV5a3w9ZZuXPHzPA8j5UrVzJvXmGlEZfQ3gI2B8d34s/d/g9YGyy9uxN/ZdQwUCNpkaQbSXPYBmsOMjpsnZU7d44dO8bQ0FDBh02IyWELfAV4RdJrwGHgW2Z2TtL3gVfxe7LfDk53DtsikUgkqK2t5Yorrih423E6bN8EHkmr9zj+HWdqmXPYFoHJyUk6OztZvXo1yWWShcRNrB0AHD16lJGRkUiGTXBCcwR4nkddXR1Lly6NpH0nNAcTExNTw2ZUOKE56O7u5uzZs8UXmqStwec1kn4sKX2i7yhhPM9j8eLFNDQ0RHaNsD3a54PPPwe+B/xNNOE44mZ8fJyurq7IbgKShBVaZeDCmGdmu4D8Fvk55gyHDx9mfHw80mETwgttB/7T/H8Ncp85oZUJiUSCxsZG6uuj9SqEEpqZPWRmHcDuwJlxV6RROWLh7NmzHD58OPJhE8LfDNwt6S1gd5Aq9AcRx+WIga6uLiYnJ1m1alXk1wo7dP4V/iugY4F/rC2yiByxkUgkWL58OYsWLYr8WmGFNhHkMUsmvJ2fqfJ0SPp4SmLko5LulTSYUrYkqBfK3u3Ij9HRUT788MPIbwKShBXaU5J+ArRIehz4z2wvZGb/bWZbzGwLvoNjB77DdkvwczJLe7cjD5Kr0OMYNiGke8PMviXpBaAd2G9m+3K9oKRVQK+ZDUlaI+lV4HV8AU3ZuyXtALal2ru5MIGyIw8SiQRNTU3U1NTEcr2MQku+EUhhErha0tVm9lSO19yKn7MW4FfwvWj/DNwD9BHe3j1TzPfjJ14u6CqecuLMmTP09PSwefPm2SsXiNl6tHUzlBuQq9DuwRfbVLJkSc8Avwb8lPD27ukDM9uGny+Xjo4OtxnpNHiex7x581i5cmVs18woNDO74FWTpBozG8n1YpKWA2fNrE/SQmDUzCbwkyXvZQZ79wwJlB054nleQVehhyHUHE3SJ4CHgdPBP/hXzOzZHK73KfxeC/xh8weSzgCdwENmNpGFvduRA4ODgxw/fpxrr413HXZYK/fX8DefGAx2LnkJyFpoZvZYyvH/MH2y5FD2bkdueJ7HggULaG1tjfW6YR9vdHF+Qn4KOBJJNI7IiWIVehhCuzeAPZIew186t0DStyV9O7rQHIXm5MmT9Pf3x/JuM52wsv67lOMnowjEET2JRILKysqCr0IPQ9gHti9LagWW4K+txMzeiTIwR+FJrkKfPz/rN4h5E/au8wlgKdAdFBnw+1EF5Sg8x44d4/Tp07G920wn7NDZmr6Bq6O08DyPmpoaVqxYUZTrh70Z+LmkOyQ1SFqSdFo4SoNk8pZVq1YVPHlLWML2aPOBR/GHTuEPnbdHFZSjsBw9epTh4eGiDZsQXmi3m9lVkUbiiAzP81i0aBHLli0rWgxh+9Hdkq6ONBJHJESdvCUsYYV2B/4mru9J2ivp3WwvJKlNUm/gpt0elOWcQNkRju7ubsbGxoo6bEL4VVBrzWxV8LnOzHJ9I/ti4Kb9WAESKDtC4Hke9fX1eW+DmC9hn6M14j83a+L8A9sv5HC92wJH7VP4lqCdQfkOfKfG+ziHbcFIrkJft24mW2F8hB06/wPfhLgJ2A9U5HCto/hW7dvwvWYd5JdAeVpcDtvzHDlyhHPnzhXl3WY6oR+qBM7VQTN7FL9nywozGzOzM8FyvWeBBPklUJ4xTpfD1ieRSNDQ0MDixYuLHUpooZ0LXK99kv4UyNrMJKku5etGfKHlk0DZkYHkKvRi3wQkCfscLel0/RPgM8B9OVxrk6Sv4+fteM3M3pSUTwJlRwYOHTrExMREyQntKjN7I8go1EgOCfzM7DngubSyR8gxgbIjM8m90Ovq6mavHANhBfNw8PlVYB/w3WjCcRSC0dFRuru750xvBuGFlqxXa2ZPAuciisdRAPbt20dXVxfLly8vdihThBXaB0E2oZ8GD1WL9y7DkZGenh5ef/11hoaGih3KBYSdo30A/DGwEpgA7o0qIEfu9PT08Pzzz7NkyRIqKnJ51Bkds6VE+HfgPWAL8B3gMTO7nVn2Y3LET1JkDQ0NbN68eWq7nbnCbEPnffh3ivXAj/CfcT0s6dcjj8wRmlSR3XXXXdTX17N+/fqSEtpPgE8Dw/hD5z7gX4C58xdc4qSLbK4NmUkyCs3M7gW+j/9u8+/x01b9Gbm963QUmFIRGYS46zSzA8CrZvZZ/B7tb/HTSzmKSCmJDMKv6/yL4PAhM+sBfhZdSI7Z6O3tLSmRQZavkszslVwvJOmmwDn7qqR/CMpcDtss6e3t5bnnnispkUG8m44dwl/ksglYKmkdLodtVpSqyCBGoZlZj5mNBl/H8R/8rgl6uG/KXzkxlcMW/yX6zakOWzN7Ez8Z3yVHKYsMirCNoqRrgUYzex8/Gd+t+C7ae3AO22kpdZFBzEIL5mHfBf4A/By2ZmbAM8BanMP2IspBZBCj0AIv24+BL5pZj6SFgZMWfMeth3PYXkC5iAzCv1QvBL8B3AA8Eixk/TLwTy6H7fSUk8gA5I9c5UlHR4ft2bOn2GFkTSmITNLbwY6HoXB7qs8xSkFkueCENodIimzJkiVlJTJwQpszpIrs7rvvLiuRQbw3A45pGBgYYN++fWzfvp329vayFBk4oRWF/v5+Ojs76ezspL+/n5GREYaGhti0aVNZigyc0GLj5MmTU+IaGBigoqKClpYWOjo6aGho4MCBA9TX1xc7zMhwQouQvr4+Ojs7OXjw4JS4WltbufHGG2lubr5g95L169cXMdLocUIrMCdOnJgS1+DgIJWVlbS2tnLTTTfR3NxclBz/cwEntAJw/PjxKXGdOnWKyspK2tra2LBhA01NTZesuFJxQsuSiYmJqcl7V1cXBw8e5PTp01RVVdHW1sbGjRtpamoqWpr1uUrJCC1w5XYA75jZg4Vse3x8nJGREUZGRhgeHp46nq7s7NmzgJ8Wqq+vjw0bNrBp0yZWrFjhxJWBkhCapOuBhWa2SdKjkm4ws7eyaePAgQOcOnVqWjGdO3dxKpHq6mpqamqoqalh4cKFNDY2UltbO1UG/kPW9vb2ObV+cq5SEkIDNnA+bdUO4Gb87RxDs3fvXs6cOTMllLq6OpYtWzb1PflTW1tLdXV1qN6ppaUl6z/kUqVUhLYY368GvvO2faaKku4H7ocLhbB169boonPMSqlMKga42Hk7LeXmsC0XSkVob+BvqgHn8906SoiSEFqwCe1osEfBpJntLnZMjuwoa4etpOP460mTNAInCniJQrY3l2Obrr1WMws9NylroaUjaU829uM425vLsRWivZIYOh2ljxOaIxYuNaFtm8PtzeXY8m7vkpqjOYrHpdajOYqEE5ojFpzQHLFQVkKTtELSO5JGg6QyOdeXtDYlQ+UPg/xtmdq6KKNlnrHN1l628WWsn0N8WV0fMyubH6AaP5/aTmBBPvWBipTjHwI3zNLWcqA6OH4CWJdnbLO1l218GevnEF9W1y+rHs3MRs2sP7VM0leDHLk/l9Q2W/2U36W6IceAI7O0NV1Gy3xim629bOO7qH6e8WVsL52yElo6QZ7cJjPbAjxAlvlvJX1S0nvAUuDyMG2lZbTMO7ZM7WUbX1r9jCn8w8SXTXtlLTRgDbBF0k7gUeAySb8T/C/90mwnm9l/mdla4EP83LkZ20rPaJlvbLO1l218afU/kW982bRXKg7bXPkA2G5mnweQVBF0+Y9nPg0kVZnZWPD1FP6caca20jNa5hvbbO3lEF96/ZE848uqvbLq0SRVSNoBXAe8gD/B7Qn+F74EfDZTfUk3pfz645JelvQysAx/d78Z2+LCjJY7JW3IJ7bZ2sshvvT62/OML2N76bhXUI5YKKsezTF3cUJzxIITmiMWnNAcseCE5ogFJzRHLDihOWKh3N8MlASS/hfYD/wy0G1mny5ySAXH9WhFRlIV0AT8pZltxN80d3WRwyo4TmjFZy3wjJklV9RPAsNFjCcSnNCKz3WAACStB86Y2dHihlR43Byt+PwqMCzpReAc8LnihhMNTmjFZy1wd4qbtixxQ2fxWVjuIgNnE3LEhOvRHLHghOaIBSc0Ryw4oTliwQnNEQtOaI5YcEJzxIITmiMW/h80F76ZEdgfHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 159.624x151.2 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ac29346a02a1473a8af59c4ad8b084bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "<zeta-acc | Q > Iteration over p:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cee058482bb84adb93750fe479865a9b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pathlib.Path(\"plots\").mkdir(exist_ok=True)\n",
    "\n",
    "env = \"rl_ltl_pac_paper\"\n",
    "reward_types = [\"zeta-reach\", \"zeta-acc\", \"zeta-discount\", \"multi-discount\", \"reward-on-acc\"]\n",
    "rl_algos = [\"Q\", \"DQ\", \"SL\"]\n",
    "# Comment out below to run for only multi-discount and Q-learning\n",
    "# reward_types = [\"multi-discount\"]\n",
    "# rl_algos = [\"Q\"]\n",
    "for rl_algo, reward_type in itertools.product(rl_algos, reward_types):\n",
    "    probs, errs = collect_pac_probs(env, rl_algo=rl_algo, reward_type=reward_type, skip_after_stable=3)\n",
    "    plot_pac_probs(env, probs, errs, rl_algo=rl_algo, reward_type=reward_type, title=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counterexample MDP Construction <a class=\"anchor\" id=\"mdp-construction\"></a>\n",
    "\n",
    "Below we implement an algorithm to construct two counterexample MDPs for any LTL formula not in `Finitary`.\n",
    "\n",
    "The algorithm is derived from the proof of existence of counterexample MDPs in Section 4.4 and Appendix \n",
    "D.\n",
    "\n",
    "Note that the construction of a counterexample MDPs for an LTL formula may not be unique, and our algorithm only returns one valid counterexample. In particular, it does not guarantee a minimal construction (e.g., the number of the states in the MDPs is minimal).\n",
    "\n",
    "In particular, the algorithm below constructs a different counterexample for the formula $\\mathsf{F} h$ than the one used in the empirical justification section (Section 5) of the paper. \n",
    "This is, again, because our algorith below does not return a minimal construction, while in our evaluation, we manually construct a minimal counterexample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "AP = Tuple[bool]\n",
    "\n",
    "@dataclasses.dataclass\n",
    "class TrapChain:\n",
    "    prefix: AP\n",
    "    cycle: AP\n",
    "        \n",
    "@dataclasses.dataclass\n",
    "class RegularState:\n",
    "    label: AP\n",
    "    trap: TrapChain\n",
    "\n",
    "@dataclasses.dataclass\n",
    "class ConstructedMDP:\n",
    "    aps: Tuple[str]\n",
    "    g_prefix: List[AP]\n",
    "    g_cycle: List[AP]\n",
    "    h_prefix: List[AP]\n",
    "    h_cycle: List[AP]\n",
    "    q_prefix: List[AP]\n",
    "    q_cycle: List[AP]\n",
    "    is_accepting_example: bool\n",
    "\n",
    "def _collect_formula_aps(formula: spot.formula):\n",
    "    \"\"\"Collects the set of atomic propositions in formula.\"\"\"\n",
    "    aps = set()\n",
    "    formula.traverse(lambda f: aps.add(f) if f.kind() == spot.op_ap else None)\n",
    "    return aps\n",
    "\n",
    "\n",
    "def _bdd_to_tuple(aps, bdd_dict, bdd):\n",
    "    \"\"\"Converts a buddy.bdd to a Python tuple.\"\"\"\n",
    "    ret = []\n",
    "    for ap in aps:\n",
    "        v = buddy.bdd_ithvar(bdd_dict.varnum(ap))\n",
    "        ret.append(bool(buddy.bdd_implies(bdd, v))) \n",
    "    return tuple(ret)\n",
    "\n",
    "\n",
    "def _rejecting_run_starting_with(aps: Tuple[str], aut, state):\n",
    "    \"\"\"Computes a trap chain by a rejecting run in `aut` starting from `state`.\n",
    "    \n",
    "    Args:\n",
    "        aps: the atomic propositions.\n",
    "        aut: the automaton to compute the rejecting run on.\n",
    "        state: the state in `aut` from which the run should start.\n",
    "    \n",
    "    Returns:\n",
    "        a rejecting run.\n",
    "    \"\"\"\n",
    "    current_init_state = aut.get_init_state()\n",
    "    aut.set_init_state(aut.state_number(state))\n",
    "    c_aut = spot.complement(aut)\n",
    "    rej_run = c_aut.accepting_run().reduce()\n",
    "    aut.set_init_state(aut.state_number(current_init_state))\n",
    "    assert rej_run is not None, f\"Must exists at least one rejecting run from state {state}\"\n",
    "    return rej_run\n",
    "\n",
    "\n",
    "def to_dra(formula: Union[str, spot.formula]):\n",
    "    \"\"\"Converts a LTL formula into a Deterministic Rabin Automaton.\"\"\"\n",
    "    if isinstance(formula, str):\n",
    "        formula = spot.formula(formula)\n",
    "    aut = formula.translate(\"deterministic\", \"complete\", \"parity\", \"sbacc\")\n",
    "    aut = spot.to_generalized_rabin(aut, False)\n",
    "    return aut\n",
    "\n",
    "\n",
    "def _construct_mdp(formula: Union[str, spot.formula], is_accepting_example=True):\n",
    "    \"\"\"Constructs an counter-example MDP for a particular LTL formula. \n",
    "    \n",
    "    One can show that the constructed MDP has no sample efficient RL algorithm.\n",
    "    \n",
    "    Args:\n",
    "        formula: the LTL formula.\n",
    "        is_accepting_example: whether try to create an accepting or rejecting \n",
    "          example.\n",
    "\n",
    "    Returns:\n",
    "        a ConstructedMDP, or None if the example does not exist.\n",
    "    \"\"\"\n",
    "    if isinstance(formula, str):\n",
    "        formula = spot.formula(formula)\n",
    "\n",
    "    if not is_accepting_example:\n",
    "      formula = spot.formula_Not(formula)\n",
    "\n",
    "    # aut = to_dra(formula)\n",
    "    aut = formula.translate(\"deterministic\", \"complete\", \"parity\", \"sbacc\")\n",
    "    terminal_aut = spot.decompose_scc(aut, \"t\")\n",
    "    if terminal_aut is None:\n",
    "        # Case pure G\n",
    "        weakly_accepting_run = aut.accepting_run()\n",
    "    else:\n",
    "        weakly_accepting_run = aut.intersecting_run(spot.complement(terminal_aut))\n",
    "    \n",
    "    if weakly_accepting_run is None:\n",
    "        return None\n",
    "    \n",
    "    weakly_accepting_run = weakly_accepting_run.reduce()\n",
    "    aps = tuple(_collect_formula_aps(formula))\n",
    "    \n",
    "    rej_run = _rejecting_run_starting_with(aps, aut, weakly_accepting_run.cycle[-1].s)\n",
    "    \n",
    "    def to_ap(x):\n",
    "        return _bdd_to_tuple(aps, aut.get_dict(), x.label)\n",
    "    \n",
    "    h_prefix = [] \n",
    "    h_cycle = [to_ap(x) for x in weakly_accepting_run.cycle]\n",
    "    q_prefix = [to_ap(x) for x in rej_run.prefix]\n",
    "    q_cycle = [to_ap(x) for x in rej_run.cycle]\n",
    "\n",
    "    if not is_accepting_example:\n",
    "        h_prefix, q_prefix = q_prefix, h_prefix\n",
    "        h_cycle, q_cycle = q_cycle, h_cycle\n",
    "    \n",
    "    return ConstructedMDP(\n",
    "        aps=tuple(a.to_str() for a in aps),\n",
    "        g_prefix=[to_ap(x) for x in weakly_accepting_run.prefix], \n",
    "        g_cycle=[to_ap(x) for x in weakly_accepting_run.cycle],\n",
    "        h_prefix=h_prefix, \n",
    "        h_cycle=h_cycle,\n",
    "        q_prefix=q_prefix,\n",
    "        q_cycle=q_cycle, \n",
    "        is_accepting_example=is_accepting_example)\n",
    "\n",
    "\n",
    "def construct_mdp(formula: Union[str, spot.formula]):\n",
    "  \"\"\"Constructs an either accepting or rejecting example MDP.\n",
    "  \n",
    "  Args: \n",
    "    formula: the LTL formula.\n",
    "  \n",
    "  Returns:\n",
    "    a ConstructedMDP, or None if no example exists (the formula is in B).\n",
    "  \"\"\"\n",
    "  mdp = _construct_mdp(formula, True)\n",
    "  if mdp is None:\n",
    "    mdp = _construct_mdp(formula, False)\n",
    "    if mdp is None:\n",
    "      assert spot.mp_class(formula) == 'B'\n",
    "  return mdp    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_graphviz(mdp: ConstructedMDP, m1_or_m2 = \"M1\"):\n",
    "    assert m1_or_m2 in (\"M1\", \"M2\")\n",
    "    \n",
    "    g = graphviz.Digraph('MDP')\n",
    "    g.attr(rankdir='LR')\n",
    "    g.attr('node', shape='circle')\n",
    "\n",
    "    def label_to_str(lab):\n",
    "        return \",\".join([f\"{ap}={int(x)}\" for ap, x in zip(mdp.aps, lab)])\n",
    "    \n",
    "    def add_single_edge(g, s1, s2, action=\"a1\", prob=\"1\"):\n",
    "        if prob == \"1\":\n",
    "            g.edge(s1, s2, label=\"\", color=\"black\")\n",
    "        else:\n",
    "            color = dict(a1=\"blue\", a2=\"orange\")[action]\n",
    "            g.edge(s1, s2, label=f\"{action}, {prob}\", color=color, fontcolor=color)\n",
    "        \n",
    "    def add_parallel_edge(g, s1, s2, prob=\"1\"):\n",
    "        if prob in (1, \"1\"):\n",
    "            add_single_edge(g, s1, s2, prob=prob)\n",
    "        else:\n",
    "            add_single_edge(g, s1, s2, action=\"a1\", prob=prob)\n",
    "            add_single_edge(g, s1, s2, action=\"a2\", prob=prob)\n",
    "\n",
    "    def add_path(prefix, cycle, name: str, init=True, color=\"green\", wrap_cycle=True):\n",
    "        states = list(itertools.chain(prefix, cycle))\n",
    "        with g.subgraph(name=name) as bg:\n",
    "#             bg.attr(rank=\"same\")\n",
    "            bg.attr('node', fillcolor=color, style='filled')\n",
    "            if init:\n",
    "                bg.node(\"init\", style=\"invis\")\n",
    "                bg.edge(\"init\", f\"{name}-0\")\n",
    "            for s, state in enumerate(states):\n",
    "                bg.node(f\"{name}-{s}\", label=label_to_str(state))\n",
    "            for s1 in range(len(states) - 1):\n",
    "                s2 = s1 + 1\n",
    "                add_parallel_edge(bg, f\"{name}-{s1}\", f\"{name}-{s2}\", prob=\"1\")\n",
    "            if wrap_cycle:\n",
    "                add_parallel_edge(bg, \n",
    "                                  f\"{name}-{len(prefix) + len(cycle) - 1}\", \n",
    "                                  f\"{name}-{len(prefix)}\", \n",
    "                                  prob=\"1\")\n",
    "\n",
    "    add_path(mdp.g_prefix + mdp.g_cycle, [], name=\"g_region\", wrap_cycle=False, color=\"grey\")\n",
    "    add_path(mdp.h_prefix, mdp.h_cycle, \n",
    "             name=\"h_region\", init=False, color=\"green\")\n",
    "    add_path(mdp.q_prefix, mdp.q_cycle, \n",
    "             name=\"q_region\", init=False, color=\"red\")\n",
    "    \n",
    "    last_state_in_g = len(mdp.g_prefix + mdp.g_cycle) - 1\n",
    "    add_parallel_edge(g, \n",
    "                      f\"g_region-{last_state_in_g}\", \n",
    "                      f\"g_region-{len(mdp.g_prefix)}\", \n",
    "                      prob=\"1-p\")\n",
    "    \n",
    "    if m1_or_m2 == \"M1\":\n",
    "        h_act, q_act = \"a1\", \"a2\"\n",
    "    else:\n",
    "        h_act, q_act = \"a2\", \"a1\"\n",
    "\n",
    "    add_single_edge(g, f\"g_region-{last_state_in_g}\", f\"h_region-0\", action=h_act, prob=\"p\")\n",
    "    add_single_edge(g, f\"g_region-{last_state_in_g}\", f\"q_region-0\", action=q_act, prob=\"p\")\n",
    "\n",
    "        \n",
    "    return g\n",
    "\n",
    "\n",
    "def visualize_constructed_MDP(formula: str):\n",
    "  \"\"\"Helper to visualize our construction.\"\"\"\n",
    "  formula_tex = spot.formula(formula).to_str('mathjax')\n",
    "  mdp = construct_mdp(formula)\n",
    "  if mdp is None:\n",
    "    msg = f\"\"\"\n",
    "    The formula ${formula_tex}$ is in `Finitary`. \n",
    "    It therefore does not have a counter-example MDP.\n",
    "    \"\"\"\n",
    "    msg = textwrap.dedent(msg)\n",
    "    display(Markdown(msg))\n",
    "    return\n",
    "\n",
    "  msg = textwrap.dedent(f\"\"\"\n",
    "  - green: staying in these states creates an **accepting** string for the formula ${formula_tex}$\n",
    "  - red: staying in these states creates an **rejecting** string for the formula ${formula_tex}$\n",
    "  \"\"\")\n",
    "  display(Markdown(\"$\\mathcal{M}_1$: \"))\n",
    "  display(to_graphviz(mdp, m1_or_m2=\"M1\"))\n",
    "  display(Markdown(\"$\\mathcal{M}_2$: \"))\n",
    "  display(to_graphviz(mdp, m1_or_m2=\"M2\"))\n",
    "  display(Markdown(msg))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may use the code cell below to visualize the construction. \n",
    "\n",
    "We use the same synatx for LTL formulas as in the Spot (https://spot.lrde.epita.fr)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$\\mathcal{M}_1$: "
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.48.0 (20210717.1556)\n",
       " -->\n",
       "<!-- Title: MDP Pages: 1 -->\n",
       "<svg width=\"340pt\" height=\"165pt\"\n",
       " viewBox=\"0.00 0.00 339.68 165.40\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 161.4)\">\n",
       "<title>MDP</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-161.4 335.68,-161.4 335.68,4 -4,4\"/>\n",
       "<!-- init -->\n",
       "<!-- g_region&#45;0 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>g_region&#45;0</title>\n",
       "<ellipse fill=\"grey\" stroke=\"black\" cx=\"107.19\" cy=\"-66.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-63\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- init&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>init&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.84,-66.7C53.93,-66.7 63.29,-66.7 72.18,-66.7\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"72.3,-70.2 82.3,-66.7 72.3,-63.2 72.3,-70.2\"/>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"blue\" d=\"M103.28,-91.27C103.06,-101.12 104.36,-109.4 107.19,-109.4 108.96,-109.4 110.14,-106.16 110.71,-101.34\"/>\n",
       "<polygon fill=\"blue\" stroke=\"blue\" points=\"114.21,-101.4 111.11,-91.27 107.22,-101.12 114.21,-101.4\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-113.2\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"blue\">a1, 1&#45;p</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"orange\" d=\"M100.29,-90.69C97.8,-108.92 100.1,-127.4 107.19,-127.4 112.96,-127.4 115.56,-115.2 114.99,-100.84\"/>\n",
       "<polygon fill=\"orange\" stroke=\"orange\" points=\"118.46,-100.34 114.1,-90.69 111.49,-100.96 118.46,-100.34\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-131.2\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"orange\">a2, 1&#45;p</text>\n",
       "</g>\n",
       "<!-- h_region&#45;0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>h_region&#45;0</title>\n",
       "<ellipse fill=\"green\" stroke=\"black\" cx=\"220.59\" cy=\"-109.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"220.59\" y=\"-106\" font-family=\"Times,serif\" font-size=\"14.00\">a=1</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;h_region&#45;0 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;h_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"blue\" d=\"M130.47,-75.3C146.94,-81.65 169.59,-90.4 187.94,-97.48\"/>\n",
       "<polygon fill=\"blue\" stroke=\"blue\" points=\"186.76,-100.78 197.35,-101.11 189.28,-94.25 186.76,-100.78\"/>\n",
       "<text text-anchor=\"middle\" x=\"163.89\" y=\"-97.5\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"blue\">a1, p</text>\n",
       "</g>\n",
       "<!-- q_region&#45;0 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>q_region&#45;0</title>\n",
       "<ellipse fill=\"red\" stroke=\"black\" cx=\"220.59\" cy=\"-24.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"220.59\" y=\"-21\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;q_region&#45;0 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;q_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"orange\" d=\"M130.47,-58.3C146.94,-52.09 169.59,-43.55 187.94,-36.63\"/>\n",
       "<polygon fill=\"orange\" stroke=\"orange\" points=\"189.23,-39.88 197.35,-33.08 186.76,-33.33 189.23,-39.88\"/>\n",
       "<text text-anchor=\"middle\" x=\"163.89\" y=\"-54.5\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"orange\">a2, p</text>\n",
       "</g>\n",
       "<!-- h_region&#45;1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>h_region&#45;1</title>\n",
       "<ellipse fill=\"green\" stroke=\"black\" cx=\"306.99\" cy=\"-109.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"306.99\" y=\"-106\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- h_region&#45;0&#45;&gt;h_region&#45;1 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>h_region&#45;0&#45;&gt;h_region&#45;1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M245.55,-109.7C253.81,-109.7 263.24,-109.7 272.14,-109.7\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"272.26,-113.2 282.26,-109.7 272.26,-106.2 272.26,-113.2\"/>\n",
       "</g>\n",
       "<!-- h_region&#45;1&#45;&gt;h_region&#45;1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>h_region&#45;1&#45;&gt;h_region&#45;1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M299.2,-133.41C298.57,-143.64 301.16,-152.4 306.99,-152.4 310.72,-152.4 313.12,-148.8 314.2,-143.53\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"317.7,-143.59 314.77,-133.41 310.71,-143.2 317.7,-143.59\"/>\n",
       "</g>\n",
       "<!-- q_region&#45;0&#45;&gt;q_region&#45;0 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>q_region&#45;0&#45;&gt;q_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M212.81,-48.41C212.17,-58.64 214.77,-67.4 220.59,-67.4 224.32,-67.4 226.73,-63.8 227.81,-58.53\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"231.31,-58.59 228.37,-48.41 224.32,-58.2 231.31,-58.59\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7fbad864d3a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$\\mathcal{M}_2$: "
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.48.0 (20210717.1556)\n",
       " -->\n",
       "<!-- Title: MDP Pages: 1 -->\n",
       "<svg width=\"340pt\" height=\"165pt\"\n",
       " viewBox=\"0.00 0.00 339.68 165.40\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 161.4)\">\n",
       "<title>MDP</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-161.4 335.68,-161.4 335.68,4 -4,4\"/>\n",
       "<!-- init -->\n",
       "<!-- g_region&#45;0 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>g_region&#45;0</title>\n",
       "<ellipse fill=\"grey\" stroke=\"black\" cx=\"107.19\" cy=\"-66.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-63\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- init&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>init&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.84,-66.7C53.93,-66.7 63.29,-66.7 72.18,-66.7\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"72.3,-70.2 82.3,-66.7 72.3,-63.2 72.3,-70.2\"/>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"blue\" d=\"M103.28,-91.27C103.06,-101.12 104.36,-109.4 107.19,-109.4 108.96,-109.4 110.14,-106.16 110.71,-101.34\"/>\n",
       "<polygon fill=\"blue\" stroke=\"blue\" points=\"114.21,-101.4 111.11,-91.27 107.22,-101.12 114.21,-101.4\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-113.2\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"blue\">a1, 1&#45;p</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;g_region&#45;0 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;g_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"orange\" d=\"M100.29,-90.69C97.8,-108.92 100.1,-127.4 107.19,-127.4 112.96,-127.4 115.56,-115.2 114.99,-100.84\"/>\n",
       "<polygon fill=\"orange\" stroke=\"orange\" points=\"118.46,-100.34 114.1,-90.69 111.49,-100.96 118.46,-100.34\"/>\n",
       "<text text-anchor=\"middle\" x=\"107.19\" y=\"-131.2\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"orange\">a2, 1&#45;p</text>\n",
       "</g>\n",
       "<!-- h_region&#45;0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>h_region&#45;0</title>\n",
       "<ellipse fill=\"green\" stroke=\"black\" cx=\"220.59\" cy=\"-109.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"220.59\" y=\"-106\" font-family=\"Times,serif\" font-size=\"14.00\">a=1</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;h_region&#45;0 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;h_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"orange\" d=\"M130.47,-75.3C146.94,-81.65 169.59,-90.4 187.94,-97.48\"/>\n",
       "<polygon fill=\"orange\" stroke=\"orange\" points=\"186.76,-100.78 197.35,-101.11 189.28,-94.25 186.76,-100.78\"/>\n",
       "<text text-anchor=\"middle\" x=\"163.89\" y=\"-97.5\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"orange\">a2, p</text>\n",
       "</g>\n",
       "<!-- q_region&#45;0 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>q_region&#45;0</title>\n",
       "<ellipse fill=\"red\" stroke=\"black\" cx=\"220.59\" cy=\"-24.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"220.59\" y=\"-21\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- g_region&#45;0&#45;&gt;q_region&#45;0 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>g_region&#45;0&#45;&gt;q_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"blue\" d=\"M130.47,-58.3C146.94,-52.09 169.59,-43.55 187.94,-36.63\"/>\n",
       "<polygon fill=\"blue\" stroke=\"blue\" points=\"189.23,-39.88 197.35,-33.08 186.76,-33.33 189.23,-39.88\"/>\n",
       "<text text-anchor=\"middle\" x=\"163.89\" y=\"-54.5\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"blue\">a1, p</text>\n",
       "</g>\n",
       "<!-- h_region&#45;1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>h_region&#45;1</title>\n",
       "<ellipse fill=\"green\" stroke=\"black\" cx=\"306.99\" cy=\"-109.7\" rx=\"24.9\" ry=\"24.9\"/>\n",
       "<text text-anchor=\"middle\" x=\"306.99\" y=\"-106\" font-family=\"Times,serif\" font-size=\"14.00\">a=0</text>\n",
       "</g>\n",
       "<!-- h_region&#45;0&#45;&gt;h_region&#45;1 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>h_region&#45;0&#45;&gt;h_region&#45;1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M245.55,-109.7C253.81,-109.7 263.24,-109.7 272.14,-109.7\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"272.26,-113.2 282.26,-109.7 272.26,-106.2 272.26,-113.2\"/>\n",
       "</g>\n",
       "<!-- h_region&#45;1&#45;&gt;h_region&#45;1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>h_region&#45;1&#45;&gt;h_region&#45;1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M299.2,-133.41C298.57,-143.64 301.16,-152.4 306.99,-152.4 310.72,-152.4 313.12,-148.8 314.2,-143.53\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"317.7,-143.59 314.77,-133.41 310.71,-143.2 317.7,-143.59\"/>\n",
       "</g>\n",
       "<!-- q_region&#45;0&#45;&gt;q_region&#45;0 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>q_region&#45;0&#45;&gt;q_region&#45;0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M212.81,-48.41C212.17,-58.64 214.77,-67.4 220.59,-67.4 224.32,-67.4 226.73,-63.8 227.81,-58.53\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"231.31,-58.59 228.37,-48.41 224.32,-58.2 231.31,-58.59\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x7fbaf8bb3820>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "\n",
       "- green: staying in these states creates an **accepting** string for the formula $\\mathsf{F} a$\n",
       "- red: staying in these states creates an **rejecting** string for the formula $\\mathsf{F} a$\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize_constructed_MDP(\"a & X a\")     # formula is in Finitary\n",
    "visualize_constructed_MDP(\"F a\")         # formula is in Guarantee\n",
    "# visualize_constructed_MDP(\"G a\")         # formula is in Safety\n",
    "# visualize_constructed_MDP(\"G a & F b\")   # formula is in Obligation\n",
    "# visualize_constructed_MDP(\"FG a\")        # formula is in Persistence\n",
    "# visualize_constructed_MDP(\"GF a\")        # formula is in Recurrence\n",
    "# visualize_constructed_MDP(\"GF a & FG b\") # formula is in Reactivity\n",
    "\n",
    "# visualize_constructed_MDP(\"!a & X G ((a -> X !a) & (!a -> X a))\") # a larger example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To support the fact that the above constructed MDPs are indeed counterexamples, we convert the constructed MDPs into PRISM format and estimate the LTL-PAC probabilities on these constructed counterexampe MDPs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rewrite_ltl_formula_for_prism(ltl_formula: str) -> str:\n",
    "    \"\"\"Rewrites ltl_formula with proper quoting so it works for prism.\"\"\"\n",
    "    ltl = spot.formula(ltl_formula)\n",
    "    relabel_map = {}\n",
    "    for ap in spot.atomic_prop_collect(ltl):\n",
    "        # we add $ here to make spot satisfied with the double quotes...\n",
    "        # this means that we assumed there are no $ in the AP names\n",
    "        relabel_map[ap] = spot.formula(f'\"${ap}$\"')\n",
    "        relabel_map = spot.relabeling_map(relabel_map)\n",
    "        ltl_transformed = spot.relabel_apply(ltl, relabel_map)\n",
    "        # format formula with parenthesis, then remove $\n",
    "        return f'{ltl_transformed:p}'.replace(\"$\", \"\")\n",
    "\n",
    "\n",
    "def to_prism(formula: str, env_title: Optional[str] = None, m1_or_m2 = \"M1\"):\n",
    "    \"\"\"Writes a PRISM model file for the counterexample for a given LTL formula.\n",
    "    \n",
    "    Args:\n",
    "        formula: an LTL formula.\n",
    "        env_title: the title to prepend to the PRISM file as a comment.\n",
    "        m1_or_m2: whether to generate M_1 or M_2 as in Figure 2 of the paper.\n",
    "    \n",
    "    Returns:\n",
    "        a file `name` for the model, written to \"models/{name}.prism\".\n",
    "    \"\"\"\n",
    "    assert m1_or_m2 in (\"M1\", \"M2\")\n",
    "\n",
    "    if env_title is None:\n",
    "        env_title = f\"Constructed MDP for {formula}\"\n",
    "    \n",
    "    mdp = construct_mdp(formula)\n",
    "\n",
    "    env = f\"constructed-mdp-#{formula}#\"\n",
    "    with pathlib.Path(f\"models/{env}.prism\").open(\"w+\") as f:\n",
    "        print(textwrap.dedent(f\"\"\"\n",
    "        // @title: {env_title}\n",
    "        \n",
    "        mdp\n",
    "        \n",
    "        const double p = {{P}};\n",
    "        \"\"\")[1:], file=f)\n",
    "        states = {s: ap for s, ap in enumerate(itertools.chain(mdp.g_prefix, mdp.g_cycle, \n",
    "                                                               mdp.h_prefix, mdp.h_cycle, \n",
    "                                                               mdp.q_prefix, mdp.q_cycle))\n",
    "                 }\n",
    "        \n",
    "        for i, ap in enumerate(mdp.aps):\n",
    "            ap_states = \"|\".join([f\"s={s}\" for s, sap in states.items() if sap[i]])\n",
    "            print(f'label \"{ap}\" = {ap_states};', file=f)\n",
    "        \n",
    "        print(\"\", file=f)\n",
    "        print(\"module grid\", file=f)\n",
    "        print(f\"  s : [0..{len(states)-1}];\", file=f)\n",
    "        print(\"\", file=f)\n",
    "\n",
    "        def print_det_chain(start, end):\n",
    "            for s1 in range(start, end-1):\n",
    "                s2 = s1 + 1\n",
    "                print(f\"  [] s={s1} -> (s'={s2}) ;\", file=f)\n",
    "        \n",
    "        g_start, g_cycle, g_end = 0, len(mdp.g_prefix), len(mdp.g_prefix + mdp.g_cycle)\n",
    "        h_start, h_cycle, h_end = g_end, g_end + len(mdp.h_prefix), g_end + len(mdp.h_prefix + mdp.h_cycle)\n",
    "        q_start, q_cycle, q_end = h_end, h_end + len(mdp.q_prefix), h_end + len(mdp.q_prefix + mdp.q_cycle)\n",
    "\n",
    "        print_det_chain(g_start, g_end)\n",
    "        print(f\"  [a1] s={g_end-1} -> 1-p: (s'={g_cycle}) + p: (s'={h_start});\", file=f)\n",
    "        print(f\"  [a2] s={g_end-1} -> 1-p: (s'={g_cycle}) + p: (s'={q_start});\", file=f)\n",
    "\n",
    "        print_det_chain(h_start, h_end)\n",
    "        print(f\"  [] s={h_end-1} -> (s'={h_cycle});\", file=f)\n",
    "\n",
    "        print_det_chain(q_start, q_end)\n",
    "        print(f\"  [] s={q_end-1} -> (s'={q_cycle});\", file=f)\n",
    "        \n",
    "        print(\"endmodule\", file=f)\n",
    "    \n",
    "    with pathlib.Path(f\"models/{env}.ltl\").open(\"w+\") as f:\n",
    "        prism_formula = rewrite_ltl_formula_for_prism(formula)\n",
    "        print(prism_formula, file=f)\n",
    "    \n",
    "    return env"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78a5bf3c16d648369548491684764293",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "<multi-discount | Q > Iteration over p:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iteration over N:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAAE1CAYAAACWU/udAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACAjklEQVR4nOy9d3hc6V3o/3mnj6Zq1LssW5Z7W3vtXW/L9vQspLIBshA2pCxJgPtcCD9IAZLL5d4EEiA3ARJICCUNSA9s9xbX9dqSLVmyei8z0vQ+7++PM7ZlW5Zkr6U5Y7+f59GjM+ecOfOdo6PzPd8upJQoFAqFQqEAQ6EFUCgUCoVCLyilqFAoFApFHqUUFQqFQqHIo5SiQqFQKBR5lFJUKBQKhSKPUooKhUKhUORRSlGhUCgUijxKKRYxQohGIcQPhBBnhRCnhRBfE0KUXKdjv00Ise46HeseIcTuq3zPbwohPrXA+n8QQviFEKZ5614WQhyc91mzQojjQohuIcSPhRA75+37rBCiUwhxQgjxlBCiZoHP+N38+fzCVX3Rq0AIMSCEsM17faXv+6wQYsMSx3qfEOJ9KyDjRP73DiHE/fPWf0YIcfsS77UJIQbyy28RQvzW9ZZvOSx27S123oQQ785fA11CiIOLHKP53LW30gghaoUQ31iNz7qZUUqxSBFCCOA/gX+TUq6TUm7Kv3Zdp494G7CgUhRCGK/yWPcAV6UUl2AcuDcvSwNw6YPAz6WUO6WU64G/Bn4mhKiet/0RKeV24CXgEwsc//3AbVLKjy9HmGs4H8XGDuC8UpRS/pGU8qXlvllK+QMp5RdXQrBlcA9Xee3lFeCngfullBuA/wF8Xwhxvf63Fvts05W2SSnHpJS/stIy3OwopVi8PADMSim/dW6FlPI/pZSTQogKIcRPhBAnhRA/F0JUwXmr48+EEMeEEK/kFQpCiI/ln4hPCCH+In9TeAvw10KIV4UQrrxl80f5p+I98y2d/BP3/8ovbxRCPJc/1kEhRBPwm8Af5I/Vlv95Mi/HT4QQZfn3vjVv3b0M7Fnku38HeHt++e3A9660o5TyJ8CPgfcssPlFYO38FUKILwItwPNCiF8RQqwTQhzIn8tvCyEc+f0uOh+XHOPTQogjQogOIcT/WeR7LJcP5M/ni0KIuuW+SQjxKSHE3wkhXhBC9Aoh7hJC/LMQ4owQ4tP5fS6ydIRmiT98yaE+AzyW//u97gr7IISoEkI8I4RoB/6/eevnXx//e9619j/y6y69ZsxLXMMb8sv3CCH+dZ7cfymEOCQ0T8AOIUQtl1x7yzx1Hwc+J6UcA5BSHgCeBh5d5vsRQrw5L8urQogvzVv/o/x13yGEeFd+XXP+//HrwLH89/q5EOKHQogeIcQfzNvvnEfkU0KIv533t30wv94hhPgPIcQpIcSXRd7aVywfpRSLl43Aq1fY9mk0a2kbmvX4J/O2RaSUtwDfQrthAPwBsDNvPX1SSnkU+AHwYSnlDillOL/fpJRyn5RyMXfRN4BP5Y/1MDAE/D/gT/PHOgP8FfBYXo5/BX5fCGEH/gLtyf4uYNMin3EIuEVoFtpb87Iuxglg/QLrXw+cnr9CSvlbwBiwV0r5DeCLwP/Jn8sB4GPzdr/S+fhLKeUeYBuwTgixdwn5liKZP59/h6agroY6tPP528B/oCmrbWhKzrnMY/wR8PX83++ZRfb7JPA9KeVW4LKbsRDCB7wD2JT/Pl/Nb7r0msmw+DV8JZxSyr1o3/F/5JXapdfectgEvHLJuuPAom7scwghKoDfAu6SUu4AHEKIN+U3/0r+ur8d+CMhhDm/fhvadbY9/3oH8D5gO/Cb5x7GLqEJuBt4J9q5B/gwcEZKuRntYbBqOTIrLqCU4o3J7cC/5Jf/Cbhj3rZzCuQ42j8VaDeAb+SfXJOLHPe7i32oEMINuM7dOKWUc/KS5rpCc0HdDvxQCPEqmmuqAWgDTuddRGng+4t+Q3gOeC8QAeaW2Fdc8vrfhRAngErgs0u8d4eU8j/zy5eeyyudj/uEEIfRHlr2sbCCX6jp8JUaEZ+zhL8L7F9U2sv5qZQyB7QDA1LKPillEk3BVy/6zqvnduDf8sv/tsD2ENrf66tCiDcDkUWumcWu4Sux0LVdCPahKblD+Wt8Pxc8Eh/LX3sH0GSsz6/vlFKemneMl6SUfillDOgHLot9Az+WUma5+PveTv66lFL+CIhdt291k3BF/7VC93SiuTgXQnDlG+w5pZfjwt//DcDr0FyRH0Sz1hYiOm85y4WHKuu89Ut1mBfAaP4J+sLKeckw8/ZbjO8CPwJ+Z4n9QLtBtc97/YiUsmsZ74OLv8+lMkUveY3QXMqfB/ZIKSeEEP+Xi8/POQJAKVp8FMAHzCxDhqvt4J/K/87NWz732oR2Pcx/OF5I1gURQrydC27SN3Dx+bns7yelzAjNNf8Q8CvAL6JZVAt9pytdw1e67mDha/ta6AR2ASfnrdsJPCWEaEazuEGzoE9yOQL4Tynl4xetFOJ1aApzr5QyIYQ4lv8OCS6/luY/nJ77PpmF9pFS5sSFWORi//uKZaAsxeLlSaBMCPHucyuEEO8UWkLJi8C59b8EvHClgwghDEC9lPJJNBfbxvymCIsn7QwCO/LvfwOAlDKE9vR/d/7YXiGEmH+s/D6z82Ig1nyspwvYJLQMOzPwyBLf/2U0K+/fF9tJCPEG4E1csDqulhN5qwa0uOQVz2UeG9pNzC+E8KK5dxfiAJqlS951/Mgix54fP112gssymQTq87EoDwtbZAteC1LK7+bdkjvyrsoXgXflN7/z0v3z7lqPlPIHwP9Es8KvdM1c6RoeRHMtgvZ3XYqlruOF+ALwiXxMEiHEHWhuzO9KKQfmfecrue0PonkL6vLvrxBalrMbCOQV4g60h7XrzUtoDxvnrv3rko1+M6EsxSIl/3T4NuCvhBB/gvbU+DKa9fQp4B+FEL+BFtv55UUOZQT+Oe/WFGjxRYBvA1/JB/nvXOB9fwL8IzACdM9b/yto7jE3muvmrrxM3xFCPIp203wv8P+EEH+e//w/llKeEUJ8HHgGzXq6NKZz6feXwP8B7UZ6yeaH8m4rB9ADPCSlvNaEg98Cvi6E+FM0xf3YEnLNCSH+CTiF5qK8Uvz1M2jn91W08/5VKeWJKx1WCHEIzUp6x1V/g8XlTQkh/gLNBdfJwpbPs2hx31fRklCuxKeBfxNamcN/LLDdBfxACGFBs2b+ML9+oWvmUyx8DX8h/xkfZIlrJM9F195y4opSyiNCS0R6Ki+rB9gspYxf4S07hRAj814/iBbb+0H+AS+JFh/8GfDB/Hk8DRxbhvxXy1+j/T8fR0sOGlqBz7ihEWqeokJR3OSVEFLKfyisJMXFcs5bXin+K9ArpfwfqyPZtZN3oxqllMm8hfsnUsp7CixWUaEsRYVCobgCUsoU8AuFluMq8AD/lVeOSeA3CixP0aGUokJR/LxaaAGKlFcLLcD1RkrpB24ptBzFjHKfKhQKhUKRR2WfKhQKhUKRp2jcp+Xl5bK5ubnQYigUCoWiyDl27NiMlLJioW1FoxSbm5s5evRoocVQKBQKRZEjhBi80jblPlUoFAqFIo9SigqFQqFQ5FFKUaFQKBSKPEopKhQKhUKRRylFhUKhUCjyKKWoUCgUCkUepRQVCoVCocizIkoxPxPvFSFEYt7wy/nbnhZCvCSEuH8lPl+hUCgUimthpYr3A8B9LDwA9vfQpnWfRJt19uQKyaBQKApMOpcmlo5hNpgxGoyYhAmjwVhosRSKK7IiSlFKmQAS2gDty9gGfFRKKYUQYSGES0oZXgk5FIpLiaeyTIYS1/z+WCpD/0yUIX+MVDZ3XWQS2TSGdGTxnRZs3K+tM+TS5ISR5uAhHOkAAKfLH2bTzM85Xf7QRb/PvWeltx8rvYe1Uz/kSccOWqd/xJmKN7A/dBAhJQJ4yXMH+0Mvc8h7B7eHXuaQ505uC76IIX+EQ959bAg9z3H3draFTzBR+lZsc9+mLufEhokfujzUx84wUtJ20e+N2XJsmPipu5K3hed4yr0OR+Tp89uFzEE8ydSUj7eWmHiypJkNmdPMGO4nZfwBWaH9TUcsa6nP9DFsXsvt8WmmrQ+Tzn6PnelKoqkoJ9NWWvwDxDJlhFsSTBibubXzJE0ZN3FDhkNNG9njSNLXF8GdGmDWVInNEmB0nYuqoSjGcA2bk0HarS5K15oZzW6nsveCDfHq2k247GP4uixsSEaZtrQSd3RRX7uWyNAQYyk7pZkp/OYqZjemSAQr2T7QQ2W6gRmDk47WOnZFY4xMzeDNdTFnuAu8R5mtKKem1wo5QX3ayumSLC2NMDixjbLwv5LDRw4PhzZ5aQlksE5H2Zac5az5NowVL+OzbyA6OMmMKYI742XcKUg2zWIb3Mr68AHc2UamjWlObb2NbSNjTIQHKM1OETTsJFF7kqSljvq+MBjK8GZsnPLF2Oqz0z9egS/+M7LUI0nz3C3N7DybxRieZFMqR495Hdb6bnY++kV2Npdf9f/Y1VCINm9GeWE0RxAoBRZUikKIx4HHARobG1dHOsUNSzKT5ZkzU3z32Aj3b6xadN9MNsd0JMlkKMFk6NzvBJFkhiq3jTVlDsaCcdZWOF+jVDms4RHCU4O4y2oW3bM066cyN0lldpKXrXfwK5G/55+c7+NTc3/Ap0v/lNujPyGUlBhLvHSm97E5eIbO9N6LfhvJAiy5/ZWhcjaXXO37u5AyRVpmeSncwsZYL6OJzbQmxhmYG2dtqpeEIUvckKM9CveGxujODWGQJxmY2o3fcpTTtjhJkYP4QQ4Zs0SjfRiTCaYm7qYhdYbanBdj1ki/KMOcHqE/Z73o94ZYUNtOAnNgkrNZB9W5C9tFNoOQMIwRSzxOP17WMEZvNI6tYoKMzOGISUI5I5OuCeraozSPBagLTmE2jpK9M4StO8TdIxlS7gzjpjADdTCdNWGdDGBPRUDkCDpKMWfmMI2lsGYD2EUa6Y4ymQlTNpfCPpHDlgphtbjIVQqG+o00BAPn/9ZzZSPETJPUjZmwp5IIYSNXM4XVIYjPTGOPWLHKMDZjhonWDK7JeszTORypUfyOTUzNTkPCgNUfxpYJYHZYiRlnmC21UButwRIYAmsTpkQIc+ks2fBObNMB0iYnWZORifgodbEqvNNz2GOTCK+BjHUcc6YCQxhKEmMY7NUYU+PM1M/QGDFgn5zGaKjAZDEwFAixNSyxT/mxyjBmmxF/6RQJWY45KJHRUXLOjWTEIGaTBULl2OYCJC0NILIMJUfYFm/BMTmFJWPE6MghvaP82//+Ojv/ZmVnPa/o6CghxLPA/VLKzPx15yZBCyF+ALxXShla6li7d++Wqvep4lrJZHM83TXFTCTFJ/69nc8+shWAnJQEoikmQwkm5inA2WiKUoeFareNKreVWq+dHQ1ebl3jo9Zjx2AQfOG/u/n4A+uvXSgp4exTEByGH30M3vgFSEUgOg3uOhg4AI374JVvwPb3aL8NJnBUwIY3afutuRO8TeBtBJOFl//+d7nt1//Paz5fV3OcUCpE31wfg6FBUrkUAKlsCoFgJDLCN09/kzpnHXaT/aKffTX7SOfSzCXnqHfVE01F8dq81DhqKLWV4rV68Vq9WIwWAH7/62/jc4/9xzV9n2wkSvTA82RDYXKJBFOf+xyeX/gFsnNzWFtbMZaX4//yl8mGQhisVuw7d+J64AFiR49CLofB68XS2EjJzh2YqqoxVZQjDAb+9dfv4N1//8I1yTSf5R4nl5PEgikiswnCgQT+kQiuMhuj3XM0bfaxfm81X3/8j/i1v/3j1yzT137jD6/Lca7nsa7XcYQQx6SUuxfcVgCl+EXgX9Biij8+pyCXQilFxbWSy0me655mPJhgOpzkC092s6vRy0RIe+20mqhy2+b9WKlwWrGZjdR67TSVlVDrtWM0LBgOuHZGjsHESej5bzjzY/C1QHhCU3o73wuBfqjcBJk42H1gMEJJGXjqwdMAjnJYOESx4uRkjtHIKH1zfUzFpy5abxAGftj7Q9aXrqfN18ZnXv4Mf3TbH130fqvRisfiwWvz4rF4KLWV4rK4MIjrn/uXCQSIHjhALpEEIDUwQODrX8e2dStGrxfbli2YqqrIzsxg9HgQFk0JG0rsmKqqMFdXY6qqwmC1XnfZlkMqniE8myAymyQ6lySXzZFJ5xDAqQNj3PL6Zho2+ihxWwoiXzGymFJcEfepEMIM/BTYDvxcCPEJNIvwCeB/A98A7MAnV+LzFYpzSCl5uc/PeDDBVDjB3x3oB6CpzMGta8qoclmxmi8kfhgEVHtsNJc5qCu1YzauUNWSv1dTiKmoZhUC7P1NMM67sTkqwGgG1yZNCXrqwVKyMvIsk1g6Rl+wj/5gP4nshdhsIpPg8MRhjk8d5ze3/yZvbHnjRQqu3ll/3vLzWD2UmFfne6THxoi+/DIykyXR3U16eBjXffcB4H372y/a11xTjamy8rwiNLrdqyLjpeSyOaLBFJHZJJHZBKm4ZlNIKRFCMH42SCKSpmVnBfc82kZ1iwfDSl2nNyErlWiTBi4ttziU3zYC3LsSn6tQXMrRwVkG/TFC8TT/8NIAD2+u5ruvjLCn2Xd+HyGg0mWlqcxBg8+O1bTC2ZGRaRh8EbIpePmvYP9HtfXnFKLVdUEJuqo1C7GASCmZiE7QG+xlIjqB5IJ3KZ6Jc2j8EEcmjtBa2sqjGx/FatQsKp/Nx46KHUzHprmt9rZVlzvZ20vs2DGQEG9vJ/zTn+L9pV+6aB9jqTdvCV5wiRaCVCJD2K+5RGOhFDJ3uQfv5NMjrL2lkuq1HkwWA3WtpbjL7QWQ9samaOYpKhRXy8mROXomIyTSWf7hpQFubfaxq6mU774yAkC500JTmYNGXwl2yyopnlQUep+GRBDGT8Cdv3NBGdbv0RSh3bs6sixBIpOgP9hPX7CPWCZ22fbnhp/j8MRh2krb+PUtv47Prj1o2Iw2tpZvpdnTDMCHdnxoNcUGIH7yJInOLgCS/f2Ef/5zSn/1VzFXVWGq0GbLet721oK5ROcjc5Kn/rGTxk2+y7bFwymGTgVo21fNhttqsNiNuHx2atd7Ma/WNXuToZSi4obkzESYjtEQmWyObx4cpLncwd3rK7CaNEvgLTtqcVpX+fLPZjSFmIpoSTOeBi1xpmk/3P17UL1ldeVZhE+99Cm2lm8lx8VlJ+cU5cayjXhtXt6/9f2U2koBMAojrd5WNpRtwGwwF0JsZC5H7PARUoPaDNlcNIqlsZGy3/gNjB4P5uoqHHfcQfmHP6wLhQjgH4vQd3z6IqU4PRzGbDHirrBTuUZz49qcZqrWuCmrfa0Zz4rFUI5oxQ1H/0yUY4Oz5KTkO8dGcFiMvGlbDWajgbvbKvjofa2rrxABBl+A6AwMvqS9bnsDlK2D8nXwut9ffXmuQCwd43s937tIIUZSEWYTs0gpGQ4PI6Vke8X28wqxzlHHg00PsrVia+EUYipF5LnnSA0OIqUk/NRTzH7722AwaAqxphrHHXcgjEYqnvhIQWS8lEwqy/SQVqMqc5LJgRBSSkxmA0aTAYNBUFpVgs1pZu3OCqUQVwFlKSpuKEbn4hzq8yOl5Mft44QTGR7b34zRILhtbRnlTutrK6O4VsZPaNmkgT6o2Q61u8BRBo2rH2tbiu7Z7vPL4VSYl8Ze4sTUCe5vup9dVbt4sPnB89s9Fg/bK7ZT5Vi87nOlycViRJ4/QDYYROZyhH/6U1LDw/h++ZcRQmgKcf9+hFFfLsepwTDZjFb7iYDgdBxfjYPSasf5fcrrXVQ2uRDXO/tZsSBKKSpuGKbCCV7smSEn4UDPDH3TER6/cy1mo4GdjV4afAXK3JwdhNFXtOWZbs1l6muBlnvAqK9/wXPuUYCf9f+Mk9Mn2V65nQ/u+CAui+v8fhaDhc1lm1nrXcsVOletGpnZWaIvvEAuFgcgOzdHdm4O3/veh8Fm061CTETTzE5EiYW02k4hBOv3XHi4MFuN1LWV4vDow817s6Cv/0iF4hqZjaZ47sw0mZzk+NAsB/v8fODutdgtRtqqnWyoLkx6PbGAVoSfy8KRv4Mdj4LVqRXl20sLI9MidM92E0honVWMwsiHdnwIp+WCy86AgRZvC5vLNp8vqi8k6clJoi++hEynkek08ePHse/ZQ+mjjwLoViECTPQGARZUep7KEmrWejCaVIRrtVFKUVH0hBNpnjkzRTor6Z4M85OOCd5/xxo8djP1pXZ2NRZI+aQTWmJNJqV1n2l7A1gc+Thia2FkWoRUNkXvXC92k5bm/0DzAxdtryqpYnvFdjxWTyHEu4xkf3++44wkl0wy9y//gqGkBPstt4DRqGuFGPLHiQaTjPXMwTxD22gyULPOg6eisPWoNzNKKSqKmngqy9NdUyTSOUZn43zn6DDv3ddEldtGmdPC7WvLCuPey+Wg7xlIhrXfEyfh9o9q5RY6jCMC9Mz2EEwFmYxNXrTeaXayrWIbdc66Akl2OfFTp0h0nAIgl0wy+41vYKqsxP3mNyMMBsy1NThuv12XClHmJJP9WmfL8non2XxjeYfHSl1bKWar/mS+mVBKUVG0nGvwHU1m8UeSfOPgAI/srKOpzIHTZuLu9RWYCtXpY+hlrWXbTI+mFPd/XOtOs/Z1uosjAqSzac7OnSWUDDEeGQfAJExsLNvI+tL1K9J+7VqQuRzxY8dI9mlxTyklwmymZO9ebFu3akk1OlaIAP6xKKl4hsn+EN7qEuxOCxtvr6F528pOf1AsD31c6QrFVZLJ5njuzDRzsTSRZIZ/eGmAezdUsqnWg8Vk4O71FdjMBbopTp7WEmrSMSjxwS2Pab91GkcE6A32Ek1HqXJUsb9uPwAPr3mYDb4N+lGI6TTRF144rxCzc3ME/u7vkPE49m3bikIhZtJZpoe0oUCJaBohwGIz8br3biiwZIpz6ONqVyiuglxOcqBnhplIimQmyz++NMC2ei9715RhNMBd68vx2AtTK0doDEYOa8tdP4a5YS3TtGytLuOIAJlchu7Zbl4ae4mnBp/CKIy8f+v7z8cW9cLYJ/6A9PgEAJnpafxf+xq2rVsxOLTyBb0rRNBKMHLZHPFIiqYtZVhsJqrWuFW5hY5QSlFRVMxv8J3NSf750BA1Hhv3b6wEYF9LGZUuW2GESwSh71ltJNTgi7D+Ya0m0eaBxtsLI9My6Av2EUqFODJxhF1Vu2h2N/PRXR8ttFgXkUskCP3wh4DmQp379rdx3Xsvjn37gOJQiIlomtnxKMl4hs4Xx5FSUuKxqv6lOkN/wQ2FYhHONfiWUvL9V0YwGgRv3VGHEILtDR6ayhxLH2QlyKS02YiZpFZ+kYqBwazFEXVYj3iObC7LmcAZjk0cY41nDRX2Cjb49OfKSw8PA5qFaCwrw/frv47Bpj38mGtrtCzTAjXzXi4TfUGklFjtJnY+2IgQguqWApUKKa6Ivq8ihWIe5xp8A/zX6UlmIknevacRo0GwrtLJ5toClQpICf3PaZZicBhe/AtYdx+YbdCwV4sn6pSB0ADhVJiD4we5o+4OGt2NqzbW6WpI5ZWi/2tfIzMzU3QKMRxIEJ1LMtEXYuDkDEIIvJVako1CX+jz8VWhuIRzDb4BXuqd4dRYkA/ctRaLyUCN18ae5gImsIwcheCINgHj6Ndh41tAGLQ4YkUBWsotk5zM0RXo4vjUcWqdtVQ7qmnztRVarMvIRaMkus4AUPqe92Cu1FzlxaIQZU4y0acV6lc2u8gkswiDON/oW6Ev9H01KRRcaPAN0D4a5PnuaR67fQ0Oqwmfw8wd68oL12pspgcmO7TlyJQ2/ql2Rz6OqM96xHMMhgaJZWIMhAa4s+5O6px1uC36u1GnhodJj44CYGlsBIpHIQIExrUSjLGeOZLRDNYSMxUNLjX6Safo/4pS3NSMB7UG3wB9MxF+8Ooov3JbM6UOCw6rkbvXV2IuVC1iZOrCxIuzT2ru07bXa71NW+7R4ok6RUrJmYBmfb1j/Tuoc9Wx0bexwFItTGpoCPu2bedfm2tri0YhZtO58yUYBpPAaDZgshgpqytQ7FuxJEXjPk2lUjz//PP4fD6qq6s5ffo0FouF7du3c+TIEQBuu+02Dh8+TDabZdu2bQwMDBAKhVi7di2xWIzx8XGqq6txOp2cPXsWp9PJ2rVrOXHiBAaDgX379vHSS9pNbvfu3XR0dJBIJNi4cSNTU1P4/X4a80+qQ0NDlJWVUVlZSWdnJzabjS1btnD06FEAbr/9dg4ePEgul2P79u309vYSiURYt24dkUiEiYkJampqKCkpobe3F7fbTXNzMydPnsRoNHLrrbfy8ssvA7Bnzx5OnDhBKpVi06ZNTExMEAgEaGpqIpvNMjIyQnl5OeXl5XR1dWG329m0aRPHjh0DYP/+/bz44osA7Nixg56eHqLRKK2trQSDQaampqitrcVms9HX14fH46GhoYGOjg5MJhO7d+/m4MGDANx6660cP36cdDrN5s2bGR0dZW5ujubmZtLpNKOjo1RUVFBaWkp3dzclJSVs2LCBV1555TJZdu3aRVdXF7FYjPXr1zM7O8v09DR1dXWYzWYGBgbonsvy9FAGXy7Ik5M23nlLPcH+k0QNgg++82GOHz1EJpNhy5YtDA8PEwwGaWlpIZFIMDY2RmVlJR6Ph56eHhwOB62trbz66quXyXLLLbdw+vRp4vE4GzZsYGZmhpmZGerr6zEajQwODl587ZFhu3WEIyc6sWfm2LFhG8fPjpOYOMO2e9/OQPcQoVCHbq896ZIMhAf49sC3ucd1D3va9nDysP6uvZ3r13Poxz/GMTJCSS7HK4EApqpKtszOrvi15/V6qaur49SpU5jNZnbu3Mnhw1q5zb59+zh69OiS115iVnDiaAemnIO165toP3OU8gYX67iP559//tquPXXfe83X3mIIKeWiO+iF3bt3y3MnXnFzkMnm+N4rI/zP77XjsZt5eEs12+u9GATc01ZJtadApRegxRGf+99a2cXRv4cNbwJXdX76xd2Fk2uZ/Hzg54RSIaZj01SUVPC6htdRbtdfR5V4xymmv/AFMBiIHjjAhtOnisJCBEjG0vS+Mk06leHUgTG23l2Hw2OjZWdFoUW76RFCHJNS7l5oW9FYioqbj8lwknAiA8D+deVsr/cCcOsaX2EVImhZpj0/17rW7PpVzVVqc0OTfusRzzESHiGYDHJo/BB7qvdQaa/UpUIESA30E+/ooPTRR7Fu2FA0ChFgoj9ELpfDaDKw/d4GAKpUCYbuKZ4rTHHTMTob458ODgJwxzrtpr21zkNLRYGnjyfDENMSf8gktRiiwQgtr9N1HPEcXYEuzgTO0D7TjkEYdFmXCNqcxNTwCEaPB/umjVT9zm8XWqRlE5lNEAkk8I9E6D6kNVh3l9vVbMQiQFmKCt1yZCBAKG8pAqwpd7C1Xgdji4IjMKv132Tr20EI3dcjnmM8Mk4gEeDA6AHurL+TMnsZVY6qpd9YANJDQ5h8Pny/9muYGxoLLc6ykVIy0aeVD5XVO/FUliAMgipVglEUKEtRoUuCsTQnhoO0VWvT3qs9Vvau0YnS8feB1Mb9YLRoccQK/dX3LURnoJPeuV4yuQxtpW26zTgFSPYPMPe970Eui6WxodDiLJvZ8RjJWJqRM7OEZhKYrUbKap1YbMoGKQaUUlTokrFgnDOTYdqqNKV4x7oKDHpompzNgL8HRrTMv2KJIwJMxabwJ/znrUSv1Uuts7bQYi1IZnqaxOnTZP1+zNU1GOzF0R80m8kxNaRZiQ6vFWuJCZPZQHljgV3+imWjlKJCl/RMhpkIJthUqylFi0knl2poFCKT0PYG7XWRxBEBOv2dDIYGiaajbCrbpNtYImi1iYn2dmxbtmBpKh7X6fRQmGw6h38siqfCjs1hprLZjbFQtbSKq0b9pRS6I5XJ8XKfnzXlDnY1+vjofToauRQcgdlBreE3FEUcEWAmPsNUfAqX2cWbWt6E2+KmwaVPl6SUUutiMzGBbdtWzPX1hRZpWSTjGQJjUXJZyWR/EJmT2BxmvFX66yWruDLKya3QHZOhBF3jYdZXuajz2vn4AzrqHzp+ApyV0HBroSW5Kjr9nQQSAUwGE82eZtp8bYVrjbcEmclJZDJF2Qc+gKW2BoO1ODI2J/uC5HI5pJRs2q+5patbPLo9z4qFUZaiQncMz8bongxza3Mpdj31h4z6YfggTHeBbw3c/XuFlmhZzCZmmYhNMBQa4uzcWUpMJTS7mwst1hVJDQ4ReeYZMlNTmBv0ac1eSmQ2STiQIDQd5/QLYwC4fDYc3uJQ6IoLKEtRoTuO9AcosZrYki/W1w3BYZg8rXWs8dTD636/0BIti05/JzmZY0flDgDWl67HIPT5PCyzWZL9/UQPHsRx2z4sReA61UowtCkYnsoSNvlsCCGoatFB+ZDiqtHnf4bipiUQTXFyJEhblYtab4G71lxKcBiqNkPTfjAVhwUQTAYZjY7yHz3/QVegC5vRxhrPmkKLdUXS4xMkOjowNzRgWbcOYdZ/EtPshFaCMdwZYGYkgtFkwFfrwGpXNkcxopSiQleMzWmlGJvr3FQ4daR40nFtIobVBWXrCi3NsukMdDITn6Ev2EeLp4XW0lZMBv3erNNDg6THxrBv2XJ+TJSeyWZyTA/mC/XrnDhLrRhNBioaXAWWTHGt6Pe/Q3FT0jUeYjqc5PaWMn0lKARHwOoGgxm8xRHnCqfCjIRHeHH0RW6tuRWH2cFaz9pCi3VFZDpNenwc90MPgdGIuaam0CItyfRwmEw6x+RAiLI6ByazkYpGF0azsjeKFfWXU+iGZCbLoYEAayucNPh0lsYe6Ifun0HlRrCXFlqaZdEV6GI2MUv3bDd7qvfQ6m3FrOOayvTYGPH2DuLt7Vjq6xEmfT+z57I5jv10ACklYX8CAIvdhK9GzUosZpRSVOiG8bkEZybCtFW7qPXqqINJLgfhcai/Vcs6LQKi6ShDoSFeGnuJXVW7cJldrCvVt9s3NTRE7PBhhMFQFG3dosEUAyf9pBJZ1t1Sicls1Eow9NB5SXHNKKWo0A3DgRg9kxFuXePDZtZRKUZ4HM4+CdmklnVaBHQFukjlUvTO9bKvZh8t3hasRh3FaC8hl0qROHuW9NgYti1bMFVXF1qkJQkHNOuw88VxQGvr5vLpLDlMcdUopajQBVJKDvUH8JaY2Vyrs2kCwREYPQp2H7j0H+eKpWMMhgYxGUx8eOeHcVvcrC/VUQOEBUiPjJAeGsa2cSPWNc1FMTcxkleK2+/THpSqVQnGDYH+rzzFTYE/mqJjNMj6Kp25TkHrYhOfg6bbtLmJOqd7tptwKsx3u7+LQNDsbsZu0tk5vYTU4CC2DRtwv/WtmBubCi3OkiSiaYZOBwAQQlBa7cDm0G+8VrF8VkwpCiG+IIQ4IIT4y0vW3y2EOCSEOCiE+M2V+nxFcXGuFGNLnZsyh6XQ4lwgPqcpwjt/B0r1H09MZBL0B/uxGC3sqd6DURhp8+l7rFUuHifZ00P4v/4LY4kdU2VFoUVaknAgQVmdllBjMBqobFIlGDcKK6IUhRC7AIeU8k7AIoTYM2/z7wDvAG4HHluJz1cUH6fHgszF0uxbo8NSjFe+AdGZoogn9sz2EElH6J3rpcndRKO7EYdZ39mQqeFh4u0d5BIJzPUN+vr7X4HgVIxkTBuAXdHoxKSndoSK18RKWYq3AU/ml58E9s3bdgrwAFYguthBhBCPCyGOCiGOTk9Pr4igisKTSGc53D/LukodlmLMDsC2d0H9brDoTLZLkFLyd+1/x5GJI3T6OxEI3VuJAOlzY6K2bi2KMVHZdI7AeJSRM7MA+GrVrMQbiZVSil4glF8OAvMLu/4D+BHQBXxrsYNIKb8qpdwtpdxdUaF/l4ri2jjnOt1Q7aLao6PsvUwKhg/D8W+Cr6XQ0ixJIBHg6eGnOTxxmP11+6lz1uG26Cxp6RKykQjJwSFkJoNt0yZMZWWFFmlJwrMJHB4rG2/Tkq50Mfxacd1YKaU4B5z7b3TnX5/jz4E7gFbgV4QQ+n78Vqw4Q4EYvdM6LMUIjYCjHHb/Gnj0Xzc3EZ0AoNHVSEVJha6HCJ8jPTSE0emk/CMfwdqs/wQbgLA/zukXx5BSsvuNzYUWR3GdWSml+DJwX375fuDgvG1ZYE5KmQJygErZuonJ5SQHe/2UO61srNGZVTM3DAf+L6Si4NS/p2IoNATAnfV3Ul1STalN/513koNDBH/0I2QyWRS9TqWURGaTVDW78VSWsPfN+vcgKK6OFVGKUspXgIQQ4gCQk1IeFkJ8Kb/5z4AnhRAvA89IKYMrIYOiOJiJJDk1HmJ9lYsaPblOpdRcpwYj1O4qtDRLEkvHODihPXtWO6rZWLaxwBItTTYYJHnmDKmzZzFVVmD0egst0pLEQilioRSl1Q7cZTq6XhXXjRVrLiil/Oglr5/I//4Z8LOV+lxFcTE6F+fMRJj33NpImZ6mYkSmIDSqKcRS/VswE7GJ81mmFfYKyu3lBZZoaVLnEmyKZCIGQCSQZOiUn/IGF5v21xZaHMUKoIr3FQWlYzRIJJnh1jW+QotyMcERsDih7Q3griu0NEsyGh6lya3F5IohlgiaUkwND2tZp0WiFMOBBOtvraZ+g09NwrhBUX9VRcGIJjMcHZxlfZWLhlKd5VvNDkDvU+CsBB1PlgDI5rKcnD7JN09/E9Dcp3onEwiQi0TxPfYYtg1tGF36L35PJTKM984xORDC5dORV0NxXVFKUVEwxoPx81MxdFWKkYzA1GnY9NaimIoxHZ+m0lHJeza8p9CiLJvU4CDRgwfJjI9jaSqOrNNIIInFbsLmNOMs1dH1qriuKKWoKBj9MzH6Z6LsW1OGxaSjSzE4ApFJ7acISjHGo+McHD+IURj5wLYPFFqcJZFSkhoYIPr88wibDUu9/jsFAYT8cWwOE+X1TtXn9AZGR3cixc1ENic51Oenym1jfbXOOoLM9kM6Dk37waazMpEFOD1zmhdGXqDeWc9Hdn6k0OIsSWZqmkRnF0avF9uGNgwOfbehA22g8PRQmJPPjOJSVuINjVKKioIwFU5wejykDRT26GiCQzYDgy/B8CHw6t+tF0wGOTlzkrXetdS6iiMbMj08pE3FKKKs02gwRYnbws4HGnCqmYk3NEopKgrCWL4UY1u9m1I9TcUIj8FEO1RuKooG4BPRCXpme1hfup4ah/5nPcpcjtTICK5776Xktn2YG/TvngYt63Sww08uK3F4VZLNjYxSioqC8OrwHKlsjj1NOizF8DZCw15wVhVamiUZDA0yGB7klspbdD8zESAzMUGis4v4q69iqanBYC0OBRP2xxECPOV21ev0BkcpRcWqE06kOT40x/oqF3V6K8UYOw6xGajbBTqf/p7Opjk2eYwaRw0t3uJoN5YaGiZ+7BgyncbcUByu00Q0TSqRpXFzGa5y/T94KF4b+v6vV9yQjM0lODOhw6kYUT9k0+CsLoqs04noBKlciu0V24uiNlFmMiT7+kj29GDfugVLvf6bIoDmOu18cRz/aASXiife8KxYmzeF4kr0z0QYDMT4rfvWYTbq6LksOAx9z0Lz/qKIJ45Hx9lUtgmb0YbPpjM39AKkx8dJDQ5iaW7GsnYdwqKjWPIihP0JNu6vwVZiwmzV0RQXxYqgozuS4mYgk81xsM9PndfO2gqddTEJjkBFG3jXgFnfFoGUkhPTJ/h+9/epdlQXxbT61NAQ1pYWvO9+N5ZG/VviAJl0lqnBEMHpOO4Knbn6FSuCUoqKVWUynKRzPExblYs6r47iM+k49D8Hc0NQ0VpoaZbEn/DjtDi5q+Gu4nCdplKk+voI/exnGCxmzLXFUT4SmU2STefIJLPKdXqToJSiYlUZnY3RPRlmR6MHT4mOuoIER2H8VTBZiyKeOB4ZZyA4QKm1lOoS/SvF1Ogo8fYOssEg5tpahKk4IjchfxxPpZ3qtR7sLh1dr4oVQylFxapyfGgOKWFng84G4M4NwlQn1O2GEv3H587OneV7Pd+j3FaOWecNywHSw8MkOjqwF1HBvpSSqf4QJ54awVlqKwoXteK1o5SiYtUIxtK8OjzH+moXdaU6cp3mchAah9ue0EoxdE4sHeP41HGa3c00uPVv1eaSSZJDQ2T8fmybN2Oq0X+TAdAGCtucZjbtr1FTMW4ilFJUrBqjc3HOTIbZWO2i2q2j+ExkAs78WJuMUQyu0+g4PbM9tJa2UuPUv4JJDw9jMFuo+NjHsKxpRui8/vMc4UCCmZEIJotRTcW4iSiOq1NxQ9A3HWFkNs7eFh8mXZVijMCau7UG4C79K5mR8Ahn586yvWI7bov+G5anhoYIP/kkuVisaFynoJVi+Eci2N0WjHqa4qJYUdRfWrEqpDJaKUajr4SWCp1NxZg6Dc/9LyhtBqO+E0AyuQyvTL6C1+Zlg29DocVZklw8jv8fv0Hs2DFMpV5MVfpvnQfaQOFUPMOG22rwVurI1a9YcZRSVKwKk6EEXRNaKUaNnqZiJILaUOF7fl/reapzpmPTTMYn2eDbUBSlGOnxCeKHDmHbtAlLc3PRJKuEAwl6jkwyOxFVrtObDH0/FituGEbypRgPb6nGY9dRtuTcMJz4F6jdCTseLbQ0SzIeHeeWqlswCROVJZWFFmdJ0uNjANqYqCKZiAEQCSRp2OTD5jSrgcI3GcpSVKwKxwZnMRkNbK/3FFqUiwn0gb9Hm4ph1ZlbdwFO+0/z1OBTVJVUYRD6/veVuRyZySkAbBs3YqqoKLBEyyOXzeEfi5DLSnzV+h+ArLi+6Pu/SnFDEIimODkSpK3KRa2epmJkUjB5CsrbtPmJOieYDJKVWZrcTUWRdZqZniF29CgA1mb9D2w+R2QuSXgmTmAsqgYK34Qo96lixRnLl2Lcv6GKKpeO6r1Co5BLw45fAq/+XXvj0XGS2SRrvWuLIp6YGR8jevAgAOYiyjqNBJKUN7gQBqEGCt+EKEtRseL0TIaZCCa4VY+lGIMvQSoGDv279gaCA3zlxFcoMZUUxUDhRHc32UAAAFOpzjoYLcLUUIiul8dxeq1qoPBNiI7uUIobkUQ6y+GBAM1lDprLdBSfkRIm2rVYYsMe0HlWZDKb5NjkMWqdtazxrim0OEuSjURJdJzC0lIcw4/PkYikEUJQ0+pVrtObFKUUFSvKRDBB13iYtmoXtV4d3WSi0xAeh0B/UcxOnIxO0j3brXWxcRRBPHFiHNvmzXh/4RcKLcpVEQ4kSMbSuH02NRXjJkUpRcWKMjIXo3sqws5GLy6bjlLbg8MQ80PT7eDWv1IcjYzSM9fDlrItRTFQODU8TOS55zD6Sin/8IcLLc6yCfnjnDk0icliVAOFb1KUUlSsKH/19FlKLEa21XsLLcrFTHVBx3e1tm4mfU+Al1JyYuoEJmFiS/mWQouzJDKbJXbsGImuLiy1tVQ88ZFCi7QsMuksiUiaXQ82UlqjI1e/YlVRSlGxYgTjaQ72BbRSDD25TlNRGHoZytZB2dpCS7Mk/oSfM3NnWF+6nlqn/ofzZqamSHZ2YW1txVSjf3nPEQkkGe4MEJqJK9fpTYwqyVCsGFOhBACbatxUunR0k5kbBpsXmu8qinjiWGSM22puI5vLUuXQf+/Q9Pg4ybNn8TzyCKaK8kKLs2zCswlcPht2t0UNFL6JUZaiYsXomYwAsLfFh1FPqe2zgzBxUpudaPcWWpolOTt3luNTx6lz1WE26P9mnR4fx/erv4r9ll1FMyZK5iRzEzFcZTbK611F06NVcf1RlqJixXipbwaARp+OuthkM1p9oqsGvPrvshJNRwkkAmRltigK9rOhEInTnQiLBUtdXaHFWTaxcIqpwRCRuSRNW8oKLY6igBTHY5yi6JiLpeie0CzFWq+OCs0jk3Dkb6HEVxSu0/HIOE6zkz3Ve6h16D8+lx6fIPrCC2SmpzHX6L905BzhQILqFg/rdlWqqRg3OVetFIUQf7oSgihuLCZCCXqnNaXosOrIIREeh+lOcFYVzUDhv3jlLzAJE06L/huWp4YGSQ0MYL9lFwa7jh6GlsA/GmGgfYYSj1UNFL7JuZa//hevuxSKG47O8TDxdLbQYlzO+Antd/VW0Hm8K5PLcGzqGG6Lmxav/jvDyHSa2NGjmCorsa1bV2hxls25gcJ2p0VlnSqWF1MUQqwDWoEeKeXZZb7nC8Bu4BUp5UfnrbcBfw2sAU5JKZ+4aqkVukZKyaE+P2vKHexr0VF8Jp2AwRe15SIYKDwVm6Ir0FU0XWzSk1NY6hvw/MIvFJ3r1GgyULXGjcunGoDf7Cz5qCyE+Azw/4C7gL8RQnx6Ge/ZBTiklHcCFiHEnnmbfwv4ZynlvUoh3pgEoil6JiOsq3Ty+6/fUGhxLhAeh6lObdmt/ySQ8eg4Z+fOstG3kXK7/ksbMuNjJPv7MVVWYCzT0cPQEoRm4hz58QAA1hL9Z/cqVpbl+I8ekFLeL6X8feAh4MFlvOc24Mn88pPAvnnb7gHeIoR4VgjxlqsRVlEcTIaS9M5EuKWxVF9TMcLjsOUXtWWLjjJir0DHTAfxTJydlTt1P1AYIN7eTuiHP8RcW1s0JQ25bI54OM2tb16Dr1b/MVvFyrOc/7R2IcS5XPBK4PQy3uMFQvnlIDB/bsxa4MfAG4E/FEJc0YUrhHhcCHFUCHF0enp6GR+r0AMdo0GyWcmOBm+hRbmY/ufh9A8KLcWymEvM0T7TTqu3lTqn/q3a7Nwc8ZPtWNetw1Kr/yzZc0TmkkwNhIiHU8p1qgAWUYpCiHYhxEk0y65fCNENDAB3LOO4c4A7v+zOvz5HEHhOShkFzgJXbNEhpfyqlHK3lHJ3RYX+590pIJeTHO4PsKbCQbWeWrslI1BSDluKY2rDWHSMreVbubvh7qKoT0yPj5Ps6cHath5Ttf7lPUckkCCbzSGEoMSjlKJiEaUopdwqpdwmpVwvpbTP+922jOO+DNyXX74fODhv20vANiGEEWgGlAl4AzETTdIzFaa10km5Q0c3mdAYHPmqNlC4CBgIDjAZm2SNZw02k44eLq5Aenwc5+teR8mePRisOvq7L0HIn6B6jYeqNW41UFgBLC/RpkII8X+FED8SQnxeCFG51HuklK8ACSHEASAnpTwshPhSfvOfAX8KvAj8nZQy9Vq+gEJfTAYT9M1EubXZp6+bTHAYdj0GjXvh7t8rtDSLkswmGQ4P0x3oLoqsU5lKkTx7FqPXi6W5udDiLJt4JMXw6QC9x6dVKYbiPMspyfhH4OvA54H9wDeAh5d60/wyjPzrJ/K/x1leso6iCDkxModRCDbXuZfeeTXpfQbO/je861uw7r6l9y8g45FxfHYfD615qCiUYnpigsjzBzD5fHh/sTjc06BNxahZ5yGbzuFUSlGRZzmJNiVSyu9IKUellN8GiqdNhWJVyeYkRwdmaalwUO3R0WUSC8BkO/hawKX/eNdoZJR/6PgHpJSU2kqXfkOBSY1p8UTb1i2YSvUv7zlmJ2NMDYZx+myYLWqgsEJjOZaiXwjxcbS44G1AYGVFUhQrM5EkZ6cibK7zUObQ0eDe8DiEJ7UkG6u+0+5zMser068STUdp8ei/iw1A4sSrCCGwbdtWaFGWTSadJeyPEwullOtUcRHLsRSfACTwy/nf711RiRRFy3gwno8nluqrTs3fB9vfDU3LSZwuLDPxGTr9nVoXG6f+XacZvx+MJtxvfnNxlWIEktgcZtZsK1cNwBUXsayYopTygRWXRFH0HB+aw242sqnWU2hRLpDLwdgxGHgBNr+t0NIsyXh0nJ7ZHh5ofoCqkmIYKDyBzGSwtK7DXKV/ec8R8sc59rNBdjzQoAYKKy5iOZbiq0KIN+ezUH1CCN+KS6UoOjLZHMcGZ2mpcFLt1tGTd3QaDGbY+k5w6d+S6Z7tJpAIsKtyFyaDjqaLXIFkXy/+L38Zo9uDsOjIZb4IMieJzqXYtL8GX7VDX14NRcFZzn/d7vzPOSRw78qIoyhWpiNJeqci7GoqpVRX8cQx6H0atr0DTDqSawEiqQivTr3KGs8aGlwNhRZnSXKJBPGjxzA3NmJt1v/A5nPEQinmpqJYS8y4ynSUEKbQBYsqRSHErcCXgU4pZfvqiKQoRkZn4/T7o3z4Xp2NDJo6A2vuhro9S+9bYMaj49Q6a6ksqSyOUozxcZLd3VhbW4tqKkbIn6Dv1Zl8PLF4Gg0oVofF2rz9JfApYCfwZSHEx1ZJJkURcmxwFrfNzPoqV6FFuUA2DX1PQ/9z4Nb/TXskPILNaKPN11YUA4Uz4+NY1q2jZOdOjB4dxZGXIOyPMzcRo3adVw0UVlzGYlfELVLKN+SnY7wOeNvqiKQoNlKZHMeH52ipcFLl0lE8MTwB02egYgM49Z0Eks6l6Qp08cO+HxaFlSilJDk0jG3zZmybNhZanGUTD6cY7ZkDUAX7igVZzH0qhRClwLkotPHcaymlqlVUnGcqnKBvOsLta8vw6GkeXXgc7F5ovgMM+i7OnopNUV5SzqMbHy0KpZidmSH63HPIVArXXXcWWpxlE5pJ4K3U4ohqKoZiIRazFDPA9+b9ZIDvA99dBbkURcTobJxBf4x9LTobLDv2KpSth6rNhZZkScYj4/y478dkcpmiGCh8fipGayumIirFmJuKkoxlADVQWLEwV7QUpZSvW01BFMXL4YEAZQ4L6yp1FE9MxyEyCYkguPVfitE+0053oJvHtz5eFAOFE909ZKanKdl9C8Kk/9IRgEQ0TWgmwXCncnQproz+//sUuiaRznJiOKjFE906ckeFxmBuCOp2QYnOLNhLCCQCdPg7aC1tpc6l/4HCuViMrN+P8957MTc2FlqcZRP2JyhxW9i0X/8PSYrCoZSi4jUxFUrSNx1hU60Ll01H7qjwBGQSmpWo8+LsscgYh8cP01raWjQDhQ1uN459+zAXUWu34FSMjudHyeVkoUVR6JjFSjJKhRDVl6yrFkJ4V1wqRdEwPBtjZDbO3jU6s8aGDmrNv8vXF1qSJTk7d5a55By3VN6C1agja/sKpEbHmP7CFxAGA0an/ktHAFKJDPFomtpWL54KVbCvuDKLWYpfBS5t6VYG/O3KiaMoNg71+6lyW1lT7ii0KBdIhGD4IESmdF+fGEwGOTl9EoBmT3NhhVkGMpcjfvQohpISLGvXFlqcZROaSZCIpPFWluAus7Pnjc2FFkmhUxZTilVSytPzV0gpTwHFk2qmWFHiqSztI6F8PFFHNV/hca0+sXo72PRdVD4aGcVu0iyXYijFyExPk+js1LrY1Opf3nOEZuL0vTpNIpLG5bNx65uLYyyXYvVZTCkqx7tiUSZDWn3i5lo3DquOMhBDY7DmLmjeX2hJlmQkPEKjW0tW8dq8hRVmGaTHxjCUlGDbsgVTRUWhxVkW6VSWeDjFlrvqqGhyYTSrVArFlVns6nheCPGnQggLgBDCKoT4LPD86oim0DuD/ijjwQT7WnQ2OKXvWUiGobS50JIsSjQd5bT/NN/q/FahRVk2mfFxHLffTsmuXQijvhsinCM8k2C0e5bZiShu1QBcsQSLKcVPA1G00VF9wCtAGPjkagim0D8H+wPUem00+nQUT4z6NZeprwVc+nbvjYRHqHPV8e62dxdalGWRjUQIP3+A8NNPF5fr1B/HVWbD5jTjKteRm1+hS66oFKWUGSnlZ6WUm6SULVLKzcD/A35t9cRT6JVoMsOp0aAO44ljcPo/wdsAlpJCS7MoI+ERnh95HrNBR6Usi5AeGyN55gymsjLM1fovHQHIpnMEp+KUuCyU1ToxW4rDulUUjiWd60IIuxDiXUKI7wH/DJSuvFgKvTMZStA3E2VrnRu7nm40oXFovhPKNxRakkWJZ+K0z7TTMdNBjbOGD27/YKFFWpL0yAip/n7sO3dicOjIO7AI4UAC/2iEwVMBNTtRsSyumB0hhHgD8B6gFvh3wCulvG+1BFPom76ZCFPhpL76neaycOanmrXo1feQ3rHIGKcDp9ng20Cts5Y76u4otEiLIjMZkr192DZuxLa+tdDiLJvQTJyqNW6qALdynSqWwWKW4veAGeAXpZR/BSRXRyRFMXCwN0BDqZ36Uh25KCNTMNkOVhe49O3eGwmP0OXvYlPZJuqc+m/tlpmawuh243nb2zAVyUDhbDbH9HCYgZMz2JxmLDYdZUgrdMtiSrEaaAe+JYT4DlAlhNB/uw3FihNKpOkcD7G2wkmlrvqdjsJMDzTuA6N+43SpbIqOmQ4kkq1lW4uii016fJzZf/kXcrEYpnL9T/EAiASSGAwCp8+msk4Vy2axRJuglPJrUso3Ah8BvgH8txDiyVWTTqFLpvLxxO31HqwmPcUTx2D/x6BmZ6ElWZT5rtN6V32hxVkW8Y5TpIeHsbauQxiKo84v5I+DgPJ6p3KdKpbNYr1PdwghfiqEOAX8NfCfUsq7gMdXTTqFLjk7FSEQTbG7WUf1iZkUdP4Q+p8Dj77dkWORMWxGG1vKt1Dr1H9D7WwwSOLkSSzr1mGu0/e5PUcuJ5mbjHH8v4YwWY1qdqJi2SzmZP9r4LeBV4H7gb8E3iql7FsFuRQ65uU+P01lJdT7dBRPDI9D/W4QBnBUFlqaK5LJZZiITrCvdh8+m48Ss47O4RVIT0yQDYexbdiAuUjiidE5zXV665vX4K3U/zlW6IfF/CBJKeUhKWVSSvljQEcTZBWFIhhL0zUeZm2FkwqnjmJhoXF49nPa7EQdu/fGo+O8MPYCL42+RL2zOFyn6bExXPfdh+OO/RhsxeGGDM3EGezwk4pnletUcVUsZimuF0J8Mb8s5r+WUv7Wikum0CWT4QR901Hu31iFxaQj5RPohbt/T/ejokYjo+yp3kMqmyqKrFOZShE7coT0yChlv/EbhRZnWUgpCQcSmG1G7C4zdqel0CIpiojFlOKjl7z+7koKoigOusZDRJIZdjXpqIdDKgZnfgypKGx7Z6GluSLZXJbuQDeDoUH21uzFadH/LML01BSJjlMY3O6iae0WC6ZIxjLUrPVQWl0cTQYU+mGx7NPn5v8ATfOWFTchUkoO9gdYU+6g1qsjl1R4HKa7oXITlOgo+ecSpmJTtM+0c2b2DHUu/VuJAKnRUZI9Pdi2bMbo0++5nU9oJs7ZY1NMD0dwlenoOlUUBVfj/3rfSgmhKA5mY2m6J8KsrXDoK54YHIF0TGvvpmNGIiN0Bjq1UowiiSemenoQZjP27TsQQhRanGUR8ido21tFzRo3JW7lOlVcHVejFHtXTApFUXCu3+ktTT5MRh3FE6e74NbHoWxdoSW5IlJKzs6dZSI6wfby7Xis+h5+DJCZnUWYLZR96ENYisV1GkoxMxwmMBbFXVFSNIpcoR+WfWeTUv4GgBDiKysnjkLPnBoLkkxn2dmooxt6fA7GT8LYK+DW7417Oj5N+3Q7a71rafY0F1qcZZEeGyPy7LPkwqGiae0W8scRBoHBKFTWqeKauJbH/ZbrLoVC9+RyksP9AVoqnNR4dNQyKzyuzU9se4PW81SnjEZGz7tOi6FgHyA1MED0xRcx19VhsBSHGzI4pc1OLKtz4fDoyMWvKBquRSnK5ewkhPiCEOKAEOIvF9gmhBAnhBDvv4bPVxSAQCxFz2SEdZVOyvQUT5wbhv5nwan/BuBGYWRb+TbK7DqaLHIFcskksaNHMTc0YG1uLrQ4yyIRTTM9HObU86O4ymwIg3KdKq6exUZH/ZDLFaAAti91UCHELsAhpbxTCPFlIcQeKeWRebu8BZi6FoEVhWEiGKdvJso7dtdj1MvNRkoYPwFtb4TSpkJLc0X8cT+JbIJ3bXgXa71rCy3OsshMTJDqH8Da2lo8rtOZOKVVJbjLbbhV1qniGlmsTvEjr+G4twHnGoc/CewD5ivF9wD/9hqOr1hl2kdDSCnZWqejeGJ0BvqeAUc5uPR74x6NjPKTvp+wtWIrd9ffXWhxlsX0X/8NrocfxmAxYyrVUU3qIsxNxRk6HaBxsw9nqY68GYqiYjGl+HdSygeu8bheLmSrBoHN5zYIIR4CngMyS3w+QojHyTcgb2xsvEZRFK+VXE5yJB9PrNZVPHEMZrqh7fVg1q9lMBYZY1/NPkptpVTYKwotzpLIdJrQD36AMJnwPvK2QouzLFLxDLFQEiklbp8dg56yoxVFxWJXzmuZCTQHuPPL7vzrc7wf+PpyDiKl/KqUcreUcndFhf5vJjcqM5EkPVMRWquclDl0lHAxNwwVG6H+1kJLckWCySB9wT4kkhZPS1GUCKRGRgBIj45iqimOpKCQP47JZKBpcxkulXWqeA0sphS3CSF+cMnPD4UQP1jGcV8G7ssv3w8cnLetFfgP4HeAjwkhNlyL4IrVYyKo1Sfe2uzDoJd4YjYDI4ehZjt49FsIPxIZ4cWxF+kL9hXN7MT00BAA1rb1mKv0O3FkPsHpOEd+PEAqnsXlU0pRce0s5r7sAZ64loNKKV8RQiSEEAeAE1LKw0KIL0kpn5BS7gAQQrwPMEkpu67lMxSrx/HhOSxGwRY9xRMjk5BNAwKcVYWW5ooMBgfpme3hoaaHqCzRv4LJJRIkBwcBKNm9G2HW/xzCdDJLIpJm10ONeKscGPXUqF5RdCymFBNSysFLV+aL9z+w1IGllB+95PUTl7z+h2XKqCggmWyOY4NaPLHKraMn8NC4Fk/c9T4wLhqaLhiRVISTMyfx2Xy0+dowCP3frNPDwxis2t/Z0qTfjN75hPxxZkYiOLwWVbCveM0s9l/6+SusV8X7NxEzkRRnp6K0Vbko1dP08rlBrVjf21BoSa7IaGSUrkAXG30bi2JMFEBqaIjYES1RvFgGCodnEiRjaXJZqRqAK14zi03J+OGVNq2QLAodMh6M0z8TZW+LTz9JIpkkDByAknJw6zcRZCg8RFegi81lm6l26Lu5AEA2EiXRc5bw008DYHS7l3hH4cmks0TmEtStL6Wi0Y3J/FryAxWKFSreV9w4HB2cxWE1saFaRzfI8DgMH4bG28Chz6zkWDpGf7Cfakc1m8o2YTLo08U7n/TwEImODmybNhE/erTQ4iyLsD/BYEcAk9nA7jc0F1ocxQ3AShXvK24A0tkcx4dmWVvhoMqto2LouWEI9MJdvwt6sV4vYSwyRqmtlEc3Plo8sxMHh0icOoX74YexbdlSaHGWRcifoHFTKdmMxF2uoxpaRdFyRaW4UJKN4uZiKpykd0pznXpLdFSfGBqFnb8MFfqt5hkJj/C97u/x5pY3U+vQr4v3HNm5ObLBIL5f/mWMZT48b3lLoUVakmwmh380QnQuSV1rKWarcp0qXjv6T4dTFIyxuTgD/ij71uiogXUyAj3/BdmkbkdFJbNJpuPTbK/YTr2rHrNRRwlKVyA1NESis5NcPI61qUk/8eNFCAcSpGIZYsGUKthXXDeUUlRckaMDAUpLLLRWOQstygXC41C5SetkY9dnT86xyBhT8Ska3Y00uPWbHTuf5OAQoZ/9DJlOYy6SUoywP4GrzEb9hlLlOlVcN5RSVCxIMpPl1eE5WiocVHl09BTu74Xj34SqTYWW5IqMRkb5bvd3mYnNFEUpRmZmhmRXF8Jsxrq+tSgagOdykunhMK/8fBCL3YTVrv9EJkVxoK4kxYJMhbR44t1tFbhtOnL/xQOw9V3g1mfLtHQuzamZU6SyKbZWbMVq1FGC0hVIDWkJNvatW4umYD8SSGCxGdm4vxZPhbISC00ulyObzRZajAUxGo0YDMu3/5RSVCzI6FyModkY+1p8hRblArEAdHxfK8PQaTxxIjLBKf8pNvg20ODSv+tUSkl6eBjX/fcjs9miUYphf4LAWBRvVQnuMqUUC0kulyMajWKz6cijNI9oNIrD4Vi2YlRKUbEgh/oCVDittFToLJ44dQpu+TWwOAotzYKMREboCnTxQNMDxeE6nZwk2ddPLpHAsWc3RqeO/t5XQOYkIX+cyYEQlU1ubE4deTJuQrLZLDabDbNO++RKKclms8tWiiqmqLiMRDrLyZEgLRUOqvXU73S6G+Kz0HR7oSVZkGwuS6e/k2AyyI6KHZSYSwot0pKkBoeIHjxIZmoKc5HMLI0Gk+Sykk37a/FW6f8cKwrL1WZSK6WouIypUJLe6Siba904rDpxJuRykJiDe/9It6OiJmOTDIYH2Vi2sSiyTmU2S6qvj+SZM9i2bMbSoH+ZQSvY73x5nJA/rrJOFdcdndzxFHpi0B9lbC7OvhYd1SdGp+HUv4OrGva8v9DSLMhoZJQt5VvYXLa5KFyn6bExEmfOYK6rw7p2LQa7/hWMlJKwP0HzljJK3FZK3DpqKqG4IVBKUXEZB/v9VHtsNJfpKG4XHofWB8FTByb93QhzMseZwBleGHuBR9Y9gsviKrRIS5IaHMLa1oalsRFLY3Ek2MRCKYLTMYwmA95q5TpVXH+U+1RxEbFUhlOjIa0+UU/xxKnT8NJfgk7dklOxKbIyS52zjnqnPt2785GpFKn+PuLHj2N0u7DU69+yBS3rNDgdJzQTV1mnNyDvfve7ede73sXevXtpamrixz/+8arLoJSi4iImQ0l6pyNsrfNgt+ikl2Q2A9EZuPN3NUtRh4xFxkjlUrSVthWF6zQ1Mkq8vYPk2bOYamoQFv1Z3wsR8sepXeelZp0Xh6c4ZFYsnxMnTtDS0sKhQ4f41re+xac//elVl0EpRcVFDMxEmQwl2aun+sTIBHR8D8ZPgLOy0NJchpSSntkevnLiK1iNVrw2b6FFWpLU0CCJ9vaiKtiPR1JM9AYZaJ/BVWpDGPTfn1WxfOLxODMzM3zyk58EYNOmTczOzq66HEopKi7i5T4/daV2Gkt1FE8MjsJMDzTtB4NOrNd5+BN+2v3trPGsodnTXGhxliQXj5Ps6yM9MYFt0ybMNfpshHApYX8Cd7mdikaXagB+A9LR0UFra+v5JgCvvPIK27ev/vhepRQV5wlEU/y0Y5y1FQ4q9TQ/cboTSpugemuhJVmQ0cgonf5ONvg2FIfrdHgYo8NJ+Uc+gmVNM8KovweNS5E5yZGf9BOPpHD5bDhLlVK80Thx4gRDQ0MkEgmi0Sif/OQn+fjHP77qciilqDhPx2iQ4UCcbfVebGad3CjTCc1S3PZucOtzLmHfXB8DoQG2lm+l3F5eaHGWJD00ROzwYYTZXDSu07mpGAMn/Ix0zeIstWFQrtMbjhMnTvDoo49yzz33sGfPHj74wQ+yf//+VZdDlWQoAJiNpuieDANw6xodTUkIj8FU3lIs0VHdZJ7ZxCwnZ07S4GxgnXddocVZkmwkQuJMN5HnnsNxx35MlfqL0V6KlJKZkQgAm/bX4lau0xuSEydO8Ld/+7f82Z/9WUHlUJaiAoCOsSADM1EA6kt1VP811QXBYVh3H+hw8O1YZIzeuV42lG2g3qX/UozUoJZgY9u8GUtTc3EME/YnSEbTAOSyEqdPKcUbkd7eXlpbWwsthrIUFTAXSzEciHN6XLMUK106uelEpmD6DIy/Ci59uk5HIiO8seWNmA1mKuwVhRZnSdL5MVHut70NS1Nx9DqdGY6cX3ZX2DEa1bP8jcjo6GihRQCUpagAOkZDBKIpOsdDAFhMOrksRl+ByKS2rMNRUaFUiNP+0/TM9tDgatC91ZWZnSUbCuN77DGsbesxlenPHX0pkdkk8UiKuakYAO4ynTywKW5YdHL3UxSKYCzNUCDGs2em9DU7MTSutXbr/qn22uYprDwLMBYZQwiBURiLootNemiIxOnT5JJJrE3NhRZnWcyMhMmks4z1BAFwKaWoWGGUUrzJ6RgLEoimOD0eYv86HWVOjr2idbGJTBVakisyFBqi0l7JBt8GKkv0n7CSHBgk9JOfQDZbFK7TeDhFdC6J0WRg852a+9ykl6xoxQ2LUoo3McG4ZiU+c2aKvWvKcOplTNTc8AVluPUdhZXlCsTSMV6dfpV/6vwnapw1GHXYVGA+melpEp2dGBwOrOtbMbrdhRZpSWaGI8RCKV797+FCi6K4iVBK8Sbm1GgQf0SLJd6xrpw15Q4+el/hs78YO67FEnt+DrU7Cy3NgvSH+nlq8CnWedcVR8H+0BCJjo6iaeuWjGcI+eOUuC1surMWo8nALW/Qv9yKhfnCf3cXWoRlo5TiTUookWYwEOOZrin2tZThsBrZUufh4w+sL6xgs4MQ84PVA2vuAruOaibzJLNJuvxdjEXH2FK+hRqH/pKA5iNzOVLDw7gfegj7rXuKYpiwfySs9TrtC2K1m/DVOtj3lrWFFktxjfzlUz1Xtf/Y2Bi7du3CZrORyWSu+/6LoRN/mWK16RgNMhNO0jkR4nceaGNtpRNHod2nUmpWYnwWjvydNhXDrT8r7LT/NEPhIQA2l23GZND3v1FmYoJkbx9kMpTcugdDiY7qUBcgncoyNxlH5gABwiAoq3UWWizFMvnnQ0PLWv9Le68c1/b5fDz11FM88sgjy/rMq91/MZSleBMSTqQZ9Md45sw0+1rKcNqMbK7VQYwp0KcpxEwCbn9Ca/5duRHu/r1CS3aeSCpC31wfnYFOgKJxncZefpnM9HRRuE79IxHSqQzWEhPVazz4ahwYzepWdTNhs9koLb3YS/SZz3yGe+65h3vvvZeBgYEl979W9P2Iq1gROkZDzISTdOWtxHWVTkosOrASx1+F2QE49g9w7x9CxQawOuF1v19Y2ebRPtNOVmY5EzgDQK1Tn00FziEzGZJ9/SR7enC/+U2Y6/VdOpJN55idiDIzHCE6l2TdLVWU1SkrsZj4xL+3L2v9J/69nYH/9cZlHbO9vZ3R0VGeffZZOjs7+dznPsdXvvKV1yzrQiileJOhWYlRnjkzxW0tZbhsJjbX6qAG0N8LiRCc+Qm0PgQmK9RsK7RUFxFIBBiJjCCE4H2b38fnj30ei1Hfg27TY2MkTp/G0tSEtWUtBp0PEw6MR8llJdUtHqSUeKvsmK36zuxVXMxnH7l8ms0n/r39svWLuU8vpbOzk2effZZ77rkHgJqaGr75zW/y93//9zz88MP83u9dP2+SUoo3GafHQkyHk3RNhDUrscpZ+IkYuZxmJWYzWgyx4Vao3ARme2HluoT26XYSmQTPDD/Dw80PF1qcZZEaHMK2eTPWlhbd1ybmsjn8YxFGu2ex2E1UNLgoq3cVWiyFDmhra+PBBx/kS1/6EgDpdBqz2cwv//IvX/fPUkrxJiKazNA/E+XpriluW1uGy25iU40OYokz3ZqV2Pmf0PqgpgyrtxRaqosYj4wzFZ/CIAw0u5sxCP3HuHKpFMnesyS6zuDYswdzrb5dvbMTMbLpHKXVDoQB3OV2rHZ1iyo2FrIAP/Hv7VdlGabTaV7/+tdz4sQJHnroIT772c9SXV3NPffcgxCC97znPTz++OOL7r93795rkl9dcTcRp8ZCTIWTnJkM87vb21hf5dKBlZiFiZOAhLJ1YHFA1WbNfaoTpJS0z7QTToUJJUNsLNtImU3/fUPTIyMk2jtIDQ3hfdvbECb9/rvLnMQ/GiE0E8fmNGOxmShvULHEG4WrrX82m808+eSTF63bu3cvf/AHf7Ds/a+VFXvcFUJ8QQhxQAjxl5es/6QQ4uX8z30r9fmKi4kmM/RNR3ima4rb12qxxI01OnBNTXdBKgqHvgLOKjCXaEpRRwyGBgmmgjwz9AxdgS4A2nxtfHD7Bwss2eKkBgeJt7dj27ZN967T4HScdDJLYCxKPJzG4bVid+o7/qlYPgWvf74KVkQpCiF2AQ4p5Z2ARQixZ97mb0gpbwNeD3xyJT5fcTmnx0NMhTQr8fa15bRVu7CaCmwlZjMw0Q7JMGx/t6YUq7eC0VxYueaRzWU55T/FVGyKnrke9tftp8JeQa2zlg/t+FChxbsiuViMZG8vmZkZ7Js3YaqqKrRIizIzEkFKSfO2cjwVdioadfDAprgpWSlL8TbgnC37JLDv3AYpZX9+MQnIFfp8xTxiKc1KfPrMFLevLcdtN7OhWgexxKnTmpX48l9DLKC5Tis2FFqqizg7d5ZYJsZTg09xR90d2Ew2tpZfnl2nN1LDwxhdbso//GEszc0Ig35joCF/nGQszatPDhMNJrG7LDg8+nGfK24uVuo/xQuE8stBYKGqyk8BixaaCCEeF0IcFUIcnZ6evq4C3kycHgsxEUzSPRnm9rVlbKh2FX5mYiYFkx3azESzHXwtULMdjPqJe6WyKboCXfQH+5mOT3NL1S3UO+sps+s/npgaHCR66BAGqxVLo75dp+eGCG++o5YSt0XFEhUFZaXujHPAOVPEnX99HiHEI0CZlPKfFzuIlPKrUsrdUsrdFRX6n2quR+KpLL3TEZ45M8X+deV47GbaqnXgmpo6BZkkJOZgwxvB5oZyfcUdOgOdJLNJnhp8insb78VisBSFlZgNh0l0nSH6wgsY3G5M5ToaCXYJ0WCSWCjJ0Ck/BpPA5rDgLtNXKY7iOvDM5wotwbJZKaX4MnAuieZ+4OC5DUKIbcCH8z+KFeb0eJDxYIKeyTC3tZSxocaF2VhoKzEJk6e1rFNfi5Z1WrMDdOTii6Vj9M71IpHsrt7N5rLNrPWuxWnRvxWTGhwi0d6OfetWrM36bus2MxxB5oMoBqNBWYk3Ks/9r6vavaOjg9tvv50777yTxx57DCkXj7QdOnTo/P4f//jHX4ukK1OSIaV8RQiREEIcAE5IKQ8LIb4kpXwC+HOgCvi5ECIopXzrSsigyFuJU1Ge6dKsRG+JmbYqHViJEx2QTYHRqiXV2DxQpq8JCKf8p0hlUwyGB9lRuQOzwczGso2FFmtZpAYHSJw6hfed79R1r9NENE1kNkE6maVxcxlmmxFPhbISi56jX1/e+t2PXfEQbW1tvPTSSwA89thjHD16lD179lxx/6amJp5++mlsNhuPPvoo7e3tbN16bV6dFQvgSCk/esnrJ/K/H1qpz1RcTOdEiLFgnLNTEd66o46NNW5MhbYS0wktwWbwRcAAFW3azEQhCivXPILJoFaGkQxyxn+GFk8LG3wbsBr1n/yRmZ0lF4ni+/Vfx1xbi9GjgxZ+V2BmOEw6meXkMyPsfriJ8joXQkfXgaJwmM0XMtCtVisNDQ185jOf4emnn8ZgMPC1r32N5ubm8/tUV1efXzaZTBiN155Zr5+sBsV1JZHOcnYywtN5K7HUYaa1UgeuqYmT2hSMqq2QS0OJD3xrCi3VRZycPkkqm8Jj9fD6ltdjN9lZ511XaLGWRXpwkHhHB+a6Ol27TlOJDMHpOGarkd2vb8JsMeKt1vdIK8UyOfNTbUA4aH2M568/R8/PITyxaLP/H/zgB3ziE59g/fr1TE9PL6sh+MmTJ5mZmWHTpk3XLL5SijconeMhxubi9E5FeGRHHZtqPIW3ElMxrVi/7zlIhmDzI1C7q7AyXcJUbIqJ2AQHRg+Qkznub7q/KGYmAmSDQab/31eIHzmC7/3v13XW6cxIhHg4xdCpAG37qvHVOTEYlJV4Q9D2eu1nPj0/v3hd2+sXdZ8CvOUtb+Etb3kLTzzxBKdPn16yIXggEOAjH/kI3/72t1+T+Pr/T1dcNYl0lp5JrS7xjnXl+JwW1unBShw/AckI9D0Dt30EnJXg1dcU+PbpdkLJEMcmj/H4tsfxWDw0u5sLLdaS5KJRIs8/T/TZZzHX12NrW4/B4Si0WAuSSWWZm4xhsZuoXuvBYDTgq9GnrIrCkEwmsVq1cIXb7WZiYmLRhuCZTIb3vve9/Pmf//lFrtRrQSnFG5CuiTCjc3F6p6N5K9GNsdBP4cmI1vg7OKJlmrqqdWclDoeGCSQDPDfyHDsrd+KxethWsU33ca5cMknk+efJxeIA2LZu1bWV6B+NkoymSUTTeCrs+GodGAtdN6u4fixkAf7oY0tahvP52c9+xuc//3kAWltb+eM//mM+97nPXbEh+He+8x2OHDnC//yf/xOAz33uc9x2223XJL5SijcYiXSW7skwT3ddsBLXVujBSnxVa+eWCMLWt4O7Ftw1hZbqPDmZo8PfwVRsijOBM3xk50eoKqmi2vHanjpXGpnJED1wgGwofH6d/ZZdmBv0ZYGfI5vJERiPEo+kmZuM4aksoaxOWYk3PHdf3bzDt771rbz1rRcXJvzBH/zBFRuCv+c97+E973nPNYs3H/V4doNxZiLM6Gycvpko+1p8bK71FD5WkwhqQ4TTMUhp3Uv0ZiX2zvUSSUd4aqh42rnJXI7oyy+T8QcAiBw4AIC1qQmDVZ+ZsoHxKNl0Fk+FnaYtZZRWOzAVelKLYuVZJKFGbyileAORzFywEu9cV06Z00pLuQ6ewsdehfgc9D0La+7W4ohO/XQoSufSdAY6SWVTGIWR3dW7aXQ1UmpbqDuhfogdOUp6bByA6EsvEX/1VQDduk5zOUlgLMpo9xyDHX6EEMpKVOgOpRRvIM5MhBnOW4l7W3xsqXUX3kqMz8JsPxgtULVFq0fUmZXYHegmkUkgELyz7Z1YjVbdW4nxEydIDQycX44ePIgvn3RgrqsroGRXZm4yRiaVpXa9l9pWL54KOxabiuAo9IVSijcIqUyOMxMXrMQKl5U1hbYSpYSRo9oEjJf/SivU963RahN1QjwTp3u2m+HwMN/t+S4Aa71rKTHrt2YuceYMia4z518by8rwvfe9GL1eAF0OE5ZS4h+JMD0cJhJIYLYaVUs3hS5RSvEGoXtSsxIHZqLsayljS62n8FmTw4fh8N9Cz39pCtFg1LrX6IjT/tNkZIZGdyNvb307FoOFjT79tnNLDQ4Sf/WEtjw0RPjJJ7HU12OqrCywZIsTmomTSmQwGATCIHCV2bCW6GdupmJl+ZtX/6bQIiwbpRRvAFKZHF0TYZ7unOTOVs1KbCorsKUz0a61c+v5OQT6Ye29WuNvm37ajoVSIfqD/RybPMbJ6ZOYjWY2+DZgMepz4nt6YoLo4cPa8uQkc//6rxf1NrVt2Uz5h/XZZ39mJEIylsZX68Dls1HRoIMevIpV48snvnxV+4+NjbFr1y5sNhuZTGbJ/XXfEFyxunRPhhnyxxjwx3j7LQ1srSuwlejv1dym57jzt8Fk0+oTdUTHdAeJTIJnh5/lvRvfS4mpRLft3DJ+P9EXX4ScJBsKMftP/4Tr4YextrYCYG1txb55M/bNmwss6eVEZhMkIml6j09Ts9ZL/YZS7C59PngoXjvf6f7Osta/Y/07rngMn8/HU089xSOPPLKszyyKhuCK1SGd1azEp7o0K7HKbaWxkFZicBQGXtCWZU77nctp7lOrfmJIM/EZRqOjvDj2Iuu866hyVLGlfAtGg/7KA7KhEJEDB5CZLAAGpxPPI49gbWkBwNLURMkufbml53NuiPDG27W61PJ6ZSUqFsdms2Gz2S5apxqCK5ZF92SYQX+UIX+Md9zSwJa6Arono36thZvMQcwPp/5DW291QvW2wsm1AJe2cyu1ltLk1l8D7Vw8TuS555HJFLlkkrnvfAfPm998XiGaa2soufXKI3UKzfRQmPbnR0hGM6zZVk5ptQNnqT5rKBXXh2eHn+X5kecvW//5o5/nt3f/Ns8OPwssbileSnt7u2oIrliaTDZH13g+47S1nGqPjQZfgazEZBjO/jdk0xCd0Uov1j8Ek+1QuREs+snmHI2MMpOYuaydm96QqRSR554nF4shMxnm/uVfMJaVYXC7ATBVlOO47TaEjoYzn0PmJGM9c8xNxRjqCLDj/gasJSaVcXoTcE/DPdzTcM9F6z7z8mf47d2/fX771dLZ2akagiuWpnsywsA8K3FrfYGsxHRCyzBNx7XuNYe+DK0PQsNebXu1fmr+pJS0z7QTSUXome3hgzs+SHVJNZUl+srelNkskRdeIBsMAhD60Y8QdjvuN74RIQRGrxfnHXfosvwik84y3DlLLJg8v85dbsfqMOMuV0OEFVdPW1ubagiuWJxMNkfXRIinOqe4c30FNV4bdd4C3HCyGTj7JCRCWl3i0b/XlOE5hQhg0oe7bDo2zYnpE/yw94fc03APH975YWxGm+6sRCkl0ZcPkpmeQUoJgGP/foylpQiDAYPTgfPuuxAW/SWrJGNphk4FSCW0jMFcNh9XRqqM05uEhdyin3n5M1flLk2n07z+9a/nxIkTPPTQQ3z2s5+lurpaNQRXXJmeqQh901GGZ2O8c3cD2wphJUqptW6LTkM2pWWd7ngvOPIt3Eqb4O7/ufpyXUIkFeHk9ElGo6MAPD/yPAZh4K76u2hyN+Gx6qdMBCB+9CjpUU3WyDPPIIxGnHffDYChxI7z7nswXJKEoAeic0mGOwNkM5oijIVTGPMzPK0lZjwVykq8Wfng9g9e1f5ms5knn3zyonV79+5dlYbgSikWIZlsjs7xEE93TXFXawV1pXZqPAW44Qy+BMFhLbEmEYSZM7Dpbdo2dx2suUerTywQ6Wya04HTnJ09S47cRdvWeddhEiY2l+mrhCHe3kGyrx+A6MGDJDo6KPv1XwdAWCw477oLo1N//UJnJ6KMnw2et2wB/KMRbPkC/bI6J6LQLQcVBeNDOz5UaBGWjVKKRcjZ6QtW4rv2NLC9EFbi6CvafESA7p9BSdkFheis0pRhgRJApJT0Bfs45T9FMpu8aFvvXC8Atc5a1nnX6aqdW7Knh8Tp04DWuSb64ov4fu3XMDgcCJMR5113YvToy6qVUjI5EMI/Ejm/brw3iN1lpmHDhXZ+pdX6Oc8KxWIopVhkZHOSzvEQT+WtxPpSO5XuVXalTZ+B8ROa+3T0GDTdcSFuWFIG6+4HY2EurYnoBCemTxBKhS5aH0wGeWrwKQbDgwBYjVY2+DYUQsQFSQ0PEzt+HNBmJJobGih7//s1JWgQOPbvx1RWVmApLyaXzTFyZpawPwFoGacIsDvNWOwX//0NRv1lyCoUC6Gu1CLj7FSE3qkoo7Mxbl3jY1u9d3UFmB2EoZe15d6noPdJMJg0pWjzaFmnptVPAAmlQhwYOcCB0QMXKcRUNsUzQ8/w1ZNfxWf38eEdWhu0jb6NmI2F770ps1mSvb1MfvZzIDXlOPM3fwPZrKYQBTj27cP8GjPqrjfpVJb+kzPnFSJAz9EpZkYieKtKKMl3rDFZ9NcMQaFYDGUpFhHnrcTOSe5aX0Gjr4QK1ypmdkamoP95zUIM9Gsxxdt/S6tBtLpg/cNgXl2rNZlNctp/mr65vsvihhPRCf6l819o9jTz+LbHL0qoWetdu6pyXkoukSDZc5Zk71lkMkXkmWewbdrE3L/+K563vhVh1hR2ye7dWBoaCirrpcQjKYZOBciktA47sVAKa4mJ5m1lmOcpwRK3hYaNPoJTzQWSVKEXpr/0V1Q88ZFCi7EslFIsIvqmI/RMRRidi/PuWxtXN+M0Pgdnn4JcBiZPgaNc62lqcYK5RLMQV7FAPydznJ07y2n/adK59EXbRsIjCCGotFfyjrZ3UO+qv+z9BlEYJ0l2bo5EdzfpoSFk9mIlPvf97+N68EGs69cDYN+29XznGr0QmokzcmZWc5XmGTkzS2WjC2/Vhb9/abWDmrUehEFw65v19R0Uq8/MX//1VSnFjo4OHn/8cYxGI+vWreNrX/vaov2cr3b/xVBKsUjI5SSn51mJTWUllDlXyUpMxaDnvyGT0Eov0jHtt6VKa/S9/kGwuVdHFrSONO3T7YTT4YvWB5NBXBYX0XQUgzBgdpovUohGYaTV23rV6eHXg/T4OIkzZ8hMTl22LZfQXJC+X/1VDHYti9jWth7bRn2NsJoZiTDZHzz/erR7lrJaJ+v3VJ1fJ4SgusWDr1Z/GbKK1WH23xbuKHPp+tJ3vfOKx2hra+Oll14C4LHHHuPo0aPs2XPldoZXu/9iKKVYJPTNROieiDA2F+c9tzauXiwxk9Lat6UimiJ84fNw+0fBbAejRbMQ7aWrIspcYo4T0yeYil+sWJLZJC+OvsixyWO8d9N7afO1XfbeOmcd28q34bQ42VqxOh12ZCZDanCQZHc32VD4su25eJz4q68SO3gQ4LxCtKxpxr5jx6rIuBxkTjLeO8fsREx7LaXWVcdsgHkP4yazgfqNPhwefTRrUBQvZvOFeL/VaqWhoWHRhuAL7X+tKKVYBORyklNjIZ7qmuTu9RWsKXfgc6xCMksuB71PQywAkWmY7Yf9H9OsQ4NJyzJ1rHxGZCKToGOmg4HQAJILbjspJSemT/DM0DOs8a7hA9s+gNt6scXqtXrZXrF9Vdu45eJxkmd7z8cLL0XmcshEgukvfhHr2rV4fvEXCfz93wNgrquj5BqfcFeCbCbH8OkA0Xkt206/MEbj5jKq11xw39ucZho2+rDY1C3lZif8zDNEn332/OuyD3wAgIlPfvL8Osc99yxqKQL84Ac/4BOf+ATr169nenp6yYbg8/cvew2Z2uoKLgL6ZqKcmQjPsxJXKZY4cADC41o88dCXofUBTSEKg1aH6Kpa8hCvBSkl3bPdnPafJiMvHjQ6Eh7hJ/0/wWQw8c62d1Lnqrtou81oY3P5Zta416zabMnM7CzJ7h5SQ4MwL+YG+SzTri6ihw9jaWzEdd99VHz0o+etQwBTZQWO2/YVdhbmPFLxDIOn/KTi2rlPRNJYHSZadlRgc154MvdU2Klt9aqyCwUArte9DtfrXnfZ+upPf/qqjvOWt7yFt7zlLTzxxBOcPn16yYbg8/f/0Y9+tOxZjJeilKLO0azEIE91TnF3WyVrK5x4S1bBShw+AoE+bbnrR9B0OzTepk2/aLkbPHWLv/81EkqFODJxhEAiwLPDz57vrD+bmMVlcZHJZbi99nY2l22+SIkYMNBa2rqqJRfpsTEtXjg1fdk2KSXkcsz81V9hcLkoufXW87HC+QoR0Bp8v4Y5cNeT6FyS4a4A2bSWDCSlpPvwJGtvqbjIPVq1xkN5vZp8obi+JJNJrFbtOnO73UxMTCzaEPzS/e32a+/wpZSizun3a1biRCjBo/saV2cSxuQpmOzQYog9/wVbflGzEEEr1C9tXrGPPmcdnvKfIiu1lP/nR57n7vq7EULw3Mhz3FJ1C82ey2Woc9SxrUKLG640MpMhNTBAsqfnsnihlJL0yAixw4cRJhOet74V32OPYXQvnIxkqtR6xQpz4esmE9E0U4MhTj4zwppt5UgpGe8NUtnkZuvr6s4/gBhNBuo3lOIs1V8PVkVhWcgtOvHJTy7pLp3Pz372Mz7/+c8D0Nrayh//8R/zuc997ooNwS/d/8EHH7xm+cX8XoV6Zvfu3fLo0aOFFmPVyOUkvdMR/vznZxj0x2irdvHuPQ3cvq58ZT840Ad9z0Euq2WZTnRA417NZdqwF6qufXjnUkRSEY5MHGEmMXN+XU7m+JODf0KZrYz3b3s/FoPlMveix+JhR+WO6xI3lFIiEwly8TgyHicXj5OLJ5CJC8u5eAyZTBF+5pkF3USBb36TrN9Pya23Yt+58zKLELQ+ppamJqxrWzB6PAWv40olMkwPhpmb0pJpnv3WGe55tA0pJUOnAtSs9ZzvUmMtMdO4yXdZ1xrFzUk6rZVEmRd5qCvk9b2QfEKIY1LK3Qvtr65qnZHLSfpmopwaCxJNZvmv05N47Gbeu6+RLStlJWaSWuwwPAEvfEErwp9oh0CvZiUC1O1aUYV4dvYs7TPt52OH2VyWzkAnL46+CMAjrY9gNV6c1Wg1WtlStoU1nuXHDXPxONm5ubyCi2sKMBYnl8grwUQClvmcGH32WVyvex2Z2VniR4+SnZ3F+8534n796zH6fAsO/zVVlGNtacHc0HCRq7RQN4xMKsv0cITZ8SiXPiC3PzvC+luradpyIWnBXWants17fvqFQrEciqVwH5RS1A1SSvpnonSMhYgkLk4quXt9BesqXbht18m9lk1rCjA8kU+kCWhdakBr7i0l1O68MBy4eivUbL8+n30JsXSMIxNHLiqzeHnsZQ6OH8Rn83FPwz3825l/o9ZZe367AQPrStexybdpWXHDXDRKamSU9MgwmRn/FS28ayH0058SP3EC+44dOO+7DwBT+cXW/KVWoR7IZnL4RyP4RyPksheUYTqZxT+qNfdu3laOxX5BcVc2ualoVDMRFTc2SikWGCklg/4YHWNBQvksv3Q2x9mpCB2jWqH0rc2lry2WmMtqLdrCY5oijM5o456uxES7lkxjMELFBqhf0Mvwmumb6+PE9AkyMkMkFaEv2Me2im1YjBbeveHd1DhqLntPraOWbRXbcFkWvzlnIxHSIyOkR0bI+AMXbTtn4S2FlBIZjyOzWYwuF7GjR0mPj5Odmzu/j33XLpz33YdhgYG/V7IKC0kuJwmMRZkZCZ9PogEYOzuHs9RKidtKPKy5m1w+LV5oMArq2kpxl6l5iIobH6UUC8iQP0b7aJBgPE0yk+XMRJhTYyF6psLUeOxsrnVzfHiO1ioXTutV/KlyOYjNQOicEpzSFONCRGfA5oXRo+DL9wPd+5ta6zZfCzRd2/TqxYilYxydPMpEdIL+UD9OsxOn2UkwqT0E3FJ1y2XvcVvcbK/YTrXjyo2xs+Ew6eFhUiMjZGfnliWLTKfJBAJk5+bIzs5iaWnB6HYT+NrXNOUnBCV79uC6/35kNoupshJrWxum0lJm/uqvMFddXJaiR6sQtAL82ckY00NhMqksUkoS0QwDJ2fYeHsNNocZk8WIyWxgzfZyBjv8AFjsJho3+bCWFD4JSKFYDZRSLADDgRgdo0HG5hJ0TYToGAvRNx2h0VfClloPb95ee14J/ujkOFvqlri5SqkV2IfHtZ/IpOYihQsW4ewAOKth6CUtmSY6Dbf8GnT9UJuDmIpcGP9kc4O3EZrvvO7ffSA4wJHJIxyfOs6h8UNIKXl4zcNUllRyZ/3ln3cujvhg04MLxg2zc3N51+gI2WDwsu1SSmQyicFmI9nfT6pXm6fo/9u/pfSXfon06Cihn/8cY2kpJq8XmpoQViuet70NY2npRUkyjr17r/i99GgVniM4HWNqMEwslMJkNtB9eJLS6hLK6p3UrNOuLV/N5W3ZnKU26jeUYjSp+KHitXH4h31F0wN3xZSiEOILwG7gFSnlR+etrwX+CbABfySlfHKlZNAbo3NxXjo7zUu9AU6NBRn0x2gpd7C5zsPbd9Vjv8KYHcdCVmJ8TlOAoTFNCUamtZFNMz3aBIvoFLQ+pNUYtj6glVk02KGkHJyV4KgERwXcmk9rXnf/hWO7a6Hlnus6JDieifPs0LP8qP9HHJs8Ro2jhgeaHqDF03Je2Rkw4La6KbWW4rP5KLWV4rF6zjf4Pkdmdva8a/RcOUQuFiMTCCBMJszV1YR++lNSQ0NkAwGMXi/lH/wgMhY7r7BcDz2EsFqxrl9PRb4B93zMtbWXrVsIa2ur7qzCc4QDCYZO+8mksiRjGbpenuCWh5to2lqG2WrEYBB4Ky808RZCYLYZseazShs3+3TTSEBR3Bz58cBVKcWxsTHe9KY3cfr0aSKRCCbT4qpK9w3BhRC7AIeU8k4hxJeFEHuklEfym38P+P+Ak8CPgBteKZ4cmeOfDw9xuC/AWDDOukoXuxpLec+eRqxm7SZttxhwWs24bCacVtP53+PBuHaQRAiCozByRHN1hkahepu2XL5eU3o1OyAZAosdvLuhpFRTegYjlK3TjuO6PE53GWvv095znXhx9EWmY9O8OPYi0XSUX9n0K1SUVOAyuyi1XVCAXqsXk+HyS/JDOz5EJhAg2dtL7JVXSA+PkPX7sba1YfR6mf7iF7X5gz4f9h07MFdXY12/HtuWLZjKyjCUaDd+2+bNsHkzkWefxdLYeFXfQZiMGFwujG639tvjpexDH6Rk187rco6uJ7MTUc4em6LEbaH78AQV+QkWOx/UvnOJ24LVbsJaYsJi136sdhMWmwlh0G4ke97YrBSi4po4dWB0Wes333nlBiA+n4+nnnpq2V1piqEh+G1cUHZPAvuAc0pxG/BRKaUUQoSFEC4p5eXdkleAWCrD337odygLj2JJdzPnfDsvbT6IK+xlZ18N1XOnCbjWIiv81JXP0dN/F/bwv5Mx1RLyNPP09rPc3b6ODdN9VCUmeaX8bcSbx0jOmSibSWBKHCNt389IYzdT5XPsOPEwtRPPMeNspq7MwS83GpFjZlxdhzHkYsy5HuLg7S/Q1L2GllFJbXKYV517aGk8TtIoeKZ/B03BH/Otf9jObEWasxuH2XrobnYEjpA1uDjlS7Bp6wa6juQw+X1Y0s8QdryJ9h3HcIb9tHVb2Bw/w1nbWlJ1JhrX1tJ9MIkr8kMypmrmPC2c2DvAtqPNrAlMU52N8Yr191j3b99l7Gw/mbMT2FPHiVn3MdA6y3RVml0vtbEj2cG4wcdUeSnrbt9B99OnKAkdwyDjhO13cuTeHtad8rJmRFI23EtP0zY8W7w8YN9A4lgNazPHCJiaMDU0UXVbGaf+81ma4h2MCBOTjvW0vWkXwwdOYxrvpSw1Skd8HW3760lMBpg5M05j7izj5btxBDx47U66HvgkbWvSBBMmhqJGGoIphuVW3MYsrlSWrl47mxrjzEWMRBNGYs1vYKTbQZk7jd2So2fMzoaGOIGwiZQ0U99kpW/MTG2LC5vXwdnOGNvub2a8N0huNkddZSmnnx4jbb+X9Sem6Xt1hs131jLaPYvBaKBmrYf250ZYv6eKdDLHSFeADbfVMHw6gKXERFWzmxNPDbPpjlri4RQT/SEyySwlHgvOUhvl9U5e+fkgO+5vJDgdJzAWYe2uSnqPT+OrLsFTVcKr/z3Eroea8I9GCPkTrNlWTs+RSYIzcWrWehjtnmXLXfXsfv0aTemVXFB+RvPSHoBicXUpbkxsNhs228WNIYq9IbgX6M0vB4HN87YZ5YWCqCBQCiyoFIUQjwOPAzRe5ZP9QhiEIGvMkvIZyYoKMsSostViwk7WmyXlXUNWmjB4PcQ8VWScEWLOGnJ4wAbN5iaEG1LCSSRnIWOK4CqvxWBIks4mSIk6yGVwlHipMjrIOKOkHC7qWxzEXFY2b2vmeHwQk9uNwI7NkGSzex3WyhJsMoPZth57PI2tag1GowlrOEUkV4q92klJVSlr6zZirhZkmh8ATNizZlx77sQjj5PwgMFUgiUuWNe4F1NEYo1nsHr3YwtncK6ppGzzVqx9z2Jo2IQxV4LD6mZLza046pJYK51YrTYyXbOUNdaSTCcIxMMYLOsxJ63U1DZT7jNiqchg8W3BlpB4XCWUNdbirOtDlNcjRAZrBraUb6GkNonVlCMWjrB1TRsNOzaDycRooAO3aze5hBVTTRXOujLKGt04rBvJ5gxkZCnOujJK11Qg3BlKSuowHM/h27uNVMqAaU8SX42Z7HgCmxWcJZJqkcRdC6ZYFlNY4ii14JM5SqwG7GYo92axW3Jk7WA2SWJAWaUJT4Ube5mT2jJB+a4qLFEDmYygrNFFvCuAr86JyWKkKjWL3WWhtKqEbDaH3WWhotHF9EgEh9dKZZNL217tQAiB3WWhstGNw2sjm85R3pDfXuPAbDVq25vdODxWTBYDmXSOmZEIvhoHdpcFu8tCdYsHu8uiTaQwgN1loazOgctno8RloWqNG7vLgrvcfv6Y5Q1OkvE0W++uZ9dDTcrSUxSE/pNa840128o59IM+tt6jjW87+tMBdr++mbPHprCWmBa1FC+lvb191RqCr0hHGyHEh4FpKeW3hRC/ANRLKb+Y3/aslPKe/PIPgPdKKUNLHfNm62ijuP7IbBaZzTLzla9Q+cQThRZHobghuLRjzELu03MdkuazHKV4zz338OSTT/L973+fP/zDP6SmRgv/1NTU8IY3vOGihuDneOKJJ7j33nvPu1710tHmZeADwLeB+4F/mLftpBDiNrSYons5ClGhuB4IoxFhNCqFqFAUGW1tbcXdEFxK+YoQIiGEOACckFIeFkJ8SUr5BPC/gW8AduCTix5IoVAoFEXFQhbgs986c1Xu0nQ6zetf/3pOnDjBQw89xGc/+1mqq6tVQ/D5KPepQqFQ6I/lNAQvZJ3i1bpPVVWuQqFQKFaUYspmVkpRoVAoFDcsV+sNVW3eFAqFQnHNGI1GolFt9JjeyoCklCSTSRyOy9sYXgllKSoUCoXimjEYDDgcDt0pRNBaFzocDgxX0bJSWYoKhUKheE0YDIarUjx65sb4FgqFQqFQXAeUUlQoFAqFIk/R1CkKIaaBwfzLc/NuYpfsdun6hfYrB2auu4ALcyU5V+r9y9l/sX2We16vZp0631feR53vq3v/Uvtf63Z1vq9t32I+301SyoqFNhSNUpyPEOIWACnlscXWL7SfEOLolYo2V0vOlXr/cvZfbJ/lnterXKfO9xX2Uef7+p7va92uzve17Xujnu+iVIqvhdW8iBXqfK826nyvLup8ry6rcb5VTFGhUCgUijw3o1L8aqEFuMlQ53t1Ued7dVHne3VZ8fN907lPFQqFQqG4EjejpahQKBQKxYIopahQKBQKRR6lFBUKhUKhyKOUokKhUCgUeZRSVCgUCoUiz02vFIUQ24QQ/yGE2FBoWW4GhBD3CyF+Xwjxl4WW5WZACLFXCPG7Qog/KbQsNwNCiLcLIb5caDluBoQQHxNC/IEQ4heu53FvWKUohKgVQrwihEgIIUzz1n9BCHHg3E1ZSnkS+I9CyXmjcBXn+0kp5ecAY8GEvQG4ivN9CBCAp1CyFjvLPddCiK1AFAgWStYbgeWeb2AWTYfZrufn37BKEQgA9wEHz60QQuwCHFLKOwGLEGJPoYS7AVn2+RZCfBz4p4JIeeOw7PMtpfxzYKQgUt4YLPdc3wu0AbuFELUFkfTGYFnnW0r5j1LKPwZ2iOs44fiGHTIspUwAiUvO1W3Ak/nlJ4F9Qogp4EFgvRDi01LK5OpKemNwFed7F7AHkEKIQ1J1j7gmruJ8NwCbgLLVlfDGYbnnWkp5zmKsllKOra6UNw5XcW1XA9uA1PW8j9ywSvEKeIHe/HIQ2CylHAR+qWAS3dh4ufx8fwb4SsEkurHxcvn5/hLw/YJJdOPi5ZJzfW6DlPL3CiHQDY6Xha/tH17vD7qR3acLMQe488vu/GvFyjGHOt+ryRzqfK8Wc6hzvZrMsUrn+2ZTii+j+aoB7meez1qxIqjzvbqo8716qHO9uqza+b5hlaIQwiyEeBLYDvxcCLFXSvkKmq/6AJCTUh4urJQ3Dup8ry7qfK8e6lyvLoU+32pKhkKhUCgUeW5YS1GhUCgUiqtFKUWFQqFQKPIopahQKBQKRR6lFBUKhUKhyKOUokKhUCgUeZRSVCgUCoUij1KKCkWBERr/KYSwCCG+vUqf+SkhxNtX47MUimJCKUWFovC0At1ozY07CiyLQnFTc7M1BFcodIUQ4s+BO4EM2uihkBAiAaxHG0Mkgd8HzgD/DFjyb32nlHJKCHEKeArYm/+dzR9nSEr5HiHEPcD/AFLAGuCPpJQ/uESGjwO/gHY/+FMp5Y+EEJ8GHgDSwJeklN9dmTOgUOgL1dFGoSgwQoh/BD4OfBl4H5qCPAjskVLmhBAG8kOZpZRpIcQHgFIp5f8SQvQDbwQ60azND0sp/0sI8RTwfqAJ+ALauC4n8IKUcosQ4lNoVulp4DNSyrcLIWz57buFEK+gjUNKCSEMUsrcKp0OhaKgKEtRoSgQQoj/v717V80qiMIw/H4IKpgQEGtvICKIWCeVV2AjIlgpSO7AQrSysLYIIghWNtHK1iiISMDGQ2MjErAUPDXCZ7En8JNKJBB/eZ9qz2aYQ7VYszZ7VoAbwAmm4LUMPAGuALeA+0l+AjeZMrY7SY4BS8DzMcyPtu/GeNvA6/H+E3B0PG+1/QV8SfItyZGZZSwDp5I8He2FJItM2eV6kgK3gbd7unnpH2VNUdonbTeB88BD4DJwr+0q8AF43PYisAlcBS4Am21XgHVg5wbW3Uc9s+2dPqeTHEiyBCy0/T7T5z3wsu3qmPtk26/Ai7aXgLvAtb3YrzQPzBSl/XUG2GKqCe78+X8ReDSytIPAGlOt8EGSs8A2U+b4pz4DG8BxdgW4tm+SvErybMzxkekIdyPJIeAwcP3vtibNH2uK0n9sfGhzru3aPi9Fmgsen0qSNJgpSpI0mClKkjQYFCVJGgyKkiQNBkVJkgaDoiRJg0FRkqThN4pTBmKlYRkqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 460.8x316.8 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAAFDCAYAAABLKBVBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAulElEQVR4nO3de3xV93nn+88jCYEEEroggwEDxpirAIEESLLj+7FjJ42dTJu0TpM05zS2OzmTaWfGJ3F7cm182rQ9E7dNJzPOdDqJk5ykTdPEcdpJ48RObOsucb9jBzDGYIEQEkjotp/zx17CG6HLltDeS3vr+369eGnt315rr2cvlvTdv7XWb21zd0RERAQywi5ARERkqlAoioiIBBSKIiIiAYWiiIhIQKEoIiISUCiKiIgEFIoiIiIBhWKaMrMlZvasmR0xs31m9j/MLHeSXvshM1sxSa91h5lVjHOZx8zsc8O0/08zO2tmWTFttWZWF7Ouc2a23cwOmdmPzWxTzLwvmtl+M9tpZj8zs+uHWcd/Crbnl8f1RsfBzI6a2ayYxyO93xfNbPUYr/U7ZvY7CajxVPCzzMzuiWn/gplVj7HsLDM7Gky/x8w+Mdn1xWO0fW+07WZmvxnsAwfMrG6U11g2uO8lmpktNLNvJGNd6U6hmIbMzIAfAt919xXuvjZ4nDdJq3gIGDYUzSxznK91BzCuUBzDm8BdQS03AEM/CPzE3Te5+0rgb4D/ZWYLYp5/r7tvBGqAPxzm9X8XqHL3P4inmAlsj1RTBlwORXf/jLvXxLuwuz/r7n+ViMLicAfj3PeCAPw8cI+7rwYeB75vZpP1uzXaurNGes7dT7r7hxNdw3SgUExP/xtwzt2/Ndjg7j9099NmVmJm/2xmu8zsJ2Y2Hy73Or5kZs1m1hIECmb2+8En4p1m9lTwR+E9wN+Y2Q4zywt6Np8JPhVvie3pBJ+4/zSYXmNmvwheq87MlgKPAX8UvNaq4N/zQR3/bGbFwbIPBr27WmDLKO/9H4BfD6Z/HfjHkWZ0938Gfgz81jBPvwLcFNtgZn8FLAd+aWYfNrMVZvZSsC3/3sxmB/NdsT2GvMbnzazRzPaY2V+M8j7i9WiwPV8xs0XxLmRmnzOz/25mL5vZq2Z2m5l928wOmtnng3mu6OlYtCf+ziEv9QXgo8H/350jzIOZzTezF8xsN/B/x7TH7h9/FrOvPR60Dd1nZoyxD68Opu8ws+/E1P2XZlZv0SMBZWa2kCH7Xpyb7g+AP3H3kwDu/hLwc+CDcS6Pmf1aUMsOM/vrmPbngv1+j5l9IGhbFvw+/h3QHLyvn5jZj8zssJn9Ucx8g0dEPmdmX4v5v703aJ9tZj8ws71m9lULevtyJYVieloD7Bjhuc8T7S1tINp7/GLMcxfcvRz4FtE/GAB/BGwKek+fdfcm4Fng4+5e5u6dwXyn3b3S3Uc7XPQN4HPBa70TOA78V+DJ4LUOAl8BPhrU8R3gCTPLAZ4i+sn+NmDtKOuoB8ot2kN7MKh1NDuBlcO03w/si21w908AJ4Ft7v4N4K+Avwi25VHg92NmH2l7/KW7bwE2ACvMbNsY9Y2lJ9ie/51oQI3HIqLb8z8APyAaVhuIhtycOF/jM8DfBf9/L4wy32eBf3T39cBVf4zNrAj4DWBt8H6eDp4aus/0M/o+PJI57r6N6Ht8PAi1oftePNYCLUPatgOjHsYeZGYlwCeA29y9DJhtZu8Onv5wsN9XA58xsxlB+wai+9nG4HEZ8DvARuCxwQ9jQywFbgfeT3TbA3wcOOju64h+GJwfT83TjUJx+qkG/r9g+pvArTHPDQbIdqK/VBD9A/CN4JNrzyiv+73RVmpm+UDe4B9Od2/3ITfeteghqGrgR2a2g+ihqRuAVcC+4BBRH/D9Ud8h/AL4beAC0D7GvDbk8T+Z2U7gOuD/GWPZMnf/YTA9dFuOtD3uNrMGoh9aKhk+4Ie7IfFINyke7Al/D7hl1Gqv9i/uHgF2A0fd/TV37yEa8AtGXXL8qoHvBtPfHeb5DqL/X0+b2a8BF0bZZ0bbh0cy3L4dhkqiIVcf7OO38PYRid8P9r2XiNa4OGjf7+57Y16jxt3PunsX8CvgqnPfwI/dfYAr3281wX7p7s8BXZP2rtLIiMeoJaXtJ3qIczjGyH9gB0Mvwtv7xgPAnUQPRf4e0d7acC7GTA/w9geumTHtY9193oA3gk/QbzfGXAwTM99ovgc8B/zHMeaD6B+o3TGP3+vuB+JYDq58P0NrujjkMRY9pPyfgS3ufsrM/l+u3D6D2oBCoudHAYqAM3HUMN67+/cGPyMx04OPs4juD7EfnIerdVhm9uu8fZj0Aa7cPlf9/7l7v0UPzd8HfBj4N0R7VMO9p5H24ZH2Oxh+356I/cBmYFdM2ybgZ2a2jGiPG6I96F1czYAfuvsjVzSa3Uk0MLe5+yUzaw7ewyWu3pdiP5wOvp/+4eZx94i9fS5ytN99CainmJ6eB4rN7DcHG8zs/Ra9oOQVYLD9YeDlkV7EzDKAxe7+PNFDbGuCpy4w+kU7x4CyYPkHANy9g+in/9uD1y4wM4t9rWCeczHnQGYG53oOAGsteoXdDOC9Y7z/WqK9vH8abSYzewB4N2/3OsZrZ9Crgeh5yRG3ZWAW0T9iZ82sgOjh3eG8RLSnS3Do+L2jvHbs+dO4L3CJ02lgcXAuai7D98iG3Rfc/XvBYcmy4FDlK8AHgqffP3T+4HDtXHd/Fvgk0V74SPvMSPvwMaKHFiH6/zqWsfbj4XwZ+MPgnCRmdivRw5jfc/ejMe95pMP2dUSPFiwKli+x6FXO+UBbEIhlRD+sTbYaoh82Bvf9SbkaPd2op5iGgk+HDwFfMbMvEv3UWEu09/Q54Otm9jGi53Y+NMpLZQLfDg5rGtHziwB/D/y34CT/O4ZZ7ovA14ETwKGY9g8TPTyWT/TQzW1BTf9gZh8k+kfzt4H/amZ/Hqz/j939oJn9AfAC0d7T0HM6Q9+/A38B0T+kQ56+LzhsNRs4DNzn7hO94OATwN+Z2ZNEg/ujY9TVbmbfBPYSPUQ50vnXLxDdvjuIbven3X3nSC9rZvVEe0m/Me53MHq9vWb2FNFDcPsZvufzItHzvjuIXoQyks8D37XoMIcfDPN8HvCsmWUT7c18Omgfbp/5HMPvw18O1vF7jLGPBK7Y9+I5r+jujRa9EOlnQa1zgXXu3j3CIpvM7ETM43uJntt7NviA10P0/OD/An4v2I77gOY46h+vvyH6+7yd6MVBxxOwjpRn+j5FkfQVhBDu/j/DrSS1xLPdglD8DvCquz+enMomLjiMmunuPUEP94vufkfIZU056imKiEyAu/cC7wu7jnGYC/xrEI49wMdCrmdKUiiKpLcdYReQonaEXcBkc/ezQHnYdUx1OnwqIiIS0NWnIiIigbQ+fDpv3jxftmxZ2GWIiMgU0tzcfMbdS4Z7Lq1DcdmyZTQ1NYVdhoiITCFmdmyk53T4VEREJKBQFBERCSgURUREAgpFERGRgEJRREQkoFAUEREJKBRFREQCCkUREZGAQlFERCSgUBQREQkoFEVEJCV0dXXR3NxMV1dXwtahUBQRkZSgUBQREUkihaKIiEhAoSgiIhJQKIqIiAQUiiIiIgGFooiISEChKCIiElAoioiIBBSKIiIiAYWiiIhIQKEoIiISUCiKiIgEFIoiIjLlRSIR6urqEnozcICshL66iIjIJDhw4ACHDh0iEokkdD0J6Sma2UIzazGzS2aWZWY3mtlLZvZLM/u2mWUG833QzGrM7Dkzyw/a7jKzWjN7wcwWB22lZvaymb1iZhsSUbOIiExNPT09NDU1cdNNN5Gfn5/QdSXq8GkbcDdQFzxuB37N3W8DfgU8YGYzgMeA24BngEeDeT8N3At8CngiaPtj4LeA9wfTIiIyTTQ3N9Pf38+mTZsSvq6EhKK7X3L3czGPz7l7e/CwHxgAVgK73b0feB6oNLNcoNvdO929HlgbLFPk7q+7+xvA3NHWbWaPmFmTmTW1trZO8jsTEZFkOnfuHHv37mXz5s3k5uYmfH1JvdDGzBYC9wD/ChQAHcFT54HC4F9HzCKZwc/YOket2d2fdvcKd68oKSmZjLJFRCQktbW1zJkzh/Xr1ydlfUkLRTObCXwd+FjQO2wHBg8O5wePz8W0AUSG/Bw6LSIiaerYsWOcOHGCqqoqMjMzx15gEiTz6tOngf/i7vuCx4eA0uCim3uAOnfvMrMcM5tD9NDp4LxtwUU3EaK9ShERSWMDAwPU1taycOFCli1blrT1JiQUg4to/gXYCPzEzL4AvA9Yamb/HvhLd/8nM/sa8BLRHuLDweJPAj8FLgEfCdo+C3wHMODjiahZRESmjj179tDZ2cm9996b1PUmJBTdvY9o7y9W3jDzPUP0ytPYtueJXngT27YLuHWSyxQRkSmoq6uLlpYW1qxZQ1FRUVLXrTvaiIjIlNLY2EhGRgYVFRVJX7dCUUREpozW1lYOHjxIeXk5s2bNSvr6FYoiIjJl1NTUUFhYyNq1a8eeOQEUiiIiMiUcOXKE06dPU1VVRUbG1fGUm5tLeXl5Qgfx64bgIiISuv7+furr61m6dCmLFy8edp7BUEwk9RRFRCR0O3fupLu7m8rKylDrUCiKiEioLly4wI4dO1i/fj1z5456e+uEUyiKiEio6uvryc7OTsq3YIxFoSgiIqE5deoUr776Klu3biU7OzvschSKIiISDnenpqaGefPmsXLlyrDLARSKIiISkoMHD3LmzBmqq6sxs7DLARSKIiISgt7eXhoaGlixYgULFiwIu5zLFIoiIpJ0LS0t9Pf3s3Xr1rBLuYJCUUREkqq9vZ09e/ZQVlbGnDlzwi7nCgpFERFJqrq6OnJzc9mwYUPYpVxFoSgiIknz+uuvc/z4cSorK8nKmnp3GlUoiohIUkQiEWpra1mwYAHLly8Pu5xhKRRFRCQp9u7dS3t7O9XV1WGXMiKFooiIJNylS5dobm5m9erVzJs3L+xyRqRQFBGRhGtsbARgy5YtIVcyOoWiiIgk1NmzZzlw4ACbN28mJycn7HJGpVAUEZGEqq2tJT8/n9LS0rBLGZNCUUREEuZXv/oVJ0+epKqqioyMqR85U79CERFJSf39/dTV1XHDDTewZMmSsMuJi0JRREQSYvfu3Vy8eJGqqqqwS4mbQlFERCbdxYsX2b59O+vWraOgoCDscuKmUBQRkUnX0NBAVlYWmzdvDruUcVEoiojIpDp9+jSHDx9my5YtzJw5M+xyxkWhKCIik8bdqampoaioiNWrV4ddzrglJBTNbKGZtZjZJTPLCtoeN7OXzexbZjYjaPugmdWY2XNmlh+03WVmtWb2gpktDtpKg2VfMbOp910jIiICwOHDh2ltbaW6uhozC7uccUtUT7ENuBuoAzCzEuBOd78V2AU8FATjY8BtwDPAo8GynwbuBT4FPBG0/THwW8D7g2kREZli+vr6aGhoYPny5SxcuDDsciYkIaHo7pfc/VxM01bgxWD6eaASWAnsdvf+wTYzywW63b3T3euBtcEyRe7+uru/Acwdbd1m9oiZNZlZU2tr6yS+KxERGc327dvp6elh27ZtYZcyYck6p1gAdATT54HCEdoKY9oAMoOfsXWOWrO7P+3uFe5eUVJScm1Vi4hIXDo6Oti1axcbN24kLy8v7HImLFmh2A7kB9P5wePh2s7FtAFEhvwcOi0iIlNAXV0dOTk5lJWVhV3KNUlWKDYCtwfT9xA913gIKDWzzME2d+8CcsxsjpltBfYFy7SZ2WIzW0i0VykiIlPEG2+8wdGjR9m2bRtZWVlhl3NNElJ9cBHNvwAbgZ8Afwj80sxeBo4DT7l7n5l9DXiJaA/x4WDxJ4GfApeAjwRtnwW+Axjw8UTULCIi4xeJRKipqWH+/PmsWLEi7HKumbl72DUkTEVFhTc1NYVdhohI2tq7dy+vvPIK733ve0mV6zjMrNndK4Z7ToP3RURkQnp6emhqamLlypUpE4hjUSiKiMiENDc3E4lE2Lp1a9ilTBqFooiIjNu5c+fYu3cvmzZtIjc3N+xyJo1CUURExq22tpa8vDzWr18fdimTSqEoIiLjcuzYMU6cOEFlZSWZmZljL5BCFIoiIhK3gYEBamtrWbRoEcuWLQu7nEmnUBQRkbjt2bOHzs5Oqqqqwi4lIRSKIiISl66uLlpaWli7di1FRUVhl5MQCkUREYlLY2MjGRkZlJeXh11KwigURURkTK2trRw8eJCKigpmzZoVdjkJo1AUEZEx1dTUUFhYyJo1a8IuJaEUiiIiMqojR45w+vRpqqurychI79hI73cnIiLXpL+/n/r6epYtW8aiRYvCLifhFIoiIjKiHTt20N3dTWVlZdilJIVCUUREhtXZ2cnOnTvZsGED+fn5YZeTFApFEREZVn19PTNnzmTTpk1hl5I0CkUREbnKm2++yWuvvcbWrVuZMWNG2OUkjUJRRESu4O7U1NRQUlLCzTffHHY5SaVQFBGRKxw8eJCzZ89SXV2NmYVdTlIpFEVE5LLe3l4aGhpYsWIF8+fPD7ucpFMoiojIZS0tLfT397Nt27awSwmFQlFERABob29nz549lJWVMXv27LDLCYVCUUREAKirqyM3N5cNGzaEXUpoFIoiIsLx48c5fvw4lZWVZGVlhV1OaBSKIiLTXCQSoa6ujuuvv57ly5eHXU6oFIoiItPc3r17aW9vp7q6OuxSQqdQFBGZxrq7u2lubmbNmjUUFxeHXU7oFIoiItNYU1MTAFu2bAm5kqlBoSgiMk2dPXuW/fv3U15ezqxZs8IuZ0pIWiiaWa6Z/djMXjSzH5rZTDN73MxeNrNvmdmMYL4PmlmNmT1nZvlB211mVmtmL5jZ4mTVLCKSzmpqaigoKGDdunVhlzJlJLOn+E6g3t3vABqA3wTudPdbgV3AQ0EwPgbcBjwDPBos+2ngXuBTwBNJrFlEJC299tprvPnmm1RVVZGRoYOGg5K5JV4FZgbTBcAy4MXg8fNAJbAS2O3u/YNtZpYLdLt7p7vXA2tHW4mZPWJmTWbW1NraOulvQkQk1fX391NXV8eSJUu44YYbwi5nSklmKB4GtpnZXqACOAJ0BM+dBwqJhuXQtsKYNoDM0Vbi7k+7e4W7V5SUlExe9SIiaWLXrl10dXVRWVkZdilTTjJD8SPAT9x9HfBjIAvID57LB9qDf0PbzsW0AUQSX6qISHq6ePEiO3bsoLS0lIKCgrDLmXKSGYoGtAXTZ4gePr09eHwPUAccAkrNLHOwzd27gBwzm2NmW4F9SaxZRCSt1NfXk5WVxebNm8MuZUpK5g3uvg1818w+BPQBHwA+ZmYvA8eBp9y9z8y+BrxEtIf4cLDsk8BPgUtEe5wiIjJOp0+f5siRI9x2221kZ2eHXc6UZO4edg0JU1FR4YMDU0VEpjN35wc/+AGRSIT3ve99mFnYJYXGzJrdvWK453QdrojINHD48GFaW1uprq6e1oE4lrhC0czeF/xca2bfNLN7EluWiIhMlt7eXurr61m+fDnXX3992OVMafH2FP9d8PM/Av8F+HxiyhERkcm2Y8cOent72bZtW9ilTHnxhmK2mWUBGe5eA/QksCYREZkkHR0d7Nq1i40bN5KXlxd2OVNevKH4PNEhE98ws1koFEVEUkJdXR05OTmUlZWFXUpKiCsU3f2zwZU6De5+yd3vT3BdIiJyjU6cOMHRo0fZtm0bWVnJHIGXuuK90OYBM2sEGswsy8z+R4LrEhGRaxCJRKitrWX+/PmsWLEi7HJSRryHT/8IuBV4K7hZ97KEVSQiItds//79nDt3jurq6rBLSSnxhuKAu/cAgyP9R70pt4iIhOfSpUs0NTWxatUq9MUI4xNvKH7fzP4BWGJmzwD/mMCaRETkGjQ3NxOJRNiyZUvYpaScuM68uvtTZvYTYB2w3933JrYsERGZiLa2Nvbt28fWrVvJzc0Nu5yUM2ooDt7JJkYEWGVmq9z9+4krS0REJqKmpoa8vDxKS0vDLiUljdVTXD9CuwMKRRGRKeTo0aOcPHmS++67j8xMXfoxEaOGortfcTs3M8tx9+7EliQiIuM1MDBAXV0dixcvZunSpWGXk7LiHaf4bjPbDfzMzHaZ2bsTXJeIiIzD7t276ezspKqqKuxSUlq8tzj4AnCru583s7nAC8BziStLRETi1dXVxfbt21m3bh2FhYVhl5PS4h2ScRToCKY7gNcTUo2IiIxbQ0MDGRkZlJeXh11Kyou3p5gNNJlZE1AOnDazvwJw908kqjgRERlda2srhw4d4pZbbmHmzJlhl5Py4g3FP4+Z/nYiChERkfF75ZVXKCwsZM2aNWGXkhbiHbz/CzNbChQBFrS1JLIwEREZ3ZEjR3jrrbd417veRUZGvGfDZDRxhaKZfQu4DjgRNDnwvyeqKBERGV1fXx91dXUsW7aMRYsWhV1O2oj38OlSd781oZWIiEjcdu7cyaVLl6isrAy7lLQSb3/752Z2t5kVm1mRmRUltCoRERlRZ2cnO3fuZMOGDeTn54ddTlqJt6eYCXyV6OFTI3r49K5EFSUiIiOrr69n5syZbNq0KexS0k68oXiXu69MaCUiIjKmkydP8tprr3HHHXcwY8aMsMtJO/EePm0ws1UJrUREREbl7tTU1FBSUsLNN98cdjlpKd6e4t3Ar5lZF9FDp+7uGxJXloiIDHXgwAHa2tp48MEHMbOwy0lL8Y5T1BdziYiEqKenh8bGRm6++Wbmz58fdjlpK95xivOIjktcxNuD93V7NxGRJGlpaaG/v5+tW7eGXUpai/ec4t8D7cA7gP3AhM7umtmHzexnZvaimS0ys8fN7GUz+5aZzQjm+aCZ1ZjZc2aWH7TdZWa1ZvaCmS2eyLpFRFJVe3s7e/fuZdOmTcyePTvsctJa3PcFcvengfPu/lWiPcZxMbNFwO3ufre73wH0AncGNwXYBTwUBONjwG3AM8CjweKfBu4FPgU8Md51i4ikstraWmbPns369evDLiXtxRuKfWaWCZw1s/8TmMjXOt8HZAY9xb8GtgIvBs89D1QCK4Hd7t4/2GZmuUC3u3e6ez2wdrSVmNkjZtZkZk2tra0TKFNEZOo4fvw4r7/+OpWVlWRlxXttpExUvKH4cPDz3xI9dPrbE1jXfCDb3e8GuoAC3v6OxvNA4QhthTFtEL2RwIjc/Wl3r3D3ipKSkgmUKSIyNUQiEWpra1m4cCE33nhj2OVMC/GG4kp3HyB6Yc68cSwX6zzwi2D658AyYPD+RPlEz1m2D9N2LqYNIDKBdYuIpJw9e/bQ0dFBVVVV2KVMG/GG258EPz8D7AW+MoF11QCDYxvLgNeB24PH9wB1wCGgNDhUew9Q5+5dQI6ZzTGzrcC+CaxbRCSldHd309LSwpo1ayguLg67nGkj3gPUg+GZ6+7fNrPfHe+K3H2HmXWb2YvAGaKHZK83s5eB48BT7t5nZl8DXiLaQxw8bPsk8FPgEvCR8a5bRCTVNDY2AlBRURFyJdNLvKF40MwagT8NrhCd0K0U3P0/DWn6UvAvdp5niF55Gtv2PNELb0RE0t6ZM2c4cOAAVVVVzJo1K+xyppV4D58eBKqIDp0YAB5KVEEiItNdTU0NBQUFrFu3LuxSpp1RQ9HMvmtmnwbuJ3rV539z94i7n09KdSIi08xrr73GqVOnqKqqIiNjItc0yrUYa4v/NvDPwFzg60QvgvkTM3tvwisTEZlm+vv7qaurY8mSJdxwww1hlzMtjRWK/wB8gOi4wt8jeuXp3wK5Ca5LRGTa2bVrF11dXRqCEaJRQ9HdHwK+RnTA/n8G1gG/zwTvfSoiIsO7cOECO3bsoLS0lLlz54ZdzrQ15gFrdz8MvOTuHyXaU/wicDbRhYmITCcNDQ1kZWWxefPmsEuZ1uI6i+vu/1cw+Vl3P+XuP0pgTSIi08qpU6c4cuQIW7duJTs7O+xyprVxXdrk7r9MVCEiItORu1NTU8O8efNYtWpV2OVMe7reV0QkRIcOHeLMmTNUV1djNqH7osgkUiiKiISkt7eXhoYGbrrpJhYsWBB2OYJCUUQkNNu3b6e3t5dt27aFXYoEFIoiIiE4f/48u3fvpqysjDlz5oRdjgQUiiIiIairqyMnJ4eNGzeGXYrEUCiKiCTZiRMnOHbsGJWVlWRlxftlRZIMCkURkSSKRCLU1tayYMECbrrpprDLkSEUiiIiSbRv3z7OnTtHdXV12KXIMBSKIiJJcunSJZqamli1ahXz5s0LuxwZhkJRRCRJmpqacHe2bNkSdikyAoWiiEgStLW1sX//fjZv3kxurr59b6pSKIqIJEFNTQ15eXmUlpaGXYqMQqEoIpJgR48e5eTJk1RVVZGZmRl2OTIKhaKISAINDAxQW1vL4sWLWbp0adjlyBgUiiIiCbR7924uXLhAVVVV2KVIHBSKIiIJ0tXVRUtLC+vWraOwsDDsciQOCkURkQRpaGggKyuL8vLysEuROCkURUQS4K233uLQoUNUVFQwc+bMsMuROCkURUQmmbtTU1NDUVERq1evDrscGQeFoojIJDty5AhvvfUW1dXVZGToz2wq0f+WiMgk6uvro76+nhtvvJGFCxeGXY6MU9JD0cz+g5m9HEw/bmYvm9m3zGxG0PZBM6sxs+fMLD9ou8vMas3sBTNbnOyaRUTitWPHDnp6eti2bVvYpcgEJDUUzWwmsDGYLgHudPdbgV3AQ0EwPgbcBjwDPBos+mngXuBTwBPJrFlEJB5dXV388pe/pKmpiQ0bNpCfnx92STIBye4p/i7w9WB6K/BiMP08UAmsBHa7e/9gm5nlAt3u3unu9cDa5JYsIjK2rq4uvv/97wNQVlYWbjEyYUkLxaAXeLu7/zxoKgA6gunzQOEIbYUxbQCj3jjQzB4xsyYza2ptbZ2c4kVExnDq1CnOnDnD5s2bmTFjRtjlyAQls6f4IeDbMY/bgcHjC/nB4+HazsW0AURGW4m7P+3uFe5eUVJScq01i4iM6cKFC9TX15Ofn8/y5cvDLkeuQVYS17UKKDOzx4B1QAXRQ6h/BtwD1AGHgFIzyxxsc/cuM8sxszlED53uS2LNIiIjunTpEtu3b2fv3r309PSwYsUKzCzssuQaJC0U3f2Tg9Nm9rK7f97MPhlciXoceMrd+8zsa8BLRHuIDweLPAn8FLgEfCRZNYuIDKevr4/du3ezc+dOADZt2sTChQv50Y9+FHJlcq2S2VO8LLjiFHf/EvClIc89Q/TK09i254leeCMiEppIJMK+fftoaWmht7eXtWvXsmnTJnJycjhz5kzY5ckkCCUURURSibtz5MgRmpqa6OzsZOXKlZSXl5OXl3d5ntzcXMrLy8nNzQ2xUrlWCkURkVEcP36choYG2traWLp0Kffddx9FRUVXzTcYipLaFIoiIsM4deoUDQ0NnDp1igULFvCe97yHBQsWhF2WJJhCUUQkRltbG42NjRw7doyioiLe+c53smTJkrDLkiRRKIqIAJ2dnTQ3N3Po0CHy8vK48847NcRiGlIoisi01t3dzfbt29m3bx/Z2dlUV1ezdu1afeXTNKVQFJFpqa+vj127drFr1y4ANm/ezPr163WLtmlOoSgi08rAwAD79++/PNZw3bp1bNq0iVmzZoVdmkwBCkURmRbcncOHD9PU1MTFixe5+eabqaioYM6cOWGXJlOIQlFE0t6xY8dobGykra2NZcuWcf/991NYWBh2WTIFKRRFJG3FjjW8/vrrefDBB5k/f37YZckUplAUkbTT1tZGQ0MDx48fp7i4mPvvv58bbrgh7LIkBSgURSRtdHZ20tTUxOHDh8nPz+euu+7ipptu0lhDiZtCUURSXuxYw5kzZ3LLLbewZs0ajTWUcVMoikjK6u3tvTzWMCMjQ2MN5ZopFEUk5QwMDLBv3z62b99OX18f69ato6ysTGMN5ZopFEUkZQwdazj4vYYaayiTRaEoIinh6NGjNDY2cu7cOY01lIRRKIrIlPbmm2/S0NDA6dOnWbhwIbfffjvXXXdd2GVJmlIoisiUdPbsWRoaGnj99dcpLi7mgQceYPHixWGXJWlOoSgiU0pHRwdNTU0cOXKE/Px87r77bpYvX66xhpIUCkURmRK6u7tpaWlh//79zJo1i1tvvZXVq1drrKEklUJRREI1dKxheXk569evJytLf54k+bTXiUgoBgYG2Lt3Lzt27NBYQ5kyFIoiklSRSITDhw/T3NzMxYsXWbVqFeXl5cyePTvs0kQUiiKSPEePHqWhoYH29nZuvPFGHnjgAQoKCsIuS+QyhaKIJNzJkydpaGjgrbfeYuHChdxxxx0aayhTkkJRRBLmzJkzNDQ0cOLECebNm6exhjLlKRRFZNJ1dHTQ2NjIq6++Sn5+Pvfccw833nijxhrKlJe0UDSzbcCXgQGgyd3/wMweBx4EjgG/4+59ZvZB4ONAG/Cwu3eY2V3Ak8Al4EPufiJZdYtI/Lq6umhpaeHAgQPMmjWLd7zjHaxatUpjDSVlJLOneAy4y90vmdm3zOwdwJ3ufquZfRJ4yMx+ADwG3Ab8G+BR4M+BTwP3AmuBJ4iGpohMEb29vezcuZPdu3eTkZFBRUUFpaWlGmsoKSdpe6y7n4p52A9sAF4MHj8PPAzsA3a7e7+ZPQ88bWa5QLe7dwL1ZvanyapZREbX399/+XsN+/v7KS0tpaysjJkzZ4ZdmsiEJP1jnJltAOYB7UQPpQKcBwqBAqBjSFthTBtA5hiv/wjwCMCSJUsmqWoRiRWJRDh06BDNzc10dXVprKGkjaSGopkVAV8B3g+UA4uCp/KJhmR7MB3bdi6mDSAy2jrc/WngaYCKigqflMJF5LJf/epXNDY20t7ezvLly6moqNBYQ0kbybzQJgv4JvC4u58ys0bg3wJ/BtwD1AGHgFIzyxxsc/cuM8sxszlEzynuS1bNIvK22LGGixYt4s4776SkpCTsskQmVTJ7ir8BbAG+FFyW/QTwSzN7GTgOPBVcffo14CWiPcSHg2WfBH5K9OrTjySxZpFpb+hYw3e9610sWrRo7AVFUpC5p+8RxoqKCm9qagq7DJGUdP78eZqamnj11VeZO3cuW7ZsYfny5WGXJXLNzKzZ3SuGe07XS4vIFbq6umhububAgQPk5ORw2223sXLlSo01lGlBoSgiQHSs4Y4dO9izZw8ZGRls3bqVdevWaayhTCva20Wmuf7+/svfazgwMEBpaSkbN27UWEOZlhSKItPU4FjDpqYmuru7Wb16NeXl5eTm5oZdmkhoFIoi09Brr71GY2Mj58+fZ/ny5WzZsoW5c+eGXZZI6BSKItPIG2+8QUNDA62trSxevJi7776befPmhV2WyJShUBRJc5FIhGPHjvHss8+SmZnJokWLePe7383ChQvDLk1kylEoiqSRnp4ezp49S1tbG2fOnKGtrY22tjY6Ojo4cOAAjz76KOXl5WGXKTJlKRRFUpC709HRwdmzZy//a2tr48KFCwBkZGRQVFREUVERN998M2ZGbm4uS5cuDblykalNoSgyxfX19dHW1nZVAPb39wOQk5NDcXExy5cvp7i4mOLiYgoKCq4YbN/V1UVPT4+uLBUZg0JRZAq5cOHCFeF39uxZOjqi35xmZhQUFFBcXMyyZcsuB2A8QZebm6vDpiJxUCiKhGBgYIBz585d1fvr6ekBIDs7m+LiYpYsWUJRURHFxcUUFhbq7jIiCabfMJEE6+7uvqr3197ezuDN+PPz8ykuLmb9+vWXAzAvLy/kqkWmJ4WiyCSJRCKcP3/+qgDs7u4GICsri6KiIhYsWMC6desuXwiTnZ0dcuUiMkihKDIBvb29V4XfuXPnGBgYAGD27NkUFxezevVqiouLKSoqYu7cuQTfJSoiU5RCUWQU7k5nZ+dVARg79KGwsJDi4mJWrFjBvHnzKCoqYtasWSFXLiIToVAUCfT39w879KGvrw+AWbNmXR76UFRUxLx5864a+iAiqU2hKNPShQsXrgrA8+fPA9GhD3PnzqW4uJilS5eOa+iDiKQ2haKktYGBAdrb2686/Dk49GHGjBkUFxezePFiNm7cePniFw19EJme9JsvaWNw6ENsD7C9vZ1IJAJAXl4excXFlJaWXr74JS8vTxe/iMhlCkVJOZFIhI6Ojss3vB4MwK6uLgAyMzMpKiriuuuuY82aNZcDUEMfRGQsCkWZ0np7e6/4xofBnuDg0Ifc3FyKi4tZuXLl5XN/+fn5uvhFRCZEoShTwuDQh6EB2NnZCUSHPgze9/Omm266HIAa+iAik0mhKJNmYGCAvr4+ent76e3tpa+v7/Lj0X52dHSwf/9+SkpKyM7OZubMmVfd9LqgoIDMzMyw36KIpDmF4jQXiUTiDq+Rfg5OD17QMpIZM2aQnZ191c+CggJWrlxJRUUFixcvZvbs2Ul69yIiV1IopiB3vxxGEwmv2MeD38k3kqysrGHDbM6cOWRnZ19+PFLgDf7MysrSVZ4iMuUpFJOov79/wuE19PFoMjIyrgqr7OxscnJyyM/PHzW8Yn/OmDFDF6yIyLSiUBxDJBK5pvCKXWbwq4JGMlzPKzs7mzlz5owZXrHL6tybiMjEKBTH8Ld/+7ejhllWVtawYZWTkzNmeMU+1h1URETClzJ/ic3sy0AF0OLu/z5Z67399ttHDL4ZM2boPJmISBpJiVA0s83AbHd/h5l91cy2uHtjMta9cuXKZKxGRESmgFS5iqIKeD6Yfh6oDLEWERFJU6kSigVARzB9HigcaUYze8TMmsysqbW1NRm1iYhImkiVUGwH8oPp/ODxsNz9aXevcPeKkpKSJJQmIiLpIlVCsRa4O5i+B6gLsRYREUlTKRGK7t4CXDKzl4CIuzeEXZOIiKSflLj6FCCZwzBERGR6SomeooiISDIoFEVERAIKRRERkYBCUUREJGBjfXNDKjOzVuDYJLzUPODMJLyOjEzbOPG0jRNP2zjxJmMbL3X3YQeyp3UoThYza3L3irDrSGfaxomnbZx42saJl+htrMOnIiIiAYWiiIhIQKEYn6fDLmAa0DZOPG3jxNM2TryEbmOdUxQREQmopygiIhJQKIqIiAQUiiIiIgGFImBmC82sxcwumdmY3xwy3vmnOzMrNbMaM3vJzP7OzGyM+bfFzP/lZNWZyiawD2sbj9ME9uNxzS8T2o8nfRsrFKPaiH6JcbxfXjze+ae7g+5e7e7vCB6PNfD2GHBXMP91ZrY+seWlhfHuk9rG4zfe/Xi888v49+NJ38YKRcDdL7n7udg2M/uMmb1oZj83s2VjzS8jc/e+mIc9wOtjbN9T7n4peNgPDCSp1JQ1gX1Y23icJrAfXzV/EspMaRPYjyd9GysUhxF8al7k7ncAHweeCLei1Gdm7zGzPcB1QAlxbF8z2wDMc/d9SSs0TcS7D2sbj8949+Mh859NZq3pIJ79eLK3sUJxeGuAO8zsReCrQL6ZfSj4tPKpcEtLTe7+rLuXAm8Aaxlj+5pZEfAV4P8Iq+YUN+Y+rG08fuPdj4fM/+6w6k5hY+7Hk72NdZHI8A4C/+ru/w7AzGYE3fRnwi0rNZnZTHfvCR52AAsYZfsGJ9i/CTzu7qfCqDkNjLoPaxuP3wT246Hzdye75jQw1n486dtYPUWiG9rMngc2Aj8BZgGngk8jLwAfHW1+M9uW9KJTyzvN7Bdm9gtgPvDXjLJ9gd8AtgBfCuapSnK9KWe8+zDaxhMx3v146Pz/muR6U84E9uNJ38a6zZuIiEhAPUUREZGAQlFERCSgUBQREQkoFEVERAIKRRERkYBCUUREJKBQFBERCeiONiJpzsx2AvuBG4AT7v6BkEsSmbLUUxRJY2Y2E1gEfNLdbyH6NVE3hVyWyJSlUBRJb6XAD9z9WPA4AnSFWI/IlKZQFElvGwEDMLNy4KK7vxluSSJTl84piqS3MqDLzH4K9AEfC7cckalNoSiS3kqBB9z9UtiFiKQCHT4VSW+zFYgi8dNXR4mIiATUUxQREQkoFEVERAIKRRERkYBCUUREJKBQFBERCSgURUREAgpFERGRgEJRREQk8P8D9lpXvkF35LAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 460.8x316.8 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "env = to_prism(\"a U b\")\n",
    "# env = to_prism(\"FG a\") # For more complicated LTL formulas\n",
    "\n",
    "reward_types = [\"zeta-reach\", \"zeta-acc\", \"zeta-discount\", \"multi-discount\", \"reward-on-acc\"]\n",
    "rl_algos = [\"Q\", \"DQ\", \"SL\"]\n",
    "# Comment out below to run for only multi-discount and Q-learning\n",
    "reward_types = [\"multi-discount\"]\n",
    "rl_algos = [\"Q\"]\n",
    "for rl_algo, reward_type in itertools.product(rl_algos, reward_types):\n",
    "    probs, errs = collect_pac_probs(\n",
    "        env,\n",
    "        # You might need to change these parameters, if the counterexample MDP is large (i.e., more than 4 states).\n",
    "        # Otherwise, for the given number of samples, the algorithm might not converge to the\n",
    "        # desired LTL-PAC probability. This will cause the intercept plot (e.g., right of Figure 3) to have\n",
    "        # infinities, which matplotlib will emit a waring when during plotting.\n",
    "        ep_numbers = default_ep_numbers,\n",
    "        ep_length = default_ep_length,\n",
    "        rl_algo=rl_algo, \n",
    "        reward_type=reward_type, \n",
    "        skip_after_stable=3\n",
    "    )\n",
    "    plot_pac_probs(env, probs, errs, \n",
    "                   rl_algo=rl_algo, \n",
    "                   ep_numbers = default_ep_numbers,\n",
    "                   ep_length = default_ep_length, \n",
    "                   reward_type=reward_type, \n",
    "                   plot_for_paper=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gridworld environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "env = \"gridworld_sadigh14\"\n",
    "\n",
    "reward_types = [\"zeta-reach\", \"zeta-acc\", \"zeta-discount\", \"multi-discount\", \"reward-on-acc\"]\n",
    "rl_algos = [\"Q\", \"DQ\", \"SL\"]\n",
    "# # Comment out below to run for only multi-discount and Q-learning\n",
    "# reward_types = [\"multi-discount\"]\n",
    "# rl_algos = [\"Q\"]\n",
    "\n",
    "ps = np.geomspace(0.9, 0.6, num=5)\n",
    "ep_numbers = np.geomspace(100, 3000, num=21, dtype=int)[1:]\n",
    "ep_length = 30\n",
    "\n",
    "for rl_algo, reward_type in itertools.product(rl_algos, reward_types):\n",
    "    probs, errs = collect_pac_probs(\n",
    "        env,\n",
    "        ps=ps,\n",
    "        ep_numbers = ep_numbers,\n",
    "        ep_length = ep_length,\n",
    "        rl_algo=rl_algo, \n",
    "        reward_type=reward_type, \n",
    "        skip_after_stable=3\n",
    "    )\n",
    "    plot_pac_probs(env, probs, errs, \n",
    "                   rl_algo=rl_algo, \n",
    "                   reward_type=reward_type, \n",
    "                   ep_numbers = ep_numbers,\n",
    "                   ep_length = ep_length, \n",
    "                   plot_for_paper=True, \n",
    "                   title=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}