notebooks/nonstat_scm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-stationary SCM and DAG (NONSTAT.) from figure 3(c) in paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from dcbo.examples.example_setups import setup_nonstat_scm\n",
    "from dcbo.utils.sem_utils.toy_sems import NonStationaryDependentSEM as NonStatSEM\n",
    "from dcbo.utils.sem_utils.sem_estimate import build_sem_hat\n",
    "from dcbo.experimental.experiments import run_methods_replicates\n",
    "from dcbo.utils.sequential_sampling import sequentially_sample_model\n",
    "from dcbo.utils.plotting import plot_outcome\n",
    "\n",
    "from dcbo.methods.dcbo import DCBO\n",
    "\n",
    "from numpy.random import seed\n",
    "seed(seed=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n"
     ]
    },
    {
     "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.42.3 (20191010.1750)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"242pt\" height=\"152pt\"\n",
       " viewBox=\"0.00 0.00 242.00 152.00\" 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 148)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-148 238,-148 238,4 -4,4\"/>\n",
       "<!-- X_0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>X_0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-126\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-122.3\" font-family=\"Times,serif\" font-size=\"14.00\">X_0</text>\n",
       "</g>\n",
       "<!-- Z_0 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Z_0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z_0</text>\n",
       "</g>\n",
       "<!-- X_0&#45;&gt;Z_0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>X_0&#45;&gt;Z_0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-107.79C27,-105.31 27,-102.83 27,-100.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5,-100.14 27,-90.14 23.5,-100.14 30.5,-100.14\"/>\n",
       "</g>\n",
       "<!-- X_1 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>X_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-126\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-122.3\" font-family=\"Times,serif\" font-size=\"14.00\">X_1</text>\n",
       "</g>\n",
       "<!-- X_0&#45;&gt;X_1 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>X_0&#45;&gt;X_1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4,-126C62.39,-126 71.31,-126 79.82,-126\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.92,-129.5 89.92,-126 79.92,-122.5 79.92,-129.5\"/>\n",
       "</g>\n",
       "<!-- Y_0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y_0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y_0</text>\n",
       "</g>\n",
       "<!-- Z_0&#45;&gt;Y_0 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Z_0&#45;&gt;Y_0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-53.79C27,-51.31 27,-48.83 27,-46.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5,-46.14 27,-36.14 23.5,-46.14 30.5,-46.14\"/>\n",
       "</g>\n",
       "<!-- Z_1 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Z_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z_1</text>\n",
       "</g>\n",
       "<!-- Z_0&#45;&gt;Z_1 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>Z_0&#45;&gt;Z_1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4,-72C62.39,-72 71.31,-72 79.82,-72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.92,-75.5 89.92,-72 79.92,-68.5 79.92,-75.5\"/>\n",
       "</g>\n",
       "<!-- Y_1 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Y_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"117\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"117\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y_1</text>\n",
       "</g>\n",
       "<!-- Y_0&#45;&gt;Y_1 -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>Y_0&#45;&gt;Y_1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M54.4,-18C62.39,-18 71.31,-18 79.82,-18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"79.92,-21.5 89.92,-18 79.92,-14.5 79.92,-21.5\"/>\n",
       "</g>\n",
       "<!-- X_1&#45;&gt;Z_1 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>X_1&#45;&gt;Z_1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M117,-107.79C117,-105.31 117,-102.83 117,-100.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"120.5,-100.14 117,-90.14 113.5,-100.14 120.5,-100.14\"/>\n",
       "</g>\n",
       "<!-- X_2 -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>X_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-126\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-122.3\" font-family=\"Times,serif\" font-size=\"14.00\">X_2</text>\n",
       "</g>\n",
       "<!-- X_1&#45;&gt;X_2 -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>X_1&#45;&gt;X_2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M144.4,-126C152.39,-126 161.31,-126 169.82,-126\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"169.92,-129.5 179.92,-126 169.92,-122.5 169.92,-129.5\"/>\n",
       "</g>\n",
       "<!-- Z_1&#45;&gt;Y_1 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Z_1&#45;&gt;Y_1</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M117,-53.79C117,-51.31 117,-48.83 117,-46.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"120.5,-46.14 117,-36.14 113.5,-46.14 120.5,-46.14\"/>\n",
       "</g>\n",
       "<!-- Z_2 -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>Z_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-72\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-68.3\" font-family=\"Times,serif\" font-size=\"14.00\">Z_2</text>\n",
       "</g>\n",
       "<!-- Z_1&#45;&gt;Z_2 -->\n",
       "<g id=\"edge11\" class=\"edge\">\n",
       "<title>Z_1&#45;&gt;Z_2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M144.4,-72C152.39,-72 161.31,-72 169.82,-72\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"169.92,-75.5 179.92,-72 169.92,-68.5 169.92,-75.5\"/>\n",
       "</g>\n",
       "<!-- Y_2 -->\n",
       "<g id=\"node9\" class=\"node\">\n",
       "<title>Y_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Y_2</text>\n",
       "</g>\n",
       "<!-- Y_1&#45;&gt;Y_2 -->\n",
       "<g id=\"edge12\" class=\"edge\">\n",
       "<title>Y_1&#45;&gt;Y_2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M144.4,-18C152.39,-18 161.31,-18 169.82,-18\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"169.92,-21.5 179.92,-18 169.92,-14.5 169.92,-21.5\"/>\n",
       "</g>\n",
       "<!-- X_2&#45;&gt;Z_2 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>X_2&#45;&gt;Z_2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207,-107.79C207,-105.31 207,-102.83 207,-100.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"210.5,-100.14 207,-90.14 203.5,-100.14 210.5,-100.14\"/>\n",
       "</g>\n",
       "<!-- Z_2&#45;&gt;Y_2 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>Z_2&#45;&gt;Y_2</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M207,-53.79C207,-51.31 207,-48.83 207,-46.34\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"210.5,-46.14 207,-36.14 203.5,-46.14 210.5,-46.14\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.sources.Source at 0x7f83d0a6cb50>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = 3 # In this example, as in the paper, we consider three time-slices\n",
    "init_sem, sem, dag_view, G, exploration_sets, intervention_domain, true_objective_values, all_causal_effects  = setup_nonstat_scm(T=T)\n",
    "dag_view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the DAG used for the STAT. example (figure 1 in paper) we need to modify the structure of this DAG to represent the non-stationary evolution of the SEM used in the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OutMultiEdgeView([('X_0', 'Z_0', 0), ('X_0', 'X_1', 0), ('X_0', 'Z_1', 0), ('Z_0', 'Y_0', 0), ('Z_0', 'Z_1', 0), ('Y_0', 'Y_1', 0), ('X_1', 'Z_1', 0), ('X_1', 'X_2', 0), ('Z_1', 'Y_1', 0), ('Z_1', 'Z_2', 0), ('Z_1', 'Y_2', 0), ('Y_1', 'Y_2', 0), ('X_2', 'Z_2', 0), ('Z_2', 'Y_2', 0)])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simple networkx syntax is used to manipulate an object\n",
    "G.add_edge('X_0','Z_1')\n",
    "G.add_edge('Z_1','Y_2')\n",
    "\n",
    "# We print out the edge set just to make sure they have been properly added\n",
    "G.edges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above DAG is the graphical structure we will be working with and is faithful to the one used in figure one in the paper. `dag` is a networkx object i.e."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "networkx.classes.multidigraph.MultiDiGraph"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and is the one we will use for exploring the optimization methods further down. `dag_view` is simply a convenient way to visualise the structure at hand and is just a string that we visulise using `pygraphviz`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "digraph { rankdir=LR; X_0 -> Z_0; Z_0 -> Y_0; X_1 -> Z_1; Z_1 -> Y_1; X_2 -> Z_2; Z_2 -> Y_2;  X_0 -> X_1; Z_0 -> Z_1; Y_0 -> Y_1; X_1 -> X_2; Z_1 -> Z_2; Y_1 -> Y_2;  { rank=same; X_0 Z_0 Y_0 } { rank=same; X_1 Z_1 Y_1 } { rank=same; X_2 Z_2 Y_2 }  }\n"
     ]
    }
   ],
   "source": [
    "print(dag_view)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**OBS. that the `dag_view` will not show our added edges as the object was generated before we modified the topology**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we make use of two types of structural equations models: `init_sem` and `sem`. The latter concerns interactions with the first time-slice in the DBN, which has no incoming edges from the previous time-slices, and is only active at $t=0$. For all other time-slices i.e. when $t>0$ the `sem` model is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other setup parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exploration sets: [('X',), ('Z',), ('X', 'Z')]\n",
      "Intervention domains: {'X': [-4, 1], 'Z': [-3, 3]}\n",
      "True optimal outcome values: ['y^*_0 = 4.317', 'y^*_1 = -8.329', 'y^*_2 = 12.387']\n"
     ]
    }
   ],
   "source": [
    "# Contains the exploration sets we will be investigating\n",
    "print(\"Exploration sets:\", exploration_sets)\n",
    "# The intervention domains for the manipulative variables\n",
    "print(\"Intervention domains:\", intervention_domain) \n",
    "# The true outcome values of Y given an optimal intervention on the three time-slices\n",
    "print(\"True optimal outcome values:\", [r\"y^*_{} = {}\".format(t,val.round(3)) for t,val in enumerate(true_objective_values)])\n",
    "# Number of trials\n",
    "N = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamic causal Bayesian optimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/stationary.py:168: RuntimeWarning:overflow encountered in true_divide\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:52: RuntimeWarning:overflow encountered in square\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:76: RuntimeWarning:invalid value encountered in multiply\n",
      "Time index:   0%|          | 0/3 [00:00<?, ?it/s] /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index:  33%|███▎      | 1/3 [00:05<00:11,  5.83s/it]reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index:  67%|██████▋   | 2/3 [00:14<00:07,  7.76s/it]reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:29<00:00,  9.84s/it]\n"
     ]
    }
   ],
   "source": [
    "D_O = sequentially_sample_model(init_sem,\n",
    "                                sem,\n",
    "                                total_timesteps=T,\n",
    "                                sample_count=10, # How many samples we take per node in each time-slice\n",
    "                                epsilon=None) # If we employ a noise model or not\n",
    "change_points = T*[False]\n",
    "change_points[1] = True # We encode where the time-series changes stationary regime\n",
    "DCBO_input_params = {\n",
    "    \"G\": G,\n",
    "    \"sem\": NonStatSEM,\n",
    "    \"base_target_variable\": \"Y\",\n",
    "    \"make_sem_estimator\": build_sem_hat,\n",
    "    \"exploration_sets\":exploration_sets,\n",
    "    \"observation_samples\": D_O,\n",
    "    \"intervention_domain\": intervention_domain,\n",
    "    \"intervention_samples\":None,\n",
    "    \"number_of_trials\": N,\n",
    "    \"sample_anchor_points\": True,\n",
    "    \"seed_anchor_points\": 1,\n",
    "    \"debug_mode\":False,\n",
    "    \"ground_truth\":all_causal_effects,\n",
    "    \"change_points\": change_points,\n",
    "}\n",
    "\n",
    "dcbo = DCBO(**DCBO_input_params)\n",
    "dcbo.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAGSCAYAAAD9+4bmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3sklEQVR4nO3deXTc9X3v/+d7RqN9tySv8oaxA5jNC2nJZocGSBOycJsF0l7S9sa/X2/gNumW3pubpcntuf21OYG25LSnzUbbQGKytIUkhqSxIaQJ8YKxDbbANrIlWbb2fZfevz9mbAkjjSR7pO93Rq/HOXMsfef7nXmNbM9L32U+H3N3REREphIJOoCIiISbikJERJJSUYiISFIqChERSUpFISIiSakoREQkKRWFiIgkpaIQEZGkVBQiIpKUikJERJJSUUjgzGzUzA6a2Qtm9ryZ/aGZRSbcf5OZPW1mNWb2nJl92czyzezDZtY8Ydtvm1l+Yhszs/9tZi+b2UtmttvMrgnuVaaWmS02s4fN7KSZ7Tezn5vZe+c5w2ozOzLF8rsv8TE/dv7vMPF9z+VklNRQUUgY9Lv7De5+DfA24O3AZyD+hgg8CnzC3Te4+43ALqAose23Jmw7BHwgsfyjwM3A9e6+Hvi/wL+bWe68vao5YmYG/CvwtLuvdffNwAeBFZOsmzXP8QBWA5MWxQzyfAzIn2YdmWdB/CMSmZK7N5nZDmCvmX2W+Bv+Q+7+8wnrfBsg/n4Zl3gDKgDaE4s+AbzF3fsS2zxpZv8JfAj4yjy8lLn0VmDI3f/+/AJ3PwX8LYCZfRi4EygEook9ja8Ca4E+YIe7H0r8fHvc/QuJ7Y4A70w85A+BZ4iXbQPwbnfvN7PNiccCeHKKfH8BXGVmB4GHiP+dTMzzGeCP3P2died9ENgHFAPLgN1m1uLu2xP3/3kiV38ix7lL+aHJpVNRyAWr//T7n52Lx639i3fM6nHd/aSZRYEqYCPxN5upfMDM3ggsBV4CHjOzYqDA3U9etO4+IBMOP10DHJhmnU3Ade7eZmZ/Czzn7u8xs7cC/wTcMM32VwJ3uftHzGwn8F+AfwG+Btzr7k+b2V9Nse2f8uoi+PBFebZNtpG7/42Z/QGw3d1bEosLgF+4+yfN7C+BjwD/Z5rskmI69CTp7lvufgOwBDgM/PG8JzC7G7PHJtzWJW4Tl92dWPehCcseSCy796J1y2f39PalxLmdvRMW/8jd2xJfvxH4ZwB3/wmwKFGmybzi7gcTX+8HVptZKVDq7k8nlv/zLGJOzDMbQ8DjE3NcwmPIZdIehVww29/854qZrQVGgSbgBWAz8G/JtnF3N7PHgPvc/S/MrNfM1l60V7EZeCrlgd0fBh6e5J47Jln3nkmWPQg8OItnfIH4b/iJzf2jZlZBfI/pvN4ZPM4Ir/5lceL5m8EJX48CebPIN5mJeZI978WGfXzSnFH0nhUI7VFIqJhZJfD3wIOJN4gHgXvM7PUT1rkzcZL7Ym8ETiS+/ivgb8wsL7HNryXun+wNPd38BMg1s9+bsCzZCeCfEj83Q+KwT4u7dwG1xA8JYWabgDXJntTdO4COxKE+zj/mJLoZv9hgMqeAq80sJ7GXcssstpUAqJ0lDPISJz5jxH/b/GfgiwDufs7MPgh8wcyqgDHgaeJXPsH4OYoIUA98OLH8b4Ey4LCZjQJnSZyQnZdXNIcSe0/vAe43sz8Bmon/xv6JKTb5LPBVMztE/GT2+b2a7wD/1cxeAJ4lfo5nOr+deCxn6pPZh4BRM3se+DrjFxicz1+XOO9xBHgFeG7C3f8A7DKzM+dPZkvwTFOhiohIMjr0JCIiSakoREQkKRWFiIgkpaIQEZGkVBQiIpKUikJERJJSUYiISFIqChERSUpFISIiSakoREQkKRWFiIgkpaIQEZGkVBQiIpKUikJERJJSUYiISFIqChERSUpFISIiSakoREQkqYyZM7u0tNTXrVsXdIxp9fb2UlBQEHSMpNIhIyhnqilnaqVLzv3797e4e2WydTKmKBYvXsy+ffuCjjGtPXv2sG3btqBjJJUOGUE5U005UytdcprZqenW0aEnERFJKmOKYtSDTiAikpkypih6hpz69r6gY4iIZJyMKQoHvnuggZqz3UFHERHJKBlTFDlRGB1zfnC4kX21bbjrWJSISCpkTFHkZRlvXl8BwE9fbmF3TRNjYyoLEZHLlTFFAbB5VTnvuG4p0YjxfF0njx06w9DIWNCxRETSWkYVBcD6xUXcuWk5ubEoJ5t7+c6BenoHR4KOJSKStjKuKABWlOXzga3VFOfFONs5wLf21tHeOxR0LBGRtJSRRQFQXpDNB7dWs7g4l87+Yb65t46Gjv6gY4mIpJ2MLQqAgpwsfmPzCtZWFjAwPMp399fz8jldPisiMhsZXRQA2VkR7rhuGdetKGFkzPn+4Ub2n2oPOpaISNrI+KIAiESMt76uijdeWYE7PP1SM3t0+ayIyIwsiKIAMDO2ri7n7dcuIRoxnjvdwfcPNzI8qstnRUSSWTBFcd7rlhTz3huXkxOLcLyph+/sr6d/aDToWCIiobXgigKgujyf92+ppig3i8bOAb659zQdfbp8VkRkMguyKAAqCnP44E0rqSrOoaMvfvlsY6cunxURudiCLQqAwsTls6sr8ukfGuU7++s53tQTdCwRkVBZ0EUBkJMV5V3XL2fj8hKGR53HD53hudO6fFZE5LwFXxQA0Yjxa1dVcfMVi3CHPTXNPP1Ss4YqFxFBRXGBmfH6tYu47ZolRMzYf6qdHxw+y4gunxWRBU5FcZGrl8Uvn83OivDSuW6+e6CBgWFdPisiC5eKYhIrF41fPtvQ0c+39tbR2TccdCwRkUCoKKZQWZTDB7ZWU1GUQ1vvEN/ce5pzXQNBxxIRmXcqiiSKcmO8b/MKVpbn0zc0yqP76jjZrMtnRWRhUVFMIzcW5T03LufqZcUMjzr//vwZDtV3BB1LRGTeqChmIBoxbr16Ma9fW447/MfRJp55uUWXz4rIghDaojCzUjP7tpkdM7OjZvarAefh5isqeNvVi4mYsbe2jV1HdPmsiGS+rKADJPHXwC53/w0zywbygw4EsHF5CYU5WXz/cCPHznbTMzjCHdcvIzcWDTqaiMicCOUehZmVAG8GvgLg7kPu3hFoqAlWVxTwvi0rKMzJor69n5376ujs1+WzIpKZQlkUwBqgGfiamT1nZl82s4KgQ01UVZTLB26qpqIwm9aeIb619zRNunxWRDKQhfGErJltAX4BvMHdnzWzvwa63P1TF623A9gBUFlZuXnnzp3znnVo1Hm2cYSWficrAluXZLGkYOr+7enpobCwcB4Tzl46ZATlTDXlTK10ybl9+/b97r4l2TphLYolwC/cfXXi+zcBf+ru75hqmw0bNnhNTc08JXy1kdExfnz0HEcbu4mYcctVVWxcXjLpunv27GHbtm3zG3CW0iEjKGeqKWdqpUtOM5u2KEJ56MndzwJ1ZrYhsegW4MUAIyWVFY1w2zVLuGlNOWPu/OjFc/zncV0+KyKZIcxXPd0HfCNxxdNJ4LcDzpOUmfGGdRUU58b4ybEmnn2lja6BEd529WKiEQs6nojIJQttUbj7QSDp7lAYXbuihMLcLL5/6AxHG7voHRzhHdct1eWzIpK2QnnoKd2tqSjgfVuqKciJcrqtj0f31dE1oMtnRSQ9qSjmyOLiXD6wZSXlBdm09Ayxc28dTd26fFZE0k9oDz1lgpL8GB/YWs2/P3+GhvZ+Ht1XT2fjCB3PNQQdLamaxhGKT7eztqKA0vzsoOOISMBUFHMsNxblzhuX8+SL56g528253jHyWnqDjpXUmZ4xnqpp5qmaZsoLslldUcDaigKWlebpxLzIAqSimAdZ0Qhv37iEG6pLeXqonptuWBZ0pKSe6a+jckkRta19tPUO0dY7xIFT7eTEIqwqL2BNRfyWl60T9CILgYpinpgZy0rzWFIQYW1luD+tebo4yrZrlzI25pzp7OeVll5eaemltWeIl85189K5bsxgaUkuayoKWVNRQEVhNmba2xDJRCoKmVIkYqwoy2dFWT5vurKSzr5hTrb08EpLL/Xt/ZzpGOBMxwA/O95CUW4WaysLWFNRSHVZHllRXSchkilUFDJjJfkxblxZxo0ryxgcGaWurY+Tzb3UtvbSPTDC83WdPF/XSSxqVJfns7aikNUV+RTlxoKOLiKXQUUhlyQnK8q6qiLWVRXh7pzrGrywt9HUNcjJ5l5ONsdP2lcV57CmooC1FYUsLs7RISqRNKOikMtmZiwpyWVJSS43X1FB98AwtS19nGzpoa6tj6auQZq6Bnn2ZBv52dF4aVQWUF2eT06WToiLhJ2KQlKuKDfGtStKuHZFCcOjY9S391Pb0suJ5h66B0Z44UwXL5zpIhoxVpTlXbiKSp/ZEAmnyy4KMytw914zK3T3nlSEkswRi0YuFMG2DZW09g7Fr6Jq7uVMZz+nWvs41drHnppmFhVms6aigNWLClhemkdEn9kQCYVU7FGUmdlvA8eBXSl4PMlQZkZFYQ4VhTlsXV1O/9Aota3xS29rW+OX37b2DLGvNv6ZjdWLCi4Uhz6zIRKcVBTFLcCHga+aWZW7N52/w8zuBt4FjAIGPObuj6TgOSUD5GVHuWppMVctLWZ0zDnTMf6ZjbbeIWrOdlNzNv6ZjWUleaypjBfHogJ9ZkNkPqWiKH4J/A5QPbEkEt7i7h88/42ZfQlQUchrRCPxS2qry/N58/pK2nuHeKU1foiqvr2fho747ZmXWyjOi7G2ooAzPWOcbA7/0c7uIU1gJentsovC3Y8mvjw0yd05ZvYOoA5YAeRd7vPJwlBWkE1ZQTabVpYxMJz4zEZLL7UtvXT1D3OwroPaxhGaDp4JOuq0Tp0eZlV9B9etKA06isglmeurnv47cCdwLfGyuHeOn08yUG4sypWLi7hycRFjY8657gFeae5loPk0ayoKgo6X1MiYU1sL/3G0iYHhMbauLtNhM0k7KS0KM/sR8Efu/jyAu/cB/5LK55CFLRIxlpbksbQkj6H6LLbduDzoSNPqrq+h0+Bnx1sYGB7lTVdWqCwkraR6QJ5PAA+Y2dfMbGmKH1skLa0pifL2jUuJRoz9p9r50YvnGBvTeQtJHyktCnc/4O7bgceBXWb2GTPTeQlZ8DYsKeJd1y8jFjVeONPF9w83MjI6FnQskRlJ+RCfFt+nrgH+DrgPeNnMfivVzyOSblZXFHDnphXkxCIcb+rhXw+eYXBkNOhYItNKaVGY2c+ABuB+YDnxz1dsA24ys39I5XOJpKNlpXm8b3M1BTlR6tr6+O6BBvqHVBYSbqm+6mkH8KK7X3wA9j4zOzrZBiILTWVRDh/YspLvHKjnbOcAj+6v4703Ltdw7BJaqT5H8cIkJXHeO1L5XCLprCQ/xvu3VlNRmE1rzxA799XT3jsUdCyRSc3bNGTufnK+nkskHRTmZPG+LdUsK82lq3+YnfvqaOoaCDqWyGtovkqRAOXGorz3xhWsWpRP39Aoj+6vp769L+hYIq+iohAJWHZWhHffsJz1i4sYGhnjewca0mIMK1k4VBQiIRCNGG/fuIRrl5cwMuY89nwjRxu7go4lAqgoREIjEjFuuaqKm9aUM+bOriNnOVjXEXQsERWFSJiYGW9YV8Gb11cAsPtYE7842crUFxOKzD0VhUgIbV5VztuuXowZ/PxEK3tealZZSGBUFCIhtXF5Ce+8Lj6Y4MHTHTzxwjlGNZigBEBFIRJi66qKeM8Ny8nOinC0sYvHD51hWIMJyjwLbVGYWa2ZHTazg2a2L+g8IkFZuSifOzctJzcW5WRzL//6XIMGE5R5FdqiSNju7je4+5agg4gEaWlJHu/bsoLCnCzq2/v59v56+oZGgo4lC0TYi0JEEioKc3j/1mpK82M0dQ3y6L56ugaGg44lC0CYi8KBJ81sv5ntCDqMSBiU5MV4/5ZqKotyaOsdYufeOto0mKDMMQvrJXdmttzdG8ysCvgRcJ+7P33ROjuID21OZWXl5p07dwaQdHZ6enooLCwMOkZS6ZARFnbOoVHnF40jtPY7OVH41WVZlOVe3u99C/nnORfSJef27dv3T3d4P7RFMZGZfRbocfcvTLXOhg0bvKamZv5CXaI9e/awbdu2oGMklQ4ZQTmHR8f4/qFGXmnpJTsrwruuX0Z1ef4lP95C/3mmWrrkNLNpiyKUh57MrMDMis5/DdwKHAk2lUi4xKIR7rh+Ga9bEh9M8F+fa+CEBhOUORDKogAWA8+Y2fPAL4Hvu/uugDOJhE40Yty+cQk3VJcyMuY8/nwjL57RYIKSWqmeCjUlEpMcXR90DpF0YGZs21BJTizCsyfbeOKFswyMjLJpZVnQ0SRDhHWPQkRmwcy4+YoK3rKhEoCnapr5zxMtGh9KUkJFIZJBNq0s47ZrlhAx49mTbeyuaVJZyGVTUYhkmKuXFfPO65eSFTGer+tk15GzGkxQLouKQiQDXVFZyHtujA8meOxsN489r8EE5dKpKEQyVHV5Pr+xeQV52VFeaenlewcaGBjWYIIyeyoKkQy2uDiX92+ppig3i4aO+GCCvYMaTFBmR0UhkuHKC7J5/9ZqyvJjNHcPsnNfHZ39GkxQZk5FIbIAFOfGeP/WaqqKc+joG2bn3jpaewaDjiVpQkUhskDkZ2fxG5tXsKIsj57BEXbuq+ds50DQsSQNqChEFpCcrCjvuXE5aysLGBge5TsH6jnd2hd0LAk5FYXIAhOLRnjndcu4amlxfDDBgw0cb+oOOpaEWCjHehKRuRWNGLdds5icWISDpzt4/FAj0fYR8k62Bh1tWqe6RmnuHmRRQTaRiAUdZ0FQUYgsUGbGtvWV5MWi/PxEK8c7xhg5Ef6iqD03StsvTpEVMSqLcqgqzqGqKJeq4hwWFeQQVXmknIpCZAEzM35l7SKqinLY1VXH9WvKg46UlDsMt5ymJC9GZ/8wjZ0DNHYOAJ0AZEWMiqIcqopyWFycS1VRDosKVR6XS0UhIqytLOTqRVHesK4i6CjTGmnIYtsb1zAwPEpT1yDnugdo6hqkqXuAjr5hznYOJK7mipdHNGJUFOawOLHnsbhY5TFbKgoRSUu5sSgrF+WzctH49K8Dw/HzF+e6BmhK/NnRN8y5rgHOdb22PKoSh64WF+eyqCCbrKiu75mMikJEMkZuLEp1ef6r5g4/Xx5NiT2Pc10DtE8sj4b4etGIsagw+8JeR1VRLhWFKg9QUYhIhpusPAZHRhOHqwZpSux9tPcNxZd1DXIkUR4RO18eiXMexTlUFOYQW2DloaIQkQUnJ2vy8ojveYyXR1vvEM3dgzR3D/JCYi7yiBnlhdksLsqhqji+95Hp5aGiEBEhXh4ryvJZUTZeHkMjYzT3xIvjXNcgzd0DtPYO0dI9SMvF5VEQoypxpVVVcS5dQ54x42mpKEREppCdFWF5aR7LS/MuLBsaGaOlZ/yEeVP3IG09Q7Qkbi8m1qs9NcwJOxVM8BRTUYiIzEJ2VoRlpXksm1Aew6PnyyO+99HcM0hrdvz8RiZQUYiIXKZYNMLSkjyWloyXx57+V9j2q6uDCzVD98xgnYwpilhXF9xxx/iC+++P//nxj48vu+suuPtuuOceaGuLL7viCnjgAXjwQXjiifF1H3oIjh+Hz39+fNlHPwq33/7q59m6FT79afjc52Dv3vHljz0Gu3bBl740vuxTnyLr4py33Qb33gsf+xicOBFfVl4ef/6HH4ZHHpn317SqogK2bZvxa2Lduvjzz/Nr2tjaCo8/Pid/T6l8TRtbW+PbBvxvb7rXtLG1Fa68MtB/ezN5TRtbW2HRokD/7c3kNV3IGeC/vRm9pplw94y4rV+/3tPB7t27g44wrXTI6K6cqaacqZUuOYF9Ps37a+ZezyUiIimhohARkaRUFCIiklTGnMw+d+4cW7ZsCTrGtHp7eykoKAg6RlLpkBGUM9WUM7XSJSewaboVMqYoFi9ezL59+4KOMa09e/awbdu2oGMklQ4ZQTlTTTlTK11ymln/dOvo0JOIiCSlohARkaQypigGRoNOICKSlpqnWyFjimJwxBkeHQs6hohIummZboWMKQoHjjZ2BR1DRCTjZExRAOw/1c7YmAcdQ0Qko4SuKMzsr8zsmJkdMrPvmVnpTLaLGHT0DXOiuWeOE4qILCyhKwrgR8BGd78OeAn4nzPZKCdqAOw71U58nCsREUmF0BWFuz/p7iOJb38BrJjJdtlRyMuOcrZzgPr2aT8/IiIiMxS6orjI7wA/nMmKBtxQXQrEz1WIiEhqWBCHaczsx8CSSe76pLv/W2KdTwJbgDt9ipBmtgPYAVBZWbn5nx/5Fk/UDjM6Bm9dGaMkx+boFVy6np4eCgsLg46RVDpkBOVMNeVMrXTJuX379v3unnygvOkmrAjiBnwY+DmQP9Ntzk9c9JNj5/yLT9b4Dw83XuZ0HnMjHSYzSYeM7sqZasqZWumSk3ScuMjMbgf+BHiXu/fNdvtNK8uImFFztpuugeHUBxQRWWBCVxTAg0AR8CMzO2hmfz+bjUvyYqxfXMiYO8+d7piTgCIiC0nohhl393WX+xibV5Vx7Gw3Rxo6ef2acnJj0VREExFZkMK4R3HZqopzWbUon6GRMQ7VdwYdR0QkrWVkUQBsWVUOwHOn2xnRYIEiIpcsY4uiujyPquIc+oZGeVGDBYqIXLKMLQozu7BXocECRUQuXcYWBcCVVYWU5MU0WKCIyGXI6KKIRIxNq8oADRYoInKpMrooAK5ZVqzBAkVELkPGF0UsGtFggSIilyHjiwLg+hWlxKLGKy29NHcPBh1HRCStLIiiyMuOcs3yEkB7FSIis7UgigI0WKCIyKVaMEWhwQJFRC7NgikKiA8WCHCkoZOB4dGA04iIpIcFVRQaLFBEZPYWVFGABgsUEZmtBVcUEwcLPNrYHXQcEZHQW3BFMXGwwH2n2jRYoIjINBZcUYAGCxQRmY0FWRQaLFBEZOYWZFGABgsUEZmp0BaFmf2hmbmZVczF48eiEa5fUQpoWA8RkWRCWRRmVg3cCpyey+e5oVqDBYqITCeURQHcD/wJMKcnDzRYoIjI9EJXFGb2bqDB3Z+fj+fTYIEiIslZEFf8mNmPgSWT3PVJ4H8Bt7p7p5nVAlvcvWWKx9kB7ACorKzcvHPnzkvKs/fsCPXdY6wrjXBtZdYlPcZM9fT0UFhYOKfPcbnSISMoZ6opZ2qlS87t27fvd/ctSVdy99DcgGuBJqA2cRshfp5iyXTbrl+/3i/Vuc5+/+KTNf7gT172/qGRS36cmdi9e/ecPn4qpENGd+VMNeVMrXTJCezzad5fQ3Xoyd0Pu3uVu69299VAPbDJ3c/O5fNqsEARkamFqiiCpMECRUQmF+qiSOxZTHp+ItU0WKCIyORCXRTzaeJggfs1WKCIyAUqignODxbYrsECRUQuUFFMoMECRUReS0VxEQ0WKCLyatMWhZkVJP4M/ydHUkCDBYqIvNpM9ijKzOxe4I1zHSYsNFigiMi4mRTFLcCHgbVmVpVsRTP7RCpCBU2DBYqIjJvJwEa/BH4HqHb3pol3mNnEwZUMuAH4/1KWLkCbVpZxqK6TmrPd3LxuEcW5saAjiYgEYtqicPejiS8PTXJ3l7v/t/PfmNnfpSpY0EryYqxfXMixs908d7qDt6yvDDqSiEggLveqpz+/6PtPXubjhcrmxKWyRxo6GRgeDTiNiEgwZlUUZvYjM7v+/Pfu/srE+929LVXBwkCDBYqIzH6P4hPAA2b2NTNbOheBwkaDBYrIQjeronD3A+6+HXgc2GVmnzGzvLmJFg4aLFBEFrpZn6MwMwNqgL8D7gNeNrPfSnWwsNBggSKy0M32HMXPgAbgfmA58c9XbANuMrN/SHW4sJg4WODJFg0WKCILy2wniN4BvOivHS3vPjM7OtkGmeD8YIG7jzWxt7adKyoLie9YiYhkvtmeo3hhkpI47x0pyBNaGixQRBaqlI0e6+4nU/VYYaTBAkVkodIw47OgwQJFZCFSUcxCXnaUa5ZpsEARWVhCWRRmdp+ZHTOzF8zsL4POM9GmVWVEzKg5203XwHDQcURE5lzoisLMtgPvBq5392uALwQc6VXODxY45s5zpzuCjiMiMudCVxTA7wF/4e6DABcPbR4GGixQRBaSMBbFeuBNZvasmT1lZluDDnQxDRYoIguJTf2xiDl8UrMfA0smueuTxIcu3w38D2Ar8C1g7WSf3zCzHcQ/BEhlZeXmnTt3XrzKnGnqG+NnDSPkZMFtq2JEIzP7AF5PTw+FheGefjwdMoJypppypla65Ny+fft+d9+SdCV3D9UN2AVsn/D9CaByuu3Wr1/v82lsbMz/5Re1/sUna/xQXceMt9u9e/fchUqRdMjorpypppyplS45gX0+zftrGA89/SuwHcDM1gPZQEuQgSajwQJFZKEIY1F8FVhrZkeAbwL3JFovdDRYoIgsBKErCncfcvffdPeN7r7J3X8SdKapnB8sEGBvbTsh7TMRkcsSuqJINxosUEQynYriMmmwQBHJdCqKFNBggSKSyVQUKaDBAkUkk6koUkSDBYpIplJRpIgGCxSRTKWiSCENFigimUhFkUIaLFBEMpGKIsXOD+vx3Ol2RkbHAk4jInL5VBQpVl2eR1VxDn1Doxxt7A46jojIZVNRpJgGCxSRTKOimAMaLFBEMomKYg5osEARySQqijkycbDAhg4NFigi6UtFMUcmDha4r1bDeohI+lJRzCENFigimUBFMYc0WKCIZAIVxRzbtFKDBYpIelNRzLGSfA0WKCLpLSvoAAvB5lVlHDvbzZGGTtaNeeg/hDemy3lFZILQFYWZ3QD8PZALjAD/3d1/GWioy3R+sMBTrX38oHaYF8deDjpSUqdODVOfW0tlYQ5VxTlUFuZSVZxDbiwadDQRCUDoigL4S+DP3P2HZvbrie+3BRvp8v3qFYs41zVIxCBiFnScpNyhtWeI1p4hjp0dH6+qOC9GZVEOVUU5F/4szMnCQv56ROTyhLEoHChOfF0CnAkwS8osLcnj97ZdwR7q2LbtyqDjJPUfP6ln45aVNHUN0tQ9QHP3IC09g3T1D9PVP8yJpvFhSfKyo1QV5VBVlHuhPErzYyoPkQxiYRtewsyuAp4AjPjJ9pvd/dQU6+4AdgBUVlZu3rlz57zlvFQ9PT0UFhYGHSOpyTKOudMzDJ2DTufgGB2DTuegMzTJ/ExZESjJMUpzjJLErTjbUr4nlQ4/S1DOVFPO1Nq+fft+d9+SbJ1AisLMfgwsmeSuTwK3AE+5+3fM7P3ADnf/tekec8OGDV5TU5PipKm3Z88etm3bFnSMpGaa0d3pGhihuXuApu5BmhO37oGR16wbjRiLCrMT5z1yqSrKoaIwh+ysS7/wLh1+lqCcqaacqWVm0xZFIIeekr3xm9k/Ab+f+PZR4MvzEkpmzcwoyYtRkhdjXVXRheV9QyM0dQ3S3DMY/7N7gPa+4fihrK5BXjjTldgeyvKzLzrvkUtetk6ai4RJGM9RnAHeAuwB3gqE+xIheY387CxWV2SxuqLgwrLBkVFaeoZo6oqf82jqHqS1Z4i23vitZsJJ86LcrAulUVkUv/KqSCfNRQITxqL4CPDXZpYFDJA4ByHpLScryvLSPJaX5l1YNjI6Rlvv0IXDVudPnHcPjNA9MMLJ5t4L6+ZlR8cv102UiD7vITI/QlcU7v4MsHm228W6uuCOO8YX3H9//M+Pf3x82V13wd13wz33QFtbfNkVV8ADD8CDD8ITT4yv+9BDcPw4fP7z48s++lG4/fZXP8/WrfDpT8PnPgd7944vf+wx2LULvvSl8WWf+hRZF+e87Ta491742MfgxIn4svLy+PM//DA88si8v6ZVFRWwbduMXxPr1sWff5avKQuouv9+qgD+OP6axoCO991N063voPlT/4em/lGas/Lpr6jk9G23c3rXnguP2d/Xx+Ge8vjrfuqp8ee/6SZYdyU8/I3xZcuXw1u2wVN7oKFhfPndH4LjL8MvJ3xU5y1vief93vfGl61bBze9Hn74Q2hP/Jzz8uC9d8LhQ3D48Pi6t98e/3PXLgD6+vp4/rlWuPY6+N53oT8x7HxZObz97fDLZ+N/L+e9972BvKa+vj6ef/RZ7M474dAhOHwIS3Sx/fqvx7/4wQ8wEguvvz5+e/RR6O+PLy9fBHe8E/vPn8PLL13Yng+8H2tthR//OP54ADffDK/bgH31qxciWXU13HorPPkk1NXFl7nDRz4CR49izzxDb28vJ7/2JNx2K1RUwjcm/Exe9zp405viP+eW1viy/Hz40IfgwH7Yf+DVP2d49c9k8ybYtDn+mH198WUVi+J/zz/9KRw7Nr7uhz4ELc3wxJPjy970RnjdVfCP/ziec+XK+P+JJ56A06fH1/3IR+DYUfjpM+PLgnhNMxC6q54ulU5mp07YMro73YMj8b2OCZfsHj52nNWrVwcdb1q1tbXKmULKmVp/cOuGcJ7MFpkNM6M4N0ZxbowrKscvN/zJyCne8pZwfyYF4KmnGtIq58RfHif+Gjnxd0pP3DPV75mTrTtx+asfd/rn48J2zjPewBvfsHrK1xEW6ZLzD2awjopC0lbEjGgk/Ce40y9nuLMWZhul+dlBx5hWuuScCY0eKyIiSakoREQkqYw59HTu3Dm2bEl6PiYUent7KSgomH7FAKVDRlDOVFPO1EqXnMCm6VbImKJYvHgx+/btCzrGtMJ2RdFk0iEjKGeqKWdqpUtOM+ufbh0dehIRkaRUFCIiklTGHHrqGHT+5j/CPSxUxKCoZ5JxuUVEgtM83QoZUxQAoyGfi3oUeLF1jJHRMbKi2pkTkVBomW6FjCmK0hzjf9wS3k+/ujvf3FtH7YjzYmMX160oDTqSiMiMZExRACH/9KuxdXU5ew/Dvtp2Ni4rIRLqvCIicTr+MY+urCqkMGZ09g/zUlP39BuIiISAimIeRSLGlWXxH/ne2nYyZeReEclsKop5Vl0UoTAni5buQV5p6Z1+AxGRgKko5lk0YmxaVQrA3to27VWISOipKAKwcXkJubEoZzoGaOiY9tPzIiKBUlEEICcryvXVJUB8r0JEJMxUFAG5sbqMWNSobemjqWsg6DgiIlNSUQQkLzvKxuXn9yraA04jIjI1FUWANq8qIxoxXm7qpr13KOg4IiKTUlEEqCg3xuuWFOEO+05pr0JEwklFEbAtq8sxg6ONXXQPDAcdR0TkNVQUASsvyGZdVSGjY86B0x1BxxEReQ0VRQjctLocgCMNnfQPab4KEQkXFUUIVBXnsmpRPkMjYxys6wg6jojIq6goQmJrYq/iYF0HQyNjAacRERmnogiJFWV5LC3JZWB4lMMNnUHHERG5QEUREmbG1jXxvYrnTreHflpXEVk4VBQhsraigIrCbLoHRjja2BV0HBERQEURKmbGlsS5in21bYxpr0JEQkBFETIbFhdRnBejvW+YE809QccREVFRhE0kYmxZVQbALzWxkYiEgIoihK5eVkx+dpSmrkFOt/UFHUdEFjgVRQjFohE2nd+reEUTG4lIsFQUIXXt8hKysyLUt/fT2KnpUkUkOCqKkMqNRbmhuhTQXoWIBCvQojCzr5pZk5kdmeS+PzQzN7OKILKFwY0rS8mKGCebe2npGQw6jogsUEHvUXwduP3ihWZWDdwKnJ7vQGGSn511YbrUfbXaqxCRYARaFO7+NDDZO+D9wJ8AC/7a0E2ryoiYUXO2h84+TWwkIvMv6D2K1zCzdwMN7v580FnCoCQvxoYlRYy5s/+09ipEZP5Z0B/oMrPVwOPuvtHM8oHdwK3u3mlmtcAWd2+ZYtsdwA6AysrKzTt37pyn1Jeup6eHwsLCWW3TNej8x+lhohG4dVWM3Cybo3Rxl5IxCMqZWsqZWumSc/v27fvdfUvSldw90BuwGjiS+PpaoAmoTdxGiJ+nWDLd46xfv97Twe7duy9pu3872OBffLLGf/pSc2oDTeJSM8435Uwt5UytdMkJ7PNp3l9DdejJ3Q+7e5W7r3b31UA9sMndzwYcLXBbV8c/gPd8fQcDw5ouVUTmT9CXxz4C/BzYYGb1Zva7QeYJs6UleVSXx6dLPVSviY1EZP5kBfnk7n7XNPevnqcoaWHr6jLq2vp47nQ7N64sJRYN1Q6hiGQovdOkkZXl+SwuzqVvaJQXzmhiIxGZHyqKNGJmF85V7D+l6VJFZH6oKNLMuqpCyguy6eofpuZsd9BxRGQBUFGkGTNjc2II8n2nNLGRiMw9FUUaumppMUW5WbT2DHGiuTfoOCKS4VQUaSgasQsTG+3TdKkiMsdUFGlq47IS8rKjNHYOUN+uiY1EZO7MuCjMrCDxZ/gHL1kAsrMiFyY22qshyEVkDs1mj6LMzO4F3jhXYWR2bqguJTsrwqnWPs51DQQdR0Qy1GyK4hbgw8BaM6uaeIeZvc3M/tHMbkh8vyNlCWVKubEo1yYmNtJehYjMldkUxS+B3wFOuXvTRff9DvDHwG+a2VuBG1ITT6azaVUZ0YhxvKmHtt6hoOOISAaacVG4+1F3P+Tu35/k7m5373D3PyI+henWlCWUpApzsrh6aTHumi5VROZGqq56ulAe7v6nwD+l6HFlBrasLsMMjjZ20zWg6VJFJLUuqSjM7Edmdv3579393ybe7+5/e7nBZOZK87NZvzg+XeqBU+1BxxGRDHOpexSfAB4ws6+Z2dJUBpJLsyUxWOCRhk76hzSxkYikziUVhbsfcPftwOPALjP7jJnlpTaazEZVUS5rKgoYHnWeq9NehYikziWfozAzA2qAvwPuA142s99KVTCZva1rygE4WNfB4Ij2KkQkNS71HMXPgAbgfmA58c9XbANuMrN/SFU4mZ3lpXksL81jcHiMIw2aLlVEUuNSp0LdAbzorx2N7j4zO3qZmeQybF1TTsNzDRw41cH1K0rJ0nSpInKZLvUcxQuTlMR577iMPHKZVi/Kp7Ioh57BEY42amIjEbl8Kf91091PpvoxZebi06XGz1XsO9XGmKZLFZHLpOMSGejKqkJK82N09A3zclNP0HFEJM2pKDJQJGJsWRXfq9iriY1E5DKpKDLUVUuLKMzJorl7kNrWvqDjiEgaU1FkqKxohE2rSgHY+4oGCxSRS6eiyGAbl5eQG4vS0NFPQ4emSxWRS6OiyGA5WVGur05MbKS9ChG5RCqKDHdjdRmxqPFKSy9N3ZouVURmT0WR4fKyo2xMTJe6r1aDBYrI7KkoFoDz06W+dK6bjj5Nlyois6OiWACKc2O8bklRYrpU7VWIyOyoKBaILavLMYMXG7voGRwJOo6IpBEVxQJRXpDNuqpCRsc0XaqIzE6gRWFmXzWzJjM7MmHZ583skJkdNLMnzWxZkBkzyfnBAg83dDIwrImNRGRmgt6j+Dpw+0XL/srdr3P3G4hPtfrp+Q6VqRYX57JqUT5DI2McrOsIOo6IpIlAi8LdnwbaLlrWNeHbAkAj2qXQ+b2Kg3UdDI2MBZxGRNKBBT2yqJmtBh53940Tlv058F+BTmC7uzdPse0O4rPtUVlZuXnnzp1zH/gy9fT0UFhYGNjzuztP14/QNuBcWxllXWn0NesEnXGmlDO1lDO10iXn9u3b97v7lqQruXugN2A1cGSK+/4n8GczeZz169d7Oti9e3fQEfzlc93+xSdr/B+fPuEjo2OvuT8MGWdCOVNLOVMrXXIC+3ya99egz1FM5xvAfwk6RKa5orKARYXZdA+McLSxa/oNRGRBC11RmNmVE759N3AsqCyZymx8YqP9p9o1XaqIJBX05bGPAD8HNphZvZn9LvAXZnbEzA4BtwK/H2TGTLVhSRHFeTHaeoc40azpUkVkallBPrm73zXJ4q/Me5AFKBoxNq8qY/exJvbWtrOuqhAzCzqWiIRQ6A49yfy5Zlkx+dlRznUNcLpN06WKyORUFAtYLBrhxpVlAOzVYIEiMgUVxQJ33YoSsrMi1LX10dip6VJF5LVUFAtcbizK9StKAe1ViMjkVBTCjStLyYoYJ5p6aOkZDDqOiISMikIoyMnimuXFgCY2EpHXUlEIAJtXlRMxo+ZsN73D+gCeiIxTUQgAJXkxNiwpZMyd4x2aq0JExgX6gTsJly2ryzna2M2prjH2n2on7J+/O94xSvHp8B8qU87USoecBpzpGaOjb4iSvFjaf5hVRSEXVBTmsLaygNpaePqlSUd2D5Xa5lG6a5QzVZQztWobRzj3s1piUWNRYQ6LCrJZVJhDZWEOiwqzyc+Opk2BqCjkVba/roq648e4ZmVp0FGmFe2IcK1ypoxypo67098coTAni57BEc52DnC2c+BV6+RlR6lIlEZFQQ4VRdmUF2STk/XaOWKCpqKQVynOjbGxIsq2DVVBR5mWNWYpZwopZ2pFzh5l25vX0j80SmvvIC09Q7T2DNLaM0RzzyD9Q6PUtfVRd9HwOcV5MSoKsy+UyKKCHMoLsolGgtv7UFGIiMyhvOwoK7LzWVGWf2GZu9M9OEJrzxAtPYO09sSLpK13iK7+Ybr6hznZ3Hth/YgZ5QUxFhXmvGovpDgva14OX6koRETmmZlRnBujODfGmoqCC8vHxpz2viFae+MFcn4vpLN/mJaeIVp6hqih+8L62VmRC+c+FhVmTzj/kdq3dhWFiEhIRCKJE9+FOaxfXHRh+dDIGG2J8mjtHaKle5DW3kF6B0dp7Byg8aLzH/kTz39MOISVnXVpn4hQUYiIhFx2VoQlJbksKcl91fL+odHEnkf83Mf5cyF9Q6Ocbut7zfQBJXmxCXse8QKZCRWFiEiaysuOUl2eT3X5q89/dA2MXDjv0dozSEvvEO29Q3T2D9N50fmPmVBRiIhkEDOjJC9GSV6MtZXjy0fPn/9IlEdzYi9kJlQUIiILQDRiVCSumoLx8x+/O4NtNdaTiIgklTF7FLGuLrjjjvEF998f//PjHx9fdtddcPfdcM890NYWX3bFFfDAA/Dgg/DEE+PrPvQQHD8On//8+LKPfhRuv/3Vz7N1K3z60/C5z8HevePLH3sMdu2CL31pfNmnPkXWxTlvuw3uvRc+9jE4cSK+rLw8/vwPPwyPPDLvr2lVRQVs2zbj18S6dfHnn+fXtLG1FR5/fE7+nlL5mja2tsa3Dfjf3nSvaWNrK1x5ZaD/9mbymja2tsKiRYH+25vJa7qQM8B/ezN6TTPh7hlxW79+vaeD3bt3Bx1hWumQ0V05U005UytdcgL7fJr3Vx16EhGRpFQUIiKSlIpCRESSsvghqvRnZt1ATdA5ZqACaAk6xDTSISMoZ6opZ2qlS84N7l6UbIWMueoJqHH3LUGHmI6Z7Qt7znTICMqZasqZWumUc7p1dOhJRESSUlGIiEhSmVQU/xB0gBlKh5zpkBGUM9WUM7UyJmfGnMwWEZG5kUl7FCIiMgfSvijM7HYzqzGz42b2p0HnmYqZfdXMmszsSNBZpmJm1Wa228xeNLMXzOz3g840GTPLNbNfmtnziZx/FnSmqZhZ1MyeM7PHg86SjJnVmtlhMzs4k6tggmBmpWb2bTM7ZmZHzexXg850MTPbkPgZnr91mdnHgs41GTP7eOL/zxEze8TMcqdcN50PPZlZFHgJeBtQD+wF7nL3FwMNNgkzezPQA/yTu28MOs9kzGwpsNTdD5hZEbAfeE/Yfp4Wn02+wN17zCwGPAP8vrv/IuBor2FmfwBsAYrd/Z1B55mKmdUCW9w9tNf9m9lDwE/d/ctmlg3ku3tHwLGmlHh/agBe7+6ngs4zkZktJ/7/5mp37zezncAP3P3rk62f7nsUNwHH3f2kuw8B3wTeHXCmSbn700Bb0DmScfdGdz+Q+LobOAosDzbVayXGMutJfBtL3EL3G4+ZrQDeAXw56CzpzsxKgDcDXwFw96Ewl0TCLcCJsJXEBFlAnpllAfnAmalWTPeiWA7UTfi+nhC+saUjM1sN3Ag8G3CUSSUO6RwEmoAfuXsYcz4A/AkwFnCOmXDgSTPbb2Y7gg4ziTVAM/C1xKG8L5tZQdChpvFB4JFp1wqAuzcAXwBOA41Ap7s/OdX66V4UMgfMrBD4DvAxd+8KOs9k3H3U3W8AVgA3mVmoDueZ2TuBJnffH3SWGXqju28C3g58NHGoNEyygE3A37n7jUAvEOZzktnAu4BHg84yGTMrI370ZQ2wDCgws9+cav10L4oGoHrC9ysSy+QSJY75fwf4hrt/N+g800kcftgN3B5wlIu9AXhX4tj/N4G3mtm/BBtpaonfMHH3JuB7xA/rhkk9UD9hz/HbxIsjrN4OHHD3c0EHmcKvAa+4e7O7DwPfBW6eauV0L4q9wJVmtibR4B8E/j3gTGkrcZL4K8BRd/9i0HmmYmaVZlaa+DqP+MUMxwINdRF3/5/uvsLdVxP/d/kTd5/yN7YgmVlB4uIFEodzbgVCdXWeu58F6sxsQ2LRLUCoLrK4yF2E9LBTwmngV8wsP/H//hbi5yQnldaDArr7iJndCzwBRIGvuvsLAcealJk9AmwDKsysHviMu38l2FSv8Qbgt4DDieP/AP/L3X8QXKRJLQUeSlxVEgF2unuoLz8NucXA9+LvF2QBD7v7rmAjTeo+4BuJXwpPAr8dcJ5JJcr2bcD/E3SWqbj7s2b2beAAMAI8R5JPaKf15bEiIjL30v3Qk4iIzDEVhYiIJKWiEBGRpFQUIiKSlIpCRESSUlGIzIKZLZowMuhZM2uY8H12Yp13TTeSsZl92MwenJ/UIpcnrT9HITLf3L0VuAHAzD4L9Lj7F87fb2ZZ7v7v6IOfkkFUFCKXycy+DgwQH0TxZ2Z2iPiQ3fea2R3A/waygVbgQxcP62Bm7wM+A4wSH5wtbOMsyQKnohBJjRXAze4+amYfnrD8GeBX3N3N7L8RH032Dy/a9tPAbe7ecH5oEpEwUVGIpMaj7j46yfIVwLcSk0JlA69Mss7PgK8nJo8J/UCMsvDoZLZIavROsfxvgQfd/VriY/+8ZrpJd/9/iR+eqgb2m9miOUspcglUFCJzq4Txoe/vmWwFM7vC3Z91908Tn5ynerL1RIKiohCZW58FHjWz/cBU81H/lZkdNrMjwH8Cz89XOJGZ0OixIiKSlPYoREQkKRWFiIgkpaIQEZGkVBQiIpKUikJERJJSUYiISFIqChERSUpFISIiSf3/jowV9SOoT1kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_outcome(T,N,outcomes=[dcbo.optimal_outcome_values_during_trials],labels=['DCBO'],true_objective_values=true_objective_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore all optimization methods with one function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the demo in `stat_scm.ipynb` here we are going to demonstrate the replicate method as used to generate the results in the paper. In all these examples we do not employ any interventional data, just observational."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Experiment count:   0%|          | 0/3 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>DCBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/stationary.py:168: RuntimeWarning:overflow encountered in true_divide\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:52: RuntimeWarning:overflow encountered in square\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:76: RuntimeWarning:invalid value encountered in multiply\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:26<00:00,  9.00s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>CBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:26<00:00,  8.85s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>ABO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:01<00:00,  2.95it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>BO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:00<00:00,  3.26it/s]\n",
      "Experiment count:  33%|███▎      | 1/3 [01:06<02:12, 66.35s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>DCBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/stationary.py:168: RuntimeWarning:overflow encountered in true_divide\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:52: RuntimeWarning:overflow encountered in square\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:76: RuntimeWarning:invalid value encountered in multiply\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:29<00:00,  9.97s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>CBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:28<00:00,  9.38s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>ABO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:00<00:00,  3.45it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>BO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:00<00:00,  3.76it/s]\n",
      "Experiment count:  67%|██████▋   | 2/3 [02:17<01:08, 68.98s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>DCBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/stationary.py:168: RuntimeWarning:overflow encountered in true_divide\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:52: RuntimeWarning:overflow encountered in square\n",
      " /Users/nd/anaconda3/envs/scibo/lib/python3.7/site-packages/GPy/kern/src/rbf.py:76: RuntimeWarning:invalid value encountered in multiply\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:27<00:00,  9.17s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>CBO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      " /Users/nd/cloud/turing/DCBO/src/utils/sem_utils/toy_sems.py:181: RuntimeWarning:invalid value encountered in double_scalars\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:27<00:00,  9.09s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>ABO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:00<00:00,  3.90it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\t>>>BO\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "reconstraining parameters GP_regression.rbf.variance\n",
      "Time index: 100%|██████████| 3/3 [00:00<00:00,  3.92it/s]\n",
      "Experiment count: 100%|██████████| 3/3 [03:25<00:00, 68.53s/it]\n"
     ]
    }
   ],
   "source": [
    "# See method for argument details\n",
    "change_points = T*[False]\n",
    "change_points[1] = True # We encode where the time-series changes stationary regime\n",
    "R = 3\n",
    "results = run_methods_replicates(G=G, \n",
    "                                sem=NonStatSEM, \n",
    "                                make_sem_estimator=build_sem_hat, \n",
    "                                base_target_variable='Y',\n",
    "                                intervention_domain = intervention_domain, \n",
    "                                methods_list = ['DCBO', 'CBO', 'ABO', 'BO'],\n",
    "                                obs_samples = None,\n",
    "                                exploration_sets = exploration_sets,\n",
    "                                total_timesteps = T,\n",
    "                                number_of_trials = N, \n",
    "                                reps = R, # Number of replicates (how many times we run each method)\n",
    "                                n_restart = 1,\n",
    "                                save_data = False,\n",
    "                                n_obs = 10, # The method samples 10 time-series for each replicate\n",
    "                                num_anchor_points = 100,\n",
    "                                sample_anchor_points = True,\n",
    "                                controlled_experiment=True,\n",
    "                                change_points=change_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyse results and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.experimental.analyse_results import get_relevant_results, elaborate\n",
    "from src.utils.plotting import plot_expected_opt_curve_paper\n",
    "from matplotlib.pyplot import rc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since we didn't save the results we cannot use the pickled file so we have to convert results to the correct format\n",
    "data = get_relevant_results(results=results,replicates=R)\n",
    "exp_optimal_outcome_values_during_trials, exp_per_trial_cost = elaborate(number_of_interventions=None, \n",
    "                                                                         n_replicates=R, \n",
    "                                                                         data=data, \n",
    "                                                                         best_objective_values=true_objective_values, \n",
    "                                                                         T=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_params = {\n",
    "    \"linewidth\": 3,\n",
    "    \"linewidth_opt\": 4,\n",
    "    \"alpha\": 0.1,\n",
    "    \"xlim_max\": N,\n",
    "    \"ncols\": 5,\n",
    "    \"loc_legend\": \"lower right\",\n",
    "    \"size_ticks\": 20,\n",
    "    \"size_labels\": 20,\n",
    "    \"xlabel\": r'$\\texttt{cost}(\\mathbf{X}_{s,t}, \\mathbf{x}_{s,t})$',\n",
    "    \"labels\": {'DCBO': 'DCBO', 'CBO': 'CBO', 'ABO': 'ABO', 'BO': 'BO', 'True': r'$\\mathbb{E} \\left [Y_t \\mid \\textrm{do}(\\mathbf{X}_{s,t}^\\star = \\mathbf{x}_{s,t}^\\star) \\right]$'},\n",
    "    \"colors\": {'DCBO': 'blue', 'CBO': 'green', 'ABO': 'orange', 'BO': 'red', 'True': 'black'},\n",
    "    \"line_styles\": {'DCBO': '-', 'CBO': '--', 'ABO': 'dashdot', 'BO': '-', 'True': ':'},\n",
    "    \"width\":10\n",
    "}\n",
    "rc('text', usetex=True)\n",
    "rc('text.latex', preamble=r'\\usepackage{amssymb}')\n",
    "rc('font', family='serif')\n",
    "rc('font', size=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABEEAAAEVCAYAAAAPRvuSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACAwklEQVR4nO3deXxbV53//5ck74vs7GkapekW09jpQhNaK5QtaeOUFqiBOmUZkqFu2BrPMPHMD76Nv99xYeb7rT0DbmEgcSEBCsRhGiilsUNThi0ypRttrLTp3ij7asn7It3fH9eSLVtObMe2JOv9fDxurHvvufeeK8XH0kfnfI7FMAwDEREREREREZEpzhrtCoiIiIiIiIiITAYFQUREREREREQkISgIIiIiIiIiIiIJQUEQEREREREREUkICoKIiIiIiIiISEJQEEREREREREREEkJStCswUnl5edGugohE2YEDB6JdBbVFIglO7ZCIRFsstEMi8SxugiCgX3iRRBZLb/rVFokkpphqh3b9PNpVEJHJ1N4JaSnk3fbpaNdEJO5pOIyIiIiIiEgsy0gDvz/atRCZEhQEERERERERiXXJydGugciUoCCIiIiIiIiIiCQEBUFEREREREREJCEoCCIiIiIywTZs2IDb7T5vOZ/PR0VFxSTUSEREJDHF1ewwIiIiIvHC4/FQW1tLdnY2jY2NrFmz5rzH3HfffZNQMxERkcSlIIiIiIjIBHA4HFRWVgKwe/fu85Z3uVw4HA48Hs9EV01ERCRhaTiMiIiISJT5fD4AFixYEOWaiIiITG0KgoiIiIhEWX19PU6nM9rVEBERmfI0HEZEREQkilwuF6tXrx7VMR/72B1h6yVrSrhrBDlHREREEp2CICIiIiJR4vP58Pl82O32UR33q1/9coJqJCIiMrVpOIyIiIhIlNTX11NUVBTtaoiIiCQMBUFEREREosDtdisPiIiIyCTTcBgRERGRKPB4POzbty9sW2NjIz6fj6qqKpYsWaJeIiIiIuNMQRARERGRKCgqKhoS5KitrWXfvn2Ul5dHqVYiIiJTm4bDiIiIiMSI5ubmaFdBRERkSlNPEBEREZEJ4PP52Lx5My0tLXg8HqqrqyksLGT58uVDcoF4PB62b9/O7t278Xq9VFVVceutt5Kfnx+l2ouIiExNCoKIiIiITAC73R4a1lJZWXnOsg6Hg/Lycg2DERERmWAaDiMiIiIiIiIiCUFBEBERERERERFJCAqCiIiIiIiIiEhCUBBERERERERERBKCgiAiCcbtdrNy5cqoXbuurg6Xy0VDQwMulysq9RCR6ItmWxS0YcMG3G53VOsgItHjfv0tVv59WVSu3fDnp6mrf4qqH/6MdV//Jg1/fjoq9RBJRJodRiTB7Nq1C4fDMenX9Xg8bN68mQcffDC0bcOGDeTk5GgKSJEEFM22qLa2luzsbBobG1mzZs2k10FEYsOuPzbimDt70q/b8Oenyb98IY6L5gDga22jeMPX8Rw9TuknPzLp9RFJNOoJIpJgGhsbWbx48aRft7a2dsiHjfXr11NdXT3pdRGR6ItWW+RwOKisrKS8vJycnJxJv76IxI7GvzWx+IpLJ/263pbWUAAEwJ6Vyca//xTVW38+6XURSUQKgogkGLfbzfLlyyf9uvX19UO+9XU4HBoSI5KgotUWiYgEuV9/i+XXLZnUa3qOHqfioYdxv/5W2Pb8yxeG6iQiE0vDYUQSQDAHh8fjAaChoYGGhgZKS0snpTu6x+PB5/MNuZbdbgfMD0MaEiMy9UW7LRIRcb2wj4Y/P43n6HHAHJrS8OenKf3E7WG9MyaK46I53P2J24cMw/G1tZv7ozA8RyTRKAgikgCcTidOp5Oqqiq8Xi+VlZWTen2fz3fO/V6vd5JqIiLRFO22SETEed0SnNctoeqHP8Pb2kblvXdPeh3K//5TQ7bt+mMj+Vdcij0rc9LrI5JoFAQRSSD79++nsLBwxOU3bNhAS0vLiMtXVlbq21wROS+1RSISbftff4vCawtGXH7Dv32blta2EZevvPfuEfcs8bW2saPhd+ys+eaIzy8iY6cgiEgCcblclJaWjrj8wJlcRETGi9oiEYk219+aRjUTy4Nf/4cJq0vZv32bmq+VTcpwHBFRYlSRhOF2uwGzO/pkC+b+GG5YjGZoEEkc0WyLRESgP/moc5KTokZS9cOfUfrJj8REXUQShXqCiCQIl8sVlnzU5/OFghPD2bBhw6iuUV5eHrELusPhwG634/V6w64ZDIooKapI4ohmWyQiAmZy1PwBU+P6WtvOm4tjw799e1TXKF9313l7dtTVP8XyvhwlIjJ5FAQRSRAulytsDH5dXd15u6OPZxf0wsJC3G532AeTpqYmfRsskmCi3RaJiLhe2BeWD6Su/qnzDo0Z7+Ewrhf2UXDlZWHBmOCMNRoWIzKxNBxGJEFkZ2ezYMECwPwQMtnBh/LycrZs2RK2bfv27WzcuHFS6yEi0RXttkhEJDsrkwV9gQbXC/smvSeG64V9+NrasWdm4Dl6PLRsr39KARCRSaCeICIJYv369WzevJmcnBzsdvukD0FxOBxs3LiR2tpaHA4HHo+HNWvWaCiMSIKJdlvk8/nYvHkzLS0teDweqqurKSwsZPny5QrIiCSI9Xd+lM07HiMnOwt7ZkZYb4yJ5mttY93/+reI+xxzZ0ecPldExpfFMAwj2pUYiby8PA4cOBDtaohIlMRKGxAr9RCRyRcrv/95eXkc2PXzaFdDRKIg79a7YqIdEolnGg4jIiIiIiIiIglBw2FEJOG53W6amppwOByhmSrULV5EREREZOpREEREEprH42Hz5s1hs09s2LCBnJwc5SsREREREZliNBxGRBJabW0ta9asCdu2fv16qquro1QjERERERGZKAqCiEhCq6+vx+FwhG1zOBy4XK4o1UhERERERCaKgiAikrA8Hg8+n29IEMRutwNmrhAREREREZk64iYnSKalF84eG5+TWZPAlha+zd8FgZ6+/SlgSwnf39sBhn+crp8C9mlgsYzP+URkTHw+3zn3e73eiNvz8vLC1r/yla9w7733jlu9RERERERkYsRNEOTbncfg7z84PidLnQk5i8O3tb0NbQfNx5kLIXNB+P7mJug+Mz7Xz1wIly6Ff/gHmDFjfM4pIpPmwIED0a6CiIiIiIiMQdwEQSwG0NM1PiczWsFyOnxbl6///O0+8A/a39lq9hYZD+0+ePt1+M53oLwcMjLG57wiIiIiIiIiMqy4CYK0JFshc5yCBclZkJUTvq2jBTpbzMdp2ZA+aL8tG3rGaTgMqWbA5VU3bN4M994LSXHzUohMGcHcHz6fL/R4oJycnCHbREREREQkfsXNJ+9/tM3j1rqmaFdjfDz9NPzgYWg+AS88Cz/5EXxuHViVp1ZkMjkcDux2O16vNywIEswVkp+fH62qiYiIiIjIBNCn7mi44Qb46MfAPhPaTsGvH4RHKqNdK5GEVFhYOGQWmKamJpxOZ5RqJCIiIiIiE0VBkGi59Va48TrwH4bOs/Doz2D3D6JdK5GEU15ezpYtW8K2bd++nY0bN0apRiIiIiIiMlHiZjjMlGOxwN+Vwut7wL0Pujrhp7+CS2+CRYuiXTuRhOFwONi4cSO1tbU4HA48Hg9r1qzRUBgRERERkSlIQZBoSkmFr/8Avn4XHAtASw88WAP/6z646KJo104kYTidTg1/ERERERFJABoOE232GfB/fgaXXAO2ZDjyFvznf0BfYkYRERERERERGR8KgsSCmTNhwwaYPd9cf/MV+FYVdLRGt14iIiIiIiIiU4iCILHikkvgC1+EnFng74BndsD//Tz4e6NdMxEREREREZEpQUGQWHL11fCpEggc7guE/AW+99Vo10pERERERERkSlAQJNbcfCusWg1WGxjArsdgxwPRrpWIiIiIiIhI3FMQJBZ9oRqW3QCWvvUffQ9+vz2qVRIRERERERGJdwqCxCKbDf7lYbj8CjAMMALwrfug6Y/RrpmIiIiIiIhI3FIQJFalZ8H/+RHMmWMGQnp6oPIL4Hkl2jUTERERERERiUsKgsSyGfPhf/8AsrLM9dY2qPg0eE9Ft14iIiIiIiIicUhBkFh36dXwL/9hDpHBAsdPwX2fgM72aNdMREREREREJK4oCBIPlq6GL/4LYIDFCm+8Bf/2GejtjXbNREREREREROKGgiDx4tb18PFPQcBvBkKefQ6+twECgWjXTERERERERCQuKAgST9Z9E5zvNQMhVivUPwF134x2rURERERERETigoIg8cRqhX/5IVx5hdkDxGqFn/4Adj8c7ZqJiIiIiIiIxDwFQeJNShpUbofZMyBgABb4/v+Dvz4R7ZqJiIiIiIiIxDQFQeJR7mz4Pz+GjDQwgJ5e+NY/wyt/iXbNRERERERERGKWgiDxamEB/H/fAqvFXG/pgof+Aw6/Ht16iYiIiIiIiMQoBUHi2dLVsH4jGEmQsQjOBOC/HoITB6NdMxEREREREZGYoyBIvLvtS3DXRiDJXH/DA1s3w9ljUa2WiIiIiIiISKxREGQquOvzcOMy8PvBYoEXXoa6reA9Fe2aiYiIiIiIiMSMpGhXYDRuvO9rIyp3UeA9XOn/aNi2A7ZHOW59/rzHZqdk8W8ln2NZ3vwx1TEqbEnwxX+E0/fBGwch0Am7tkFmO3ykDHJmRruGIiIiIiIiIlEXV0GQI76jWLCct1xX22l6mlvCth2adoqzGecfInIEuPfH/8VT91WSmR5HT096FvzDP0PFveA5AATg0e2QPQNuXguZuVGuoIiIiIiIiEh0xdVwGMOY+GskBbJIeftWdtSdP9gSc2bNgw3/BFlp5nqvAb/8HTz7P9Dui27dRERERERERKIsjro6QN0X/nVE5VJsaWQk2cK2tfXcQU/g1nMe9/t9r/PME/l0tqfy610dXHdtJtdeG2fBkKuc8KX/BTX/BrYF0J0Mj9ZDaipc8z6zx4iIiIjEtZdfT4l2FSQCiyXyN3aWMbydjLN3oAnPYoGLZveQlRntmojI+cRVEMSZf8kFHH3+D/+Fiy9hc7PB7t+10dXby0Pf6+Tb1elkZ1/AZaNh+Z3gDcDPfwkGcPI0PLYbkpOhYDmkZkS7hiIiInIBvvHdOdGugogMkpke4NMfOY1zaSfJcfUpSySx6NdzAIsFPvtZC6++ls6Bt1o5erKLmu/18PWN2VitcRSPtyXBLcXQ3AxP/M68sbcPwe6nIMkK7ypUIERERCSO2eJqQLNIYujssvKDX8zi7UNePrLSR04OxNNHCJFEoSDIIJmZUHq3jfv/LY0jvU1se+cRMh/7KP94xweiXbXRSc2A2++Es83Q+AL426Hxl5B2ygySLLoBUtKiXUsREREZA8dF3dGuwpQwlnxzxjADVcaUu24s15+EHHkyes0tNjo6zejkU405HD6eQsnq08y7yCA9PcqVE5EwCoJEcNVVcPWKl2l84b8IGAH+648/5+Zrr6Lg0jjrepo7G4rvhDNH4dkGIAC/ewqmzQIskHcjJGtMsYiISLxZX3I82lWIeyPN0zHWoMNo8oCM9hpjyTEiE6ulzcpPH5+J52gqAK+8mc53fzaXu247ySUX9zJ9OgxKWSgiUaIgyDA2rMnniTcv5p1mD93+br70g8387l/vIyne+p/OuwLuLIGTL8E7h8xtv/yFGSCxWGHRexQIERERiTMLHNGuwdQxkT0r1GkjkQT453tO8Mivcnn6RTOh4KmzSWyum8vHbzlN/hUdTJsGWZkKYolEW5x9op88GWnJVN11N2nJyQC87XuLr217LMq1GgOrFa68Ef5uI8yYZm4L+OGnm+Ho6/Das+DvjW4dRURERKLEYpm4xaoloZYcO/z9J5v5+C2nSU42Q2Dd3RZ+/puZPOnK5cwZOHIUOrui/J9eJMEpCHION161gM8svQOLxYLVYuEXL/6Gp154PdrVGr3kFLh6BXxmA2T2DUrs7IJt34LmY/D6cwqEiIiIiIhcoPR0uOV97dzzyePMnNb//voPf81m2y9n095p5fhxOHUa/IEoVlQkgSkIch5fv3MVV816FwHDAIvBP/1sC80tndGu1uilZ8HSIvjUPZDUNwqq2Qc/qoZ2L7z5PATUEouIiIiIXIj0dLh6cQ/33HmMqy7vCG1/42AqD/54LsfPpNDVCUeOQGubkt2KTDYFQc4jyWblwbWfJyslA8MwONt1kq9s+Vm0qzU20+bC9avhk5+GYGbzw0fhZ/8JHa3w1t8UCBEREZkAGzZswO12R9zX0NBAbW0ttbW1VFRUUFdXN8m1E5Hxlp4Olyww+NRtp7jlvd5QHhBvi43v/WwOz7qzSE6GM2fg6FHo7olufUUSSUwkRm1oaMDj8QDgcDiw2+04nc4o16rflfNnUPbBT/Pvv30YgD++8yd+tOdqPrdyaZRrNgbzroSlq8F7Gnb9xtz28gF44ofw4c/D2y/BwqvNXCIiIiIyZh6Ph9raWrKzs2lsbGTNmjVDyjQ0NJCfn09RUVFoW3FxMT6fj9LS0smsroiMs/Q0mDsX3rfUx8Vzutn+xAw6Oq34/fDo7mkcPJJC8aozWDADIfZsyMk184uIyMSJ+ifd2tpavF4vpaWloT/2FRUVUa7VUKWrC3nP/KUEDAOb1cK/PfEjDh73Rbtao2e1wmXXwU3FcOON/dv37oWnfwMtZ+CgWz1CRERELpDD4aCyspLy8nJycnIiltmyZQvbt28P21ZYWKjeICJTRGqKGQi5YkEn9372GPPm9Hf5eGZfJt95ZC7eNhvp6dDebg6R6eg4xwlF5IJFNQji8XhwuVyUlJSEtjmdTiorK6NYq8isVgvfuftz5KZOw/AnMfPsar77PX98juFLToHLr4eiv4N3Lerf/tgv4fVnwXsSDr0cvfqJiIgkiI0bN3LrrbeGbWtpacHh0By4IlNFSrIZCMnN9vPFNce5vqAttO/wsWS+vW0uB95KIzkFkpPh5Ck4dhx6NG+ByISIahCkoqIirPsnEHNDYQaaPS2Tfy/+Inkn/5lZ7e/Dvc/Ko4+3RLtaY5Nhh4VL4JNfgnlz+zYa8JNaOPkOnD4Khw5EtYoiIiJTndPpJD8/P7Tu8/mor69n48aNUayViIy35CSYMweSbAafXHWGO245i81m7mvvsPLwjlk81WgHIC0NAn6zV0izV4lTRcZbVHOCNDU1hf7Iu1wuHA5HzH/z8eHCKzl+oIdHf91OWkoSj/ysh3df3cVlC1OjXbXRmznfnBnmc1+F//oGeH3Q2wu11fDVf4OTh8CaBPMuj3ZNRUREprRgfrTm5mZ27tx53vdDH/vYHWHrJWtKuCtCzhERiR3JSTBnLhw/BjcsaeXi2d38+Fez8LVaMQyo/0MOB4+ksObDp0lPM0hKhpYWaGuF6dPNZKsicuGi1hPE5/Ph85k5NRoaGigoKMDn81FRURHaPlheXl7Y8tBDD01mlUM+86lkFi5IprvXj8VI4t+/7aWzO077q82/CqY7YN0/QlpfIKetHb5/PyRZ4fibcOyt6NZRRERkiisqKqK0tJQlS5ZQW1s77HuhoF/96pdhiwIgIvEhyWYGQgzg4jnd/OO6o1zm6Artd7+WzoM/nsuxk8kApKaCLQlOnDQXEblwUQuCBGeD8Xg8FBUVYbfbyc/Px+l0ct9990U85sCBA2HLvffeO5lVDklJga9uSMNqhY6ko/yp+0HKNv93VOpywaxWuPxamLUA/u5L/bPCnDwNtZWQkQ1HX4eTB6NaTRERkURQVFSEw+Fg7dq10a6KiEyQJJuZI8QCpKUEWL/mBDct7R9if/JMEg/+ZA5/ezkDMN+ep6dDj6bRFRkXUZ8dxm63h607nU5279593m9Aou3yy6zcfMcxXp7+AF0pR9j92pPs/FNTtKs1NilpcNm1cHEefPJT/dvfOgK//CVkTzfzg5w6FLUqioiIJAqn04nb7cbtdke7KiIyQWxWs0eI1WLm//joymY+/ZHTJCebCUC6uy088tgMfv1UbmjSxuTkKFZYZAqJWhAkGPwYPOY1uL2pKfYDCl+48wouyV1AIGBgsRr8719t49iZ1mhXa2yyp8OCq6DgA7BiJVjTIWMRvPg6/GEvZE0zp849eyzaNRUREZkSPB4Py5YtC/WODQq+Fxq8XUSmllAgxAY93XDd4nY2fPY4M6f1D7P/4zPZbN4+m9a2qH93LTJlRO23KRj8GK7HR6wnSAVIsln53ufXk2xNAaCl9wxf+v4j9PoDUa7ZGM1aALmzYOWnYfnHoe++2PMH2LcfsmfAmy9C84no1lNERGQK8Pl8FBQUkJOTE7Y9GPwYOGuMiExNVos5a4wtCbq74aLZPWz43DEWX9ERKvPGwVS+tW0ubx9OiWJNRaaOqIYUV61aNaTHh9vtxm63x0UQBCB/4Wy+dNNd+AMGFgs8d+xpHvrV3mhXa+wWLDZnhLn1FrhsYf/2nb+Bd94B+wx460XwnopaFUVERKaCYC60wWpraykpKYmb90IicmGsFpgz28w72NUFGWkG6z5+iqL3ebFYzDLeFhvf+9mc6FZUZIqI6hS55eXlrFu3jpKSktC2zZs3h6bNjRf/cMf72eP+G+5TL2K1Wvj+n3fw/iWX8+5F86JdtdGzJcGVS2H/Xii5Ax7+CZw8BV3NsPl/Q1kFzLkM3ngBrrzeHEYjIiIiQ/h8PjZv3kxLSwsej4fq6moKCwtZvnx5KPhRWlpKXV0dBw8eJDc3F4/Hg9PppLS0NMq1F5HJZLHArJlw6jR0dpqzwqx0+pg/t5uf/noGHZ1W/P5o11JkarAYhmFEswIej4ft27eTm5tLc3MzS5YsoaioaEi5vLw8Dhw4EIUajszBE15u+ff76DLaALg06yp+9bUN2DNTo1yzMfKegjeeg0AyPFQNp18BDMi1w1cfMHOEtLfAomWQYT/v6UQuVKy0AbFSDxGZfLHy+5+Xl8eBXT+PdjVEZAIYBpw5C+3tZiAE4HSzjR/9chZHjifz9MvvjYl2SCSeRT3DjsPhoLy8nNLSUsrLyyMGQOLBgtk5fG315/AHzJjSmy0vU/nThvjND5IzE+YtAlsAPvtJcy4vAG8nPFIHAcOcPvfVZ6AjTpPBioiIiIjEEIsFpk+DzAyzRwjAjFw/X/nMca4vaItu5USmiPMGQV5++eXJqMeU8LlbluKcv9ycLcYCv9q/i/qnX412tcZu7qVgnw5zL4fPlIItw5wx5qgXfvGYmco6PRMOPA1vvwQnD0JbM3R3RrvmIiIiIiJxyWKB6dMhO6s/EJKSbLDmw2eiWzGRKeK8QZCsrCx+8IMfhNYbGxsntELx7vtf+gw5yTPBsDCr9YP8ansuh062RLtaY7fwavNn/k3wyXv7Z4zZfwAanoKkFMieBh1tcOxNePVZeNkFL/4OXvkLvOOG00eg3WfO/SUiIiIiIuc1bRpkZ0NH30QxwSSpInJhzpsY1eFwUFdXF/r5r//6r+cs39raSlZW1rhVMN7kZqXx758s5T8ebGNaYBFnW+BHP2/mH+9JIyMtOdrVG71gotSXG+GGpXCmGRqfMfftfRpyM6DQCSmp5jJQIABtXvCeAH8AMMxhNSkZkDXdHE6TlgUpaeZ1REREREQkZFouWK3Q3Azp6dGujcjUMOwnz8bGRvbu3YvL5eLgwYNs2rSJp556KmKA47e//S0Oh4OrrrqK+vp6PvnJT4adx+FwMH/+/Im5gxj04RsX8VbTWX7b4Cc5ycZzLju78o7xsZUXk2SLehqW0UvNgIVLzKlxV6+EZi+8/Cr0+qCuBtI74NoVQ4+zWiE1zVwG8vdC83E45QEMMICkZEjPgsxpkJljBkaS08xziIiIiIgkqBy72QvkjEbDiIyLYYMghYWFFBYWAnDzzTezatUq7rvvPr797W8PKWsYBt///vfZv38/FouFgwcP4nQ6Q+dobGxMqCAIwN2ftfPSvmZOHbNhsVj59c5MFl5yjKXvisNpcwFyZ5tT4558Bz75Mfj+t+CNNwADfrIF9j4JOTmQMw1yZ0DOLHOZNsfs9TEwmGFL6uv5kRF+jd5uMzBy7A2w9JVPSoFMO2TmmrPQJKeZARIRERERkQRhz9Z3gyLj5bxjEFpaWli1ahUbN26kpaWF3/72t9xyyy1hZVatWsWqVasAqK6upqCggO3bt1NRUYHdbqegoCAUUEkUaak2/uErafyvih56LG28bv0p33zUzve++BnmTo/T4UJzL4X2s9DZDn/3OfjPl8Hrg95eePW14Y+z2iA702y97XbIyQX7gGDJtLkw3QGpKWbQY7CONmg5Y17HghkgSc2A9Gxzqt70LDM4khzhWBERERGRKSArM9o1EJkazhsEaWpqYv369QBkZ2cPCYAMduutt7J48eJQUKSlpYXs7OxxqGr8WXxlJituf41vNdYQsLXzt5Ow5Ykr2PjJ98dnfhCrFS65Gg40msGH9f8Cm/+fGQg5l4DfLOP1AYcjnDcdMq8ygyDZWWYGqPQeSPWas9PMuARmXGYGUbKzICmpP9/I2ePmhOoYkJwMKenmkJrsaWZgJDVDYXMREREREREBRhAEGdyD4xe/+AUOh4Mbb7wxYvnFixeHrSdqACTo7o9fyqP7HbzdegCA7S88zruvmM+tN+RhtcZhiufkFLj8ejNR6kWLYNP34O0m8J6C5lPgPQ2+ZvB6oaUFWlqh6zyzwlj6AkJd3dB1Bk6dga7j0N0XMEmeDWkDhlNlpEPyGUjvhmy7OQQnZwZMm2kGZ7KnQ1q2OeTGMJRvRERERERERIARBEEG27dvH5s2bcJisbB48WKcTidLlizB6XQm9Kwww0lLTaJm3ef4xEPfoId2OgItfLv+N8zMySQzfeKHb1w6dxpZ432d9Cy4dAm8vQ/sM8zZY86ls9VMhHr2BHhPmsES31nwNoPPBx0p0GMDv7//GKOn/7FlUK+Z9g5oPwz+c/RAsVjMPoPZWWY2qSw7ZGWbP3NmgH0mzJgPsxxm0ET5RkRERERERKa8UQdBKisrqaysxO1209jYiMvlora2FovFgsPhoLCwkLvuuot3vetdE1HfuHT15XMpXf4JvvvnH2OxwBu+/fzsD3/lPVdeFiqTmzwTy4DJv/2BXny9Z0d0fgsWclNmhm3rCXTTTRut7V0U5i8YnxsZaNpcczjK6cOQlXvusmlZMDcL5l4+fBnDgI5Os+eIrwVOvQan3zaDJd3ToCPZ3NfaCgEjPEgy3PlaWs3lyLHhy1mtMO1SmH55f94SWxJ0vgH0mMGU9CvAmmo+tljBakD7ATPXSfAcwVwlWPp7mFhtfRO6922z9j0OnsdiNbenXwrJGf3nIgAdb5nrtmTIurLvmL7zBDqg69iA+ljA0nctW9956Tt3UhKkZprDhNIyzOFBKRmQlmlOUWxL7q9L6Pr0101ERERERGQKGXUQJCg/P5/8/HzuvvtuANxuN1u2bGHfvn3U1dWRn5/Pj370I/UO6VP2sZv40yv72Xf6WQAeP/AEjx/o33/N8WpsRmpovdN2nP2zHhrRuW1GKtccrw7b1pLyGq9Nf5Ar7AU89vV7JyYHybwrocNrJkpNyzh/+XOxWMxhLhnpMGcWDAgQhQkEoK0NTr0F3hNwtm8Ijvcs+PqG4PhazIDKSAQC0NYL3afg2Kn+7W37IdB3jsyzYB3QQyTQC60vnefExjlXw2S8ywyyhMr6oW1f34oVMpeEn8Df1hekGcH5LX37LAPKDByFZbGC/VpISTYDJjYrWLuh921IToKUHMh+l7k/ORmSbGC0AGfMoVHJKeYUyCmpZi+alL48LCnp5pKabgZcUjMgJRPSM80eN7YkM8hiS+qvx8AAzHDbRERERERELsCYgyCD5efnU1NTw8MPP8zOnTupra3lc5/7HI8++uh4XSKupSTbqF5Xwsf/801aA5Mxybf5ifc1bxPbf/8if190niErY2G1wsJr4cBfoLU5/NJGoL93hEFfLwjCe0FYrH09Jfr2D+wdca5rZmdD9tXnrltPp5k0tfk4+E6ZwRLfWTNYMjBfSXfP0OE20JdsdRgWS3ggIXKhc66GsdnAOuBX0RjQmwSbmdMkTNc49tKwAVbo9psLQG8LdJzs290CGYMCXN3HoWtAgtuwp8oYJiAzcKMFkm1mwCU1B7KvMoMrNpsZaDHOAifNoEz6XMhYYO4XERERERG5QKMOguzfv5+6urrQtLfz588P25+TkwNAaWkpTqeT//iP/+Cf/umfxqe2cW7R/Bl87cOf4aE9j9HpbwvbN216ABuB0HqnxTpkiMtwrEYq02YEwrYdsT4NmJ/XH/nz//CZFdeRkjwBHySTU+AqJ/R2m70qjL56GIH+dSNgTm9rBMDfawYY/L1mj4pAwOz54A/u73scXDdP1nc+oz84EQyeGAPuOzj0JBiESc82E6RetKgvyBIMxAwIwnS2QmsHtHVCaxu0dZh1ap8P/r6eIGmXgyWlrx4GEIDW6X3XMcLrYRjm8cEP/YbRXx9j4LYB95K+CGxp/dsCfmgJns8CWYvDj+1tgTZr+HUDwedpwHkDAfNcPb3m89/T0/+4t297xOjMwP9LEYItgwNElkErEQM+gzb6A31LD9Aevm9gUtwUH6SG/66IiIiIiIiM1aiDINu3b8dut/Pzn/+cioqKsOSoAHv37uWTn/wkYPYOaWxsHN8ax7nimwqYMz2LnIz0QbPDdPctQSnAP4zizC1ha68cvJH/9fhfAXjT9zJ/bnqHD103zBCTC2VL6h/WMBmCwZXAgIDL4KDLwHKDgzABf1+AoBcys2G6MSDwEgwiLDtPJd7bf43RMkZyzKrRn/dc1zP8Q4NIgeDzYUBXT19wpMfsRdN1ynwcsIF1uhk0CQTMn+3HoOOoub+7u++Y7r4gS9/SO1zgpWdQEGUMQZgJ4vP5qK6uprKyctKuKSIiIiIik2vUn1wLCgq48847AfB4PGzfvp29e/dSV1eHw+HgG9/4BgA/+MEPWLVqFQsWTEBSzjiWlpLEVQvmsO+toyRPYBf/uTOymZ9xBZ6217FY4Tv1v+UD16yPz2l5B7P2Jf7UCImxCQY/IgWPgsGh3r6eOr09/QEif4+5ze83l0AvYb10BgoGeiwWM6YR/H9nYAZkervMwIk/AKSbganevqXrjDkcp6cLDDtYss3gyfNNE/q03HfffRN6fhERERERib5RB0EKCwtDAQ6Hw0F5eXnEcnv37mXXrl2sX7/+gis51cybkU1GajKBwR8cx9n6DxXx9V8/hM1i4YXjz/HyOyfJv3T2hF5T4oDVCtYJmJ45MExQxRgUdOnp7u+R4+8Ff3d/UKW3F/xz+4IowwRYJoDL5cLhcODxeCb8WiIiIiIiEj2jDoI4HA4+//nP09jYOCQfyEA//OEPL6hiU11uVtr5C12gT33oOh7cczEnOo5gsfj59q/3UFv2qQm/riSoieyhM5ZhRyPk8/kAWLBggYIgIiIiIiJT3JgH3BcWFo5nPWQCWK0W1i5fjdGXg+GPbzVy7HRrlGslMgbjNhvOUPX19Tidzgk7v0ySCQyUiYiIiMjUMXlZByUq7i66kezkXAwDuow2vvP4H6NdJZGY4XK5WL169aiPy8vLC1seeuihCaid4O9LptvVBe3t4PPBmTNw/DgcOgRvvw1vvQWvvw5vvgkej5msV0RERERkGJM4pYdEQ0qyjeJrV7Ltmf/GZrXweNMfuK/rFtJS9dJLYvP5fPh8Pux2+6iPPXDgwATUKEH4g7MzBfqT7A6eWcgwzO1gJtc1jP4prq1WsNnMnymDctv09sI770BODsyYYZYTERERERlAn4QTwD989ENsf/430JNBxtnl7Pmzj9tWTI92tUSiqr6+npKSkmhXY2oIBjQiBTaCUyWPNbAxGklJ5tLWZvYamTEDcnP7ryMiIiIiCU9BkAQw3Z5Oaf4/8fv6GaQlp/DEE53c+kFjakyXKwnN4/FQVlY24vKFhYWUl5fjdrvHnAck19oLvlNjOjaMxQq2QQmSA70Q6BvOYUkC26CAgL8bjN4Lv3bw/Jn24fOtDAxsBAL9wYzzBTaCj8crsDEWqanmzzNnoLkZZs+GzMzJrYOIiIiIxCQFQRLE54uvYO9vfQQMA8/BJJ7b186ya/ShQOKbw+Fg586doz7O4/Gwb9++sG2NjY34fD6qqqpYsmQJRUVFEY/9afIh+Mw4JIa2ZYB9Sfi2rpPQ/qb5OGUWZF4Wvr/9beg6fuHXDp5/VgFc/x5YtgxmzuwPbAQC/T01YiWwMRbp6Wb9jx6F5GSYO7c/QCIiIiIiCUlBkAQxfTq8770p/O4PXSQn2fjvX3YoCCIJq6ioaEiQo7a2ln379lFeXn7+E/SOQ/LNgBXavIPO2zrg3B3AoP3dbeNz7eD5Tx6F3/4annoSLr0cbrgBrr7aDB5MFRYLZGSYPVYOHoTsbDPgk6Q/fyIiIrFkw4YNNDY2kpOTE9q2ePFiysvLcTgcEY9paGigrq6OpqYmAAoKCqisrMThcFBRUUF9fX3Y+ZxOJxs3boyYE662tpb6+no8Hg8FBQVs3LiR/Pz8cb5LiQV6F5hAij+Syu//2I0FC8++epS/vdrLtYvmRrtaIjGhubl5ROW6rUBK8oVf0JYCqSlDt9F37uTkob0tLKlgHYdrA9iSAMM8X6AT3nwF3nkLHn8crrsO3vMeuPji8blWLLDZzCExnZ3mrDLTp5v5QiZw+mUREREZndWrV1NZWRlar62tpbi4mJ07dw4JhFRVVbFjxw5qampwOp34fD7q6+vDytx5551hX3BVVFRQXFzMnj17wspt2LCBlpYWampqcDgcuFwu1q5dy/333z9s72CJXwqCJJCFCy1csuQwTx15lPa0N6l67L38tPzvo10tkajyeDxs376d3bt34/V6qaqq4tZbbx028n9Hx0IO/PTpSa7lBDAMeP5Z+O0T8NZh6O2Alleg5xJo7IDGRpg/3+wdcs01kJZ2/nPGg5QUc2luDs8XouSpIiIiMae0tJTm5mbKysrChkC73W4efvhh9uzZEwqO2O328ya9r6ysZOXKlbhcrlB+uIaGBvbv3x8WGHE6nWzbto21a9fidDrHNJugxC4FQRLM+z/Uw69/8QZWi4W/HHqaQyc/wfxZ+qWWxOVwOCgvLx/ZMJigzNwJq8+ket/NUHAVNP0VnvkxvBCA5gPgnw8pc+DQIXP5zW/g2mvN3iHz50+NgEFamhkIOnbMHBozd+7UCfSIiIhMIeXl5eTl5eF2u0NfUm3evJmSkpJhh8mci9frDRsiU11dzcaNG4eUy8/Pp6CggLq6OkpLS8d+AxJz1A84wRQvX8LcrHkEAga9RjfVv9xz/oNEZOqaPh/yFsJiC3z2KvjElXBNAVgDZo6SQA90dcHTT8NDD0FNDbhc0NER7ZpfuGC+EJsNPB44csRMDisiIiIxJT8/P5T3A2D//v2jnunP7XZTXFzM6tWrw3r8ejyeYXsAL168eEgyfYl/6gmSYKxWC5+9sYiqPT/EarWw++U/0NpxG1npcTDTg4hMjDnvhvd8E16ogsUfgJtugvYOeOEleOZZOHkarGlgTTIDBb/6FTzxhDlM5oYbYMGC+O4dEswX0tMD77xj5gqZPl35QkREJO49+aSZ7qura/KvnZoKt98ON9984efKycnh4MGDoXWPxzOiISoPP/wwO3bsAMDn81FSUhKWc8Tj8YTOH0lubi67d+++kKpLDNI7vAR0d1EhOam5BAIG7X4f3/vN3mhXSUSibcZiKPy/kH61+U4pIx2W3wBlX4LP3wUFl4GlG4yAWb6nB559Fr77XfjWt+DPf4b29ujew4VKTjZ7hvh88NZb5s/gFMEiIiJx6MknoxMAAfO6Tz45Pufyer3k5uaG1h0OBz6f77zH3X333TzzzDM888wzHDhwAIfDwbJly0LBj+BwmuD6YB6PZ0xDbiS2KQiSgNJSkvjY1SsAs2fI9r8+SW9vIMq1EpGoy5wDV14LVj90dJrb/K2Q9ATcPAPKvwy33AizsoEBbcaxY/DrX8M3vgHbt8Obb8Z38CAtzVxOnjRnkon34I6IiCSsm282e2REQ2rq+PQCAcLygYA5TMXlco36PKWlpaxevZqqqqrQtuBsMJE0NTWxePHi0VdYYpqGwySor9z2QepeeIIufxenuo7yiz+9xF0fvDba1RKRaEubBpcWgGc/tLaA71Hwt4HXBWlHoPBjUHiDmUPjhVfB/Vp/Ho3eXnj+eXOZNcscKnP99eZQk3hjsUB6OgQC5hCgtDRzJpnB0xaLiIjEsJtvHr9ARLTU1taSn58flgNk/fr1FBcXU1paOuqeGg6HI6znR2lpKdXV1UOSn7pcLtxuN9u2bbug+kvsUU+QBDV7WiYfuuK9BAwDi8XC5t/Vn/8gEUkMqbmwIB+yLdBj6e/V0fk2HP4B9ByHhZfBRz4IX/0s3L4KLr44/BwnT5qzynzzm/DTn8Lrr8dn7xCr1RwiEwiY+UJOnAC/P9q1EhERSQh1dXVs2bKF+++/P2x7fn4+d999N+vWrcPtdoe2NzQ00NDQMOz5XC4XW7ZsCQt4lJSUUFBQQHFxcSg44nK5KCsro7KyUtPjTkHqCZLAym5bxe7XfkcgYPCW91X+8OKbvP+ay6JdLRGJBSk5MP8aSE6Hd34HHU+bvSP8LXD0JzB9JdiXQnoKXJcG774STnTBC/vgb3+Dzr7hNL298OKL5jJjRn/vkOzsqN7eqCUlmUtbm5krZNYssNvjOyGsSJQYhkGAAH4jQMDQcFwR6VdfX09TU1Mo38fixYvZuXNnxN4e5eXlLFmyhE2bNuHxePD5fKxatYry8vJQmR07dtDY2Bg6n91up6amZsjMMlu3bqW2tpaysjI8Hg8FBQURy8nUYDGM+PhqLi8vjwMHDkS7GlNOSdV/8RfPM1gsFt4/93Z+9LU7ol0lkYhipQ2IlXpMmm4vnHkbDr4Ivnqgu39fZj7MvBWsKWbCVH8XWJPBYgf3AXNa3QGZ3EOsVsjPh/e8BxYtis9AQkeHeR+zZ8fncB8Zk1j5/c/Ly+PArp9HuxphBgc2DMOg199Lb6CXXsP8GTD8ZrmA+dYz+KsfF29ERWLEzXf9c0y0QyLxTD1BEtyXV63m9e9Ymdn6Ac4euYgTJ/3MnmWLdrVEJFak5MD0heY0sodnQfPj0Hvc3Nfmhu4TMOfjkDwDktLB8IP/NBQ44Ppr4dhJeOYZeO65/t4hgQDs22cu06aZwZClS2GY6eliUjBfyNGjZua3OXOUL0SmnHMFNvxGLz0DAhvmYgY2gkENq8WCxWLF2rfYrHp/ISIi0aeeIELlN3p5+vk2rFYLt6628sXPZ0W7SiJDxEobECv1mHTdXmg9AkfPQPOT0LGvf58lBWbdBplX9W8L9EKgG5JzzUBKr98Mejz9tDnjymBWK7zrXeZwmbw8cz1e9Paa8wBmZ8PMmeawGZmSYuX33zKo99TB534VsdxPd+7ma9/8Xmj9rjtu5v/d9+WIZW/99FdpeuXN0PpvHqnm6quuCAU2wAxuNL38Bh/+zD+Fyl191RU8WfediOf8p3+t4ZFH+3OOVVVs4O8+cWvEsnOuLgpbP/5S5DH9P/7vXZRXPhha/8zHV/Mf/7ssYtmbS77CSy+/Hlr/7faHuGbxlUPKvbj/NW5Zc29oXfeke4okVu5p0aJFMdEOicQzvVMTPnJbEs+/ZAUs7Pmdn8/dZZCREYfd00Vk4qTkQJYF5lkg+cPgnQe+p8DoBaMbTuwE+w0w/YNgsYE1yVx6W6DXC8nT4d3vNvOBHD9uBkOef75/+tlAAPbvN5fcXLNnyLJlZk+RWBfMF9LRYQZ4Zswwe7XEUyBH4pp146cjbrcc6Qlbz3jWxUX3vYKvp40OfxcArSRz2Mim69Cx8HN++z6s2TZImcv07EWh7TlnvWHlDM+b9JR9POL1A690ha1nPvoTpj/zS053++gN9AJw2kjnFOlDjg2eMznjSrLTLwptT3szfIhdwPUkPWV/jHh9w9MRtm5/8BtMz0nheOeZ0LbDRhZvtAyq94B7yspZRkpSf/1sz/4trKy/7vv0/OkHEa8/2PR/KaUn0MuZbl9o2xtGLkeOhE/DHbwnw5LEjOmD8hEcPhy22ltdTk/20B42vS3hCZxth99h+r+U0uHvwtfTZpbByhtGLmcPnAkrG7yngC2bmbnXnfOehnvt/YP+76X+9Y9M/5emiP/3Ojzh//dC9zTo/16G/u9FvJaIjJ6CIMLVV8Pll6Rx4M12ev2w67edfOJjQ/8oiEiCS+nLjm45CknLwDYbfI9Db7O53fc0dB+BWXdAUl/iU1uq+bPnrBkMSZlhDh35yEdg9WpoajIDIm/2fwNNczPs2QNPPWX2CnnPe+Cqq8whObEste9ez5wx72HWLMhSzzoRERGRWKKvqQSrFT56ezKGrZMz9t/xb3+ootevbO0iEkGKHTIughmZMONyyP0UpA/o4tvpgSM/gM5BCVFtaWBJhs6T0HYI/J2QnAzXXQdf+AKUl8P73x8eNDAMeOUV+PGP4d//HRoazABDrEtPN/ODHD9uTqvb1XX+Y0RERERkUigniADQ1WVQuOlrnO48jgF8/ZZ7+MJthdGulkhIrLQBsVKPqOv2Qedx8HXD6bPQ8zw0/5H+lIgWmP4hc4jM4NlfDD/4u81EqqkzzBllgnp7zSExTz8Nr70W+dpXXgk33mj2Don1/Bt+v5kQNivL7BkS6/WVc4qV3/+8vDye/NkDE34di8XMPzIwuengfCQiMrlu+sSGmGiHROKZ3o0JAKmpFlYXOHnk2V9itVj40d567rn1RqxWvdkRkQiCQ2M4Dsmz4fi7YfZFcOoxCHQABpzpyxmS+97wYy02MwAS6IU2DyTbIWUaWG1mkODqq83l9GlzZplnnoGWAQOnX3vNXLKyzNwh73mPmZA0Ftls5hS6XV3w1lswfbqZ50T5QuQCzcqaNfEXMQx6/L34DX//NLeB/mluB84EY8EMmFgttr6fCpiIiEhsUhBEQr784Q/y6Iv1dPR0cqTNwxNPv8zthYujXS0RiVWhQMgxmDcXjgBz18Hpx6DrMNiyIfscSfVCyVPbzQSqydPMBKzBD04zZkBREdx8M7z8Mvz1r3DggDlMBqC1FX7/e3O5/HJ473th8eKhPU9iQUqKuXi95jJzpjmbTCzWVeJCelLa5FwnOfJ2wzDwY06dGzACBPx+egNmsMRv+EMBk0AgEDZ1bvB//MCAiYiIyGRSEERC5s3I5gOXL6f+lacA+N6TDQqCiEQwGaMI4+aDwcBAiONiOHwMZpZAy5/MKXNtmec/hy3F/NnTDL0+M3lqUkZ/gMBmg4ICczlzBp591uwd4h0wU8Abb5jLnDlmbpHrrovNRKppaWYQ5+RJ817mzjW3icQZi8VCEjazZxfAOX7dBgZM/AE/RiAQFjAxiIuR2SIiMkUoCCJhvrL6Zva8/nt6ev24T+3j2VcPs3TRxdGulkhMeePMGxN6fgODednzyEwZQQAhFgQDIZ3HwDEPjhyDrA9AasrQsu2vQ/pCsET482PrCxB0njB7iKTNMrcNNH063HILrFxpJk3961/NXiLBwNTx47BjB/z2t3DTTeZQmeCsLbHCYjGTpwYC4PGYw2VmzTITxYpMQWEBE6v+n4vI+TU0NFBXV0dTUxMABQUFVFZW4nA4qKiooL6+npycnFB5p9PJxo0bsdvtQ85VW1tLfX09Ho+HgoICNm7cSH5+/qTdi8QeBUEkTMGlc3j3vOt4+uCzADz4xG/58aJ1Ua6VSGzJSMmY0PMbhsHRlqMsyF1Aii1CICEWpdjND/cdR2H+PDh2HDq7IG1AAKL9dTheB6kXw+xiSBr6RgWLxcwXYvih/Uhf8tSZQz84Wa3m0JfFi83paPfuhb/8pX8mluZmePxxc5pdpxOWLzeDDbHEajXr1NMDb79t5gqZNi02e7CIiIhMkqqqKnbs2EFNTQ1OpxOfz0d9fX1YmTvvvJPy8vLQekVFBcXFxezZsyes3IYNG2hpaaGmpgaHw4HL5WLt2rXcf//9FBUVTcr9SOxRZjYZovSDt5CUZL4J//PbLg6d9EW5RiKJxWKxkJqUyuGWwwSMOJquOjkb0i8CowMumgsZ6dDeae7rbYWTj5mPuw7DmT3Dnwf6kqdmQMBvJk/tPGU+jiQ3Fz78Yfj6180cIgOn2W1vhz17zCl2H3ssNqfYTU42gyE+nxkM8Xr7e7aIiIgkELfbzcMPP8zOnTtxOp0A2O12SkpKcDgcwx5XWVkJgMvlCm1raGhg//79bN26NXSs0+lk27ZtbNq0CZ9Pn3ESlYIgMsSK667gstzLAOg1evnWY09GuUYiicdmtYEBx1qORbsqoxMMhAQ6YPYsyLVDW4eZGyR3OWAxE6bOWDWy81mTzGBIbzu0e6CrGYYLDKWnw4c+BF/7Gtxxhzl0Jqi72+wt8sADUFcHx2LweU1LM5dTp+DNN83l7bfNITPHjpmz5Xi9ZmCns9O8J79fARMREZkyNm/efN6Ax3C8Xm/YEJnq6mo2btw4pFx+fj4FBQXU1dVdUF0lfmk4jAxhtVr4u+Wr+D+/eYNef4An9v8P/9pxO1npcdItX2SKSE1Kpb2nnbMdZ5mWPi3a1Rm55GzzZ8cxmD7N7Olw4hTYb4DUeYBtZAlTBwpLnuqFlJnhyVPDrp8MhYVwww3w0kvwP/8DR4+a+wIBeO45c1m8GD7wAVi4cGz3ORGC+UIGMgwz4NHZadY/EOi/7+C0G8ElKclckpP7F6vVXGw282e8JN0VEZFx8fiBx/nNq78ZUdmbLrmJz1z9mbBtj7z0CH96508jOv62Rbdxe97to65j0P79+yMGLs7F7XazadMmVq9eHZbrw+PxDJv7Y/Hixezbt2/M9ZT4piCIRHTXB67ju/8zm5NnO8jxvo/f/aGHjxQpCCIy2TKSMzjZdpJUW+qE5yIZV8nZQF+OkOws8wP40WOQPt/8ID6Y9xlIc0Dq3HOfdyTJU4OsVrj2WrjmGnj1VTMY8uab/fv37zeXhQvNYMhVV41PgMAwAKP/p7lxwLoxoFzfvkCgr0wACGD2mEkBi9UcGmSzmcGNkVy7p8fMjRIMmEDf/KR91wsGQgYHTJKS+gMlwaCJAiYiIjKJPB5PxOSmgz388MPs2LEDAJ/PR0lJSWhITPA8QFjPkIFyc3PZvXv3ONRY4pGCIBJRks3KxveW8bMfpdPVBbuegNtuifzZRUQmVkZKBkdaj7AwdyFJ1jhqtpOzgIvMQEhGRt8UukchOSn8A337a3Dmt+aMMTOKIPuac583lDw10J88NWWGGTQwBgcZ+pYrLjWXgwfhD38Et7vvZAa8+Tq8+RrMmQ03OeHqfLOxM/oCEwwMUEQIZkQKahiE99YAsPT9Y/SVsQQPtfTvCxU0oHvgMJdgjw+r+TxZUswg0MBAicVmPh5JsCRYr97e/oCJYQyoa1/QJBgwsdn6AyUDe5gEgyYKmIiIyDhwOBwjytVx9913hyVGra2tZdmyZezcuROHwxEaTjNcbxCPxzOmITcyNcTRu2mZbLd/cC6/f8LP20daOXSkl+eeD7BsqaIgItv+tu28ZS62X8zNl90ctu3FYy/ywrEXznusPdXOTQtuYlbmLACsFispthQO+w6zIGcBlnj6sDkwEJISDIQcg0APpCSbs8CcajDLGr1w6jdm4tTpN59/Kk2LtT95avuhAR/C+z7ADwkyADNt8PEPwvuuhr1/hRfd4O8re9QDO+pgtx2W3wDXXwMpqYMCFFawRvn5NwIQ6DTzgZw3UJIMtuShgRLoD26MZDaaYMBkYB6SwQETi6U/IDK4h0lqqqYAFhGJktvzbr+gISqfufozQ4bITJTFixfjcrlGPXNLaWkpHo+HqqoqHnzwQYDQbDCRgiBNTU0UFhaOS50l/igIIsNKTYUVH7Lxk+1JdPb08Mtfd7Fsafr5DxSZ4vaf3H/eMpFmdTnZfnJExwK8dPwlNjo3kpViznSSZE2is6eTE20nmJM1Z3QVjraBgZCkDHDMgyNHoasbUlPgok/B8Ueh56RZvuUF6DpqTqObPIJcKFYbWEc5VOii+fCJ+XDzStj7NDz7AnT3mPt8bVD/O/h9IxQuhRuXmTPdxAqLtS+QcY4/4aFASRt09w2voe+HxWI+GEmgJHTNAQGT8wUzDMMMlPT0hAdMrFYz30lWlpkAVkEREREZZP369RQXF1NaWjrqnhoOhyM0DAbMwEh1dTWlpaVh5VwuF263m23bto1HlSUO6Wt9OacPfhCm283x9s+8+jaul45GuUYiieGDCz8YCoAEpSWn0dLVgq8rDqd0S86C9HnQ22Z+GL54nhkAae+E5Bkwby1kDvimpvsYHPkhtL8+sfXKscOtN8PGe2Hl+8ODHR0d8Ls/QdVD8MRvwRtHz7vF2jdcJtUcLpSUZi62NLCmgrUvx1Og00w223Ec2o9B22FoOwht75g/249Axyno8kJvB/i7zF4855q6eeDwmbQ0M/CRkWE+7u2FEyfgnXfM/CyHD5tTA3d3a5YbEREhPz+fu+++m3Xr1uEODV01p7ttaGgY9jiXy8WWLVvCAh4lJSUUFBRQXFwcCo64XC7KysqorKwcUe4RmZrUE0TOKTcX8m88ytPPPII35RX+Y9cynFd/OdrVEomqz137ufOWyUrOGrLtmjnXcFH2Rec87mSb2Rvi/QvfH3F/RkoGx1uPk2pLJTUpdQS1jSHJmZBxMbQfhqRMmDsHTp0GbwtkpsOsj0LafDj9JNDXk+H4Dsh9r7kM7p0wnjLS4QPvBecN8NzfzN4hzV5zX08PND4DTz8L1xTATYXm9L/xbkw9SiA0vMhiAYK9R4I9SgbmKLEOfc1stvDZb/x+OHmyfzhNWprZUyQ93QyixNPQLxERGRfl5eUsWbKETZs24fF48Pl8rFq1KiwHyI4dO2hsbAzlD7Hb7dTU1OB0OsPOtXXrVmpraykrK8Pj8VBQUBCxnCQWi2HEx1cveXl5HDhwINrVSEjPv3aEuzb/bzq6e7BZbDzx1f/L4kumwAcAiSux0gbk5eXxP8/+T1TrEDACdPu7uSTnEmzWEeRziDU9bf2BEIsFzjbDqTNmIMRigc5DcGIn+Fv6j0m/FGZ9DGyTNEOO3w/79sOfGuH4yaH733UlvM8JC+ZPTn1imRHoW/zhPUSC0/YSXEYgOIwmAFgMc1xmViakpUJKStSDInnXr4qZdujArp9HuxoiEgV5t94VE+2QSDzTcBg5r3dfOY+C2UuwWqz4DT81v9F0UiKTqcffw7a/beP1M+bQEKvFitVi5WjrUeIkjh0u2COkt83sATAtF+bOhrYOc5aStPlw8echbWH/MR1vweEfQNeRyamjzQbXLoGvlMLflcDCQeOSX3kNtvwIHv4xvPp6Yg/lGDL0pm8ZOPTGmjyyJTkNMrIhKxsy7ea25lY4fAreOQpHToKvA3oMsCT3nXsSFxEREYl7CoLIiHz+A7eQkmTDMOB/3vgzJ862RbtKIgmhx9/Dj1/8MftP7mfrC1s5cMr89ifFlkJXbxenO05HuYZjNDgQkp0FF8+Fzi6zN4AtE+beBTnL+4/x++DIj8H3/OQFHSwWWHQF3P13UPp3Zg+Qgd72wI/r4LsPw4tNZhBHxo/Vag6RyUyH9DTz9WhuBs8RePNtOHTY7EnU1ZXYgSgREREZMQVBZESKlubhsC/AarHQ2dvFQ7/5XbSrJJIQWrtbOd52HICeQA8/evFHoRlm0pPTOdtxlrbuOA1KDg6EZPRNodvdaybQtFhh+gdgzp1mjwIA/HC6Hk49bibonEyXOOAzd8KGe+C6JeFT5R47Ab94DP7zv+Avz/bPNCPjy2o1h8hkppt5XCwW8HrNoMhb74CnLyjS2amAlIiIiESkIIiMiNVq4VM3BnuDGDz20u9o79SbfJGJNi19Gl9Y+gWmpZtTxfYGevnJSz+h6UQTABnJGRxpOUK3vzua1Ry7YCDE324GQlJSYMHFEDDMKXQBMq6EeZ+HlAFTA7fugyNboefM5Nd59iz4+Efgq1+GwmXhU702e+E3u6H6Ifj9n80ZZmTiDAyKpKeBzWrO4nPoqBkUOXgITp9RUERERERCFASREfvsh5YyPX0GVosFb1czW3/7dLSrJJIQpqdPZ/3165mePh0Af8DPIy89wovHXsRisZCWlMbhlsMEzjVtaSxLzjSnzw0GQpKSYP48SE4yh8cAJE+Diz4HWVf3H9dz0swT0halBHG5OfDhW6D8XvjQTeGznrR3wJ4/mNPr7noyvqbXjWcWizn1crCnSJINWlrh8DEzKPKOxwyKdHQoKCIiIpKgFASREUtJtvGxa1aQnGTDHzB45OndBAIagy2JpzfQO+HLYNPTp/PFpV9kVoY5M1PACPCzpp/x/NHnzRliDDjWcmyyn4rxk5wJ6QN6hNhsMO8i89v99k6zjDUZZt0OM28F+mbFMbrNYMhARmBy80NkpMOH3gcbv2IGRXLs/fu6e8D1V/iP78LOx+Hkqcmrl/QHRTLSzNcpOak/KPL2QTMocuq0GRTx+6NdWxEREZkESdGugMSX9UU3sf3539Dd6+No62F+ubeJj9+0JNrVEplUqbbU8xe6AAEjQGt3K2lJaWFT4Oak5bB+6Xpqn6vleNtxDMOgzl2H3/CzbN4y2nvaOdtxNjR0Ju4kZwAXQ8dhcypciwXmzDa/uT/rhay+6XGzr4OUuXDiUUieFZ48FeDs78H7V/McuTeB/brw/W2vmPlEbOlmGWs6WNPM5UKmYE1NMYfHvOfd8FLf9Lon+gI0gQA8/xK8sA+uWgQ3FZr5T2RyBYMiwV9hw4C2dmj2mbP4JiWZwZKsTHNoli0Op6AWERGRc1IQREZluj2dlVe+j1+5G+jq7uW//7hfQRBJOHOy5py/0AXq9nebuT56uklP7h9mYU+1hwIhwSlyf+H+BcnWZK6dey0n206SakslIyVjwus4ISIFQmbOMPNunDhlDnOwWCD1IjNPCAwNXPjbAT/4WyJfo9kF3Ucj7LCYARFbhhkgsaaDdcDjYMBkcPDEMqhTpc1mJk69tgAOvG4GQ97xmPsMA/YfMJeFC+B9TrjysgsLvsjYWSyQkmwuYL4+7R3gawUMsCWZ/+cyM83gyQXYsGED69evJz8/f8i+hoYGvF4vBw8eZP/+/ZSUlFBUVHRB1xMREZHIFASRUfvy6hU803QW64kb6D51Oe+8A5dcEu1aiUwtKbYULsm5hLMdZznVcYr0pPRQr5CslCzWL13Plue2cKTlCLMyZnH5tMsByEjJ4EjrERbmLiTJGqdNfKRASI7d/Jb+yFHzm3qr1QxGRBLo6n9sS4uwf7hkpQYE2s1lpHmfZ95q9kwZ6OwfIdBt1u/Ka+Fdf2dOpfunRnjlVbOMxWIOx3j7IMydbQZDCq4y70uiZ3BQBMyhMt4Ws6fIKHk8Hmpra8nOzqaxsZE1a9YMKdPQ0EB+fj4OhwMAn89HcXExHo+H0tLSsd6JiIiIDCNO3yFLNF02bxqfKyhl9/EOfO2d7KpP44tfUJdhkfFmsViYnjGdzJRMjrQeoaenh7Rk80N9RnIG669fz86Xd3LbotvITs0GwGqxkmJL4bDvMAtyFmCJ1x4GkQIhmX1T6B4+an4rP9xQhTkfN4e7BNrN3hqDZV4FvV4zGOLv6PvZDsYYZryKdP7Wl8zzA2QuBlsmLHSYy9HD8Kv/Awd8gA0sNnjzdXjrrzDNDoVL4NqrIDXb7GVi6xumY003pwke3OtEJl5ycvgMQKPgcDiorKwEYPfu3RHLeL3eUAAEwG63s3HjRsrKyhQEEZGE5fP5WLZsGSUlJaF2dKCKigrq6+vJyckJbXM6nWzcuBG73T6kfG1tLfX19Xg8HgoKCti4cWPEnnmSGBQEkTG5+WZobEzF29bBU40n6J77+oiGTqdZM8izXx+27Uz3cTztr47oullJuVyeFT785ninh2Odb4/o+GnJs1mQmRe27XD7G5zqPjKi42elzmde+qVh295uexlvz8iSHc5Lv4xZqeF5AF5veZE2/8hmjrgk4ypyU2aGbXvZ9wzdgc4RHX9F1jVkJoX/YdjXvJcAI5sl4V3Zy0gd8M16b6AXt69xRMcCLMlZjnXAh7gOfxuvtjw/4uMTVWpSKgtzFnK64zRnO86SnpyO1WIlPTmdT1/96SHlk6xJdPZ0cqLtxKQM3ZkwoUDIITOQYLFAWhosmG9OgRoIDP/h1JoM1pzI+6Z/KPJ2ozc8KBIKkrQPs73DDNAM5h/Q08Q6qCfKzEwoWgjOLnj+JOw7Bf6+37+TZ+DXb8OTDTAvE5KtkGTt/5lkheQUSE2HlDRISYe5RZCW3fdBPQlsFuh9x9yWmg1pF53jCZZo83g8VFRUUFBQEPZmPPjY7XbrTbqIJKT6+nry8/Opr6+PGAQBuPPOOykvLw+tV1RUUFxczJ49e8LKbdiwgZaWFmpqanA4HLhcLtauXcv999+voYcJSkEQGZMFCyA/30rzX9M41vkqNb/fNqLj0vxzWXzqqrBtp9Ne5p2cn4zo+OzuRVx59vKwbccyn+NI1uMjOn565zIWeh1h2zzZezmZ8YcRHT+nbSUXt4Z/qHwzdw/NqX8b0fHzW4qZ3R6etPLA9MdpS35rRMdf2vz3TOsK/9DVNPMXdNtOj+j4RWf+gaye8OfvhTk/xmDobCSR5J+6hFR/fxCm19LOS7O3jehYgGuPX42V/g+tHUnHeHnGyI4f2/ewU4fFYmFmxkyyUrI42mLmAgn2ChnMfdLNsZZj3HjxjaQnp2NPHfqNSNxIzgDmhwdCkpNhQV+PkK7uC87VEGJJgqRsIHvs5zAMmLGyP0gyOAgS6AQsYE+FD8yHG+bC307Ciyehs+/3sKMH3mge2fWy2szeJKHzd0Gbu+9+UiD32v4ASUoy0GoGSVJS+5Z0SO0LqKSkQ2oGpGQM+JkJaVlmktDgeYI9I5L7ho3Ea2+jGOBwOLj77rvDeoKA+Q1ocL+ISCJqaGgI9YpzuVw4nc7zHlNZWcnKlSvDyjc0NLB///6wwIjT6WTbtm2sXbsWp9MZseeITG1xEwR59dVXw7p1G8NMf7hlyxbWr18fWi8tLWXLli0Ry15//fU8/3z/t9DPPvss119//ZByzz33HEuXLg2tv/vd7+a5556LeM577rmH2tra0PrmzZu55557IpYd3E093u5p1SrY15TKiXYLgb66v/Pg1rBjL9mwLmw9EDDo9Qd4+80f8uJzZaHtWfmLmLFi0AwPfY7+/Nd0n+z/kD9tpZPcaf1j8P0Bg4Bh0HXiFMe29wdDUmbN4KK7PhJ2Ln/f9f/27L2889a20PbpH3KSXRDeQyRo4D29xQ/56CfDky0G+q7f0nSAM79znfOegtf//ZM34W3+W2j73DW3kzo7vIcHMOSeTk7/MytWhL9OAcO8/umn9tLq7u9RE+meegMBev0BHvtF+Ie8wa9T0OB76rmsneuv7/+/0GsJhF77wa9TpHvq9QfwnX2OP+x5X2hbpNcpaOA9LVq0KGKZRJOWlMYluZdwqu0U3i5vqFdI0CunXuGRlx7BH/DTG+hluWM5qbZUUpMmdkabCRUpEGKzwfx5cOw4dHSaU+nGAotlaI6QgVLnwsKvmcGQYA+ThR2wygvPN0FjE/hazF4phr9v6XvM4ClcLQyZ6d4YENC0JEFvr7kEO6f0nIbOSElhz3lThIbuWGzmeS02c3hO6jwzV0swQJJkgC1g9lZJzTR7rGSkmzldcuyQkwO5fY/HOMRkqhn4LWbQrl27yM/P1xtzEUlIPp+PpqYmnE4nhYWFNDQ0jCgIAuYQw4FDZKqrq9m4ceOQcvn5+RQUFFBXV6ehhwkoboIgEnsWL4YvfsHGL/+Yi7vDDKi8Q3gQ5JpZS8PWU8lm0cXtdHd38+KAz/Iz0mcNKRvUnLyHbvo/XF++qJOLLm4Prc82pjPbWEpz90GO0R8wSE/OGHLOXBw4Frbzzlu9vDOg88X8rEtYOMz1B99TwbXtYev2wCV4SeLtrE7O0B8wiHRPF82exkxLO0/vDeBt7t++KHcxubMWDL33wfeUZh1y/eRAPt2082L6q7TSHwSJdE+XzU4i09LOY78Iv85wz/3ge5o+3Qi7fq/Rhc0wjx38OkW6p/zZnRw/0skfBvRSjPQ6BQ2+JzFZLVZmZ80mOzWbo63mB9q0JDMI8JdDf8EfMD8sP/XWU/QEelhhXcHC3IVh0+3GnVAgZECOEKsVLpoLJ0+Zs3lkDpMoNdZYLH0zzKT3d3HKAG5eCh/sBc9hc9rWnl7o6elfuruhqx26O/p/Yu8v190NnUlgnQ09XRCIEGQIjKzXVzgD6O0LxgzYbM0wgyADAy3dp6DroLk/aQakD8qa3fE29Db33XMyZKeYS1YK2FP617NTIPM8vUzsy2D6B8O3nfy1OQUywKzbzfwvAx3a0n/9GOXz+dixYwc7d+48Z7mPfeyOsPWSNSXcFSHxqogIhx6HI78Zn3PlXg2Lvjz8+efdBvNvv6BL1NfXU1hYCMCtt97Kpk2bhh0SE+R2u9m0aROrV68OG0bo8XiGHVa4ePFi9u3bd0F1lfikIIhckGXLYNmyPMDscWDZFN4oPnbflyIet2WLnfrH+tdvuW4xW4Ype/0vf8DzRzyh9S/9/XSuv372gBKzgfebvVu+//9CWy+fO2vY6594O52/Pdu//uUPf2D4HjuD7ukbX589qMRH++4pifW//vmI7qn+0SSODUhD8p93f3L4HjsD7skxNyPC9dcCcM/Bv1H77N4R3dM3/1f4+vCvU/g93XCdPcL1zWMHv07D39MRfvhf/evnep0G35OES09OZ2HuQk60naClq4WM5Aw+veTTPPLSI7x86mUA/vjOH+n2d/PRvI8y3z4/fhOlghkIscyH9kNmAMFiNT8kz55l9io42wwBA9LOkTQ11iUlwaXjNN2WYUB3T3ggpbMFOs9AZxt0t5nBlK62AYGVTvNxd1ff4y7o6YaegLn09i09ATAiDEMyBvRWiTg7kQHBHETtXeZyfJj6Wy1mcCQ7BbKTBwVJkiG1a+gxRm9/glsjQq6jgftjVFlZWWjc+rn86le/nKQaiYhMroaGBkpKSgAoKioadkjMww8/zI4dOwAzgDw4iarHY74vHdgzZKDc3Nxhk1bL1KYgiIjIKA0zcm1cjTRWYbVYmZs1F3uKnWNtx7Bg4bPXfJaf7fsZTSeaALN3SI+/h89c/RlmZc6awFpPgqR0yBgUCAGYlgu5OdDZaQZD2jvMJzE1JXGnnQ3ef1jOlFnAZaM7jxHoG77Tt/g7zJ/WdEi/NLwnypkXoPkF6OqAtMWQepXZq8XrhWYfHPbCqTZo7Tn/L1LAAF+XuUSScgZy3H3Da/qG2RhvmNuzUyClFTICcfX6V1VVUVpaOuJu3yIiU43P58PlcrF1a39PbKfTGXFIzN133x02pLC2tpZly5axc+dOHA5HKJg8XG8Qj8ej3EsJymIMl4gixuTl5XHgwIFoV0NEoiRW2oC8vDyeeGLi6zFjBkybdv5yA/kDfk60naC1u5VUWyrb3dt56fhLof1Xz7maLy79Ymg63bjW2zE0EDKQ329++D7bbA7VsNnGL4GqjF2gBzDMmX1a28zAiLdvGfy4o+M8J7MMfe2NAP1jdqxmACQ7ywyQ5djBntmfkyT4My11xFHHvNs+N+Z2aOXKlVRWVg4b4Kirq8PhcIwoAJKXl8eBXT8/bzkRmXrybr0rJt4PTZS6ujqqq6uH9N7wer0888wzofWKigqys7OH5FWqqKigubmZBx98EDDb3pKSkoh5P4qLiyksLIyYm0mmNvUEEREZpYwIs6KOt9Onzc/uM2eOvFeIzWrjouyLaOtu41jrMT6x+BPYLDZeOPYCAC8df4nv/vW7/EPhP4RyiMSt4XqEBNlsYM82l54eM9mor9UMjqQkKSlntFj7nncrkJsGuTOGL9vVDb5gcKSvJ8ngdf+gZLFDgiJG32sfntA6TErygIStg3/azf9DSRP7dsnlcg2ZJjfYjVvfUopIIqmrq4s4dW1eXt6IZolxOByh9hPMCSWqq6uHBEFcLhdut5tt27aNW90lfigIIiISgzIyoKXFDITMnTu6WUgzUzJZmLuQ463HuX3R7dgsNp49aibB2XdiH9/+y7cpd5bHd6JUOH8gJCg5GWZMh+nToKsLmr3mcJmAAempcTVcIqGkpsCsmeYSiWGYvUm8PvM1Df70tfSvt7ad/zrdPWZy3ZOnIu+3WCAr0wyITACXy4XP5xvyxn379u36dlJEEorP58Ptdg8JgIA5JGb79u3nDIK4XC62bNlCTU1NaFtJSQkNDQ0UFxeH8i25XC7KysqorKzULFwJSkEQEZEYlZZmjgg4fBjmzRvdZ3Wb1cY8+zxaulr48KIPgwWePfIsFouFS3Iu4XjrcebZ501c5SfLwECIJcnsaTBcMMRiMZ/UuWnmcIz2DmhuNn9aLZA68mEREgMsFnOoS3aWOV1yJMFeQAODJN6W/p4kXp9Z5lwMA1pazWWUfD4fmzdvpqWlBY/HQ3V1NYWFhSxfvhyn04nP52PdushTlDscDgVBRCSh1NXVsWrVqoj7SkpKKCsrC9u2Y8cOGhsb8fl8ANjtdmpqaoYESrZu3UptbS1lZWV4PB4KCgoilpPEoZwgIhIXYqUNyMvL4667JrYemZmwfDksXGiud3WZn/fmzx9br/zeQC9HW47y3y//NxdlXsT7Fr6P9p52pqVNIyM5gzp3HQBJ1iRsFhs2q40kSxJWq9Vc79tms9pY7lhOenJ62LmfP/o8VotZ1mqxYrOaPyNtc9gdYTPU+AN+Wrtbw8oEj7FarCOfzcbfDb2tfUk7u6FvmmAsFjM4YkkaPsDR29ufn6K3x+w5kqLhMgnBMKCjc1CQZMBPr88MfvS9Vcp7vilm2iHlBBFJTFM9J4jIZFBPEBGRUXr99Ym/xksvwYc+BCtXmh0UenrA4zEDIaNNZ5FkTcKR4+Dvr/17TradpDfQS0ZyBs2dzTR3NvPkG08SMAIYDI2JW7AMXGF6+nSmpU0zc1Jiob2nnZq/1GCxWLBgCQUthntcU1RDkjUpVP5wy2H+/U//TvAyA69nwdIfTLHayE7J5psrvhlWv9dOv8ZPXvqJGTAZWFcMc7rWgN+cEjXQy9zM6axf8imwGIANrMk8f8LN468/2X+Yv9fMRdHVbX7wTbKFuuBcme3gU5euDrv+7489yx9OPD+i1+Hd09/F7fPfF7btN4f+yHNnXhnR8e+f/W4+MHdp2Lafv9XAqy0HR3T8hy9+L0tnLA7b9sibuzjWeYokSxLJ1iRSrEkkWZNI7lsfuFw3LY856eE5PF5v8WAYRnjZQcdah+uZEwssFshIN5d5cyOX8fvN3iReH3zpXya3fiIiIjLuFAQREYlBhgFPPQUHDsBdd8GsWeZnsXfeAYfDDIyMVk5aDhnJGRxrPUZHTwfpyekYhoHFYsFmGVl+kPSkdFKT+i/e0dMRCnAEjIAZexgQTAk97vtxvPW4+cBibjvSeoSW7pawAMbA4y1mtAULFjp7O3m7+e1QUMWKlYPeg7x19q0R9RjpsKZywpIGgR4s/i4I+DjsfYc3zwajWsHZRsxrEvBDazcWvx8DC6nYONNxJuycR9tP8FbLoRE9d7NTcoYcf7jt+IiPvzLz4iHHe9qOjvj4Y20nOJMR/kH/gPctPO3HI5S2DPjXfJRlTSPDmhrqUGPBwvcO7KC5e9AwkQivhQX4l8WfY37mnNA2wzDY8GzViOoO8M1rv4w9OTO0frbbR8WL3x/x8Q8t++ew9Xdaj1L98k9GdGxuSjb3X/PFEV9LREREYpeCICIioxRhlrVxEwjA//wPvPmmuX7oENTUwG23wQ03mCktDh6Eiy8e2yw1ybZk5tvn4+vycbL9JIZhcNuVt+E3/BiGQcAIhB77DT8BI2BuC/gxMAgEArR3t4fO1+PvYfHMxea+CMcGHwePD+us0RfcyErJCl174DkCRsAs1xdYsWDB3zfMxTAMDAzautvo9neP6N47ejo40zlwlpBUvJYMOi0p5tSqRgBLoNfsQYJhxmlSDDCsWHoNOrvb8Z09gWG1hYbLtHW109078ut7O3xh29q7O0Z8fHt3+5DjO7o7R3x8W3fbkONbulvpCh1vEIpOhbFgsUBbdytnOs709eLpP76zt3NE1z/dfppUo79XiGEYdPScbxrcfqdaT9GV3P9/z9vdOuLjLVg42XpySH1Genwy1iHHi4iISHyKqSCIz+ejurqaysrKaFdFRGRYV145sedftAj++EdoaDB7f3R3w86d8Mor8IlPmDlDDh+GOXNgLEnNLRYLOWk52FPtGBhcPv3yyOU4f+8KgCVzlpz3esO5YvoV3HTJTRH3BQMd/kB/MGZgPhIAR46DZRcviziUZ7AkaxKzM2cPOf5Dl30o0sXB6IVAL/R2QKCTVKuVGdb0/iSZASsll6zk1kvef95rA2QkpZGbkh22bV3WDEp6Vw9zRLispHTsKVlh276SeRdd/vMk9uyTk5JFZlL487exYC3tvZ10Gd10B3rp9vfQHeihJ9BDtz+43k1XoJtLsx1kpmQQGPCaLMyaR3tvBz2BXnoCfnoMP719j/1GLz19AS0RERGRWBFTQZD77rsv2lUQEYk6iwXe/3644grYvh2O941W2L8f/vM/4ZOfhKuuMrf39sL06WO9jmXEgY5oCNbPahs+p0RaUhoXZV805mtkpmSSmZJ5/oJBAT/k9kJvJ7Q3k33mONmtrWYXnvQ0sCaZywjzYOSmZA8JjIzGjNTcMR8LhA1PGYuvX93fLSoYtApghPXsGRgEGRgQMwyD2sKKUV1v4PEzjZmjOn5wMG5m5kxqZ+WP+XgRERGJTzETBHG5XDgcDjweT7SrIiISEy6+GDZsgF27YO9ec1trK2zdCoWF8OEPw9mzZiBk1izN7joprDbABrZUSM2BaZeY3XXafHDmJLR6gU5I6QuCGEbfDDXJfcdOXaGgFRDDsTURERFJcDERBAnO7bxgwQIFQUQk5nWPLAXDBUlONj87JyfDRz8K73oX7NgBLX0pLRobzVlq7rrLXO/thblzQxOZyGSy2cA+zVx6e81I1dmz0N1l/pW1WcDfCYFOM99IKDBi61v0oomIiIhMlpgIgtTX11NSUkJdXV20qyIicl7ZYx+9MCKBgBnsMAwzEarFAnl58NWvwqOPQlOTWe7kSfjOd+CWW8yeIYcPw7x55mdyiZKkJMjNNZeuLvD5zMXIhLTpZpRqUK4RAl19OUgG9BqxWNW1R0RERGQCRD0I4nK5WL16ZEnh8vLywta/8pWvcO+9905EtUREhjVjxuRco6UFTp/uD4ZkZsJnPwvPPgu//rX5GTsQMBOoHjgAxcXmyIz5883P4hJlqanmOKWZM6GzE86cgY4OM7iRmgppA5KUBvxg9IC/B/wd5uNA3+w4Q5K+9gVMQvssA34O2Gb0/bRE2IbRH2QxBhwbirtYBmyzRNg2sFzfTwVtREREJA5E9W2yz+fD5/NhH+H0BgcOHJjgGomIxAarFXJyzF4nPp8ZDLFYzGDIsmVw2WVm0tR33jHLv/UWfPe75lS6fj84HJCSEt17kD4WC6Snm0le/H5oazMDIr29Zred1NQBuUbSgFF2NQoLiATXR7qN/m3B3ij0LUaAIQEXwwAC/WUGbjPoO2bANQaeX0RE5BwqKiqor68nJycntM3pdLJx48aInxdra2upr6/H4/FQUFDAxo0byc8fecJrSVxRDYIEh8GIiESDy+XC7XYD0NzczIIFC2KuTbJazZEVdjt4veZnZ4vF7CnyxS/C734He/aYnQY6O+G//9vsFVJUZE61m55+3kvIZLLZzBfTboeeHjPC5fWaL2BysrmMVlhPDUjIpKQDAzoT/VhERCbMnXfeSXl5eWi9oqKC4uJi9uzZE1Zuw4YNtLS0UFNTg8PhwOVysXbtWu6//36Kioomu9oSZ6IWBHG73TidzmhdXkQSnMvlAqC0tDRsW21tbdi2WGG1wrRpZu8Qr9fsGWK1wsqVZrBj+3Y4dcosu2+f2UPkIx+Bm24yh9FIDEpONqNZ06ebY5uam81eIoYRHgwZ/AF8tOuR9oUNhelbP9e+4dbHWjbS44HlBz8eKFIZiyXyY6v1/GVG+1hERCZNZWUlK1euxOVyhT47NjQ0sH///rDAiNPpZNu2baxduxan0znikQaSmKIWBPF4POzbty9sW2NjIz6fj6qqKpYsWaIonohMmO3bt/Pggw+GbXM6nWzfvj1KNRqZYDDEbjc/M589C7NnQ1kZ/OY38PTTZjmfDx55BF59FT796cnJYyJjFBznNHeu2SOkvd1cgvsGfhAf7udIygz8Ofjxha6PtawCCyIich5erzdsiEx1dTUbN24cUi4/P5+CggLq6upi8gstiR1RC4IUFRUNCXLU1tayb9++sC5QIiIT4dChQ2HfKsQbm80MbOTm9gdDbrvNnEr3v//b7FAA8Ne/whtvQGkpFBREs8YyIlYrZGWZi4iIyIV68kl4/HGzx+FkS02F22+Hm28e0+Fut5tNmzaxevXqsFwfHo9n2NwfixcvHvJFu8hg1mhXYKDm5uZoV0FEEsQ999zDunXrwqbmrq2tZf369VGs1egFgyGXXmomUV24EL78ZTMYEnT6NPzf/2sGR0ITjoiIiMjU9+ST0QmAgHndJ58c1SEPP/wwy5YtY9myZRQXF1NQUEBlZWVov8fjAQjrGTJQbm4u+/fvH3udJSHERBDE4/FQVVXF7t27aWxspKqqKpSsUERkIhQVFVFZWRlKuFVbW0tJScmIsorn5eWFLQ899NAk1PjcbDZzJtZLLzVHVZSUmD1DBqaWePRRuP9+M7mqiIiIJICbbzZ7ZERDauqoe4HcfffdPPPMMzzzzDMcOHAAh8PBsmXLQsEPh8MB9AdDBvN4PKEyIsOJ6uwwQQ6Hg/Lycg2DEZFJ5XQ6ufvuu2lsbKS6uhq73T6i2WFiebrupCQzR8j06WYSVYcDfvlLOHbMDJQcOAD//M+wfr051a6IiIhMYTffPObhKLGgtLQ09IV5MJdbcDaYSF9cNTU1UVhYONnVlDgTE0EQEZHJFpweNxh8rauro6KiArfbHdbtMl4NDIZccgk89hi4XGbKifZ2+Na34MYbYcECM1iSm9v/0243jxcRERGJNofDEdbzo7S0lOrq6iHJT4Pv7bZt2zbJNZR4o7e5IhK3PB4PZWVlIy5fWFgYCnoMnh2mpKQEp9NJcXHxiIfFxIOkJJg/H+65B667Dn70I2hpMSflcLngL38xAyNWq9lTJDijqN1uzkITDI4EAySDgyXWmBhUKSIiIlORy+Viy5Yt1NTUhLaVlJTQ0NBAcXExNTU1oZ4hZWVlVFZWanpcOS+9fRWRuOVwONi5c+eIl2AAxOVysWTJkojnu+eee4YdZwrw6quvYrFYQstwtmzZElbunnvuGbbs9ddfH1b2ueeei1juueeeCyt3/fXXD3vOe+65J6zs1q1bcDqhqgqWLgW/vz/48dOfWvjJTyxs22bhhz+00NoKR4/C/v3Q2Ai//S3s3Alf+9oWrr3WwqWXWpg2zUJ+/j2Ul8M3vwnf/a45Je/jj8Of/gSLF4ff05///BxtbQxZ/vzn8Hu67rrrI5Zra4N168Lv6aGHtgxbdmA5i8UybLmHHgp/ndatu2fYstddp3uK1Xvq7Jyc3ycREZlYO3bsoLi4mJUrV7Jy5Uqqq6upqakZMpvf1q1bWb16NWVlZSxbtoza2lpqampGNKxZRD1BRCThOBwO9u7dG3Gf3W6f0gm17Hb46lfhtdfg5ZfNJKk//Wl4GasVDMN8bBjQ2ws9PdDdHV6uuxsOHzYDJjZbeE+SkyfDy/7nf8KsWUPrM7jcoUNm/SL5y1/C13/xC2hqOvf9Bg13zpdfHnqN4coeOhS+rnuKfM5o3dO73232eBIRkfhUWVk5qiHJpaWlQ4bEiIyEgiAiknCCY0vdbnfYsBefz8fBgwcT4luEK680F4B168L3/ehH0NwMZ8+aS/Dx44/DwC/VLRYzSBIMkAQDJxaLGTgZqL3dHIYzWHt7+LrfH7kcmNcZqLNz+LKDDVduYA+C4DWGK+v3h69P1j0Fn9fBfL7I2zs6wte7u4cvO/ie2toiTyLQ1jb0uOHOOThY1tExfNnBLuSenn4aPvWpkV1HREREEpeCICKSkB588EFqa2vZtWsXubm5gNkLRLNUmdPqzpo1tEeAzwc/+1n/+gc+YA6vaW42e5QEgyXNzWa+kebm/rKZmWYvlMG6usLXbbbI5YL1Gigtbfiygw1XLi1t6DWGK2uzha9Pxj0NFwCB4c+Znh6+npLSX3bw+QbfU0YGZGcPPefgYJHNFrlc8HoDpaUNX3aw4coNfp1SUsLLWq3w/veP/DoiIiKSuCyGca63WLEjLy8vpqelFJGJFSttQKzUQ0QmX6z8/ufl5XFg18+jXQ0RiYK8W++KiXZIJJ4pMaqIiIiIiIiIJAQFQUREREREREQkISgIIiIiIiIiIiIJQUEQEREREREREUkICoKIiIiIiIiISEJQEEREREREREREEoKCICIiIiIiIiKSEJKiXQERERERERGRsaqoqADA4XBQWloa5drIRKirq8PtdgOwceNG7Hb7mM+lIIiIiIiIiIjEtcrKymhXQSZQSUkJYAZDvF7vBQVBNBxGREREREREoqqhoYGVK1dSVVVFcXExFRUVVFRUUFxczMqVK0O9PSaL2+3G5XJN+HV8Ph8NDQ2ha1ZUVLBs2TLy8vIoLi4O7aurq2PlypWsXLmSDRs2hHpFBE1GXS/EaOpXV1c3gTVREERERERERESirKioiFWrVlFeXk5hYSGVlZVUVlayc+dOSktL8Xg8Iz6Xx+OhoqKCvLw8qqqqhuyvqqpi5cqV1NbWnvM8Pp9vxNd0uVwUFxef95yD1dXVUVRUBEB+fj6VlZXcf//9oesH95WUlOBwONi6dSsPPvgg+fn5YXXcu3fvqOs8GcZSv4KCggkNhCgIEkceeuihaFchqnT/iX3/EhsS/f+h7j+x719iw8+3b492FaJK95/Y9z8WFoslbBnOli1bwsrdc889w5a9/vrrw8o+99xz41LX3NzcIdvcbncoADBSDoeDyspK7HY7S5YsiXidYHBlvDidTlavXj2qY+rq6nA6nUO2FxUV4XQ68Xg8oSBObW0tRUVFQ56HYMBnx44dbNiwIeZ6hIylfvn5+bjd7gkL6CgIEke+853vRLsKUaX7T+z7l9iQ6P8Pdf+Jff8SG+q2T2w36Vin+0/s+09EmzdvBsaW82P16tVDehTU1dVRUlJyQTklxsvevXtDPToGq6mpAeDhhx+moaGBffv2hfJiDJSfn09RURGFhYUsWbIk1HMkVoy1fiUlJdTX109InRQEERERERERkZji8/kiDmUZjZKSElwuV6hHQUNDA06nc9wCIME6NjQ00NDQMGTIjs/no7a2FpfLFTa7CZg9JCL1fAmy2+2hwE9ZWRnl5eXDlvV4PDz44IOj6i0zHurq6iguLg4bdlRbW8uyZcvYsGFD6PkYS/3y8/NDQ2jGm2aHERERERERkZjR0tJCdXU1zc3N5wwUnE9+fj4Oh4O6urrQ4/EMFKxdu5Zt27aFgiqDh3qsXbuWnTt3htY3bNhAeXk5DocDt9s9bC+QoJKSklDwpKGhYdjhO8EeImPtBTKSYNOaNWuGPHclJSWsXr2aFStWsH//fgCam5upqakJG+Yz1vodOnRoVOVHSkEQERERERGRKcowjBGVu+eee86ZB2Sg8coBMpzs7GzKy8tD+SQuRElJCVu2bKGmpiYs6ODz+di8efM5e1ici8vlIicnJ6xXycAgQUNDA/Pnzw87Zvny5dTW1lJZWYnH4zlvEMTtdpOTkwNAdXV1xJwg42GszwGYPVa2bdtGcXEx69atw+l0RsxzMhY5OTn4fL5xH7oUV0GQvLy8aFch6hL9OdD9J/b9x4pEfx10/7p/ib68W++KdhWiKE33r/uPdiVkkjgcjrA8GGP5QGy32ykoKBjywdzlckVMmjpSbreb7OzsYffv27dvSC+WnJyc0BCR5ubmUIBjONXV1dTU1FBXV0d1dTVlZWVhPUtG6kIDPueTn5/P3XffzcMPPzymZLPD1S87Oxuv15u4QZADBw5EuwoiImqLRCTq1A6JSCIJDqFwuVx4PJ6IyUHPZe/evUMCIG63my1btnDPPfeMuadBfn7+ORN3LliwYEhOC6/XG+rJkZubi9frHfb4qqoqSktLsdvtlJaWUl9fj9vtpra2dtSBhvMFfMY6HCbI5/Oxf/9+nE4nZWVlPPXUU6N6ToerX0tLy4T0fImbIIiIiIiIiIhMXc3NzRG3u91uqquruf/++0d9zsbGRtavXx+2LT8/n5ycnAuaSSUYWBkYRNm3b1/oQ3swn8dADQ0NoWSndrt9SCLVILfbTUtLS1jwpqamhpUrV0YcFuPz+aivr8fhcFBQUBAWgBhJwOdCe4jcd999odlsiouLh+RCudD6jTcFQURERERERCSqGhoa2L17N2AGLoK9E/bv309TUxOFhYWj6hXgdrvZtWsXPp8Pl8uFw+EIfcD2+Xzj0sMgOFQlmNsjNzeX3bt3s3z5cpxOJzU1NdTW1uJwOPB4PJSWloau63Q6qa2tHfIc7Nq1i8bGRhwOBw0NDaFATUNDQ+g8wdwbGzduxG63s3nzZtavX4/X6x0yfGQ8Aj7Dqa2tpb6+Pix44XA4cLlcrFu3jpKSEoqKiqJWv+FYjJFmyhGRcVdVVRVqJEVEokHtkIhEm9ohuVAVFRWhHhYjERxa43Q6hyQ3DXK73Xg8ngn9cL5hwwYefPDBCz6P2+1m8+bN5ObmDnkefD4f1dXVo3p+xttY6ud2u2lqahoy/Kmurg6n03lBQayY6AkSvEGHwxGKIo2kERzrcbGmoaEBr9fLwYMH2b9/fyhidi4ej4eysjJKSkooKCjA4/HgcrkoKSk5b5bhWHIh9zEVXv/9+/fz8MMPR9zncDjYs2dPxH3x/vpv2LCB9evXR6xrtNoDtUNqh9QODaV2SO3QZFI7pHYoErVDidkeTLRgj4rx6hEyVsuXLx/RVLnnEpyW98EHH4yY26OpqYn8/Hw8Hs+wAZ+JNNb61dXVsXHjxomplBFlBw8eNO69996wbffee6/R1NQ0IcfFmvr6euPgwYOhda/Xa6xYscLYsmXLOY87ePCgsWLFCmPp0qXGokWLjLVr18bdvRvG2O9jqrz+mzZtMpqamoyDBw+GLdu3bz/nvcTj63/w4EFj06ZNxgMPPGAsXbrU2Lt3b8Qy0WgP1A6pHVI7pHZoYBm1Q5NP7ZDaIbVD4WUSuT0Yi02bNo37OZuamoz6+vpxP+9gDzzwwAUdf/DgQWPv3r1GU1NT2P+nFStWhPaf73dpIo2lfud67rdv3x7292Isoh4E2bRp05Bf/qamJmPt2rUTclys2b59+5Bt9fX1xqJFi8553MGDB6dEgzbW+5gqr/+5frnPJd5f/xUrVkT8ox+t9kDtkNohtUNDqR3qp3Zo4qkdUjsUidqhfonUHozFRAVBIr0+483r9U5IsOVCAwUT7Vz1O9fv/ngEQawT079k5IJZYgcKJlOZiONiicfjoaKiArfbHbY92B1q8HbpNxVefyBiN9/a2tpRT/01VUSrPVA7pHZoLKbC6w9qhwZTOzT51A6N3VR4/UHt0GCJ3B6MlcvloqKiYkii0QuRn58/KUOJ7Hb7hOQd8fl8437O8XSu+kX63a+rq6OiomLIjDtjEdUgyHDjsILjgIb7ozfW42KNw+Hg7rvvHnIfwf8QIxmf5vP5Qkl74tlo7mOqvP6RuFyuUTW2U+X1h+i1B2qH1A4FqR0yqR1SOzTZ1A71UztkUjuUuO3BWO3Zs4fKykpKS0ujXZWYEet5cUZbv5KSEiorK9m5c+cF53GJahDkfNEpr9c7rsfFovLy8iHJaXbt2kV+fv55k9bs2rUrLPHRhg0bYj7iF8lo72Mqvf6DeTyeETcIU+X1D4pWe6B2SO0QqB0aSO3Q8NQOTRy1Q2qHBlI7NLxEaA9EJlrUh8NIOJ/Px44dO6ipqTlnOYfDwfr163E6ndjtdvLz81m+fDn33XffJNV0fEyV+xgPbrd7xFFNPW8ykdQOxed9jAe1QxIr1A7F532MB7VDIjLRFASJMWVlZdTU1Iyo8R/8zYjT6WT37t1xF/2eKvdxoaqrqykoKBhxeT1vMlHUDsXvfVwotUMSK9QOxe99XCi1QyIy0aIaBAk2WsM1Ujk5OeN6XKyrqqqitLR0RGMgIyWECd53PI2HHMt9TMXX3+fz0dTUNOJ5u6fK6z9QtNoDtUPh1A6Z1A6d31R5/QdSOxQb1A6Z1A6d31R5/QdSeyAy8aIaBHE4HNjt9iFj1IK/vMONBRzrcbGsrq6O5cuXj+gPfjCL+nCN+4UmipksY72Pqfj6u1yuEb9uU+X1Hyxa7YHaoX5qh4ZSOxTZVHn9B1M7FH1qh4ZSOxTZVHn9B1N7IDLxoj4cprCwcEi24qampvP+8RvrcbHI5XJRUFAQVnePx3PORr2ysnJI415fXz+iBGKx4kLuYyq9/gD79u0bcdmp8vpHEq32QO2Q2iG1Q2qHgtQORY/aIbVDaodMag9EJpgRZQcPHjTuuOOOsG333nuv0dTUFFr3er3GHXfcYdTX14/quHiwd+9eo76+3jh48GDY8sADD4TKRLr/4DGDy8Tb/Y/kPqby6x907733GmvXro24byq+/itWrDD27t07ZHu02gO1Q2qH1A6pHQpSOxQdaofUDhmG2qGgRG8PRCZaUrSDMA6Hg40bN1JbW4vD4cDj8bBmzZohXbY8Hk9Y966RHhfLfD4f69ati7jP4XBQXl4eWh98/0VFRTQ0NNDQ0BDaP9IEYrFkpPcxFV//gRwOBy0tLcPunwqvv8/nY/PmzbS0tODxeKiurqawsDCs23O02gO1Q2qH1A6pHVI7FD1qh9QOBakdim47JJIoLIZhGNGuhIiIiIiIiIjIRIt6ThARERERERERkcmgIIiIiIiIiIiIJAQFQUREREREREQkISgIIiIiIiIiIiIJQUEQEREREREREUkICoKIiIiIiIiISEJQEEQmRV1dXbSrENP0/IhMPP2enZueH5GJp9+zc9PzIyKTQUEQmXBVVVWsXr0aAJ/Px8qVK1m5ciV5eXmhZeXKlaE/fMEywX3Lli1j5cqVuFyuaN7GhFq9ejVVVVXRrobIlKV26PzUDolMLLVD56d2SEQmg8UwDCPalZCpy+Vy4fP5KCoqGrJv3bp1oT/kW7duxel0hu2vqKgAoLKycuIrGgPq6upwOBxDngcRuTBqh0ZO7ZDIxFA7NHJqh0RkoqkniIzaaL6B2L59e8Q/+AA1NTXY7XYAysrK8Pl8YddwuVwJ8wcfoKSkhO3bt0e7GiJxQe3QxFA7JDJyaocmhtohEZloCoLIqA38xuJc3G43S5YsGXa/3W7n/vvvB8wun/fdd1/ocVlZGVu3bh2fCseRJUuW4Ha7o10NkZindmjiqB0SGRm1QxNH7ZCITCQFQWTC7Nq1a9hvPYKKiopYtWoVALt376ahoYGysjI2btyIw+GYjGrGFKfTqaRgIuNI7dDoqR0SGV9qh0ZP7ZCITCQFQeKYy+Viw4YNLFu2LJRMq6GhYUi5qqoqVq5cybJly9iwYUNYN8ughoYGiouLQ+cqLi6muLg49A1HQ0NDKDEXmN9+DEzk5fF4hpyzsbFxRH+4v/GNb4R1A83OzqakpGRUz8V4qq2tDUtEtnLlytC+DRs2hCUni/RcXoj8/HyamprG9ZwiE0nt0MRQOyQycmqHJobaIRGZqhQEiVMNDQ2sW7cOh8NBTU0Ne/bswW63s2XLlrBy69atY/fu3VRWVlJTU0NLSwsrVqwI+2PV0NDApk2bKCkpYdu2bezcuZOCggLcbneoXFFREXv27GHPnj2AmZwruL5nz54L+pbCbrdTU1MTWj906NCYzzUeSktL2bNnT+gbmYH3tmbNGvLz83nqqadCz/l4G+83EiITRe3QxFE7JDIyaocmjtohEZmyDIk7Xq/XWLp0qbF9+/Yh271eb2h9+/btxtKlS8O2GYZh3HHHHcamTZtC6w888IBx7733DrnO0qVLjb179w7ZvmjRoojbB9dl7dq1I7qfYPmlS5caixYtMhYtWmRs2bJlxMdOpDvuuMNYtGiR8cADDxgHDx40VqxYMeT5HG9r1641Dh48OKHXELlQaocmj9ohkcjUDk0etUMiMpWoJ0gcCnbJHNxF0m63h0Xi9+7dy+rVq4dE5++55x7q6+tD62vWrKGxsZF169ZRW1sbmsbtmWeeGfP0ZF6vl+zs7BGXv++++7jnnntC3zJUV1dH7FI6Fg0NDRG7xY7Etm3bsNvtPPzww6xbt46tW7de8LcdF1IfkVihdmh01A6JjD+1Q6OjdkhExKQgSBzyeDwj6m65f//+iOXsdntYF0OHw8FTTz2F0+nE5XJRVlZ2zvGyIzGa4+rq6mhpaaG0tDSsG+i6devGdO3Bdu3aNebuqQO7po7Xm5ALqY9IrFA7NDpqh0TGn9qh0VE7JCJiUhAkDjkcjhFNG7Z48eKIf6g8Hk9Y9D44B31paSlbt27lmWeeYc+ePbS0tISmaTsXj8czZIo4u91OS0vLiI6trq4O/WHNz89n48aNoX1VVVXnPcf57N+/n/z8/DEfv3379tA3QOvWrbvgMaojqU9OTs4FXUNkoqkdGh21QyLjT+3Q6KgdEhExKQgSh4J/gGpra8O2ezyesG8Lli9fTl1d3ZA/UrW1taxevTq0HkwENpDD4aC0tJTGxsaIdRh4ztra2iF1cTgceL3e897LunXrqKmpCXsTUlpaGvqj+PDDDw95QxEUfLPidrsjvglyu92hNw2RnoeRqKqqwuFwsHXrVkpKSvB4PJSVlU1ofQa/KROJRWqHTGqHRKJH7ZBJ7ZCIyOgoCBKH7HY7lZWVVFdXU1VVhdvtxuVyDekuWVJSQn5+fmhqN7fbzbp16/B6vaFvF4J8Ph/FxcU0NDTg8XhoaGigoqIi7M1BUH5+PnV1dbjdbhoaGqivrx/TWNmKigqcTmfEYwd2Ay0rKxvyB9Ln84W+kXA4HOzatStiPZcsWcKqVasoKSkJ+0MavN9zzUFfW1vL/v37KS8vB8wM8Pn5+bhcriHfyFxofQbSH3yJB2qH1A6JRJvaIbVDIiJjEu3MrDJ2e/fuNdauXWssWrTIWLFihfHAAw9ELPfAAw8YK1asMJYuXWrce++9Q7J5b9myxdi0aZOxadOmUEbyFStWDJuR/ODBg6Es4UuXLh223KZNm4Zk9a6vrw/VJZj5fMWKFWFlmpqaIpZZsWJFKAt7MHv6vffea9TX1w/7HN17770RM7cPPD7S8zXw2sF7GJyxfcWKFaGM9Bdan4H3HiuZ4EVGQu2Q2iGRaFM7pHZIRGQ0FASRCbN3794J/wO2d+9e44477hj2D+ngNxSDDZ5WL9r12bJli9HU1DSudRJJZGqHRl8ftUMi40vt0Ojro3ZIRCaShsPIhHE6nezbt2/cz+vxeNiwYUPoGiUlJaHEWQMTn/l8vlBXyuHG0cZaffbt23dBSctEJJzaIbVDItGmdkjtkIjEFgVBZEKtWbNm3OeAdzgc3HrrraGxvzk5OaE/lFVVVaHxqXa7nYKCAlwuV8Qp2EY6td5k1aehoYE1a9ZccH1EJJzaIbVDItGmdkjtkIjEDothGEa0KyFTW1VVFevXr5+0BFcul2tEiclqa2spLS2Nifr4fD42b94cSjomIuNL7ZDaIZFoUzukdkhEYoN6gsiEKy8vp76+ftKuN9Kp3ybjDz6MrD719fX6gy8ygdQOqR0SiTa1Q2qHRCQ2qCeITCkul4uCgoKYmVYt1uojIhMv1n7vY60+IjLxYu33PtbqIyKJTUEQEREREREREUkIGg4jIiIiIiIiIglBQRARERERERERSQgKgoiIiIiIiIhIQlAQREREREREREQSgoIgIiIiIiIiIpIQFAQRERERERERkYTw/wPi8Ps+8dn8wQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we demonstrate how to use the plot functionality as used in the paper. Each frame corresponds to one time-slice. As we have only run each model for ten trials, we do not expect stellar results (see final frame).\n",
    "plot_expected_opt_curve_paper(T,\n",
    "    true_objective_values,\n",
    "    exp_per_trial_cost,\n",
    "    exp_optimal_outcome_values_during_trials,\n",
    "    plot_params,  \n",
    "    fig_size = (15,4))"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aaa96eb265c32c2b50eaffc5e0bc3c80016276efa552bb99d04c94d0e8842f70"
  },
  "kernelspec": {
   "display_name": "Python 3.7.11 64-bit ('scibo': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}