MSafeOpt_Memory_Efficient.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfaf034f-15b2-47f7-94d6-31c8456d4b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install trieste==1.2.0, plotly==5.17.0, matplotlib==3.8.0     # using a virtual environment is recommended"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f290f329-83dd-4d7a-82e8-368a3fdc0e65",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' \n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"0\"\n",
    "\n",
    "import tensorflow as tf\n",
    "\n",
    "gpus = tf.config.list_physical_devices('GPU')\n",
    "if gpus:\n",
    "  try:\n",
    "    for gpu in gpus:\n",
    "      tf.config.experimental.set_memory_growth(gpu, True)\n",
    "  except RuntimeError as e:\n",
    "    print(e)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ffb84b01-0335-4a33-b3e8-035fcacbc6ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import gpflow\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import trieste\n",
    "from trieste.acquisition.rule import EfficientGlobalOptimization\n",
    "from trieste.data import Dataset\n",
    "from trieste.models import TrainableModelStack\n",
    "from trieste.models.gpflow import build_gpr, GaussianProcessRegression\n",
    "from trieste.space import Box, SearchSpace\n",
    "from trieste.objectives.utils import mk_observer\n",
    "from trieste.types import TensorType"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aefaf37-004e-41a6-9c7d-c4bc7d96a183",
   "metadata": {},
   "source": [
    "### Objective and safety functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2ba9c98a-1084-4ed5-978f-076157ab6f36",
   "metadata": {},
   "outputs": [],
   "source": [
    "from trieste.objectives import ScaledBranin, Hartmann3\n",
    "hartmann_3 = Hartmann3.objective\n",
    "\n",
    "def obj_hartmann_3(x:TensorType) -> TensorType:\n",
    "    obj_tensor = hartmann_3(x)\n",
    "    obj_tensor = tf.squeeze(obj_tensor, -1)\n",
    "    return obj_tensor\n",
    "\n",
    "def g_syn_3D(x: TensorType) -> TensorType:\n",
    "    x0 = x[..., :1]\n",
    "    x1 = x[..., 1:2]\n",
    "    x2 = x[..., 2:]\n",
    "    safety_tensor = -( x0 + x1**2 + x2**3)\n",
    "    return tf.squeeze(safety_tensor, -1)\n",
    "    \n",
    "obj = obj_hartmann_3\n",
    "safety = g_syn_3D\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3234fc56-b15d-411e-a559-7307ee4c81c2",
   "metadata": {},
   "source": [
    "### Search space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b44c6dda-7305-4e1d-809f-cf4b31bf1fc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def comb_funcs(x: TensorType) -> TensorType:\n",
    "    return tf.stack([obj(x), safety(x)], axis = -1)\n",
    "\n",
    "def comb_funcs_plot(x:TensorType) -> TensorType:\n",
    "    return -comb_funcs(x)\n",
    "\n",
    "observer = mk_observer(\n",
    "    comb_funcs\n",
    ")\n",
    "\n",
    "num_objective = 2\n",
    "\n",
    "data_dim = 3\n",
    "data_num_per_dim = tf.constant(50)\n",
    "\n",
    "mins = [0, 0, 0]\n",
    "maxs = [1, 1, 1]\n",
    "\n",
    "search_space = Box(mins, maxs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c87baa42-b9ee-4d88-a5a7-e81efac7cbea",
   "metadata": {},
   "source": [
    "### Build Gaussian Process models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a41b003d-4989-4c3b-a9c2-3bb6bcc0f100",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow_probability as tfp\n",
    "\n",
    "def build_model(data, variance):\n",
    "    variance = tf.math.reduce_variance(data.observations)\n",
    "    kernel = gpflow.kernels.Matern52(variance=variance, lengthscales=[.2, .2, .2])\n",
    "    prior_scale = tf.cast(1.0, dtype=tf.float64)\n",
    "    kernel.variance.prior = tfp.distributions.LogNormal(\n",
    "        tf.cast(variance, dtype=tf.float64), prior_scale\n",
    "    )\n",
    "    \n",
    "    kernel.lengthscales.prior = tfp.distributions.LogNormal(\n",
    "        tf.math.log(kernel.lengthscales), prior_scale\n",
    "    )\n",
    "    \n",
    "    gpr = gpflow.models.GPR(data.astuple(), kernel, noise_variance=1e-5)\n",
    "    gpflow.set_trainable(gpr.likelihood, False)\n",
    "    return GaussianProcessRegression(gpr, num_kernel_samples=10)\n",
    "\n",
    "def build_stacked_independent_objectives_model(\n",
    "    data: Dataset, num_output: int, search_space: SearchSpace, variances = [1, 1, 1]\n",
    ") -> TrainableModelStack:\n",
    "    gprs = []\n",
    "    for idx in range(num_output):\n",
    "        single_obj_data = Dataset(\n",
    "            data.query_points, tf.gather(data.observations, [idx], axis=1)\n",
    "        )\n",
    "        gprs.append((build_model(single_obj_data, variances[idx]), 1))\n",
    "\n",
    "    return TrainableModelStack(*gprs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03aba379-c0dd-4b4c-a917-0dc499506221",
   "metadata": {},
   "source": [
    "### M-SafeOpt (for finding global safe optimum)\n",
    "\n",
    "[Memory efficient implementation]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f9ba6a28-9ab2-40b2-b800-312ee0ca8f89",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = tf.constant(2.0, tf.float64)   #safety threshold\n",
    "\n",
    "beta_f = tf.constant(3.0, tf.float64)   #scaling factor for confidence interval of GP_f\n",
    "beta_g = tf.constant(3.0, tf.float64)   #scaling factor for confidence interval of GP_g\n",
    "\n",
    "from trieste.acquisition import (\n",
    "    SingleModelAcquisitionBuilder,\n",
    "    ExpectedImprovement,\n",
    "    Product,\n",
    ")\n",
    "\n",
    "from trieste.data import Dataset\n",
    "from trieste.models import ProbabilisticModel\n",
    "from trieste.acquisition.rule import AcquisitionRule\n",
    "from trieste.types import TensorType\n",
    "from typing import Optional\n",
    "from collections.abc import Mapping\n",
    "import itertools\n",
    "\n",
    "class MSafeOPTMemEff(AcquisitionRule[TensorType, SearchSpace, ProbabilisticModel]):\n",
    "    def __init__(self, num_query_points: int = 1):\n",
    "        if num_query_points <= 0:\n",
    "            raise ValueError(\n",
    "                f\"Number of query points must be greater than 0, got {num_query_points}\"\n",
    "            )\n",
    "        self._num_query_points = num_query_points\n",
    "\n",
    "    def acquire(\n",
    "        self,\n",
    "        search_space: SearchSpace,\n",
    "        models: Mapping[str, ProbabilisticModel],\n",
    "        datasets: Optional[Mapping[str, Dataset]] = None,\n",
    "    ) -> TensorType:\n",
    "        \n",
    "        def create_grid_generator(shape, box_lower, box_upper):\n",
    "            for i in itertools.product(*[box_lower[x]+range(shape[x])*(box_upper[x]-box_lower[x])/(shape[x]-1) for x in range(1, len(shape))]):\n",
    "                yield [[x/(shape[0]-1)]+list(i) for x in range(shape[0])]\n",
    "                \n",
    "        @tf.function(reduce_retracing=True)\n",
    "        def acq_query_point():\n",
    "            box_lower = search_space.lower\n",
    "            box_upper = search_space.upper   \n",
    "            \n",
    "            num_dims = box_lower.shape[0]\n",
    "            num_per_dim = 50\n",
    "                                        \n",
    "            best_sample = tf.constant(box_lower, tf.float64)\n",
    "            max_interval = tf.constant(0, tf.float64)\n",
    "            max_safe_LCB = tf.constant(-np.inf, tf.float64)\n",
    "            \n",
    "            shape = [num_per_dim]*num_dims\n",
    "            grid_generator = tf.data.Dataset.from_generator(create_grid_generator, args = (shape, box_lower, box_upper), \n",
    "                                                            output_signature = (tf.TensorSpec(shape=(num_per_dim, len(shape),), dtype=tf.float64)))\n",
    "            \n",
    "            grid_generator = grid_generator.prefetch(tf.data.AUTOTUNE).cache()\n",
    "            \n",
    "            for x in grid_generator:\n",
    "                mean, var = model.predict(x)\n",
    "                stddev = tf.math.sqrt(var)\n",
    "\n",
    "                mean_f = tf.reshape(mean[:,0], [num_per_dim, 1])\n",
    "                mean_g = tf.reshape(mean[:,1], [num_per_dim, 1])\n",
    "                                    \n",
    "                stddev_f = tf.reshape(stddev[:,0], [num_per_dim, 1])\n",
    "                stddev_g = tf.reshape(stddev[:,1], [num_per_dim, 1])\n",
    "                \n",
    "                UCB_g = - mean_g + beta_g * stddev_g  # upper confidence bound\n",
    "                LCB_f = - mean_f - beta_f * stddev_f  # lower confidence bound\n",
    "                \n",
    "                mask_less_g = tf.reshape(tf.math.less_equal(UCB_g, h), [num_per_dim])  # mask for all points with UCB <= h\n",
    "                mask_more_g = tf.reshape(tf.math.greater(UCB_g, h), [num_per_dim, 1])  # mask for all points with UCB > h\n",
    "                                \n",
    "                mask_zero_g = tf.reshape(tf.one_hot(indices=[0], depth=num_per_dim, on_value=1.0, off_value=0.0), [num_per_dim, 1])  \n",
    "                mask_zero_g = tf.cast(mask_zero_g, tf.bool)\n",
    "\n",
    "                mask_zero_g = tf.reshape(tf.math.logical_and(mask_more_g, mask_zero_g),\n",
    "                                   [num_per_dim, 1])  # mask for all points with s=0 and UCB > h\n",
    "                                      \n",
    "                mask_safe_g = tf.reshape(mask_less_g,  [num_per_dim, 1])\n",
    "                mask_safe_g = tf.math.logical_or(mask_safe_g, mask_zero_g)\n",
    "\n",
    "                LCB_f_safe =  tf.math.multiply(LCB_f+tf.constant(100000, mean.dtype), tf.cast(mask_safe_g, mean.dtype))\n",
    "                LCB_f_safe_max = tf.reduce_max(LCB_f_safe) - tf.constant(100000, mean.dtype)\n",
    "\n",
    "                max_safe_LCB = tf.cond(max_safe_LCB < LCB_f_safe_max, lambda: LCB_f_safe_max, lambda: max_safe_LCB)\n",
    "                            \n",
    "            for x in grid_generator:\n",
    "                mean, var = model.predict(x)\n",
    "                stddev = tf.math.sqrt(var)\n",
    "                            \n",
    "                mean_f = tf.reshape(mean[:,0], [mean.shape[0], 1])\n",
    "                mean_g = tf.reshape(mean[:,1], [mean.shape[0], 1])\n",
    "                                    \n",
    "                stddev_f = tf.reshape(stddev[:,0], [mean.shape[0], 1])\n",
    "                stddev_g = tf.reshape(stddev[:,1], [mean.shape[0], 1])\n",
    "                \n",
    "                UCB_f = - mean_f + beta_f * stddev_f  # upper confidence bound\n",
    "                UCB_g = - mean_g + beta_g * stddev_g  # upper confidence bound\n",
    "                LCB_f = - mean_f - beta_f * stddev_f  # lower confidence bound\n",
    "                LCB_g = - mean_g - beta_g * stddev_g  # lower confidence bound\n",
    "                \n",
    "                mask_less_g = tf.reshape(tf.math.less_equal(UCB_g, h), [mean.shape[0]])  # mask for all points with UCB <= h\n",
    "                mask_more_g = tf.reshape(tf.math.greater(UCB_g, h), [mean.shape[0], 1])  # mask for all points with UCB > h\n",
    "                   \n",
    "                mask_zero_g = tf.reshape(tf.one_hot(indices=[0], depth=num_per_dim, on_value=1, off_value=0), [num_per_dim, 1])  \n",
    "                mask_zero_g = tf.cast(mask_zero_g, tf.bool)\n",
    "                \n",
    "                mask_zero_g = tf.reshape(tf.math.logical_and(mask_more_g, mask_zero_g), [mean.shape[0], 1])  # mask for all points with s=0 and UCB > h\n",
    "                                      \n",
    "                mask_safe_g = tf.reshape(mask_less_g,  [mean.shape[0], 1])\n",
    "                mask_safe_g = tf.math.logical_or(mask_safe_g, mask_zero_g)\n",
    "        \n",
    "                mask_less_g = tf.reshape(mask_less_g, [mean.shape[0], 1])\n",
    "        \n",
    "                mask_combined_g = tf.reshape(tf.math.logical_or(mask_less_g, mask_zero_g), [mean.shape[0]])\n",
    "                mask_combined_g = tf.reshape(mask_combined_g, [1, num_per_dim])\n",
    "                true_indices_g = tf.math.argmax(tf.reverse(mask_combined_g, [-1]), axis=-1)\n",
    "                const_g = (num_per_dim - 1) * tf.ones(1, tf.int64)\n",
    "                true_indices_g = tf.constant(num_per_dim-1, tf.int64) - true_indices_g\n",
    "                mask_final_g = tf.zeros([1, num_per_dim], dtype=tf.int64)\n",
    "                mask_final_g = tf.tensor_scatter_nd_add(mask_final_g,\n",
    "                                                      tf.stack([tf.range(1, dtype=tf.int64), true_indices_g], axis=1),\n",
    "                                                      tf.ones(1, dtype=tf.int64))\n",
    "                \n",
    "                safe_UCB_f = tf.math.multiply(tf.cast(mask_safe_g, tf.float64), UCB_f + tf.constant(100000, mean.dtype))\n",
    "                safe_UCB_f = tf.reshape(safe_UCB_f, [1, num_per_dim])\n",
    "                max_values = tf.reduce_max(safe_UCB_f, axis = -1, keepdims=True)\n",
    "                mask_max_f = tf.greater_equal(safe_UCB_f, max_values)\n",
    "                max_values = max_values - tf.constant(100000, mean.dtype)\n",
    "        \n",
    "                mask_final_g = tf.reshape(mask_final_g, [mean.shape[0], 1])\n",
    "                \n",
    "                # Eliminate suboptimal x's, and Limit expansion if suboptimal\n",
    "                input_s = tf.reshape(x[:, 0], [mean.shape[0], 1])\n",
    "                input_s_boundary = tf.math.multiply(input_s, tf.cast(mask_final_g, mean.dtype))\n",
    "                input_s_boundary = tf.reshape(input_s_boundary, [1, num_per_dim])\n",
    "                input_s_boundary = tf.reduce_max(input_s_boundary, axis = -1, keepdims=True)\n",
    "                \n",
    "                input_s_boundary_tiled = tf.tile(input_s_boundary,  [1 , num_per_dim])\n",
    "                input_s_diff = tf.reshape(input_s, [1, num_per_dim]) - input_s_boundary_tiled\n",
    "                \n",
    "                LCB_g_boundary = tf.math.multiply(LCB_g + tf.constant(100000, mean.dtype), tf.cast(mask_final_g, mean.dtype))\n",
    "                LCB_g_boundary = tf.reshape(LCB_g_boundary, [1, num_per_dim])\n",
    "                LCB_g_boundary = tf.reduce_max(LCB_g_boundary, axis = -1 , keepdims=True) - tf.constant(100000, mean.dtype)\n",
    "                LCB_g_boundary = tf.tile(LCB_g_boundary,  [1,num_per_dim])\n",
    "                pessimistic_LCB_g = LCB_g_boundary + l_g*input_s_diff\n",
    "                pessimistic_mask_g = tf.less_equal(pessimistic_LCB_g, h)\n",
    "                \n",
    "                UCB_f_boundary =  tf.math.multiply(UCB_f + tf.constant(100000, mean.dtype), tf.cast(mask_final_g, mean.dtype))\n",
    "                UCB_f_boundary = tf.reshape(UCB_f_boundary, [1, num_per_dim])\n",
    "                UCB_f_boundary = tf.reduce_max(UCB_f_boundary, axis = -1, keepdims=True) - tf.constant(100000, mean.dtype)\n",
    "                UCB_f_boundary = tf.tile(UCB_f_boundary,  [1,num_per_dim])\n",
    "                optimistic_UCB_f = UCB_f_boundary + L_f*input_s_diff\n",
    "                optimistic_UCB_f_masked = tf.math.multiply(optimistic_UCB_f+ tf.constant(100000, mean.dtype), tf.cast(pessimistic_mask_g, mean.dtype))\n",
    "                optimistic_UCB_f_max = tf.reduce_max(optimistic_UCB_f_masked, axis = -1, keepdims=True)-tf.constant(100000, mean.dtype)\n",
    "                \n",
    "                LCB_f_safe_max = max_safe_LCB\n",
    "                \n",
    "                expd_mask = tf.greater(optimistic_UCB_f_max, LCB_f_safe_max)\n",
    "                expd_mask = tf.tile(expd_mask,[1, num_per_dim])\n",
    "                expd_mask = tf.reshape(expd_mask, [mean.shape[0], 1])\n",
    "                            \n",
    "                elim_mask = tf.greater(max_values, LCB_f_safe_max)\n",
    "                elim_mask = tf.tile(elim_mask, [1, num_per_dim])\n",
    "                elim_mask = tf.reshape(elim_mask, [mean.shape[0], 1])\n",
    "        \n",
    "                elim_mask = tf.math.logical_or(elim_mask, expd_mask)       \n",
    "                \n",
    "                mask_final_g = tf.math.logical_and(tf.cast(mask_final_g, tf.bool), elim_mask)\n",
    "                mask_final_g = tf.math.logical_and(tf.cast(mask_final_g, tf.bool), expd_mask)        \n",
    "\n",
    "                mask_max_f = tf.reshape(mask_max_f, [mean.shape[0], 1])\n",
    "                mask_max_f = tf.math.logical_and(mask_max_f, elim_mask)\n",
    "                mask_final_f = tf.math.logical_or(tf.cast(mask_max_f, tf.bool), tf.cast(mask_final_g, tf.bool))\n",
    "                                \n",
    "                max_interval_f = beta_f * tf.math.multiply(tf.reshape(stddev_f, [stddev.shape[0], 1]),\n",
    "                                        tf.cast(mask_final_f, mean.dtype))\n",
    "                \n",
    "                max_interval_g = beta_g * tf.math.multiply(tf.reshape(stddev_g, [stddev.shape[0], 1]),\n",
    "                                        tf.cast(tf.reshape(mask_final_g, [stddev.shape[0], 1]), mean.dtype))\n",
    "\n",
    "                curr_max_f = tf.reduce_max(max_interval_f)\n",
    "                max_interval = tf.cond(curr_max_f >= max_interval, lambda: curr_max_f, lambda: max_interval)\n",
    "                best_sample = tf.cond(curr_max_f >= max_interval, lambda: tf.gather(x, tf.where(max_interval_f == curr_max_f)[0,0]), lambda: best_sample)\n",
    "\n",
    "                curr_max_g = tf.reduce_max(max_interval_g)\n",
    "                max_interval = tf.cond(curr_max_g >= max_interval, lambda: curr_max_g, lambda: max_interval)\n",
    "                best_sample = tf.cond(curr_max_g >= max_interval, lambda: tf.gather(x, tf.where(max_interval_g == curr_max_g)[0,0]), lambda: best_sample)\n",
    "\n",
    "            return tf.reshape(best_sample, (1, num_dims))\n",
    "        \n",
    "        return acq_query_point()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b5106b1-9e34-40ce-94f4-56b9caf0dd54",
   "metadata": {},
   "source": [
    "### Run safe Bayesian optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "27cedf03-c6df-49d1-90bc-abb9ac06a385",
   "metadata": {},
   "outputs": [],
   "source": [
    "rule = MSafeOPTMemEff()\n",
    "\n",
    "L_f = tf.constant(1.9, tf.float64)  #hartmann_3\n",
    "l_g = tf.constant(1, tf.float64)  #g_syn_3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4b084c82-8ac6-49a7-8f99-d70f84fa8cbb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run #1: Step #1\n",
      "Run #1: Step #2\n",
      "Run #1: Step #3\n",
      "WARNING:tensorflow:5 out of the last 5 calls to <function MSafeOPTMemEff.acquire.<locals>.acq_query_point at 0x7f77a87c1f80> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.\n",
      "Run #1: Step #4\n",
      "WARNING:tensorflow:6 out of the last 6 calls to <function MSafeOPTMemEff.acquire.<locals>.acq_query_point at 0x7f77a1f49760> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.\n",
      "Run #1: Step #5\n",
      "Run #1: Step #6\n",
      "Run #1: Step #7\n",
      "Run #1: Step #8\n",
      "Run #1: Step #9\n",
      "Run #1: Step #10\n",
      "Run #1: Step #11\n",
      "Run #1: Step #12\n",
      "Run #1: Step #13\n",
      "Run #1: Step #14\n",
      "Run #1: Step #15\n",
      "Run #1: Step #16\n",
      "Run #1: Step #17\n",
      "Run #1: Step #18\n",
      "Run #1: Step #19\n",
      "Run #1: Step #20\n",
      "Run #1: Step #21\n",
      "Run #1: Step #22\n",
      "Run #1: Step #23\n",
      "Run #1: Step #24\n",
      "Run #1: Step #25\n",
      "Run #1: Step #26\n",
      "Run #1: Step #27\n",
      "Run #1: Step #28\n",
      "Run #1: Step #29\n",
      "Run #1: Step #30\n",
      "Run #1: Step #31\n",
      "Run #1: Step #32\n",
      "Run #1: Step #33\n",
      "Run #1: Step #34\n",
      "Run #1: Step #35\n",
      "Run #1: Step #36\n",
      "Run #1: Step #37\n",
      "Run #1: Step #38\n",
      "Run #1: Step #39\n",
      "Run #1: Step #40\n",
      "Run #1: Step #41\n",
      "Run #1: Step #42\n",
      "Run #1: Step #43\n",
      "Run #1: Step #44\n",
      "Run #1: Step #45\n",
      "Run #1: Step #46\n",
      "Run #1: Step #47\n",
      "Run #1: Step #48\n",
      "Run #1: Step #49\n",
      "Run #1: Step #50\n"
     ]
    }
   ],
   "source": [
    "from trieste.ask_tell_optimization import AskTellOptimizer\n",
    "\n",
    "num_runs = 1\n",
    "num_steps = 50\n",
    "\n",
    "for num_run in range(1, num_runs+1):\n",
    "\n",
    "    num_initial_points = 2\n",
    "    initial_query_points = tf.constant([[0.0, maxs[0]*random.random(), maxs[1]*random.random()], [0.0, maxs[0]*random.random(), maxs[1]*random.random()]], dtype=tf.float64)\n",
    "    initial_data = observer(initial_query_points)\n",
    "    \n",
    "    model = build_stacked_independent_objectives_model(\n",
    "        initial_data, num_objective, search_space\n",
    "    )\n",
    "\n",
    "    ask_tell = AskTellOptimizer(search_space, initial_data, model, acquisition_rule = rule)\n",
    "    pred_means = []\n",
    "    pred_vars = []    \n",
    "    dataset = initial_data\n",
    "\n",
    "    for step in range(1, num_steps + 1):\n",
    "        print(f\"Run #{num_run}: Step #{step}\")\n",
    "        new_point = ask_tell.ask()\n",
    "        new_data_point = observer(new_point)\n",
    "        dataset = dataset + new_data_point\n",
    "        ask_tell.tell(dataset)       \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8264ae9d-d4d1-4c2a-abd2-2b330095ca8f",
   "metadata": {},
   "source": [
    "### Plot cumulative regret: $R_t/t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "221cca1b-8408-459a-ba78-fbbafa988714",
   "metadata": {},
   "outputs": [],
   "source": [
    "query_points = dataset.query_points.numpy()[2:]\n",
    "observations = dataset.observations.numpy()[2:]\n",
    "\n",
    "max_obj_value = 3.86278     # https://www.sfu.ca/~ssurjano/hart3.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "84443ed5-4cce-4c7a-bb1c-afa57330d7eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$R_t /t$')"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjrklEQVR4nO3deVhU1f8H8PeAMgLCuIGgDAhpiJm45AJmaGqiZthqaoJmlqVlmWmUS5tS+eurWe4bZZmmpaaVZSaUgKYi7mugYLJo6rCoqMz5/XGagZFFloE7zLxfzzOPzL137pwZwHlz7uecoxJCCBARERHZIDulG0BERESkFAYhIiIislkMQkRERGSzGISIiIjIZjEIERERkc1iECIiIiKbxSBERERENotBiIiIiGxWHaUbYCn0ej3Onz8PFxcXqFQqpZtDRERE5SCEQE5ODpo1awY7u4r37zAI/ef8+fPQarVKN4OIiIgqIS0tDV5eXhV+nMUHoQ8//BCRkZGYMGEC5s6dW+px69atw7Rp03DmzBm0atUKH330EQYMGFDu53FxcQEg30hXV9eqNpuIiIhqQHZ2NrRarfFzvKIsOgjt2bMHixcvRrt27co8Lj4+HkOHDkVUVBQefvhhrF69GoMHD0ZiYiLatm1brucyXA5zdXVlECIiIqplKlvWYrHF0rm5uRg+fDiWLl2Khg0blnnsp59+itDQULzxxhsICAjA+++/j44dO+Lzzz+vodYSERFRbWSxQWjcuHEYOHAg+vTpc8djExISih3Xr18/JCQklPqY/Px8ZGdnm9yIiIjItljkpbE1a9YgMTERe/bsKdfxGRkZaNq0qcm2pk2bIiMjo9THREVF4d13361SO4mIiKh2s7gglJaWhgkTJmDbtm2oV69etT1PZGQkJk6caLxvKLYiIipLQUEBbt68qXQziGyKg4NDpYbGl4fFBaF9+/YhKysLHTt2NG4rKCjAH3/8gc8//xz5+fmwt7c3eYyHhwcyMzNNtmVmZsLDw6PU51Gr1VCr1eZtPBFZLSEEMjIycOXKFaWbQmRz7Ozs4OvrCwcHB7Of2+KCUO/evXHo0CGTbaNGjULr1q0xZcqUYiEIAIKCgrB9+3a8+uqrxm3btm1DUFBQdTeXiGyEIQS5u7vDycmJE68S1RDDhMfp6enw9vY2+++exQUhFxeXYkPenZ2d0bhxY+P28PBwNG/eHFFRUQCACRMmICQkBJ988gkGDhyINWvWYO/evViyZEmNt5+IrE9BQYExBDVu3Fjp5hDZHDc3N5w/fx63bt1C3bp1zXpuix01VpbU1FSkp6cb7wcHB2P16tVYsmQJAgMDsX79emzcuLHccwgREZXFUBPk5OSkcEuIbJPhklhBQYHZz60SQgizn7UWys7OhkajgU6n44SKRGTi+vXrSElJga+vb7UO4iCikpX1O1jVz+9a2SNEytPpgHPnSt537pzcrwRLbJcltomIiCQGIaownQ4IDQVCQoC0NNN9aWlye2hozX/AW2K7LLFNRLXZAw88gNWrVyvdjFpp5MiRGDx4sMWcx+DGjRto0aIF9u7da7ZzVgSDEFVYTg6QlQUkJwM9exZ+wKelyfvJyXJ/Tg7bZYltItsycuRIqFQqqFQq1K1bF76+vpg8eTKuX7+udNMAANHR0WjQoEG5jv3hhx+QmZmJp59+uti+qKgo2NvbY/bs2WZuYfkJIbBkyRJ07doV9evXR4MGDXDfffdh7ty5uHr1qmLtqqwzZ85ApVIhKSnJZPunn36K6Ohosz2Pg4MDJk2ahClTppjtnBXBIEQV5uUFxMQAfn6FH/Dx8YUf7H5+cr+XF9tliW0i2xMaGor09HQkJydjzpw5WLx4MWbMmFGtz1lQUAC9Xm/Wc86bNw+jRo0qcWK9FStWYPLkyVixYoVZn7MiRowYgVdffRVhYWHYsWMHkpKSMG3aNGzatAm//vqrYu0yN41GU+7wWl7Dhw/Hzp07ceTIEbOet1wECSGE0Ol0AoDQ6XRKN6XWSE0Vws9PCKDw5ucnt7Ndlt8mKr9r166Jo0ePimvXrindlAqLiIgQYWFhJtsee+wx0aFDB+P9goICMWvWLNGiRQtRr1490a5dO7Fu3TqTx2zatEm0bNlSqNVq0bNnTxEdHS0AiMuXLwshhFi5cqXQaDRi06ZNIiAgQNjb24uUlBRx/fp18frrr4tmzZoJJycn0aVLF7Fjxw4hhBA7duwQAExuM2bMKPF1ZGVlCZVKJQ4fPlxsX0xMjGjevLm4ceOGaNasmYiLizO+rubNm4sFCxaYHJ+YmChUKpU4c+aMEEKIY8eOie7duwu1Wi0CAgLEtm3bBACxYcOGcr7LQqxdu1YAEBs3biy2T6/XiytXrgghhAgJCRETJkww2R8WFiYiIiKM9318fMT7778vRowYIZydnYW3t7fYtGmTyMrKEo888ohwdnYW9957r9izZ4/xMTNmzBCBgYEm550zZ47w8fEx3r/9Z+Hnn38W3bt3FxqNRjRq1EgMHDhQnD592rj/9u9NSEhIsfMsXrxYeHp6ioKCApPnfuSRR8SoUaOM9zdu3Cg6dOgg1Gq18PX1Fe+88464efOmyWN69eolpk6dWuz9E6Ls38Gqfn6zR4gqTasFVq0y3bZqldyuJEtslyW2yZxssSBcCCAvT5lbVcb6Hj58GPHx8SYz9EZFReHLL7/EokWLcOTIEbz22mt45plnEBsbCwBISUnBE088gcGDB+PAgQN44YUX8Pbbbxc799WrV/HRRx9h2bJlOHLkCNzd3TF+/HgkJCRgzZo1OHjwIJ588kmEhobi1KlTCA4Oxty5c+Hq6or09HSkp6dj0qRJJbZ7586dcHJyQkBAQLF9y5cvx9ChQ1G3bl0MHToUy5cvByBnIx46dGixmqKvv/4a3bt3h4+PDwoKCjB48GA4OTlh9+7dWLJkSYmv7U6+/vpr+Pv7IywsrNg+lUoFjUZTofPNmTMH3bt3x/79+zFw4ECMGDEC4eHheOaZZ5CYmIi77roL4eHhEFX4YcjLy8PEiROxd+9ebN++HXZ2dnj00UeNPXl//fUXAOC3335Deno6vv/++2LnePLJJ/Hvv/9ix44dxm2XLl3C1q1bMXz4cADAn3/+ifDwcEyYMAFHjx7F4sWLER0djZkzZ5qcq0uXLvjzzz8r/XoqrVLxyQrZQo/QlStCpKWVvC8tTe6vCEvt5bDEdllim8zlyhUhunUr+fUYXne3bhX/+bIkJf01mptr+v2syVtubvnbHhERIezt7YWzs7NQq9UCgLCzsxPr168XQghx/fp14eTkJOLj400eN3r0aDF06FAhhBBTpkwRbdu2Ndn/9ttvF+sRAiCSkpKMx5w9e1bY29uLf/75x+SxvXv3FpGRkcbHaTSaO76OOXPmCD8/v2LbdTqdcHR0ND7v/v37Rf369UVOTo7xvkqlEmfPnhVCFPYSLVy4UAghe0Xq1Kkj0tPTjeesTI9QQECAeOSRR+54XHl7hJ555hnj/fT0dAFATJs2zbgtISFBADC2uzI9Qre7cOGCACAOHTokhBAiJSVFABD79+83Oe7284SFhYlnn33WeH/x4sWiWbNmxl6i3r17i1mzZpmcY9WqVcLT09Nk26effipatGhRYtvYI0RVZu7RS0WLff38gLg40zqY25+jplhiuyyxTebEgnDL16tXLyQlJWH37t2IiIjAqFGj8PjjjwMATp8+jatXr6Jv376oX7++8fbll1/i77//BgCcOHECnTt3Njlnly5dij2Pg4MD2rVrZ7x/6NAhFBQU4O677zY5d2xsrPHc5XXt2rUS53D65ptvcNdddyEwMBAA0L59e/j4+GDt2rXG+wEBAcZeodjYWGRlZeHJJ580vjatVmuyNmVJr+1OhJmn5Cv6PjZt2hQAcO+99xbblpWVVennOHXqFIYOHQo/Pz+4urqiRYsWAOSkxRUxfPhwfPfdd8jPzwcge8eefvppYy3XgQMH8N5775n8DIwZMwbp6ekmReSOjo6KFJVb3BIbVD1u/7CKiZGXZYp+WBmOu1MP7rlzxYt9tVr5r2F7z55AbGzNFgFbYrsssU3mZigIL/p6Vq0CRoyw7oJwJycgN1e5564IZ2dntGzZEoAsKg4MDMTy5csxevRo5P73In788Uc0b97c5HEVXZja0dHRZB2o3Nxc2NvbY9++fcXWiaxfv36Fzt2kSRNcvny52Pbly5fjyJEjqFOn8ONMr9djxYoVGD16NAD5Qb169Wq8+eabWL16NUJDQ82+VMrdd9+N48eP3/E4Ozu7YqHJMHN5UUWXkTC8pyVtM1zGKu95ixo0aBB8fHywdOlSNGvWDHq9Hm3btsWNGzfu+DpuP48QAj/++CM6d+6MP//8E3PmzDHuz83NxbvvvovHHnus2GOLhttLly7Bzc2tQs9tDgxCNsKcH1YuLoC7u/za8MEOmH7Au7vL42qSJbbLEttUHW4Pd927y+1Fw5+1UakAZ2elW1FxdnZ2eOuttzBx4kQMGzYMbdq0gVqtRmpqKkJCQkp8jL+/P3766SeTbXv27Lnjc3Xo0AEFBQXIyspCjx49SjzGwcGhXMsmdOjQARkZGbh8+TIaNmwIQPY47d27FzExMWjUqJHx2EuXLqFnz544fvw4WrdujWHDhmHq1KnYt28f1q9fj0WLFpm8trS0NGRmZhp7Wcrz2m43bNgwPP3009i0aVOxOiEhhHH2Yzc3N5MlogoKCnD48GH06tWrws9ZlJubGzIyMiCEMIak24e9F/Xvv//ixIkTWLp0qfF7s3PnTpNjyrusRb169fDYY4/h66+/xunTp+Hv74+OHTsa93fs2BEnTpwwhvHSHD58GB06dCjzmGpRqQtqVsgWaoSEMF+tirnrjczFEttliW2qLnFxpj9b/w3eqfWsbdTYzZs3RfPmzcXs2bOFELLep3HjxiI6OlqcPn1a7Nu3T8ybN09ER0cLIYRITk4WdevWFZMnTxYnTpwQa9euFV5eXgKAcTRUabU+w4cPFy1atBDfffedSE5OFrt37xazZs0SW7ZsEUIIERcXJwCI3377TVy4cEHk5eWV+Dpu3bol3NzcxObNm43bJkyYILp27Vri8V26dBGTJk0y3u/evbsIDAwULi4u4urVqybn9ff3F/369RMHDhwQO3fuFN26dSs2AuzBBx8Un332WWlvs9Dr9WLIkCHC0dFRzJw5U+zZs0ecOXNGbN68WTz44IPGeqNFixYJJycnsWXLFnHs2DExZswY4erqWqxGaM6cOSbnx201S7fX7xw9elSoVCrx4YcfitOnT4vPP/9cNGzYsNQaoYKCAtG4cWPxzDPPiFOnTont27eLzp07mzzPzZs3haOjo/jggw9ERkaG8Xtd0s/Utm3bhFqtFv7+/uL999832bd161ZRp04d8c4774jDhw+Lo0ePim+++Ua8/fbbJsf5+PiIL7/8ssT3tzprhBiE/mMrQUgI6/2wImVZc0G4tQUhIYSIiooSbm5uIjc3V+j1ejF37lzh7+8v6tatK9zc3ES/fv1EbGys8fjbh88vXLhQADC+J6UFoRs3bojp06eLFi1aiLp16wpPT0/x6KOPioMHDxqPGTt2rGjcuHGZw+eFEGLy5Mni6aefFkIIkZ+fLxo3biw+/vjjEo/96KOPhLu7u7hx44YQQogFCxYIACI8PLzYsYbh8w4ODqJ169Zi8+bNAoDYunWr8RgfH58y2yaEDBcLFy4UnTt3Fk5OTsLV1VV06tRJfPrpp8bwdePGDfHiiy+KRo0aCXd3dxEVFVVisXRFg5AQQixcuFBotVrh7OwswsPDxcyZM8sslt62bZsICAgQarVatGvXTsTExBR7nqVLlwqtVivs7OxKHD5f9LV7enoKAOLvv/8u9t5s3bpVBAcHC0dHR+Hq6iq6dOkilixZYtwfHx8vGjRoYBJSi6rOIMRFV/9jK4uu3l4TBFj35QuqGbcXhJd02bU2/3xx0dXiZs6ciUWLFiGtBqv9MzIycM899yAxMRE+Pj7V9jxxcXG4//77cfr0adx1113V9jxUaMiQIQgMDMRbb71V4n4uukpmYe2jl0gZJRWEBwcXn1G7tHmGqHZYsGAB9uzZg+TkZKxatQqzZ89GREREjbbBw8MDy5cvr/CopjvZsGEDtm3bhjNnzuC3337D888/j+7duzME1ZAbN27g3nvvxWuvvabI87NY2kbYwuglUoatFITbulOnTuGDDz7ApUuX4O3tjddffx2RkZE13g5zLvZpkJOTgylTpiA1NRVNmjRBnz598Mknn5j9eahkDg4OmDp1qmLPz0tj/7H2S2OGeYSysopfpjD0FLm7A1u33nn4vK3Q6eR0AiUFw3Pn5Ac73yvJ2t8rXhojUlZ1Xhpjj5CN0GhkyCnpw0qrlT1Btf3DypwYHCtGoyn9fWAPIxFZMtYI2RCNpvQPJS8vfqAXxdmSiYhsA4MQUQkME1AWLfaNjy9eZ8XeDiKi2o2XxohKYYuzJRMR2Rr2CBGVQauVc+IUtWoVQxARkbVgECIqQ1qanBiwqBEjOOcSEZG1YBAiKgUnoCSyPA888ABWr16tdDPoNosWLcKgQYOUbkalMAgRlYCzJZO1GDlyJFQqFVQqFerWrQtfX19MnjwZ169fV7ppAIDo6Gg0aNCgXMf+8MMPyMzMxNNPP23c1qJFC8ydO9csbTlz5gxUKlWZq7ZXVkxMDFQqFa5cuWL2c1uCZ599FomJifjzzz+VbkqFMQgRlcAwW/LthdGGAmo/P86WTBWj05UenM+dk/urS2hoKNLT05GcnIw5c+Zg8eLFmDFjRvU9IYCCggLo9XqznnPevHkYNWoU7Oz40WVpHBwcMGzYMMybN0/pplQYf5qISmCYgDI2tnhhtGECSk6mSOVlmKAzJKT4JdW0NLk9NLT6wpBarYaHhwe0Wi0GDx6MPn36YNu2bcb9er0eUVFR8PX1haOjIwIDA7F+/XqTc/zwww9o1aoV6tWrh169euGLL74w6eEw9Oz88MMPaNOmDdRqNVJTU5Gfn49JkyahefPmcHZ2RteuXRETEwNA9pKMGjUKOp3O2Gv1zjvvlPgaLly4gN9///2Ol19UKhWWLVuGRx99FE5OTmjVqhV++OEH4/7Lly9j+PDhcHNzg6OjI1q1aoWVK1cCAHx9fQEAHTp0gEqlQs+ePQEAe/bsQd++fdGkSRNoNBqEhIQgMTGx3M975swZ9OrVCwDQsGFDqFQqjBw5EgCwdetW3H///WjQoAEaN26Mhx9+GH///bfxvIZequ+//x69evWCk5MTAgMDkZCQYPL8O3fuRI8ePeDo6AitVotXXnkFeXl5Jq87PDwcDRs2hJOTE/r3749Tp04Z97/zzjto3769yTnnzp2LFi1aGO/HxMSgS5cucHZ2RoMGDdC9e3ecPXvWuH/QoEH44YcfcO3atTK/R5aGQaiaKPnXH5mHJU5AyZ+r2smSJug8fPgw4uPj4eDgYNwWFRWFL7/8EosWLcKRI0fw2muv4ZlnnkFsbCwAICUlBU888QQGDx6MAwcO4IUXXsDbb79d7NxXr17FRx99hGXLluHIkSNwd3fH+PHjkZCQgDVr1uDgwYN48sknERoailOnTiE4OBhz586Fq6sr0tPTkZ6ejkmTJpXY7p07d8LJyQkBAQF3fI3vvvsunnrqKRw8eBADBgzA8OHDcenSJQDAtGnTcPToUfz88884duwYFi5ciCZNmgAA/vrrLwDAb7/9hvT0dHz//fcA5FpkERER2LlzJ3bt2oVWrVphwIAByLntG1ba82q1Wnz33XcAgBMnTiA9PR2ffvopACAvLw8TJ07E3r17sX37dtjZ2eHRRx8t1pv29ttvY9KkSUhKSsLdd9+NoUOH4tatWwCAv//+G6GhoXj88cdx8OBBrF27Fjt37sT48eONjx85ciT27t2LH374AQkJCRBCYMCAAbh58+Yd308AuHXrFgYPHoyQkBAcPHgQCQkJeP7556FSqYzH3Hfffbh16xZ2795drnNaDEFCCCF0Op0AIHQ6XZXPdeWKEN26CeHnJ0Rqqum+1FS5vVs3eRxRefHnSjnXrl0TR48eFdeuXav0OQzfI0D+Gxdnev/276m5RERECHt7e+Hs7CzUarUAIOzs7MT69euFEEJcv35dODk5ifj4eJPHjR49WgwdOlQIIcSUKVNE27ZtTfa//fbbAoC4fPmyEEKIlStXCgAiKSnJeMzZs2eFvb29+Oeff0we27t3bxEZGWl8nEajuePrmDNnjvDz8yu23cfHR8yZM8d4H4CYOnWq8X5ubq4AIH7++WchhBCDBg0So0aNKvE5UlJSBACxf//+MttSUFAgXFxcxObNm8v9vDt27DB5v0pz4cIFAUAcOnTIpE3Lli0zHnPkyBEBQBw7dkwIIb9Xzz//vMl5/vzzT2FnZyeuXbsmTp48KQCIuLg44/6LFy8KR0dH8e233wohhJgxY4YIDAw0OcecOXOEj4+PEEKIf//9VwAQMTExZba/YcOGIjo6usxjKqOs38Gqfn5zQsVqcPtff4Yak6J//RmO46UVKi/+XNVuSk7Q2atXLyxcuBB5eXmYM2cO6tSpg8cffxwAcPr0aVy9ehV9+/Y1ecyNGzfQoUMHALIXo3Pnzib7u3TpUux5HBwc0K5dO+P9Q4cOoaCgAHfffbfJcfn5+WjcuHGFXsO1a9fKveBt0TY4OzvD1dUVWVlZAIAXX3wRjz/+OBITE/HQQw9h8ODBCA4OLvN8mZmZmDp1KmJiYpCVlYWCggJcvXoVqamp5X7e0pw6dQrTp0/H7t27cfHiRWNPUGpqKtq2bVviuT09PQEAWVlZaN26NQ4cOICDBw/i66+/Nh4jhIBer0dKSgpOnTqFOnXqoGvXrsb9jRs3hr+/P44dO1Zm+wwaNWqEkSNHol+/fujbty/69OmDp556ytgWA0dHR1y9erVc57QUDELVwLA8g+E/vJAQYMECYNy4yi3PYO0re1P53P5z1bOnnNxxxAgu+1FbGCboNIQgoGYm6HR2dkbLli0BACtWrEBgYCCWL1+O0aNHIzc3FwDw448/onnz5iaPU6vVFXoeR0dHk0slubm5sLe3x759+2Bvb29ybP369St07iZNmuDy5cvlOrZu3bom91UqlTFg9O/fH2fPnsVPP/2Ebdu2oXfv3hg3bhz+7//+r9TzRURE4N9//8Wnn34KHx8fqNVqBAUF4caNG+V+3tIMGjQIPj4+WLp0KZo1awa9Xo+2bduWeW7De2w4d25uLl544QW88sorxc7v7e1tUgtUGjs7OwghTLbdftls5cqVeOWVV7B161asXbsWU6dOxbZt29CtWzfjMZcuXYKbm9sdn8+SMAhVE8Nff126ACkpQP/+cntF//rjKuhUFJf9qN1Km6CzJr93dnZ2eOuttzBx4kQMGzbMpLA5JCSkxMf4+/vjp59+Mtm2Z8+eOz5Xhw4dUFBQgKysLPTo0aPEYxwcHFBQUFCuc2VkZODy5cto2LDhHY8vi5ubGyIiIhAREYEePXrgjTfewP/93/8Z66Zub09cXBwWLFiAAQMGAADS0tJw8eLFCj1nSef+999/ceLECSxdutT4/uzcubPCr6djx444evSoMezeLiAgwFi7Y+j9Mjx3mzZtAMj3JCMjA0IIY9AqaRqBDh06oEOHDoiMjERQUBBWr15tDEJ///03rl+/buxJrC1YLF2NtFrgf/8z3VbRv/4sqciSLAOX/aidLGmCzieffBL29vaYP38+XFxcMGnSJLz22mv44osv8PfffyMxMRGfffYZvvjiCwDACy+8gOPHj2PKlCk4efIkvv32W0RHRwOASQ/Q7e6++24MHz4c4eHh+P7775GSkoK//voLUVFR+PHHHwHIeYByc3Oxfft2XLx4sdTLKh06dECTJk0QFxdXpdc+ffp0bNq0CadPn8aRI0ewZcsWYwG2u7s7HB0dsXXrVmRmZkL33+iDVq1aYdWqVTh27Bh2796N4cOHw9HRsULP6+PjA5VKhS1btuDChQvIzc1Fw4YN0bhxYyxZsgSnT5/G77//jokTJ1b4NU2ZMgXx8fEYP348kpKScOrUKWzatMlYLN2qVSuEhYVhzJgx2LlzJw4cOIBnnnkGzZs3R1hYGACgZ8+euHDhAj7++GP8/fffmD9/Pn7++Wfjc6SkpCAyMhIJCQk4e/Ysfv31V5w6dcqkeP3PP/+En58f7rrrrgq/BkVVqXrJipizWNogNVUIX19ZDGm4VaYoUqkiS7JMRX8eqvJzReVX1WLptLSSf2dv/91OSzNjo/8TEREhwsLCim2PiooSbm5uIjc3V+j1ejF37lzh7+8v6tatK9zc3ES/fv1EbGys8fhNmzaJli1bCrVaLXr27CkWLlwoABjfk9KKnm/cuCGmT58uWrRoIerWrSs8PT3Fo48+Kg4ePGg8ZuzYsaJx48YCgJgxY0apr2Xy5Mni6aefNtlWUrH0hg0bTI7RaDRi5cqVQggh3n//fREQECAcHR1Fo0aNRFhYmEhOTjYeu3TpUqHVaoWdnZ0ICQkRQgiRmJgo7rvvPlGvXj3RqlUrsW7dugo/rxBCvPfee8LDw0OoVCoREREhhBBi27ZtIiAgQKjVatGuXTsRExNjcq6SCrgvX74sAIgdO3YYt/3111+ib9++on79+sLZ2Vm0a9dOzJw507j/0qVLYsSIEUKj0QhHR0fRr18/cfLkSZP2Lly4UGi1WuHs7CzCw8PFzJkzjcXSGRkZYvDgwcLT01M4ODgIHx8fMX36dFFQUGB8/EMPPSSioqJEdajOYmmLDEILFiwQ9957r3BxcREuLi6iW7du4qeffir1eMNohaI3tVpdoec0dxAq+h+cWi3/dXMzTxjih5/tYihWRlWDkDWO+Pvggw+El5dXjT5nenq6aNSokThz5kyNPi/d2eHDh4W7u7u4Uk0/xNUZhCzy0piXlxc+/PBD7Nu3D3v37sWDDz6IsLAwHDlypNTHFJ2HIj093WSSp5p2+/IMQ4fK7WFhlV+egZdDiMt+1F7WMEHnggULsGfPHiQnJ2PVqlWYPXs2IiIiarQNHh4eWL58ebHRWqS89PR0fPnll9BY8g9xKSyyWPr2mUNnzpyJhQsXYteuXbjnnntKfIxKpYKHh0e5nyM/Px/5+fnG+9nZ2ZVrbAkMyzMA8kNq+3YgOho4daqw0LWiyzNYQpElKev2n6vbl/2ozM8V1RyNpvSgUxtG+p06dQoffPABLl26BG9vb7z++uuIjIys8XYMHjy4xp+T7qxPnz5KN6HSLDIIFVVQUIB169YhLy8PQUFBpR6Xm5sLHx8f6PV6dOzYEbNmzSo1NAFyJtV33323Opps/OvPMOS9Y0e5ff9+oHlz+ddfRYa8315kWXTIdNH5ZMi63f5zVZShV4FTKVB1mTNnDubMmaN0M4jMziIvjQFyIq769etDrVZj7Nix2LBhg3GY3+38/f2xYsUKbNq0CV999RX0ej2Cg4NxroxrBJGRkdDpdMZbmpmHbBRdniEgAFCrgexsGV4qsjwDL4dQUZa47AcRUW1msUHI398fSUlJ2L17N1588UVERETg6NGjJR4bFBSE8PBwtG/fHiEhIfj+++/h5uaGxYsXl3p+tVoNV1dXk1t1qVsXMEwKets6fXfEVdCpunDdsooTt004R0Q1ozp/9yw2CDk4OKBly5bo1KkToqKiEBgYaFyk7k7q1q2LDh064PTp09XcyvIzzC+1f3/FHmcNRZYFBcDOncCcOUB6utKtIUD51dBrG8OsvrVt6QAia2GYafv2GcrNweJrhAz0er1JcXNZCgoKcOjQIeMsoJbAUCdU0R4hoHYWWV67BmzbBmzaBGzeDFy4ILfHxQHr1yvbNuK6ZRVlb2+PBg0aGNeNcnJyKnMiQSIyH71ejwsXLsDJyQl16pg/tlhkEIqMjET//v3h7e2NnJwcrF69GjExMfjll18AAOHh4WjevDmioqIAAO+99x66deuGli1b4sqVK5g9ezbOnj2L5557TsmXYaJoEBICsMb/Q//9F9iyBdi4Efj1V6DoH88ajexd2LwZuHIFaNBAoUYSAK5bVhmGUal3WkSTiMzPzs4O3t7e1fIHiEUGoaysLISHhyM9PR0ajQbt2rXDL7/8YlwdOTU1FXZ2hVf1Ll++jDFjxiAjIwMNGzZEp06dEB8fX2pxtRLuvRewtwcuXpT1F0qP8tLrC9tRlZ8rIYBvvwUWLgT+/FOe18DbGxg8WM6f1KOHvDx45Ajw/ffAs89W+SVQFXHdsopRqVTw9PSEu7t7scUoiah6OTg4mHzum5NKsPoPgJxHSKPRQKfTVVvhdGAgcPCg7DH5b3kXReh0wBNPAL/9BjzyCDB/fuX+8r9wARg7VgYbg/bt5WsbPFi+3qIhKyoKeOstoHdv+dxkGeLjTVdDj4uTIxOJiGqDqn5+W2yxtDWqSp2QuZw7J3tnDEHkhx+ANm1kj07R3pw72bQJaNtWhqA6dYBp04CUFFkM/s47MhDd3tP09NPy399/B86fN8eroaoqbaLOmlwAlIhISQxCNUjpIHTwINCtG3DoEODhAaxdC3TtKgtiX3oJeOAB4Nixss+h0wEjR8oen6wsGYb++gt47z2gRYuyH+vrK3sahJDPTcqypNXQiYiUwiBUg4rOMF3TfvtN9gT984+c4HHXLuCpp+SH36efAs7O8uv27YH33wf+G6loYvt2Wev0xReyt2fyZGDv3sKpAcpj+HD57+rVZnlZVEmcqJOISGIQqkGGmpl//gEyM2vueb/8EujfX85sHRIiA4+Pj9xnbw+88oosYu7fXwag6dNlaNu1Sx5z9Srw8stAnz6yl+Cuu2Rh9EcfyRmzK+LJJ+Vz7t0LnDxp3tdJ5ceJOomIJAahGlS/PnD33fLrmugVEgL44AMgIgK4dUvW6PzyC9CwYfFjfXyAH3+UPTVNmshgFBwMPP+87CX6/HN53EsvAUlJpsW1FeHmBjz0kPyavULKsYaJOomIzIFBqIbVVJ3QzZsyxEybJu9PmQJ8/XXZPTgqFTB0qKwTCg+XQWrpUuDUKblY7C+/yBFm9etXrW3Dhsl/V6+Wz0HK4LplREQMQjWuJoJQbq4cFr9sGWBnJ8PLhx/Kr8ujSRNZB/TLL0CnTnLOn0OHCntyqmrwYMDRUQasffvMc04iIqLKYBCqYdUdhK5eBR58UF7WcHQENmyQl7Mq46GHZC3P8uUlX06rrPr1C+dR+vpr853Xlty8KYNkTdaaERFZIwahGmYYYZWSAly+bN5zCwG8+CKwZw/QuLEsen3kEfM+h7kYLo+tWSMXZaXi9Hrg7Fk54m/BAuC114CBA2WdmZOT/Ld5czl1wa1bSreWiKh24szS/6mJmaUN/PxkENq+XfbemMvixXKmZ3t7+eHZs6f5zm1uN24Anp7ApUuyrb17K90i5d26BSQkyPXafv0VOH4cuH699OPVasCwDnG3bsBXX8kRfUREtoQzS9dC1TGf0J49chg8IJeysOQQBAAODnIoPWDbl8cuXpQBZuhQOVz9gQeAjz+WI/OuXwfq1gVatwYGDQJefx1YtEjOzJ2WJi+DfvUV4Ooqpzpo3x5YsYIF6EREFcEeof/UZI/QzJnA1Kny8pA5QsDFi7KoOTUVePRR4Lvvasfq9n/8Iec1cnWVtS716indouonhJzh+8cf5W3XLtOlTRo1kvM5DRwoZ/329pZLmJTl7Fk5yu+PP+T9xx6TvYNNmlTf6yAishRV/fy2yNXnrZ05C6YLCuRszampQKtWwMqVtSMEAcD998s5a9LSgJ9+kh/g1uriRTmx5bJlxZcxaddOBp+BA+UlLnv7ip3bx0f2En3yiQzY338vL7GtXAn062e+10BEZI3YI/SfmuwRysyUa32pVHK256rMyzN9ulwSw9ER2L1bLoFRm0yZIi8FPfaY7MmyJno9sGOHnItpw4bCZUscHWVN1MCBwIABstfHXPbvl8HYELZeflnOAO7oWHiMELI2KzOz8HbhAnDffUBQkPnaQkRUE6r6+c0g9J+aDEKAHO1z/jywc2flZ2n+8Ufg4Yfl1199VbiOV21y4ICsbVGrgYwMoEEDpVtUdRkZsjdm+XLg778Lt3fsKCe5HDpUXg6sLteuyYD52Wfy/t13ywJ9Q+jJyip5lJlaDcTHF/ZYEhHVBiyWrqWqenksJQUYMUJ+PW5c7QxBgLws1KaNHP20YYPSrak8IeQElI89JmdlfustGYJcXORIvn375O2FF6o3BAGy92fePODnn2XP48mTcl6p/ftl+DaEoAYNAH9/WaDdtq38Hjz+uPmndSAismSsEVJIx45ymHRlgtC1a4UfWF27ytqQ2kqlkkXjU6fKwvFRo5RuUcXcuAF88w3wf/8HHD5cuD0oCBgzBnjqKcDZWZm2hYbKGcHXr5ej9Jo2Lby5u5sut3Lliiy4T06WhdebNpV/JnIiotqMl8b+U9OXxjZtkktNBAbKodIV8dxz8rJLkyYySN2+aGZtk5IiL92oVMA//8j5hSydTidHZn36qexlAWSt17PPystf99yjbPsqY/9+GeDy84FZs4DISKVbRER0Z7w0VksZZpg+cqTsSfNut3y5vNnZyZ6I2h6CAMDXV34ACwGsXat0a8qWlgZMmiTf9ylTZAjy9JRruaWlyWBUG0MQIH8mP/9cfj11qiz0JiKydgxCCtFq5TIYt26ZXlIpy4EDsh4IkCPF+vSpvvbVNMOSG5Y6ueLBg/KSkZ+fvBSZkyMDz8qVskdryhTrKPQePRoYOVKOeHv6adlDR0RkzRiEFKJSVaxguqBAXhLLz5cjxd58s3rbV9OeekrOn7N3ryzutQRXr8q5f3r0kJcwV62SwbVnTzli79AhGRqK1trUdioVMH++LGLPygKGDJELvBIRWSsGIQVVJAgtWCBDgkYDLFlifYWs7u5A377y62++UbYtSUmy561ZMyAiQk5xYG8vQ8GePfKS0YABtWfiyopycpIF1q6uQFyc9YVuIqKirOzjtHYpbxD65x/g7bfl11FRtaOYuDIMl8eio2VvRE3KyZEBs3NnWSuzYIEsiPb1BT74QM7cvWaNnHTQFhhmKQeA//2vYpNd6nQlz1NERGSJOGrsPzU9agwATp+WHzhqtfwgrlu35OOeeEJ+EHXtKie8s7beIIOcHFmDc/GinItnw4bKB4/Tp4H33gP+/Vf2cDg6Fr8Ztu/fL0NOXp58bN26cs22MWOABx+03ve7PCZNkjVRLi5yHqRWrYofI4SsX/vhB3nbt09OGdCli1wyJChI3rj2GRFVB84sbSZKBCG9HmjYUC6zcfBgyctjbNkiVx63t5c9R+3a1UjTFHP0qAwhJ0/KgLhokazDKS+9Xta4TJki51uqCH9/GX7CwwE3t4o91lrdvCnD4M6d8udz1y4ZIG/cAGJjC8NPauqdz9WyZWEoCgqSxealhX8iovJiEDITJYIQIAtvY2Pl5aCICNN9eXly1uXUVOCNN+SaXLZAp5OzZm/eLO+PGycvzzg4lP24lBQ5j09MjLz/4IPyctv16zIUXbsmC6ANXxtuDRvK57v/fuut+6mK8+fl5cKsLLk+mrOznKk6O7vwGEdH4KGHgEceAfr3lz1xCQmFt+PHSz63g4M8n7OzDFiGrw23Bg3k4ICHH2ZoIqKSMQiZiVJB6LXXgLlzgVdekXPQFDV5MjB7tlxd/MgR5WYoVoJeL6cIeOcdef/++4F16+SSEbcTQk5uOGmSDI9OTvJ9GzvWti9rmdOOHXK6Br2+cFvTpjL4PPKIXES26MKut7t0Cfjrr8JgtHu3aZC6k6ZN5R8Ko0fLtdOIiAwYhMxEqSC0apW8FHP//cCffxZuP3BALnlQUCB7RgyLq9qazZuBZ56RH5rNmslaqW7dCvefPSunFfjtN3n/gQeAFSuAu+5Spr3WbOlS2XPZs6cMP507Vz5o6vVyiZi8PNlLl5dX8i05WS4onJlZ+NgePWQgeuIJ2/rjgIhKxiBkJkoFoSNH5IKXzs7yw97OTn5IBAfLv5off1wOZbZlJ07IuqFjx+TlkfnzC5cZmThRFlk7OsrZncePZy+Qtbl5E/jpJ2DZMvmvoVfK1RUYOlT+LHTqxMuaRLaKQchMlApCt27J/9CvXZN1FP7+wMKFwEsvyZE6x44BzZvXWHMsVk6OLJr+/nt5399fBiRAhsbo6JJHNJF1+ecf4IsvZAhOTi7c3qJFYRF2t25yAsw71ZQRkXVgEDITpYIQIP/z3rULWL1aXnYICJAFw/PmAS+/XKNNsWhCyHmUpk6VX6vVcnHQCRPkqDqq/XQ6GXq9vIrvO3dO/nGg0cheodhY2Uv03XdyxvWi6tWTvUSG4fvduvEPCiJrxSBkJkoGoXHj5AR+b7whR4itXSvnz9m1ix/wJfntN9kz9MorQOvWSreGzEWnA0JD5ei0mBjTBYXT0uQfCe7ucsSaRlO4Lzu7sBB71y55u3Sp+Pl9feX0CGPGcE4jImtilavPL1y4EO3atYOrqytcXV0RFBSEn3/+uczHrFu3Dq1bt0a9evVw77334qeffqqh1ladYYbpL7+UIcjOTo6CYggqWZ8+MjgyBFmXnBwZgpKTZehJS5PbDSEoOVnuz8kxfZyrq/yZmDZNrgF38aK8bBodLUcOBgbK36mUFOCtt2Rv0+jRckACEZFFBiEvLy98+OGH2LdvH/bu3YsHH3wQYWFhOHLkSInHx8fHY+jQoRg9ejT279+PwYMHY/DgwThc3mXdFdahg/zXMDJmwoTCcERkK7y8ZE+Qn19hGIqPLwxBfn5yf0mXzYpSqeQQ+4gIWW+XlCR7m774Qv5e5efLkYXt28tzf/89lwQhsmW15tJYo0aNMHv2bIwePbrYviFDhiAvLw9btmwxbuvWrRvat2+PRYsWlev8Sl4ay8+XtQ83b8r/5I8dA+rXr9EmEFmMoj1ABoYQVPRyWWUIIS+hzZsnR2MWFMjt3t7yEvVzzwGNGlXtOYioZlnlpbGiCgoKsGbNGuTl5SEoKKjEYxISEtCnTx+Tbf369UNCQkKp583Pz0d2drbJTSlqtZxHCAA+/5whiGybVivn1ypq1aqqhyBA9hYFB8u15c6ckZfKGjeWtXlTpsg/RCZPLn75jYisl8UGoUOHDqF+/fpQq9UYO3YsNmzYgDZt2pR4bEZGBpo2bWqyrWnTpsjIyCj1/FFRUdBoNMab1hz/y1bBt9/KmoWwMEWbQaS4tDS55ElRI0YU1gyZi5cXMHOmPO+KFbKW6No1OSt5QID8nawd/eVEVBUWG4T8/f2RlJSE3bt348UXX0RERASOHj1qtvNHRkZCp9MZb2nm/l+2gpo0sf4FVYnupOhlMT8/IC7OtGaoOn5NHR2BUaOA/fvlTOa+vnK+oiFDgL59S18njYisg8UGIQcHB7Rs2RKdOnVCVFQUAgMD8enti3H9x8PDA5lF5+AHkJmZCY+SFqb6j1qtNo5KM9yISDnnzhUvjA4OLl5Afe5c9Ty/SiWXsjlyRK5xp1YD27fLP1DefFMu+UFE1sdig9Dt9Ho98m+fNe0/QUFB2L59u8m2bdu2lVpTRESWx8VFzhN0e2G0VlsYhtzd5XHVydERmDEDOHoUGDhQDmL46CN5uWz9el4uI7I2dZRuQEkiIyPRv39/eHt7IycnB6tXr0ZMTAx++eUXAEB4eDiaN2+OqKgoAMCECRMQEhKCTz75BAMHDsSaNWuwd+9eLFmyRMmXQUQVoNHIyRJLmllaq5UzSRtmlq4Jfn7Ali3yctkrr8ji6ieflJfLPvtMLvNCRLWfRfYIZWVlITw8HP7+/ujduzf27NmDX375BX379gUApKamIj093Xh8cHAwVq9ejSVLliAwMBDr16/Hxo0b0bZtW6VeAhFVgkZT+jxBXl41F4KKGjRI9g5Nny4vl23bBtx7rxxlxtFlRLVfrZlHqLopOY8QEdUOp0/LCU8NE9c3awZ8/DEwbJisMSKimmf18wgREVmKli0LL5fddRdw/jzwzDPAAw/IGayJqPZhECIiqgDD6LLDh+U8RE5OwM6dcrX7l14C/v1X6RYSUUUwCBERVUK9enJm6uPH5ZxDer1c2+zuu4FFiwqX7yAiy8YgRERUBVqtXLJjxw5ZRH3pEvDii0Dnzlzhnqg2YBAiIjKDnj2BxES5oGuDBnKm6qCg4uumEZFlYRAiIjKTOnWAl18GTp4EQkPl2mXh4XJl+xs3lG4dEZWEQYiIyMzc3OToshkz5P0FC4CQkOpbHoSIKo9BiIioGtjbyzXLtmyRl8p27QI6dpS1RERkORiEiIiq0cCBwL59QPv2wIULQJ8+chJGTmVLZBkYhIiIqpmfHxAfD0REyGH2U6YATzwBZGcr3TIiYhAiIqoBjo7AypVyjiEHB+D77+UQ+0OHlG4ZkW1jECIiqiEqFfDCC8Cff8r5h06elHVDkyYBOp3SrSOyTQxCREQ1rEsXWTf0yCPArVvAJ5/IGalXrpSXzoio5jAIEREpwM0N2LQJ+PlnwN8fyMoCnn0W6NYN2L1b6dYR2Q4GISIiBYWGAgcPAv/3f4CLC7BnjwxDI0cC6elKt47I+jEIEREpzMEBeP11WTM0apTc9sUX8nLZ7NmclZqoOjEIERFZCA8PYMUKeWmsa1cgNxeYPBlo2xbYsIFzDxFVBwYhIiIL06WLnHcoOhpo2hQ4dQp47DGgRw8gIUHp1hFZFwYhIiILZGcnJ2A8eRKYOlXOQxQXBwQHy8kYT51SuoVE1oFBiIjIgrm6Au+/L4PP6NEyIH33HdCmjVzp/sIFpVtIVLsxCBER1QLNmwPLlgEHDgADBsj5hz7/HLjrLmDWLODqVaVbSFQ7MQgREdUibdsCP/4IbN8uZ6XOyQHefhto1QpYvJgjzIgqikGIiKgWevBBOefQ118DPj7A+fPA2LFyyP3SpQxEROXFIEREVEvZ2QHDhgEnTgCffiqH3589Czz/vJytetky4OZNpVtJZNkYhIiIajm1GnjlFSA5GZg7VwaiM2eAMWNkDxEDEVHpGISIiKyEoyMwYYIMRHPmyDmIigai5csZiIhuxyBERGRlHB2BV1+Vgeh//ysMRM89B9xzD/D770q3kMhyMAgREVkpJyfgtddkIPrkE8DdXc5H1Lu3XOn+33+VbiGR8hiEiIisnJMTMHGiDEEvvQSoVMDKlUBAAPDNN1zDjGwbgxARkY1wdQXmzwd27pQzU1+4IEedDRwoR5sR2SIGISIiGxMcDOzfD7z3HuDgAPz8swxGc+YABQVKt46oZjEIERHZIAcHYNo0uWRHjx5yiY6JE4Fu3eQ2IlvBIEREZMNatwZiYuTyHBoNsHcv0KmTrCXKyFC6dUTVzyKDUFRUFDp37gwXFxe4u7tj8ODBOHHiRJmPiY6OhkqlMrnVq1evhlpMRFR72dnJ2aiPHQMef1xeHlu4UC7oOnUqoNMp3UKi6mORQSg2Nhbjxo3Drl27sG3bNty8eRMPPfQQ8vLyynycq6sr0tPTjbezrP4jslk6HXDuXMn7zp3jh3tJPD2B9etlD1HXrvJy2cyZMhD973/A9etKt5DI/FRCWP7AyQsXLsDd3R2xsbF44IEHSjwmOjoar776Kq5cuVKp58jOzoZGo4FOp4Orq2sVWktEStPpgNBQICtLfqhrtYX70tKAnj3lnDpbt8rLQVScEMDGjcBbbwHHj8tt3t6ywPqZZwB7e0WbR2RU1c9vi+wRup3uvz/dGjVqVOZxubm58PHxgVarRVhYGI4cOVLqsfn5+cjOzja5EZF1yMmRISg5WYaetDS53RCCkpPl/pwcJVtp2VQq4NFHgUOH5Gr2zZsDqanAyJFA+/bA5s2cf4isg8UHIb1ej1dffRXdu3dH27ZtSz3O398fK1aswKZNm/DVV19Br9cjODgY50rpG4+KioJGozHetEX/ZCSiWs3LS/YE+fkVhqH4+MIQ5Ocn93t5KdvO2qBOHbk0x6lTwMcfAw0bAocPA488Ikeb/fGH0i0kqhqLvzT24osv4ueff8bOnTvhVYH/tW7evImAgAAMHToU77//frH9+fn5yM/PN97Pzs6GVqvlpTEiK1K0B8jAEIL4t0/lXL4sA9GnnwLXrsltoaHArFlAhw7Kto1sk1VfGhs/fjy2bNmCHTt2VCgEAUDdunXRoUMHnD59usT9arUarq6uJjcisi5aLbBqlem2VasYgqqiYUMgKgo4fRoYO1b2GG3dCnTsCAwZApw8qXQLiSrGIoOQEALjx4/Hhg0b8Pvvv8PX17fC5ygoKMChQ4fg6elZDS0kotogLQ0YMcJ024gRhTVDVHnNmskh9seOyWU6VCrg22/lDNVjxpQ+Yo/I0lhkEBo3bhy++uorrF69Gi4uLsjIyEBGRgauGfphAYSHhyMyMtJ4/7333sOvv/6K5ORkJCYm4plnnsHZs2fx3HPPKfESiEhhRS+L+fkBcXGmNUMMQ+bRsiXw9ddAUhLw8MNyDqJly+T2118HLl5UuoVEZbPIILRw4ULodDr07NkTnp6extvatWuNx6SmpiI9Pd14//LlyxgzZgwCAgIwYMAAZGdnIz4+Hm3atFHiJRCRgs6dK14YHRxcvICavRbm066dHEm2cyfwwANAfr6ce+iuu4DPPuMaZmS5LL5YuqZwHiEi68F5hJQlBPDrr0BkpFzcFZDLdixeLP8lMqeqfn4zCP2HQYjIuuh0cp6gksZZnDsHuLgwBFU3vR5YsgR48035/bCzA8aPB95/H+B/s2QuVj1qjIiosjSa0ucJ8vJiCKoJdnZyZNnx47KgWq8H5s0DAgLkUh78M5wsAYMQERFVKw8PWVD966+yiPr8eeDJJ2VxdUqK0q0jW8cgRERENaJvX7lkx7RpQN26wE8/AffcA3z4IXDjhtKtI1vFIERERDWmXj25cOvBg7Jo/do1WVTdti2waRMvl1HNYxAiIqIa17o18PvvwBdfyBF8p04BgwcDvXoBiYlKt45sCYMQEREpQqUCwsNlCIqMBNRqIDYWuO8+ucr9P/8o3UKyBQxCRESkKFdXuWjryZPA8OHy8tgXXwCtWgEzZgC5uUq3kKwZgxAREVkEb2/gq6+A3buB+++X9UPvvQfcfTewYgVnp6bqwSBEREQWpUsX4I8/5FxDfn5AejowerRc4f7nn1lQTebFIERERBZHpQIefxw4ehT45BOgQQM50mzAAODBB4G//lK6hWQtGISIiMhiqdXAxInA338DkybJ+zExQNeuwBNPACdOKN1Cqu0YhIiIyOI1agTMni0LqkeOlD1G330nJ2QcO1ZePiOqDAYhIiKqNby9gZUr5WWyQYNkAfXixcBddwFvvy0XdyWqCAYhIiKqddq2BX74QRZVBwXJEWazZsni6tmzgatXlW4h1RYMQkREVGv16AHExQEbN8pV7S9dAiZPloHos8+A/HylW0iWjkGIiIhqNZUKCAuTl8uiowFfXyAzE3jlFbna/ZIlwM2bSreSLBWDEBERWYU6dYCICOD4cWDRIqB5c+DcOeCFF+TaZl9+yUkZqTgGISIisioODjL8nD4NzJ0rF3VNTpYhqW1bYO1aQK9XupVkKRiEiIjIKtWrB0yYIEPQhx/KIfjHjwNPPw0EBgLffstARAxCRERk5ZydgSlTgJQU4N135SKvhw8DQ4YA994LfPMNL5nZMgYhIiKyCa6uwPTpwJkzclV7jUYu4TFsmJyY8auvgFu3lG4l1TQGISIisikNGwLvvAOcPStXt2/YUC7VMWIE0KaNLKpmILIdDEJERGSTNBpg2jTZQzRrFtC4MXDqlCyqbt0aWLECuHFD6VZSdWMQIiIim+bqCkRGyhqijz4CmjSRi7yOHi3nIfrsM85Ubc0YhIiIiAC4uMhZqc+ckct0eHgAaWlyYsYWLYCoKK5lZo0YhIiIiIpwdgYmTZI9RAsXyhB04QLw1lty0de33gKyspRuJZlLpYPQ9OnTkZKSYs62EBERWYx69YCxY2Xd0FdfyZFl2dmyZ8jHB3j5ZVlwTbVbpYPQBx98gISEBHO2hYiIyOLUqQMMHy7XMtu4EejSBbh+Hfj8c1lDFB4OHDigdCupsqrt0tjq1avh6+tbXacnIiKqUXZ2cnHXXbuA7duB3r3lMPtVq4D27YGHHgJ+/RUQQumWUkVUKAjNmzcPjzzyCGbNmgUAyMvLK/XYgoICpKamVq11REREFkalAh58EPjtN2DPHjlDtZ0dsG0b0K+fXL7jyy859L62qFAQcnV1xa5duzB16lSoVCqMHTsWHh4e6N+/P95++2189913SElJQUFBARISEtCkSZPqajcREZHi7rsPWLNGDrefMEEWWh86JOci8vUFPv4YuHJF6VZSWSoUhEaOHImsrCycPHkSQggMHDgQ3bp1w5EjRxAVFYUnn3wSLVu2hIODAxYvXoz+/ftXqlFRUVHo3LkzXFxc4O7ujsGDB+PEiRN3fNy6devQunVr1KtXD/feey9++umnSj0/ERFRRbRoIVe6T0uTxdSensD583KNM60WePVVufgrWR6VEJW7mvnoo4/ihRdeQGhoKADg4sWLSExMxP79+5GcnAw/Pz+88sorcHR0rPC5Q0ND8fTTT6Nz5864desW3nrrLRw+fBhHjx6Fs7NziY+Jj4/HAw88gKioKDz88MNYvXo1PvroIyQmJqJt27Z3fM7s7GxoNBrodDq4urpWuM1EREQG+flyMdf/+z/gyBG5TaUCHnlE9hz17CnvU9VV9fO70kGoJl24cAHu7u6IjY3FAw88UOIxQ4YMQV5eHrZs2WLc1q1bN7Rv3x6LFi2643MwCBERkbkJAfzyi+wt+uWXwu3t2slANGyYHKZPlVfVz+9aMaGi7r+pPBs1alTqMQkJCejTp4/Jtn79+pU6xD8/Px/Z2dkmNyIiInNSqYDQUGDrVrnS/YsvAk5Ocij+6NHystm0afIyGimjQkFIq9Xi5Zdfxvbt21FQUFBdbTKh1+vx6quvonv37mVe4srIyEDTpk1NtjVt2hQZGRklHh8VFQWNRmO8abVas7abiIioqIAAYMEC4Nw5WUTt7Q1cvAh88IGcoHH4cCAhgcPva1qFglBYWBg2btyIvn37wt3dHSNGjMCGDRtwtRpXoxs3bhwOHz6MNWvWmPW8kZGR0Ol0xltaWppZz09ERFSShg2BN96QI83WrQPuv1/OR7R6NRAcDHTqBCxfzoVea0qFgtDnn3+OtLQ07Nq1C88//zz27t2Lxx9/HG5ubggLC0N0dDT+/fdfszVu/Pjx2LJlC3bs2AEvL68yj/Xw8EBmZqbJtszMTHh4eJR4vFqthqurq8mNiIioptSpAzzxBPDnn8DevXLIvVoN7N8PPPcc4OUFvP46cPq00i21blUulj5+/Dg2bNiAjRs3Yu/evbCzs8P999+PRx99FIMHD4a3t3eFzymEwMsvv4wNGzYgJiYGrVq1uuNjhgwZgqtXr2Lz5s3GbcHBwWjXrh2LpYmIqFa4eBFYuVIu9lp0Oc/QUOCll4ABAwB7e+XaZ4ksatTYP//8YwxFf/zxBwoKChAYGIhZs2YZh9mXx0svvYTVq1dj06ZN8Pf3N27XaDTG4fjh4eFo3rw5oqKiAMjh8yEhIfjwww8xcOBArFmzBrNmzeLweSIiqnUKCmSB9fz58l/DJ3WLFsDzzwOjRgGlXPCwORY1aqx58+Zo3749fvvtN2RmZmLlypVo0aIFDh8+XKHzLFy4EDqdDj179oSnp6fxtnbtWuMxqampSE9PN94PDg7G6tWrsWTJEgQGBmL9+vXYuHFjuUIQERGRJbG3BwYOBH76CTh1Sl4ia9gQOHMGeOstOdrsiSfk2mZ6vdKtrd3MPo+Qq6triUPRDZMsWir2CBERkSW7dg1YuxZYvFgu/Grg6ytrip591jZ7iSymR2j16tU4evRoqfvbt29vrqciIiKyOY6OwMiRcoj9gQPA+PGARiNrid5+W/YSPf64nLiRvUTlZ7YeoaCgIBw6dAjXr19Hx44dTW4A0LdvX1y+fNkcT1Ut2CNERES1zdWrcgj+kiVAfHzhdh8fWUcUESHriqyZRRVL6/V6uLi4YMWKFUhMTERiYiKSkpKg0+kwadIkzJo1y1xPZXYMQkREVJsdPgwsXQp8+aXpive9e8tQ9NhjslfJ2lhUEALk3D23z/AshIDKwleXYxAiIiJrcO0asGEDsGIFsH174XaNBnj6aVlL1Lmz9Sz6anFBqLZiECIiImtz5gzwxRdybqKzZwu3t2kjA9Hw4bW/wLrGi6U3btyIJ598Eg899BDeffdd5OTklHjczp078frrr1e4QURERGQeLVoAM2YAycmyd2j4cLna/dGjwKRJcvbqAQOAb76x3SU9KtQjtHbtWgwbNgxFHxIQEID4+HhoNBrs2LED33zzDX744QdcuHABAGpscdaqYo8QERHZgitX5DD8lSuB3bsLt7u4yLmJwsOBBx4A7Mw602D1qdFLY127dsW5c+fw9ddf46677sLmzZvxxhtv4JVXXkFCQgL+/PNPCCHQrFkzDBw4EA8//DAGDRpU4UYpgUGIiIhszcmTwFdfAatWyctoBt7esvdoxAggIECx5pVLjQahBg0a4M0338Sbb75p3PbRRx9h2rRp0Ov1GD16NJ5//nl06tSpwg1RGoMQERHZKr0eiIuTI86+/RYoOi9yp07AsGHAkCFA8+bKtbE0NVojlJ2dDa1Wa7ItLCwMt27dwuTJk7F48eJaGYKIiIhsmZ0d0KOHHH6fkSEvnQ0cKJf62LdPLvGh1QK9esljLl1SusXmU+ErgLcPg3dzcwMA9OjRwzwtIiIiIsU4OgJPPQVs2QKcPw98/jkQHCwXfo2JkYu+engAjzwii6zz8pRucdVUOAj9+OOP2LJlC/755x+T7Wq12myNIiIiIuW5uwPjxsnLZikpQFQU0K4dcPMmsHmzvGTm7i7/3bQJuH5d6RZXXIVqhOz+KyE39Aq5ubmhTZs2iI2NxUcffYRhw4ahWbNm1dPSasYaISIiovI5ckT2Bq1eLQOSgYsLEBYme5QeegioiT6SGi2WzsnJMS6dsW/fPiQmJuLkyZPQ6/XGcNSkSRMEBgaiffv26NChA4YOHVrhRimBQYiIiKhihAD++gtYs0aueVb0YpGrKzB4sAxFffsCDg7V0wbFZ5bOy8tDUlKSMRjt27cPx48fR0FBAVQqFecRIiIisgF6PbBrlxx1tm6drC8y0GgKQ1H//uZd3kPxIFSSa9eu4cCBA0hMTMRLL71k7tNXCwYhIiIi89Drgfh4GYrWrwfS0+X2u+8Gjh+3gSBUGzEIERERmV9BgSy2XrcO8PMDXnvNvOev6ud3HfM2h4iIiKiQvb1csuOBB5RuSclqyUoiRERERObHIEREREQ2i0GIiIiIbBaDEBEREdksBiEiIiKyWQxCREREZLMYhIiIiMhmMQgRERGRzWIQIiIiIpvFIEREREQ2i0GIiIiIbBaDEBEREdksBiEiIiKyWQxCREREZLMYhIiIiMhmWWQQ+uOPPzBo0CA0a9YMKpUKGzduLPP4mJgYqFSqYreMjIyaaTARERHVShYZhPLy8hAYGIj58+dX6HEnTpxAenq68ebu7l5NLSQiIiJrUEfpBpSkf//+6N+/f4Uf5+7ujgYNGpi/QURERGSVLLJHqLLat28PT09P9O3bF3FxcWUem5+fj+zsbJMbERER2RarCEKenp5YtGgRvvvuO3z33XfQarXo2bMnEhMTS31MVFQUNBqN8abVamuwxURERGQJVEIIoXQjyqJSqbBhwwYMHjy4Qo8LCQmBt7c3Vq1aVeL+/Px85OfnG+9nZ2dDq9VCp9PB1dW1Kk0mIiKiGpKdnQ2NRlPpz2+LrBEyhy5dumDnzp2l7ler1VCr1TXYIiIiIrI0VnFprCRJSUnw9PRUuhlERERkwSyyRyg3NxenT5823k9JSUFSUhIaNWoEb29vREZG4p9//sGXX34JAJg7dy58fX1xzz334Pr161i2bBl+//13/Prrr0q9BCIiIqoFLDII7d27F7169TLenzhxIgAgIiIC0dHRSE9PR2pqqnH/jRs38Prrr+Off/6Bk5MT2rVrh99++83kHERERES3s/hi6ZpS1WIrIiIiqnlV/fy22hohIiIiojthECIiIiKbxSBERERENotBiIiIiGwWgxARERHZLAYhIiIislkMQkRERGSzGISIiIjIZjEIERERkc1iECIiIiKbxSBERHQHOh1w7lzJ+86dk/uJqHZiECIiKoNOB4SGAiEhQFqa6b60NLk9NJRhiKi2YhAiIipDTg6QlQUkJwM9exaGobQ0eT85We7PyVGylURUWQxCRERl8PICYmIAP7/CMBQfXxiC/Pzkfi8vZdtJRJVTR+kGEBFZOq1Whh1D+OneXW43hCCtVsHGEVGVsEeIiKgctFpg1SrTbatWMQQR1XYMQkRE5ZCWBowYYbptxIjiBdREVLswCBER3UHRwmg/PyAuzrRmiGGIqPZiECIiKsO5c8ULo4ODixdQlzbPEBFZNhZLExGVwcUFcHeXXxctjC5aQO3uLo8jotqHQYiIqAwaDbB1q5wn6PYh8lotEBsrQ5BGo0z7iKhqGISIiO5Aoyk96HD+IKLajTVCREREZLMYhIiIiMhmMQgRERGRzWIQIiIiIpvFIEREREQ2i0GIiIiIbBaDEBEREdksBiEiIiKyWQxCREREZLMYhIiIiMhmWWQQ+uOPPzBo0CA0a9YMKpUKGzduvONjYmJi0LFjR6jVarRs2RLR0dHV3k4iIiKq3SwyCOXl5SEwMBDz588v1/EpKSkYOHAgevXqhaSkJLz66qt47rnn8Msvv1RzS4mIiKg2s8hFV/v374/+/fuX+/hFixbB19cXn3zyCQAgICAAO3fuxJw5c9CvX7/qaiYRERHVchbZI1RRCQkJ6NOnj8m2fv36ISEhodTH5OfnIzs72+RGREREtsUqglBGRgaaNm1qsq1p06bIzs7GtWvXSnxMVFQUNBqN8abVamuiqURERGRBrCIIVUZkZCR0Op3xlpaWpnSTiIiIqIZZZI1QRXl4eCAzM9NkW2ZmJlxdXeHo6FjiY9RqNdRqdU00j4iIiCyUVfQIBQUFYfv27Sbbtm3bhqCgIIVaRERERLWBRQah3NxcJCUlISkpCYAcHp+UlITU1FQA8rJWeHi48fixY8ciOTkZkydPxvHjx7FgwQJ8++23eO2115RoPhEREdUSFhmE9u7diw4dOqBDhw4AgIkTJ6JDhw6YPn06ACA9Pd0YigDA19cXP/74I7Zt24bAwEB88sknWLZsGYfOE5FV0+mAc+dK3nfunNxPRGVTCSGE0o2wBNnZ2dBoNNDpdHB1dVW6OUREZdLpgNBQICsLiIkBig58TUsDevYE3N2BrVsBjUapVhJVv6p+fltkjxAREZUtJ0eGoORkGXoMA18NISg5We7PyVGylUSWj0GIiKgW8vKSPUF+foVhKD6+MAT5+cn9Xl7KtpPI0lnF8HkiIluk1cqwYwg/3bvL7YYQxHliie6MPUJERLWYVgusWmW6bdUqhiCi8mIQIiKqxdLSgBEjTLeNGFFYM0REZWMQIiKqpYoWRvv5AXFxpjVDDENEd8YgRERUC507V7wwOji4eAF1afMMEZHEYmkiolrIxUXOEwSYFkYXLaB2d5fHEVHpGISIiGohjUZOlpiTU3yIvFYLxMbKEMTJFInKxiBERFRLaTSlBx3OH0RUPqwRIiIiIpvFIEREREQ2i0GIiIiIbBaDEBEREdksBiFSnE5X+lwn587J/URERNWBQYgUpdMBoaFASEjxWXDT0uT20FCGISIiqh4MQqSonBwgK6v4kgBFlw7IypLHERERmRuDECnKy6v4kgDx8cWXDuCcKEREVB04oSIpruiSAMnJQPfucrshBBmWDiAiIjI39giRRdBqgVWrTLetWsUQRERE1YtBiCxCWhowYoTpthEjihdQExERmRODECmuaGG0nx8QF2daM8QwRERE1YVBiBR17lzxwujg4OIF1KXNM0RERFQVLJYmRbm4AO7u8uuihdFFC6jd3eVxRERE5sYgRIrSaICtW+U8QbcPkddqgdhYGYI0GmXaR0RE1o1BiBSn0ZQedDh/EBERVSfWCBEREZHNYhAiIiIim8UgRERERDaLQYiIiIhsFoMQERER2SwGISIiIrJZDEJERGQ2Ol3pM8GfOyf3E1kSiw5C8+fPR4sWLVCvXj107doVf/31V6nHRkdHQ6VSmdzq1atXg60lIrJtOh0QGgqEhBRfIzAtTW4PDWUYIstisUFo7dq1mDhxImbMmIHExEQEBgaiX79+yMrKKvUxrq6uSE9PN97Onj1bgy0mIrJtOTlAVlbxBZOLLqyclSWPI7IUFhuE/ve//2HMmDEYNWoU2rRpg0WLFsHJyQkrVqwo9TEqlQoeHh7GW9OmTWuwxUREts3Lq/iCyfHxxRdW5ozxZEksMgjduHED+/btQ58+fYzb7Ozs0KdPHyQkJJT6uNzcXPj4+ECr1SIsLAxHjhwp9dj8/HxkZ2eb3IiIqGoMCyYbwlD37qYhyLCwMpGlsMggdPHiRRQUFBTr0WnatCkyMjJKfIy/vz9WrFiBTZs24auvvoJer0dwcDDOlVK1FxUVBY1GY7xp+dtJRGQWWi2wapXptlWrGILIMllkEKqMoKAghIeHo3379ggJCcH3338PNzc3LF68uMTjIyMjodPpjLe02yv7iIiqgS2MqkpLA0aMMN02YkTxAmoiS2CRQahJkyawt7dHZmamyfbMzEx4eHiU6xx169ZFhw4dcPr06RL3q9VquLq6mtyIiKqTLYyqKloY7ecHxMWZ1gwxDJGlscgg5ODggE6dOmH79u3GbXq9Htu3b0dQUFC5zlFQUIBDhw7B09OzuppJRFQh1j6q6ty54oXRwcHFC6hL6xEjUoJFBiEAmDhxIpYuXYovvvgCx44dw4svvoi8vDyMGjUKABAeHo7IyEjj8e+99x5+/fVXJCcnIzExEc888wzOnj2L5557TqmXQERkwtpHVbm4AO7uxQujixZQu7vL44gsRR2lG1CaIUOG4MKFC5g+fToyMjLQvn17bN261VhAnZqaCju7whx3+fJljBkzBhkZGWjYsCE6deqE+Ph4tGnTRqmXQERUjCEUGMJP9+5yuzWMqtJogK1bZY/W7WFOqwViY2UI0miUaR9RSVRCCKF0IyxBdnY2NBoNdDod64WIqNrFxxeGIEDW0gQHK9ceotqqqp/fFntpjIjIWnFUFZHlYBAiIqpBljiqyhaG9BOVhkGIiKiGWOKoKlsY0k9UFgYhIqIaYomjqqx9SD/RnbBY+j8sliaimqDTlTyqCpA9QUqMqrr9ct2qVbJmiWuEUW1Q1c9vBqH/MAgRkS0rGoYMGIKoNuCoMSIiqjIulEq2ikGIiIg4pJ9sFoMQEZGNs8Qh/UQ1hUGIiMiGWeKQfqKaZLFrjRERUfUzDOkHSh7S37MnF0ol68YgRERkw7hQKtk6BiEiIhun0ZQedEqa74jImrBGiIiIiGwWgxARERHZLAYhIiIislkMQkRERGSzGISIiIjIZjEIERERkc1iECIiIiKbxSBERERENotBiIiILJJOV/oaZ+fOyf1EVcUgREREFkenA0JDgZAQIC3NdF9amtweGsowRFXHIERERBYnJwfIygKSk+XCr4YwlJYm7ycny/05OUq2kqwBgxAREVkcLy8gJgbw8ysMQ/HxhSHIz0/u51poVFUMQkREZJG0WtMw1L27aQjSast/LtYbUWkYhIiIyGJptcCqVabbVq2qeAhivRGVhkGIiIgsVloaMGKE6bYRI4oHmrKw3ojKwiBEREQWqWhQ8fMD4uJMa4bKG4ZYb0RlYRAiIiKLc+5c8aASHFw80JRW93M71htRaRiEiIjI4ri4AO7uxYNK0UDj7i6PKy9LqzdioLIMDEJERGRxNBpg61YgNrZ4UNFq5fatW+Vx5WVJ9UbmLuA2Z6iy1HNVF4sOQvPnz0eLFi1Qr149dO3aFX/99VeZx69btw6tW7dGvXr1cO+99+Knn36qoZYSEZG5aTSl1+14eVU8BFlSvZE5C7jN3UtlieeqVsJCrVmzRjg4OIgVK1aII0eOiDFjxogGDRqIzMzMEo+Pi4sT9vb24uOPPxZHjx4VU6dOFXXr1hWHDh0q1/PpdDoBQOh0OnO+DCIiUlhamhB+fkIA8t/UVLk9NdV0e1pa+c9Z9LGGW9FzV/Qcfn5CxMWV3M6afH2Weq6yVPXz22KDUJcuXcS4ceOM9wsKCkSzZs1EVFRUicc/9dRTYuDAgSbbunbtKl544YVyPR+DEBGRdbpyRYhu3UoOF4YP5W7d5HEVERdnGoTi4ireNnMEqtvPU5VQZcnnKk1VP78t8tLYjRs3sG/fPvTp08e4zc7ODn369EFCQkKJj0lISDA5HgD69etX6vH5+fnIzs42uRERkfWx1Hojw/NXtYDbcB5zjYqz1HNVF4sMQhcvXkRBQQGaNm1qsr1p06bIyMgo8TEZGRkVOj4qKgoajcZ401rCd4OIiKqFJdYbGc5ljkAFmC9UWfK5qoNFBqGaEBkZCZ1OZ7ylVeanjoiIbIo55zcyZ6AynM9cocpSz1UdLDIINWnSBPb29sjMzDTZnpmZCQ8PjxIf4+HhUaHj1Wo1XF1dTW5ERERlMdf8RuaeMNLcvVSWeK7qYpFByMHBAZ06dcL27duN2/R6PbZv346goKASHxMUFGRyPABs27at1OOJiIgqylz1RuacMNKcocpSz1Wd6ij79KWbOHEiIiIicN9996FLly6YO3cu8vLyMGrUKABAeHg4mjdvjqioKADAhAkTEBISgk8++QQDBw7EmjVrsHfvXixZskTJl0FERFZGoyk96JR3vTJDoMrJKf4YQ6BycSlf7ZIhVAElh6qePcsfqiz1XNVJJYQQyjahdJ9//jlmz56NjIwMtG/fHvPmzUPXrl0BAD179kSLFi0QHR1tPH7dunWYOnUqzpw5g1atWuHjjz/GgAEDyvVc2dnZ0Gg00Ol0vExGRES1ik5XcqgCZI9LeUOVJZ+rNFX9/LboIFSTGISIiIhqn6p+fltkjRARERFRTWAQIiIiIpvFIEREREQ2i0GIiIiIbBaDEBEREdksBiEiIiKyWQxCREREZLMYhIiIiMhmMQgRERGRzbLYtcZqmmGC7ezsbIVbQkREROVl+Nyu7EIZDEL/ycnJAQBob19OmIiIiCxeTk4ONJVYuIxrjf1Hr9fj/PnzcHFxgUqlMuu5s7OzodVqkZaWxnXMahDfd2XwfVcG33dl8H1XRtH33cXFBTk5OWjWrBns7Cpe8cMeof/Y2dnBq6Tlcc3I1dWVvygK4PuuDL7vyuD7rgy+78owvO+V6QkyYLE0ERER2SwGISIiIrJZDEI1QK1WY8aMGVCr1Uo3xabwfVcG33dl8H1XBt93ZZjzfWexNBEREdks9ggRERGRzWIQIiIiIpvFIEREREQ2i0GIiIiIbBaDUDWbP38+WrRogXr16qFr167466+/lG6S1fnjjz8waNAgNGvWDCqVChs3bjTZL4TA9OnT4enpCUdHR/Tp0wenTp1SprFWIioqCp07d4aLiwvc3d0xePBgnDhxwuSY69evY9y4cWjcuDHq16+Pxx9/HJmZmQq12DosXLgQ7dq1M04iFxQUhJ9//tm4n+95zfjwww+hUqnw6quvGrfxvTe/d955ByqVyuTWunVr435zvecMQtVo7dq1mDhxImbMmIHExEQEBgaiX79+yMrKUrppViUvLw+BgYGYP39+ifs//vhjzJs3D4sWLcLu3bvh7OyMfv364fr16zXcUusRGxuLcePGYdeuXdi2bRtu3ryJhx56CHl5ecZjXnvtNWzevBnr1q1DbGwszp8/j8cee0zBVtd+Xl5e+PDDD7Fv3z7s3bsXDz74IMLCwnDkyBEAfM9rwp49e7B48WK0a9fOZDvf++pxzz33ID093XjbuXOncZ/Z3nNB1aZLly5i3LhxxvsFBQWiWbNmIioqSsFWWTcAYsOGDcb7er1eeHh4iNmzZxu3XblyRajVavHNN98o0ELrlJWVJQCI2NhYIYR8j+vWrSvWrVtnPObYsWMCgEhISFCqmVapYcOGYtmyZXzPa0BOTo5o1aqV2LZtmwgJCRETJkwQQvDnvbrMmDFDBAYGlrjPnO85e4SqyY0bN7Bv3z706dPHuM3Ozg59+vRBQkKCgi2zLSkpKcjIyDD5Pmg0GnTt2pXfBzPS6XQAgEaNGgEA9u3bh5s3b5q8761bt4a3tzffdzMpKCjAmjVrkJeXh6CgIL7nNWDcuHEYOHCgyXsM8Oe9Op06dQrNmjWDn58fhg8fjtTUVADmfc+56Go1uXjxIgoKCtC0aVOT7U2bNsXx48cVapXtycjIAIASvw+GfVQ1er0er776Krp37462bdsCkO+7g4MDGjRoYHIs3/eqO3ToEIKCgnD9+nXUr18fGzZsQJs2bZCUlMT3vBqtWbMGiYmJ2LNnT7F9/HmvHl27dkV0dDT8/f2Rnp6Od999Fz169MDhw4fN+p4zCBFRlYwbNw6HDx82uXZP1cff3x9JSUnQ6XRYv349IiIiEBsbq3SzrFpaWhomTJiAbdu2oV69eko3x2b079/f+HW7du3QtWtX+Pj44Ntvv4Wjo6PZnoeXxqpJkyZNYG9vX6yCPTMzEx4eHgq1yvYY3mt+H6rH+PHjsWXLFuzYsQNeXl7G7R4eHrhx4wauXLlicjzf96pzcHBAy5Yt0alTJ0RFRSEwMBCffvop3/NqtG/fPmRlZaFjx46oU6cO6tSpg9jYWMybNw916tRB06ZN+d7XgAYNGuDuu+/G6dOnzfrzziBUTRwcHNCpUyds377duE2v12P79u0ICgpSsGW2xdfXFx4eHibfh+zsbOzevZvfhyoQQmD8+PHYsGEDfv/9d/j6+prs79SpE+rWrWvyvp84cQKpqal8381Mr9cjPz+f73k16t27Nw4dOoSkpCTj7b777sPw4cONX/O9r365ubn4+++/4enpad6f9yoUdNMdrFmzRqjVahEdHS2OHj0qnn/+edGgQQORkZGhdNOsSk5Ojti/f7/Yv3+/ACD+97//if3794uzZ88KIYT48MMPRYMGDcSmTZvEwYMHRVhYmPD19RXXrl1TuOW114svvig0Go2IiYkR6enpxtvVq1eNx4wdO1Z4e3uL33//Xezdu1cEBQWJoKAgBVtd+7355psiNjZWpKSkiIMHD4o333xTqFQq8euvvwoh+J7XpKKjxoTge18dXn/9dRETEyNSUlJEXFyc6NOnj2jSpInIysoSQpjvPWcQqmafffaZ8Pb2Fg4ODqJLly5i165dSjfJ6uzYsUMAKHaLiIgQQsgh9NOmTRNNmzYVarVa9O7dW5w4cULZRtdyJb3fAMTKlSuNx1y7dk289NJLomHDhsLJyUk8+uijIj09XblGW4Fnn31W+Pj4CAcHB+Hm5iZ69+5tDEFC8D2vSbcHIb735jdkyBDh6ekpHBwcRPPmzcWQIUPE6dOnjfvN9Z6rhBDCDD1WRERERLUOa4SIiIjIZjEIERERkc1iECIiIiKbxSBERERENotBiIiIiGwWgxARERHZLAYhIiIislkMQkRERGSzGISIiIjIZjEIEZFVS0pKwjvvvIMzZ84o3RQiskAMQkRk1ZKSkvDuu+8yCBFRiRiEiIiIyGYxCBGR1XrnnXcwatQoAECvXr2gUqmgUqkwcuRIZRtGRBajjtINICKqLo899hjS09OxZMkSvPXWWwgICAAA3HXXXQq3jIgshUoIIZRuBBFRdYmOjsaoUaOwY8cO9OzZU+nmEJGF4aUxIiIislkMQkRERGSzGISIiIjIZjEIEZFVU6lUSjeBiCwYgxARWbX69esDAC5duqRwS4jIEnH4PBFZtc6dO8POzg4zZ87E5cuX4ezsDF9fX3Tt2lXpphGRBWCPEBFZNW9vb6xYsQLXrl3Diy++iKFDh2LhwoVKN4uILATnESIiIiKbxR4hIiIislkMQkRERGSzGISIiIjIZjEIERERkc1iECIiIiKbxSBERERENotBiIiIiGwWgxARERHZLAYhIiIislkMQkRERGSzGISIiIjIZjEIERERkc36fzcta+ZBrJE/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_cumulative_regret(\n",
    "    obs_values,\n",
    "    ax,\n",
    "    label,\n",
    "    obs_every = 2,\n",
    "    show_obs=False,\n",
    "    c = \"blue\",\n",
    "):\n",
    "    cum_values = np.cumsum(obs_values)\n",
    "    cum_values_t = [cum_values[i]/(i+1) for i in range(len(obs_values))]\n",
    "    ax.plot(cum_values_t, color=c, label=label + \" (Avg. Cumulative)\")\n",
    "        \n",
    "    if show_obs:\n",
    "        obs_values = np.array(obs_values)\n",
    "        ax.scatter(range(0,obs_values.shape[0], obs_every), obs_values[0::obs_every], c=c, marker='x', label=label+\" (Instantaneous)\")\n",
    "\n",
    "suboptimality = []\n",
    "obs_obj_values = observations[:, 0]\n",
    "for i in range(obs_obj_values.shape[0]-1):\n",
    "    curr_suboptimality = max_obj_value + obs_obj_values[i]\n",
    "    suboptimality += [curr_suboptimality]\n",
    "    \n",
    "_, ax = plt.subplots(1, 1)\n",
    "\n",
    "plot_cumulative_regret(\n",
    "    suboptimality, ax, label = \"Regret\", show_obs = True\n",
    ")\n",
    "\n",
    "ax.legend(loc='upper right')\n",
    "ax.set_xlabel(\"t\",fontsize = 13)\n",
    "ax.set_ylabel(r\"$R_t/t$\",fontsize = 13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dd614d7-63ff-445e-a38a-47349c557039",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}