diff --git a/notebooks_community/FuRBO/FuRBO.ipynb b/notebooks_community/FuRBO/FuRBO.ipynb new file mode 100644 index 0000000000..f5580dd09c --- /dev/null +++ b/notebooks_community/FuRBO/FuRBO.ipynb @@ -0,0 +1,1011 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b5831947-283e-4682-aae4-bd19bcce03e0", + "metadata": {}, + "source": [ + "# Feasibility-driven trust Region Bayesian Optimization (FuRBO)\n", + "\n", + "- Contributors: paoloascia, elenaraponi\n", + "- Last update 17 October 2025\n", + "- BoTorch version: 0.12.0\n", + "\n", + "This tutorial shows how to implement Feasibility-driven trust Region Bayesian Optimization (FuRBO) with restarts in a closed loop [1].\n", + "\n", + "In this tutorial, we optimize the 10D Ackley function on the domain $[−5,10]^{10}$ subject to two constraint functions $c_1$ and $c_2$. The problem maximizes the Ackley function while the constraints are fulfilled when $c_1(x) \\leq 0$ and $c_2(x) \\leq 0$.\n", + "\n", + "[1] [Ascia, Paolo, Elena Raponi, Thomas Bäck and Fabian Duddeck. \"Feasibility-Driven Trust Region Bayesian Optimization.\" In AutoML 2025 Methods Track.](https://doi.org/10.48550/arXiv.2506.14619)\n", + "\n", + "Since FuRBO is based on Scalable Constrained Bayesian Optimization (SCBO), this tutorial shares part of the same code as the SCBO Tutorial (https://botorch.org/docs/tutorials/scalable_constrained_bo/)\n" + ] + }, + { + "cell_type": "markdown", + "id": "762be478-50e4-4af4-aa6d-d3b566649c5e", + "metadata": {}, + "source": [ + "### Objective function\n", + "\n", + "Start by defining the 10D Ackley function for evaluation during the optimization loop." + ] + }, + { + "cell_type": "code", + "execution_count": 170, + "id": "890f1a54-b6cf-4af4-9bfb-1835b2f737a6", + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "from botorch.test_functions import Ackley\n", + "from botorch.utils.transforms import unnormalize\n", + "\n", + "class ack():\n", + " \n", + " def __init__(self, dim, negate, **tkwargs):\n", + " \n", + " self.fun = Ackley(dim = dim, negate = negate).to(**tkwargs)\n", + " self.fun.bounds[0, :].fill_(-5)\n", + " self.fun.bounds[1, :].fill_(10)\n", + " self.dim = self.fun.dim\n", + " self.lb, self.ub = self.fun.bounds\n", + " \n", + " def eval_(self, x):\n", + " \"\"\"This is a helper function we use to unnormalize and evalaute a point\"\"\"\n", + " return self.fun(unnormalize(x, [self.lb, self.ub]))" + ] + }, + { + "cell_type": "markdown", + "id": "6b710672-51d0-4fc5-a3e2-ffb3da6f5649", + "metadata": {}, + "source": [ + "### Constraint functions\n", + "\n", + "Define two constraint functions." + ] + }, + { + "cell_type": "markdown", + "id": "64122b23-fc01-4e1a-94de-2fba4839fe72", + "metadata": {}, + "source": [ + "\n", + "1. Constraint $c_1$: enforce the $\\sum_{i=1}^{10} x_i \\leq T$. We will specify $T=0$ later. " + ] + }, + { + "cell_type": "code", + "execution_count": 171, + "id": "228d816a-6452-4078-b3fd-6a42569237c3", + "metadata": {}, + "outputs": [], + "source": [ + "class sum_():\n", + " def __init__(self, threshold, lb, ub):\n", + " \n", + " self.lb = lb\n", + " self.ub = ub\n", + " self.threshold = threshold\n", + " return \n", + " \n", + " def c(self, x):\n", + " \"\"\"This is a helper function we use to unnormalize and evaluate a point\"\"\"\n", + " return x.sum() - self.threshold\n", + " \n", + " def eval_(self, x):\n", + " return self.c(unnormalize(x, [self.lb, self.ub]))" + ] + }, + { + "cell_type": "markdown", + "id": "08472c8a-67de-4206-bf0f-871e40edfe7c", + "metadata": {}, + "source": [ + "2. Constraint $c_2$: enforce the $l_2$ norm $\\| \\mathbb{x}\\|_2 \\leq T$. We will specify $T=0.5$ later." + ] + }, + { + "cell_type": "code", + "execution_count": 172, + "id": "1e0f4e8d-657f-49d8-bb59-eb8c2ff4a9f5", + "metadata": {}, + "outputs": [], + "source": [ + "class norm_():\n", + " def __init__(self, threshold, lb, ub):\n", + " \n", + " self.lb = lb\n", + " self.ub = ub\n", + " self.threshold = threshold\n", + " return \n", + " \n", + " def c(self, x):\n", + " return torch.norm(x, p=2) - self.threshold\n", + " \n", + " def eval_(self, x):\n", + " \"\"\"This is a helper function we use to unnormalize and evaluate a point\"\"\"\n", + " return self.c(unnormalize(x, [self.lb, self.ub]))" + ] + }, + { + "cell_type": "markdown", + "id": "e3c6e1cd-15de-4985-ac6c-bf96c6eafc24", + "metadata": {}, + "source": [ + "### Define FuRBO Class\n", + "Define a class to hold the information needed for the optimization loop. \n", + "\n", + "The state is updated with the samples evaluated at each iteration. \n", + "\n", + "Prior to the class, two utility functions are defined. The first one identifies the current best sample, while the second one fits a GPR model to the current dataset. \n", + "\n", + "The ```Furbo_state``` class features a function to reset the status when restarting. Notice that the state is emptied when restarting. Therefore the samples previously evaluated are extracted and saved (see main optimization loop)." + ] + }, + { + "cell_type": "code", + "execution_count": 173, + "id": "020338e2-9eaf-49cf-9e5c-fd00e1a46835", + "metadata": {}, + "outputs": [], + "source": [ + "import gpytorch\n", + "import numpy as np\n", + "\n", + "from botorch.fit import fit_gpytorch_mll\n", + "from botorch.models import SingleTaskGP\n", + "from botorch.models.transforms.outcome import Standardize\n", + "\n", + "from gpytorch.constraints import Interval\n", + "from gpytorch.kernels import MaternKernel, ScaleKernel\n", + "from gpytorch.likelihoods import GaussianLikelihood\n", + "from gpytorch.mlls import ExactMarginalLogLikelihood\n", + "\n", + "from scipy.stats import invgauss\n", + "from scipy.stats import ecdf\n", + "\n", + "from torch import Tensor\n", + "\n", + "def get_best_index_for_batch(n_tr, Y: Tensor, C: Tensor):\n", + " \"\"\"Return the index for the best point. One for each trust region.\n", + " For reference, see https://botorch.org/docs/tutorials/scalable_constrained_bo/\"\"\"\n", + " is_feas = (C <= 0).all(dim=-1)\n", + " if is_feas.any(): # Choose best feasible candidate\n", + " score = Y.clone()\n", + " score[~is_feas] = -float(\"inf\")\n", + " return torch.topk(score.reshape(-1), k=n_tr).indices\n", + " return torch.topk(C.clamp(min=0).sum(dim=-1), k=n_tr, largest=False).indices # Return smallest violation\n", + "\n", + "def get_fitted_model(X,\n", + " Y,\n", + " dim,\n", + " max_cholesky_size):\n", + " '''Function to fit a GPR to a given set of data.\n", + " For reference, see https://botorch.org/docs/tutorials/scalable_constrained_bo/'''\n", + " likelihood = GaussianLikelihood(noise_constraint=Interval(1e-8, 1e-3))\n", + " covar_module = ScaleKernel( # Use the same lengthscale prior as in the TuRBO paper\n", + " MaternKernel(nu=2.5, ard_num_dims=dim, lengthscale_constraint=Interval(0.005, 4.0))\n", + " )\n", + " model = SingleTaskGP(\n", + " X,\n", + " Y,\n", + " covar_module=covar_module,\n", + " likelihood=likelihood,\n", + " outcome_transform=Standardize(m=1),\n", + " )\n", + " mll = ExactMarginalLogLikelihood(model.likelihood, model)\n", + "\n", + " with gpytorch.settings.max_cholesky_size(max_cholesky_size):\n", + " fit_gpytorch_mll(mll, \n", + " optimizer_kwargs={'method': 'L-BFGS-B'})\n", + "\n", + " return model\n", + "\n", + "from botorch.models.model_list_gp_regression import ModelListGP\n", + "from torch.quasirandom import SobolEngine\n", + "\n", + "class Furbo_state():\n", + " '''Class to track optimization status with restart'''\n", + " # Initialization of the status\n", + " def __init__(self, \n", + " obj, # Objective function\n", + " cons, # Constraints function\n", + " batch_size, # Batch size of each iteration\n", + " n_init, # Number of initial points to evaluate\n", + " n_iteration, # Number of total iterations\n", + " **tkwargs):\n", + " \n", + " # Objective function handle\n", + " self.obj = obj\n", + " \n", + " # Constraints function handle\n", + " self.cons = cons\n", + " \n", + " # Domain bounds\n", + " self.lb = obj.lb\n", + " self.ub = obj.ub\n", + " \n", + " # Problem dimensions\n", + " self.batch_size: int = batch_size # Dimension of the batch at each iteration\n", + " self.n_init: int = n_init # Number of initial samples\n", + " self.dim: int = obj.dim # Dimension of the problem\n", + " \n", + " # Trust regions information\n", + " self.tr_ub: float = torch.ones((1, self.dim), **tkwargs) # Upper bounds of trust region\n", + " self.tr_lb: float = torch.zeros((1, self.dim), **tkwargs) # Lower bounds of trust region\n", + " self.tr_vol: float = torch.prod(self.tr_ub - self.tr_lb, dim=1) # Volume of trust region\n", + " self.radius: float = 1.0 # Percentage around which the trust region is built\n", + " self.radius_min: float = 0.5**7 # Minimum percentage for trust region\n", + "\n", + " # Trust region updating \n", + " self.failure_counter: int = 0 # Counter of failure points to asses how algorithm is going\n", + " self.success_counter: int = 0 # Counter of success points to asses how algorithm is going\n", + " self.success_tolerance: int = 2 # Success tolerance for \n", + " self.failure_tolerance: int = 3 # Failure tolerance for\n", + " \n", + " # Tensor to save current batch information\n", + " self.batch_X: Tensor # Current batch to evaluate: X values\n", + " self.batch_Y: Tensor # Current batch to evaluate: Y value\n", + " self.batch_C: Tensor # Current batch to evaluate: C values\n", + " \n", + " # Stopping criteria information\n", + " self.n_iteration: int = n_iteration # Maximum number of iterations allowed\n", + " self.it_counter: int = 0 # Counter of iterations for stopping\n", + " self.finish_trigger: bool = False # Trigger to stop optimization\n", + " self.failed_GP : bool = False # Flag to pass to failed_GP in FuRBORestart\n", + " \n", + " # Restart criteria information\n", + " self.restart_trigger: bool = False\n", + " \n", + " # Sobol sampler engine\n", + " self.sobol = SobolEngine(dimension=self.dim, scramble=True, seed=1)\n", + " \n", + " # Update the status\n", + " def update(self,\n", + " X_next, # Samples X (input values) to update the status\n", + " Y_next, # Samples Y (objective value) to update the status\n", + " C_next, # Samples C (constraints values) to update the status\n", + " **tkwargs):\n", + " \n", + " '''Function to update optimization status'''\n", + " \n", + " # Merge current batch with previously evaluated samples\n", + " if not hasattr(self, 'X'):\n", + " # If there are no previous samples, declare the Tensors\n", + " self.X = X_next\n", + " self.Y = Y_next\n", + " self.C = C_next\n", + " else:\n", + " # Else, concatenate the new batch to the previous samples\n", + " self.X = torch.cat((self.X, X_next), dim=0)\n", + " self.Y = torch.cat((self.Y, Y_next), dim=0)\n", + " self.C = torch.cat((self.C, C_next), dim=0)\n", + "\n", + " # update GPR surrogates\n", + " try:\n", + " self.Y_model = get_fitted_model(self.X, self.Y, self.dim, max_cholesky_size = float(\"inf\"))\n", + " self.C_model = ModelListGP(*[get_fitted_model(self.X, C.reshape([C.shape[0],1]), self.dim, max_cholesky_size = float(\"inf\")) for C in self.C.t()])\n", + " except:\n", + " # If update fail, flag to stop entire optimization\n", + " self.failed_GP = True\n", + " \n", + " # Update batch information \n", + " self.batch_X = X_next\n", + " self.batch_Y = Y_next\n", + " self.batch_C = C_next\n", + " \n", + " # Update best value\n", + " # Find the best value among the candidates\n", + " best_id = get_best_index_for_batch(n_tr=1, Y=self.Y, C=self.C)\n", + " \n", + " # Update success and failure counters for trust region update\n", + " # If attribute 'best_X' does not exist, DoE was just evaluated -> no update on counters\n", + " if hasattr(self, 'best_X'):\n", + " if (self.C[best_id] <= 0).all():\n", + " # At least one new candidate is feasible\n", + " if (self.Y[best_id] > self.best_Y).any() or (self.best_C > 0).any():\n", + " self.success_counter += 1\n", + " self.failure_counter = 0 \n", + " else:\n", + " self.success_counter = 0\n", + " self.failure_counter += 1\n", + " else:\n", + " # No new candidate is feasible\n", + " total_violation_next = self.C[best_id].clamp(min=0).sum(dim=-1)\n", + " total_violation_center = self.best_C.clamp(min=0).sum(dim=-1)\n", + " if total_violation_next < total_violation_center:\n", + " self.success_counter += 1\n", + " self.failure_counter = 0\n", + " else:\n", + " self.success_counter = 0\n", + " self.failure_counter += 1\n", + " \n", + " # Update best values\n", + " self.best_X = self.X[best_id]\n", + " self.best_Y = self.Y[best_id]\n", + " self.best_C = self.C[best_id]\n", + " \n", + " # Update iteration counter\n", + " self.it_counter += 1\n", + " \n", + " def reset_status(self,\n", + " **tkwargs):\n", + " '''Function to reset the status for the restart'''\n", + " \n", + " # Reset trust regions size\n", + " self.tr_ub: float = torch.ones((1, self.dim), **tkwargs) # Upper bounds of trust region\n", + " self.tr_lb: float = torch.zeros((1, self.dim), **tkwargs) # Lower bounds of trust region\n", + " self.tr_vol: float = torch.prod(self.tr_ub - self.tr_lb, dim=1) # Volume of trust region\n", + " self.radius: float = 1.0 # Percentage around which the trust region is built\n", + " self.radius_min: float = 0.5**7 # Minimum percentage for trust region\n", + "\n", + " # Reset counters to change trust region size \n", + " self.failure_counter: int = 0 # Counter of failure points to asses how algorithm is going\n", + " self.success_counter: int = 0 # Counter of success points to asses how algorithm is going\n", + " \n", + " # Reset restart criteria trigger\n", + " self.restart_trigger: bool = False # Trigger to restart optimization\n", + " self.failed_GP: bool = False # Reset GPR failure trigger\n", + " \n", + " # Delete tensors with samples for training GPRs\n", + " if hasattr(self, 'X'):\n", + " del self.X\n", + " del self.Y\n", + " del self.C\n", + " \n", + " # Delete tensors with best value so far\n", + " if hasattr(self, 'best_X'):\n", + " del self.best_X\n", + " del self.best_Y\n", + " del self.best_C\n", + " \n", + " # Clear GPU memory\n", + " if tkwargs[\"device\"] == \"cuda\":\n", + " torch.cuda.empty_cache() " + ] + }, + { + "cell_type": "markdown", + "id": "58fee9fb-ca01-4024-8859-1fc3d9d4aea4", + "metadata": {}, + "source": [ + "### Define trust region\n", + "\n", + "Define a set of functions to evaluate the trust region. First sample according to a Multinormal distribution the GPR surrogates (both objective and constraints). Rank the samples according to both the objective and violation estimation. Take the top $10\\%$ of the samples according to the rank. The trust region is defined as a hyperbox enclosing the picked samples. " + ] + }, + { + "cell_type": "code", + "execution_count": 174, + "id": "3d36cf01-c5be-44ee-96bb-12564225bd7b", + "metadata": {}, + "outputs": [], + "source": [ + "def multivariate_circular(centre, # Centre of the multivariate distribution\n", + " radius, # Radius of the multivariate distribution\n", + " n_samples, # Number of samples to evaluate\n", + " lb = None, # Domain lower bound\n", + " ub = None, # Domain upper bound\n", + " **tkwargs):\n", + " '''Function to generate multivariate distribution of given radius and centre within a given domain.'''\n", + " # Dimension of the design domain\n", + " dim = centre.shape[0]\n", + " \n", + " # Generate a multivariate normal distribution centered at 0\n", + " multivariate_normal = torch.distributions.multivariate_normal.MultivariateNormal(torch.zeros(dim, **tkwargs), 0.025*torch.eye(dim, **tkwargs))\n", + " \n", + " # Draw samples torch.distributions.multivariate_normal import MultivariateNormal\n", + " samples = multivariate_normal.sample(sample_shape=torch.Size([n_samples]))\n", + " \n", + " # Normalize each sample to have unit norm, then scale by the radius\n", + " norms = torch.norm(samples, dim=1, keepdim=True) # Euclidean norms\n", + " normalized_samples = samples / norms # Normalize to unit hypersphere\n", + " scaled_samples = normalized_samples * torch.rand(n_samples, 1, **tkwargs) * radius # Scale by random factor within radius\n", + " \n", + " # Translate samples to be centered at centre\n", + " samples = scaled_samples + centre\n", + " \n", + " \n", + " # Trim samples outside domain\n", + " for dim in range(len(lb)):\n", + " samples = samples[torch.where(samples[:,dim]>=lb[dim])]\n", + " samples = samples[torch.where(samples[:,dim]<=ub[dim])]\n", + " \n", + " return samples\n", + "\n", + "def update_tr(state, # FuRBO state\n", + " percentage = 0.1, # Percentage to define trust region (default 10%)\n", + " **tkwargs):\n", + " '''Function to sample Multinormal Distribution of GPRs and define trust region'''\n", + " # Update the trust regions based on the feasible region\n", + " n_samples = 1000 * state.dim\n", + " lb = torch.zeros(state.dim, **tkwargs)\n", + " ub = torch.ones(state.dim, **tkwargs)\n", + " \n", + " # Update radius dimension\n", + " if state.success_counter == state.success_tolerance: # Expand trust region\n", + " state.radius = min(2.0 * state.radius, 1.0)\n", + " state.success_counter = 0\n", + " elif state.failure_counter == state.failure_tolerance: # Shrink trust region\n", + " state.radius /= 2.0\n", + " state.failure_counter = 0\n", + " \n", + " for ind, x_candidate in enumerate(state.best_X):\n", + " # Generate the samples to evaluathe the feasible area on\n", + " radius = state.radius\n", + " samples = multivariate_circular(x_candidate, radius, n_samples, lb=lb, ub=ub, **tkwargs)\n", + " \n", + " # Evaluate samples on the models of the objective -> yy Tensor\n", + " state.Y_model.eval()\n", + " with torch.no_grad():\n", + " posterior = state.Y_model.posterior(samples)\n", + " samples_yy = posterior.mean.squeeze()\n", + " \n", + " # Evaluate samples on the models of the constraints -> yy Tensor\n", + " state.C_model.eval()\n", + " with torch.no_grad():\n", + " posterior = state.C_model.posterior(samples)\n", + " samples_cc = posterior.mean\n", + " \n", + " # Combine the constraints values\n", + " # Normalize\n", + " samples_cc /= torch.abs(samples_cc).max(dim=0).values\n", + " samples_cc = torch.max(samples_cc, dim=1).values\n", + " \n", + " # Take the best X% of the drawn samples to define the trust region\n", + " n_samples_tr = int(n_samples * percentage)\n", + " \n", + " # Order the samples for feasibility and for best objective\n", + " if torch.any(samples_cc < 0):\n", + " \n", + " feasible_samples_id = torch.where(samples_cc <= 0)[0]\n", + " infeasible_samples_id = torch.where(samples_cc > 0)[0]\n", + " \n", + " feasible_cc = -1 * samples_yy[feasible_samples_id]\n", + " infeasible_cc = samples_cc[infeasible_samples_id]\n", + " \n", + " feasible_sorted, feasible_sorted_id = torch.sort(feasible_cc)\n", + " infeasible_sorted, infeasible_sorted_id = torch.sort(infeasible_cc)\n", + " \n", + " original_feasible_sorted_indices = feasible_samples_id[feasible_sorted_id]\n", + " original_infeasible_sorted_indices = infeasible_samples_id[infeasible_sorted_id]\n", + " \n", + " top_indices = torch.cat((original_feasible_sorted_indices, original_infeasible_sorted_indices))[:n_samples_tr]\n", + " \n", + " else:\n", + " \n", + " if n_samples_tr > len(samples_cc):\n", + " n_samples_tr = len(samples_cc)\n", + " \n", + " if n_samples_tr < 4:\n", + " n_samples_tr = 4\n", + " \n", + " top_values, top_indices = torch.topk(samples_cc, n_samples_tr, largest=False)\n", + " \n", + " # Set the box around the selected samples\n", + " state.tr_lb[ind] = torch.min(samples[top_indices], dim=0).values\n", + " state.tr_ub[ind] = torch.max(samples[top_indices], dim=0).values\n", + " \n", + " # Update volume of trust region\n", + " state.tr_vol[ind] = torch.prod(state.tr_ub[ind] - state.tr_lb[ind])\n", + " \n", + " # return updated status with new trust regions\n", + " return state" + ] + }, + { + "cell_type": "markdown", + "id": "f49690e5-6505-47df-89f7-1a56f9b087b2", + "metadata": {}, + "source": [ + "### Sampling strategies\n", + "\n", + "Define a function to generate an initial experimental design using Sobol sampling strategy, similarly to SCBO (https://botorch.org/docs/tutorials/scalable_constrained_bo/)." + ] + }, + { + "cell_type": "code", + "execution_count": 175, + "id": "0116ca79-7555-4da3-bfd4-69941926eb11", + "metadata": {}, + "outputs": [], + "source": [ + "def get_initial_points(state,\n", + " **tkwargs):\n", + " '''Function to generate the initial experimental design'''\n", + " X_init = state.sobol.draw(n=state.n_init).to(**tkwargs)\n", + " return X_init" + ] + }, + { + "cell_type": "markdown", + "id": "55726979-92f6-488f-869c-2fa6140f6b85", + "metadata": {}, + "source": [ + "Define a function to identify the best next candidate point, similar to SCBO(https://botorch.org/docs/tutorials/scalable_constrained_bo/)." + ] + }, + { + "cell_type": "code", + "execution_count": 176, + "id": "2d969ea7-2f1e-4433-b3b4-53413546c2f2", + "metadata": {}, + "outputs": [], + "source": [ + "from botorch.generation.sampling import ConstrainedMaxPosteriorSampling\n", + "\n", + "def generate_batch(state,\n", + " n_candidates,\n", + " **tkwargs):\n", + " '''Function to find net candidate optimum'''\n", + " assert state.X.min() >= 0.0 and state.X.max() <= 1.0 and torch.all(torch.isfinite(state.Y))\n", + "\n", + " # Initialize tensor with samples to evaluate\n", + " X_next = torch.ones((state.batch_size, state.dim), **tkwargs)\n", + " \n", + " # Iterate over the several trust regions\n", + "\n", + " tr_lb = state.tr_lb[0]\n", + " tr_ub = state.tr_ub[0]\n", + "\n", + " # Thompson Sampling w/ Constraints (like SCBO)\n", + " pert = state.sobol.draw(n_candidates).to(**tkwargs)\n", + " pert = tr_lb + (tr_ub - tr_lb) * pert\n", + "\n", + " # Create a perturbation mask\n", + " prob_perturb = min(20.0 / state.dim, 1.0)\n", + " mask = torch.rand(n_candidates, state.dim, **tkwargs) <= prob_perturb\n", + " ind = torch.where(mask.sum(dim=1) == 0)[0]\n", + " mask[ind, torch.randint(0, state.dim - 1, size=(len(ind),), device=tkwargs['device'])] = 1\n", + "\n", + " # Create candidate points from the perturbations and the mask\n", + " X_cand = state.best_X[0].expand(n_candidates, state.dim).clone()\n", + " X_cand[mask] = pert[mask]\n", + " \n", + " # Sample on the candidate points using Constrained Max Posterior Sampling\n", + " constrained_thompson_sampling = ConstrainedMaxPosteriorSampling(\n", + " model=state.Y_model, constraint_model=state.C_model, replacement=False\n", + " )\n", + " with torch.no_grad():\n", + " X_next[0*state.batch_size:0*state.batch_size+state.batch_size, :] = constrained_thompson_sampling(X_cand, num_samples=state.batch_size)\n", + " \n", + " return X_next" + ] + }, + { + "cell_type": "markdown", + "id": "024e10f4-4de9-433b-8827-c58a7177784f", + "metadata": {}, + "source": [ + "### Stopping criterion\n", + "\n", + "Define a function to detect when the maximum number of iterations is met." + ] + }, + { + "cell_type": "code", + "execution_count": 177, + "id": "189aae72-f033-49db-bcc8-0e791774bd76", + "metadata": {}, + "outputs": [], + "source": [ + "def stopping_criterion(state, n_iteration):\n", + " '''Function to evaluate if the maximum number of allowed iterations is reached.'''\n", + " if state.it_counter <= n_iteration:\n", + " return False\n", + " return True" + ] + }, + { + "cell_type": "markdown", + "id": "44f15711-8503-4226-8ee0-171f26f7b8f4", + "metadata": {}, + "source": [ + "### Restart criterion\n", + "\n", + "Detect when the GPR fitting process fails to stop the optimization.curve" + ] + }, + { + "cell_type": "code", + "execution_count": 178, + "id": "f0d2b342-601b-4d07-9fe3-459f28ecead2", + "metadata": {}, + "outputs": [], + "source": [ + "def GP_restart_criterion(state):\n", + " '''Function to evaluate if a GPR failed during the optimization.'''\n", + " if state.failed_GP:\n", + " print(\"GPR failed.\")\n", + " return True\n", + " return False" + ] + }, + { + "cell_type": "markdown", + "id": "476bada4", + "metadata": {}, + "source": [ + "Detect when the radius becomes too small." + ] + }, + { + "cell_type": "code", + "execution_count": 179, + "id": "81f9bd26", + "metadata": {}, + "outputs": [], + "source": [ + "def restart_criterion(state, radius_min):\n", + " '''Function to evaluate if MND radius is smaller than the minimum allowed radius'''\n", + " if state.radius < radius_min:\n", + " return True\n", + " return False" + ] + }, + { + "cell_type": "markdown", + "id": "a80c7b75-a62d-46c2-aa80-28f6f29501be", + "metadata": {}, + "source": [ + "### Main optimization loop" + ] + }, + { + "cell_type": "code", + "execution_count": 180, + "id": "b9d52417-aaa5-40b6-bf07-12c074460283", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0) No feasible point yet! Smallest total violation: 5.72e-01, MND radius: 1.0\n", + "1) Best value: -1.72e+00, MND radius: 1.0\n", + "2) Best value: -1.72e+00, MND radius: 1.0\n", + "3) Best value: -7.71e-01, MND radius: 1.0\n", + "4) Best value: -7.71e-01, MND radius: 1.0\n", + "5) Best value: -7.71e-01, MND radius: 1.0\n", + "6) Best value: -2.39e-01, MND radius: 1.0\n", + "7) Best value: -2.37e-01, MND radius: 1.0\n", + "8) Best value: -2.37e-01, MND radius: 1.0\n", + "9) Best value: -8.69e-02, MND radius: 1.0\n", + "10) Best value: -8.32e-02, MND radius: 1.0\n", + "11) Best value: -8.32e-02, MND radius: 1.0\n", + "12) Best value: -8.32e-02, MND radius: 1.0\n", + "13) Best value: -8.32e-02, MND radius: 1.0\n", + "14) Best value: -8.32e-02, MND radius: 0.5\n", + "15) Best value: -8.32e-02, MND radius: 0.5\n", + "16) Best value: -2.84e-02, MND radius: 0.5\n", + "17) Best value: -2.84e-02, MND radius: 0.5\n", + "18) Best value: -2.84e-02, MND radius: 0.5\n", + "19) Best value: -2.84e-02, MND radius: 0.5\n", + "20) Best value: -2.84e-02, MND radius: 0.25\n", + "21) Best value: -2.84e-02, MND radius: 0.25\n", + "22) Best value: -2.84e-02, MND radius: 0.25\n", + "23) Best value: -2.21e-02, MND radius: 0.125\n", + "24) Best value: -1.57e-02, MND radius: 0.125\n", + "25) Best value: -1.57e-02, MND radius: 0.25\n", + "26) Best value: -1.57e-02, MND radius: 0.25\n", + "27) Best value: -1.57e-02, MND radius: 0.25\n", + "28) Best value: -1.57e-02, MND radius: 0.125\n", + "29) Best value: -1.57e-02, MND radius: 0.125\n", + "30) Best value: -1.56e-02, MND radius: 0.125\n", + "31) Best value: -1.56e-02, MND radius: 0.125\n", + "32) Best value: -1.56e-02, MND radius: 0.125\n", + "33) Best value: -1.56e-02, MND radius: 0.125\n", + "34) Best value: -1.46e-02, MND radius: 0.0625\n", + "35) Best value: -1.46e-02, MND radius: 0.0625\n", + "36) Best value: -1.39e-02, MND radius: 0.0625\n", + "37) Best value: -1.16e-02, MND radius: 0.0625\n", + "38) Best value: -1.16e-02, MND radius: 0.125\n", + "39) Best value: -1.16e-02, MND radius: 0.125\n", + "40) Best value: -1.16e-02, MND radius: 0.125\n", + "41) Best value: -1.16e-02, MND radius: 0.0625\n", + "42) Best value: -1.16e-02, MND radius: 0.0625\n", + "43) Best value: -7.18e-03, MND radius: 0.0625\n", + "44) Best value: -7.18e-03, MND radius: 0.0625\n", + "45) Best value: -7.18e-03, MND radius: 0.0625\n", + "46) Best value: -7.18e-03, MND radius: 0.0625\n", + "47) Best value: -3.63e-03, MND radius: 0.03125\n", + "48) Best value: -3.63e-03, MND radius: 0.03125\n", + "GPR failed.\n", + "49) No feasible point yet! Smallest total violation: 4.53e-01, MND radius: 1.0\n", + "50) Best value: -2.94e+00, MND radius: 1.0\n", + "51) Best value: -2.94e+00, MND radius: 1.0\n", + "52) Best value: -2.94e+00, MND radius: 1.0\n", + "53) Best value: -2.94e+00, MND radius: 1.0\n", + "54) Best value: -1.44e+00, MND radius: 0.5\n", + "55) Best value: -2.58e-01, MND radius: 0.5\n", + "56) Best value: -2.58e-01, MND radius: 1.0\n", + "57) Best value: -7.49e-02, MND radius: 1.0\n", + "58) Best value: -7.49e-02, MND radius: 1.0\n", + "59) Best value: -6.21e-02, MND radius: 1.0\n", + "60) Best value: -6.21e-02, MND radius: 1.0\n", + "61) Best value: -6.21e-02, MND radius: 1.0\n", + "62) Best value: -6.21e-02, MND radius: 1.0\n", + "63) Best value: -6.21e-02, MND radius: 0.5\n", + "64) Best value: -6.21e-02, MND radius: 0.5\n", + "65) Best value: -2.03e-02, MND radius: 0.5\n", + "66) Best value: -2.03e-02, MND radius: 0.5\n", + "67) Best value: -2.03e-02, MND radius: 0.5\n", + "68) Best value: -2.03e-02, MND radius: 0.5\n", + "69) Best value: -2.03e-02, MND radius: 0.25\n", + "70) Best value: -2.03e-02, MND radius: 0.25\n", + "71) Best value: -2.03e-02, MND radius: 0.25\n", + "72) Best value: -1.33e-02, MND radius: 0.125\n", + "73) Best value: -1.33e-02, MND radius: 0.125\n", + "74) Best value: -1.33e-02, MND radius: 0.125\n", + "75) Best value: -1.33e-02, MND radius: 0.125\n", + "76) Best value: -3.34e-03, MND radius: 0.125\n", + "77) Best value: -3.34e-03, MND radius: 0.125\n", + "78) Best value: -3.34e-03, MND radius: 0.125\n", + "79) Best value: -3.34e-03, MND radius: 0.125\n", + "80) Best value: -1.99e-03, MND radius: 0.0625\n", + "81) Best value: -1.99e-03, MND radius: 0.0625\n", + "82) Best value: -1.99e-03, MND radius: 0.0625\n", + "83) Best value: -1.99e-03, MND radius: 0.0625\n", + "84) Best value: -1.99e-03, MND radius: 0.03125\n", + "85) Best value: -1.99e-03, MND radius: 0.03125\n", + "86) Best value: -1.99e-03, MND radius: 0.03125\n", + "87) Best value: -1.99e-03, MND radius: 0.015625\n", + "88) Best value: -1.99e-03, MND radius: 0.015625\n", + "89) Best value: -1.99e-03, MND radius: 0.015625\n", + "90) Best value: -1.99e-03, MND radius: 0.0078125\n", + "91) Best value: -1.99e-03, MND radius: 0.0078125\n", + "92) Best value: -1.99e-03, MND radius: 0.0078125\n", + "93) No feasible point yet! Smallest total violation: 2.07e+00, MND radius: 1.0\n", + "94) Best value: -3.47e+00, MND radius: 1.0\n", + "95) Best value: -3.47e+00, MND radius: 1.0\n", + "96) Best value: -2.16e+00, MND radius: 1.0\n", + "97) Best value: -2.16e+00, MND radius: 1.0\n", + "98) Best value: -2.16e+00, MND radius: 1.0\n", + "99) Best value: -2.16e+00, MND radius: 1.0\n", + "100) Best value: -4.88e-01, MND radius: 0.5\n" + ] + } + ], + "source": [ + "import warnings\n", + "\n", + "warnings.filterwarnings(\"ignore\")\n", + " \n", + "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n", + "dtype = torch.double\n", + "tkwargs = {\"device\": device, \"dtype\": dtype}\n", + "\n", + "# Initialize FuRBO\n", + "obj = ack(dim = 2,\n", + " negate=True,\n", + " **tkwargs)\n", + "cons = list([sum_(threshold = 0,\n", + " lb = obj.lb,\n", + " ub = obj.ub), \n", + " norm_(threshold = 0.5, \n", + " lb = obj.lb, \n", + " ub = obj.ub)])\n", + "batch_size = int(1)#3 * obj.dim)\n", + "n_init = int(10)# * obj.dim)\n", + "n_iteration = int(100)# * obj.dim)\n", + "N_CANDIDATES = 2000\n", + "\n", + "# FuRBO state initialization\n", + "FuRBO_status = Furbo_state(obj = obj, # Objective function\n", + " cons = cons, # Constraints function\n", + " batch_size = batch_size, # Batch size of each iteration\n", + " n_init = n_init, # Number of initial points to evaluate\n", + " n_iteration = n_iteration, # Number of iterations\n", + " **tkwargs)\n", + "\n", + "# Initiate lists to save samples over the restarts\n", + "X_best, Y_best, C_best = [], [], []\n", + "X_all, Y_all, C_all = [], [], []\n", + "\n", + "# Continue optimization the stopping criterions isn't triggered\n", + "while not FuRBO_status.finish_trigger: \n", + " \n", + " # Reset status for restarting\n", + " FuRBO_status.reset_status(**tkwargs)\n", + " \n", + " # generate intial batch of X\n", + " X_next = get_initial_points(FuRBO_status, **tkwargs)\n", + " \n", + " # Reset and restart optimization\n", + " while not FuRBO_status.restart_trigger and not FuRBO_status.finish_trigger:\n", + " \n", + " # Evaluate current batch (samples in X_next)\n", + " Y_next = []\n", + " C_next = []\n", + " for x in X_next:\n", + " # Evaluate batch on obj ...\n", + " Y_next.append(FuRBO_status.obj.eval_(x))\n", + " # ... and constraints\n", + " C_next.append([c.eval_(x) for c in FuRBO_status.cons])\n", + " \n", + " # process vector for PyTorch\n", + " Y_next = torch.tensor(Y_next).unsqueeze(-1).to(**tkwargs)\n", + " C_next = torch.tensor(C_next).to(**tkwargs)\n", + " \n", + " # Update FuRBO status with newly evaluated batch\n", + " FuRBO_status.update(X_next, Y_next, C_next, **tkwargs) \n", + " \n", + " # Printing current best\n", + " # If a feasible has been evaluated -> print current optimum (feasible sample with best objective value)\n", + " if (FuRBO_status.best_C <= 0).all():\n", + " best = FuRBO_status.best_Y.amax()\n", + " print(f\"{FuRBO_status.it_counter-1}) Best value: {best:.2e},\"\n", + " f\" MND radius: {FuRBO_status.radius}\")\n", + " \n", + " # Else, if no feasible has been evaluated -> print smallest violation (the sample that violatest the least all constraints)\n", + " else:\n", + " violation = FuRBO_status.best_C.clamp(min=0).sum()\n", + " print(f\"{FuRBO_status.it_counter-1}) No feasible point yet! Smallest total violation: \"\n", + " f\"{violation:.2e}, MND radius: {FuRBO_status.radius}\")\n", + " \n", + " # Update Trust regions\n", + " FuRBO_status = update_tr(FuRBO_status,\n", + " **tkwargs)\n", + " \n", + " # generate next batch to evaluate \n", + " X_next = generate_batch(FuRBO_status, N_CANDIDATES, **tkwargs)\n", + " \n", + " # Check if stopping criterion is met (budget exhausted and if GP failed)\n", + " FuRBO_status.finish_trigger = stopping_criterion(FuRBO_status, n_iteration) \n", + " \n", + " # Check if restart criterion is met\n", + " FuRBO_status.restart_trigger = (restart_criterion(FuRBO_status, FuRBO_status.radius_min) or GP_restart_criterion(FuRBO_status))\n", + "\n", + " # Save samples evaluated before resetting the status\n", + " X_all.append(FuRBO_status.X)\n", + " Y_all.append(FuRBO_status.Y)\n", + " C_all.append(FuRBO_status.C)\n", + "\n", + " # Save best sample of this run\n", + " X_best.append(FuRBO_status.best_X)\n", + " Y_best.append(FuRBO_status.best_Y)\n", + " C_best.append(FuRBO_status.best_C)" + ] + }, + { + "cell_type": "markdown", + "id": "a3d28f6c-5940-4d56-a3b5-4635980db76f", + "metadata": {}, + "source": [ + "### Printing result and plotting convergence curve\n", + "\n", + "Print the best-evaluated sample (over all restarts) and the objective value (if a feasible sample was found) or the smallest violation (if no feasible sample was found)." + ] + }, + { + "cell_type": "code", + "execution_count": 181, + "id": "027a9ec9-930a-48d0-b481-ffe9e9ca0d90", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Optimization finished \n", + "\t Optimum: -1.99e-03, \n", + "\t X: tensor([[0.3333, 0.3333]], dtype=torch.float64)\n" + ] + } + ], + "source": [ + "# Print best value found so far\n", + "# Ri-elaborate for processing\n", + "X_best = torch.stack(X_best).to(**tkwargs)\n", + "Y_best = torch.stack(Y_best).to(**tkwargs)\n", + "C_best = torch.stack(C_best).to(**tkwargs)\n", + "\n", + "# If a feasible has been evaluated -> print current optimum sample and yielded value\n", + "if (C_best <= 0).any():\n", + " best = Y_best.amax()\n", + " bext = X_best[Y_best.argmax()]\n", + " print(\"Optimization finished \\n\"\n", + " f\"\\t Optimum: {best:.2e}, \\n\"\n", + " f\"\\t X: {bext}\")\n", + " \n", + "# Else, if no feasible has been evaluated -> print sample with smallest violation and the violation value\n", + "else:\n", + " violation = C_best.sum(dim=2).amin()\n", + " violaxion = X_best[C_best.sum(dim=2).argmin()]\n", + " \n", + " print(\"Optimization failed \\n\"\n", + " f\"\\t Smallest violation: {violation:.2e}, \\n\"\n", + " f\"\\t X: {violaxion}\")" + ] + }, + { + "cell_type": "markdown", + "id": "db5ae248-3f04-4173-a98a-6c0738e38510", + "metadata": {}, + "source": [ + "Plot the monotonic convergence curve" + ] + }, + { + "cell_type": "code", + "execution_count": 189, + "id": "671ec5e5", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAI1CAYAAADLpwyxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+6UlEQVR4nO3dd3QU1d8G8Gc2ZdMrkJCEEgQNvXeEUAQEERBpIk2KAipFpIkUFRELKoiCSlNBKQoKKhhCQMAAUkLv0iGElkLqJnvfP/Lu/GZrNpuySfb5nJMDmblz5+7c2cl3794iCSEEiIiIiIgIAKCydwGIiIiIiEoSBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyER2cOXKFUiSBEmSsGrVKpvy0B0/Z86cQi1bWbFr1y75Gu3atatAeUVGRkKSJERGRhZK2YiIrDFnzhz5OUbFiwGyg0lISMDWrVsxa9YsPP300yhXrpz85hs2bFi+8/vzzz/Ru3dvhIWFQa1WIywsDL1798aff/6Z57FVq1aVz6388fLyQmhoKOrXr48hQ4Zg8eLFuHnzpg2v1nrDhw+Xz//EE08U6bmILLly5QoWL16MPn36oEaNGvDw8ICbmxvCwsLQq1cv/PTTT8jOzrZ3MYmIyjQGyA4mKCgIPXr0wLvvvott27bh/v37NuWj1WoxcuRIdOvWDZs3b8bNmzeRlZWFmzdvYvPmzejWrRtGjRoFrVab77xTU1Nx69YtHD9+HN9//z1ef/11VKlSBc899xyuXLliU3ktSUtLw8aNG+Xfz58/j/379xf6eahkKYktM2+//TaqVauG119/Hb/88gsuXryI9PR0ZGZm4ubNm/j1118xcOBAtGrVCteuXbN3cQtFYXybQlQQusYaWxqJygK+B01jgOzAKleujM6dO9t07FtvvYXly5cDABo2bIgff/wRBw8exI8//oiGDRsCAL799lvMnDkzz7xCQkJw4sQJ+efIkSOIjo7GihUrMHLkSPj6+iInJwebNm1C/fr1sWXLFpvKbM4vv/yCR48eAQA8PT0BAN99912hnoOKX2RkJIQQEEKUmq4Rt2/fhhACnp6eePHFF7Fy5Urs3bsXhw4dwvfff4+mTZsCAP7991906tRJvm+JqGyaM2eO/Byj4sUA2cHMmjULW7ZsQXx8PK5evYply5blO4/z58/j448/BgA0adIE+/btw4ABA9C0aVMMGDAAe/fuRZMmTQAAH330ES5evGgxPxcXF9SpU0f+adiwITp06IDhw4fjm2++wY0bNzBp0iQAQHJyMvr3749Dhw7lu9zm6ILhxo0bY+jQoQCAdevWISsrq9DOQWSNwMBALFiwALdv38b333+PYcOGoXXr1mjcuDFefPFFxMbGol+/fgCACxcuYOHChXYuMRFR2cQA2cHMnTsXzzzzDIKCgmzO47PPPpP7QC5evBju7u56+z08PLB48WIAQHZ2Nj799FPbCwzAy8sLn3zyCT744AMAQHp6OkaOHFmgPHVu3bqF6OhoAMCgQYPw4osvAgAePHiArVu3Fso5iKy1YMECTJkyBd7e3ib3Ozk54csvv4SrqysA6HUNIiKiwsMAmfJFCIFff/0VABAREYEWLVqYTNeiRQt5sNuvv/5aKF8PTZkyBc2bNwcAHDt2DH/88UeB8/zhhx+g1Wrh5OSEgQMHomXLlqhevTqA/HWzOHnyJF577TXUrVsX/v7+cHFxQXBwMDp16oQPP/wQt2/fznfZMjIy0LNnT7lv2IIFC/KdBwAcOXIEr7zyCp544gl4eXnB09MTTzzxBMaMGYPz58+bPKZRo0aQJAk1a9bMM//79+9DrVZDkiSMHTvW6nJt3LhRfm1nz541mUY5kHPz5s0m03Tt2hWSJBndi+ZmsVi1ahUkScLcuXPlbaYGi1rq737z5k1MmjQJ1atXh7u7OwIDA9GlSxerBqcWVGBgIOrVqwcAuHTpUoHzu3v3LmbOnImGDRvCz88Pbm5uqFq1KgYPHoy9e/daPNbamVTMzQIiSRLCw8Pl35WDZfPKuyDvub1792Lw4MGoWrUq3Nzc4Ofnh4YNG2LmzJm4e/eu2eMM7ykhBJYvX442bdogMDAQPj4+aNasGb7//nu947KysrB06VK0aNECAQEB8Pb2RuvWrbF+/XqL100nPj4eb731Fpo0aYKAgACo1WpUqlQJ/fr1w44dO6zKwxpXrlzB1KlT0bhxYwQGBsLFxQXlypXDk08+iTlz5uC///4ze+yJEycwevRoeWCpt7c3ateujYkTJ1p8L5nq/xoVFYUePXogODgYarUa4eHhGDNmDG7cuGGx/Ldu3cK0adPQqFEj+Pr6wsXFBUFBQahbty4GDhyIVatWITk5WU6vuy+vXr0KAFi9erXR/ae8Z02V9ZdffkG3bt0QEhICZ2dno3t8//79mDlzJiIjIxEcHAxXV1f4+PigVq1aGDNmDE6fPm3xNeU1VsKw//S5c+cwatQoVK1aFWq1GkFBQejdu7fZcTW2vAcTExMxb948tGzZUn7vlS9fHrVq1ULv3r3x1Vdf4c6dOxZfV6kgyKFdvnxZABAAxNChQ/NMf+nSJTn9yy+/bDHt6NGj5bT//fef0f4qVaoIAKJKlSpWl3ft2rVynqNGjbL6OHPq1KkjAIjOnTvL22bPni0ACBcXF3Hv3j2Lx2dnZ4uJEycKSZLkcpn6Mby2yuu+cuVKo3yTk5NFZGSkACBUKpX4+uuvjdLojp89e7bJsuXk5ORZNmdnZ7Fs2TKjY5csWSKniY2NtXgNPv/8czntv//+azGtUkJCgnzcV199ZbT/ypUremUdP368URqNRiO8vLwEADF16lS9fTExMfKxMTEx8vaVK1darCvdz+XLl+Vj2rVrJwCIdu3aib1794py5cqZPe6jjz6y+hrYqm7dugKA8Pb2LlA+27dvFz4+Phavw7hx40ROTo7J4/O6B3WU18/U8ZZ+DPO29T0nRO57Yty4cRaP8/X1FX/99ZfJ16G8p/766y/Ro0cPs/m8/vrrQgghHjx4INq2bWs23bx58yxeux9++EF4enpaLPOIESOERqOxmE9ePvroI+Hi4mLxPIb1p/P+++8LlUpl9ji1Wi1Wr15t8ljDZ+G0adPM5lO+fHlx+vRpk/n8/fffed7LAMSWLVvkY3T3pbWvWVnWFStWiMGDB1tMb82zxsnJSSxZssRsvej+HgGmwzXd39GhQ4eKX375RXh4eJg9z08//WR0fH7fg6dPnxYhISF5HrN48WKzr6m0YIDs4PIbIG/ZskVO/+mnn1pMu3DhQjnt77//brTflgBZGVRFRERYfZwphw8flvP67rvv5O0XLlyw+k3+0ksvyWkrVqwo5s2bJ2JiYsSRI0fE9u3bxbvvvivq16+frwA5ISFBNG7cWAAQrq6uYt26dSbPberhpTR27Fg5Tdu2bcWKFSvErl27xMGDB8U333wjateuLe//9ddf9Y5NTEwU7u7uAoAYPXq0xWvQoEEDAUDUq1fPYjpTatWqJQCI/v37G+1btWqV3gO3fv36Rmn2798v7//jjz/09pkLkB8+fChOnDghxowZI+8/ceKE0U9WVpZ8jO4P6eOPPy7KlSsnKlSoID744AOxd+9ecfDgQbFw4ULh5+cngNwPHSdPnsz3tbDWnTt3hLOzswAgmjVrZnM+R48eFa6urgLI/TA4ceJEERMTIw4ePCiWLVsmwsPD5eszZcoUk3nkdQ/qmAuQT5w4IbZv3y7n89577xnVw507d/SOsfU9J4QQb775pnxseHi4WLp0qTh48KCIiYkREydOlANEV1dXERcXZ3S88p5q3ry5ACAGDRokfv/9d3H48GHx448/iieeeEJOExUVJZ599lnh7OwsxowZI/766y9x+PBhsXz5cjnIcHJyMnu/rFu3Tv4gUK1aNbFw4UKxbds2cfjwYfHzzz+Lbt26yeeaOHGixTqw5J133pHz8fPzEzNmzBBRUVHiyJEjYufOneLjjz8WrVq1EpGRkUbHKj9Mly9fXnz88cciNjZW7N27V8yZM0cO7iVJMvl3QPksbNWqlXyfrF27Vhw6dEjs2LFDDBkyRE7TokULozwyMjLk6+nt7S2mTJki/vzzT3H48GERGxsr1q5dK1599VURGhqqFyD/999/4sSJE/KxPXv2NLr/lI07yrLWq1dPABBPPvmkXlm//fZbOf0333wj/P39xbBhw8SKFSvEnj17xJEjR8TWrVvFO++8I3/QliRJREdHm6wbawPkRo0aCTc3NxEeHi6++OILsX//fhEbGyvmzJkj3NzcBADh4+MjEhIS9I7P73tQ97fJxcVFjB07VmzZskX8+++/4sCBA+Lnn38Wb775pqhevToDZCr98hsgf/XVV3L6DRs2WEy7YcMGOe3SpUuN9tsSIAshRFhYmByIFMT48eMFAOHh4SFSUlL09un++DVt2tTs8b/++qv8+lq2bCkePnxoNu21a9f0fjcXIF+7dk3+A+vh4SG2bdtmNk9Lwclff/0l71c+sJXS09NFhw4d5DowbIHStY74+vqKtLQ0k3kcOXLE6g9MpuiC1ODgYKN9w4cPFwDkVjpJksT9+/f10ixYsEAOMpKTk/X2mQuQdfL6w6OkbGmqUqWKuHHjhlGaPXv2yMGMrvWwKEyePFkuy8cff2xzPk2bNpWv3fbt2432P3jwQP4Ao1KpTAZxBQ2Qhcj72xSlgrznjh8/Lrdy1qlTx+Sxf/75p5zG1IcP5T0FQHz22WdGaW7fvi28vb3lgFGSJLFp0yajdMeOHZPPZep+uXv3rvD19RUAxEsvvWS2hXjGjBlyHZ09e9bM1TDvyJEjcjkef/xxcf36dbNpDa9pQkKC3GIZEhJitF+Xvy5IDg0N1fvgKYR+/QO53wxqtVqjfEaOHCmnOXLkiN6+6OhoeZ8yADak0WhEUlKS0XZlK6wlhmUdMmSIybLq3LhxQ6Smpprdn5iYKAfabdq0MZnG2gAZgGjcuLHJ1/fDDz/IaRYuXGjxdVl6Dyq/QbYUAGu1WvHgwQOz+0sLBsgOLr8B8ocffiin//PPPy2m/eOPPyz+Ibc1QK5fv76cr6mHgTU0Go2oUKGCACBeeOEFo/1ffPGFfI4zZ86YzKNly5ZyIHvz5s18nd/UA+ns2bOicuXKcivOvn37LOZhKTjRBb59+vSxmMfp06flfAy/Vt69e7e874cffjB5/GuvvSaA3Ba3u3fvWjyXKevWrTN7nXUtmBs2bJD/bxhoPP3002Y/yBRVgPzbb7+ZTdeiRQsBQDRs2DDPPG2xf/9+ufU4LCzM4h9fSw4cOCC/nldeecVsur1798rpxo4da7S/uAPkgrznlN8Y7N+/32w6ZSB28OBBvX2GLcjmKFs8TX07oqPremHqftG16oaGhoqMjAyzeWg0GhEaGioAiBkzZphNZ87AgQPlD6CGgWdedB9QAZj8+l7nvffek9OtX79eb5+y/itWrGj2tZ49e1ZO9/nnn+vtW7NmTYH+JtgSIPv5+Rl9KLfF5s2b5TxNdenLT4B87Ngxk2m0Wq3cSt67d2+j/da+B/ft25fnucoSDtKjfMnIyJD/rxtJb45arZb/n56eXmhl8PLykv+fkpJiUx7btm1DQkICAMgzVyj1798fLi4uAGA04AbIHZimG/TQv39/hISE2FQOnSNHjuDJJ5/EtWvXEBwcjN27d6NVq1Y25ZWcnCwPSnv++ectpq1ZsybKlSsHAIiNjdXb17ZtWzz++OMAgJUrVxodm5WVhbVr1wIAevToIeeTH+3atZP/rxxId/36dVy+fBmSJKFdu3bywBdlmpycHOzbtw8Aim2eYz8/P3Tv3t3s/saNGwOAxcFMtrpz5w6ef/55ZGdnQ5IkrF69Gh4eHjblpRzYNWLECLPpWrduLQ/ULMzBYLYo6HtOV/7atWvLg31NGTVqlNExpgwYMMDsvvr16+crnan75bfffgMAPPPMM3rPUkPOzs5o2bIlAOP3cF60Wq08sDQyMlKew95auuvj5+eH5557zmw65axDlq7p888/b/a16gYZA8bXq2LFivL/TT2rikKPHj3MzjZjTmpqKq5cuYJTp07h5MmTOHnypPx3BsgdfG6runXryoN3DUmSJNdtQZ5NyuvsCAuKMECmfHFzc5P/n9c8wZmZmfL/DaeCKwhlUOzj42NTHqtXrwYAVKhQAU899ZTR/nLlyqFLly4Acme6EAazcMTFxcnbnnzySZvKoLNnzx60b98ed+/eRdWqVbF3716zDzprHD16VF7BcODAgSZnaFD+3Lt3D0DuSHlDuuBp586d8khvnd9++01eifGll16yqaxBQUGIiIgAoB/86v5fq1YtlC9f3mSAfOTIEXlEujLQLko1atSASmX+sRkQEADA9g9u5qSkpKB79+7yKP4PPvgAHTp0sDm/kydPAsj9kNugQQOLaXXB5IULF+w6N3hB3nOZmZm4cOECAFgMjoHchY90QYvuOpmi+/Boip+fX77SGd4vOTk5iIuLAwAsW7Ysz/ewbro/U+9hSy5fvozExEQAtj3HdNenUaNGeoGeoaCgIFStWlXvGFN0zwJz/P39ARhfrzZt2qBatWoAgAkTJqBZs2aYP38+9u3bV2T3rLXP6Hv37mHGjBl44okn4O3tjfDwcNSpUwd169ZF3bp19T5w657Ftsjr2hXGsyk8PFy+Tz799FPUrl0bs2bNws6dO5GWlmZzviUVA2TKF+Un5rxW8UpNTZX/r2z1LSjdQ8TZ2Tnfn+CB3ClqdKvxDRgwAM7OzibTDR48GABw7do1vcBMWQZA/1O1LVasWCEHeuvWrcNjjz1WoPx0LeP5ZeoBN3ToULi4uEAIIX+o0FmxYgUAIDQ0VP4wYQtdcLt79255m+7/usBY9+/x48fx4MEDvTROTk4F/pBirbxabHXBsy1LrJujm+7v8OHDAIDJkydjypQpBcpTdw0DAgLM3v86wcHBAAAhBB4+fFig8xZEQd5zynJXqFDBYloXFxcEBgYC+N91MsXSvaD8EGVNOsP75cGDB/Jc8/mR3yCloM8x3fXJ65oC/7uPbL2mwP+uV05Ojt52FxcXbNmyRf62499//8WMGTPQpk0b+Pn5oWvXrli7dq3RcQWhC9YtOXz4MCIiIjB//nycP38+z+lOC/JNq63XLr9+/PFH+RuL06dP491330XHjh3h5+eHtm3bYunSpXrfNJdmDJApX8LCwuT/5zUn5fXr1+X/V6pUqVDOn5CQgFu3bgGAPM9yfq1bt05u3V60aJHZVpn+/fvLxxTl0tM9e/aEk5MTgNyg3JY5k5WUD8Bly5bpLeNt6WfevHlGeQUFBeGZZ54BkPuVmu4Bf+vWLfz1118AgCFDhsjlt4Uu+I2Pj5fnQzYMkCtXroyqVatCCIG///5bL02DBg1s/iahpMvOzka/fv0QExMDIPer6o8++qjQ8jc3t2pZVhpes/I9PHLkSKvfw7r3ZHErCde0Vq1aOHHiBDZt2oSXXnpJns8+PT0d27dvx6BBg9C8eXObGxAM5fXMy8rKQr9+/XD//n24uLhg0qRJ2L17N27fvo2MjAx5+WjlXOZ5BdAlQWhoKP755x/s2LEDY8eORe3atSFJEjQaDfbs2YMxY8agTp06ZufYL00sNx0QGahVq5b8f3OLO5jab82CE9aIioqS/9+mTRub8rAl2P3555+xZMkS+VO6sr9tQQPaXr16YcCAAXjxxRdx/vx5dOjQAbt27bJ5tUNd6xeQ26pQp06dApVv5MiR2LRpEy5fvozdu3cjMjIS3333nfxHfPjw4QXK37Afsre3Ny5evCj3P9aJjIzEqlWrsGvXLjz77LPYs2ePvL0s0mq1GDx4sPxtR//+/W1aGt4U3det9+/fR3Z2tsVWZN3X9pIkGbWaSZIEIUSeLebKb5NsVZD3nLLceS1gkJ2dLXcd0l2n4qY8rxCiwO9hcwr6HAsICMDt27etWhRCdx8V5TV1cnJCr1690KtXLwC5r2nbtm1YsmQJDh8+jMOHD+Pll1/Gpk2biqwMOjt37pT7+3755ZdmV3+11KJeknXs2BEdO3YEkPsc2bFjB77++mvs3LkTly5dQv/+/XH06FE7l7Jg2IJM+RIeHi4PjlF+JW6KrqUvNDRU7n9WEEIILFq0SP69d+/e+c7j0qVL+OeffwDkdq/48ccfLf7oWlVTUlL0HqoNGzaUW010r7MgBgwYgNWrV0OlUuHs2bPo0KGDzS0dDRo0kMumG8RWEF27dpW/OdANgNH9++STT6JGjRoFyr9ixYpyP81du3bJ95Wu/7GOsh9yXFwckpKSANje/7gktHpZ8vLLL+Onn34CkDsg6IcffrDY/zk/dAFXVlaW3NfVnIMHDwLI7X9tODBX18XJUtcLIQQuXrxodr+19VCQ95xarZbv0wMHDlhMe/ToUWg0GgAossA0L66urqhduzaAwnkPmxMeHi73g7blOaa7PkeOHLHYJSQhIUEew1Cc17RixYoYPnw4YmNj0ahRIwDA1q1bjboyFMWz4NSpU/L/ld9GGjp06FChnzu/Cvr6AwMD0b9/f0RHR+PZZ58FkDtmQNfvv7RigEz5IkkSevbsCSC3hdjc8pX79++XW5B1yyUX1IIFC+Q/1o0aNbKp36uy9Xjy5MkYMGCAxZ8pU6bILbLKYwMCAuRZJtavXy93+yiIQYMGYeXKlVCpVDh9+jQ6duxo06CN8uXLy8sur1271uLSudZQqVRyK/HGjRuxbds2+eszWwfnGVL2Q9b19zZsGVb2Q/7ll1/kstna/1g54FQ5oLQkmDRpEr799lsAuS01GzZsyLOvcH506tRJ/r+uL7kpsbGx8lK4ymN0dEvUWvoj/+eff8oDwUyxth4K+p7Tlf/UqVPyc8QU3XVXHmMPukDj7Nmz2L59e5GcQ6VSyYPEdu/ene8WP931SUxMlN+TpixfvlzuPmCPa+ri4iI/Y7Kzs43uR909WJjPAeUHBnPfoGi1WnzzzTeFdk5bFeazUNeqDBRs0GFJwACZ8m3ChAly/6vXXnvN6NN4eno6XnvtNQC5A+kmTJhQoPM9evQIkydPxvTp0wHkdhtQ/hGzlhACP/zwA4Dc9et1U3JZ4uzsLH9dFx0drfc15NSpUwHkDozp27ev3KJpSl79tXWGDBmCb7/9FpIk4eTJk+jYsaP8dW9+zJw5E0DulG/PP/+8xQAlMzMTS5YssTiw4qWXXoIkSUhLS5ODZW9vb/Tt2zffZTNF2Q95/fr1ett0qlSpgipVqkAIgS+++AJA7jRZyhkD8kM5KEnZD9De5syZg08//RQA0KpVK/z6668Wp/myRbNmzdCkSRMAwDfffIPo6GijNElJSXj55ZcB5AZSY8aMMUqjCzoOHDhgsqUzPj5efhaYExgYKLdM51UPBXnPjRkzRm6BHz16tDwwVumvv/7C8uXLAeReo6ZNm1osT1EaP368PLh5+PDhei2Spvz+++84fvx4vs8zefJkqFQqCCEwYMAAi88qw33Dhw+Xu5298cYbuHnzptExx44dw/vvvw8g99tE3fO0MO3Zs8fitxRZWVnyN1NeXl5630wB/3sWFOZzQPnNmrkp0aZPn44jR44U2jltZe17MC4uzuI3TkIIeRo/SZIK5Ztje2IfZAezd+9evQeJ8hPexYsXjd7Iw4YNM8rj8ccfx5tvvokPPvgAhw4dQuvWrTF16lQ89thjuHTpEhYsWCC3RLz55pt5fgWv0Wj0pv7RaDRITEzElStX8M8//2Djxo1ygOfr64s1a9bke75OIPe16/qE9enTx+rj+vTpg+XLlyMnJwc//PAD3nzzTQC5X3uPGDECy5cvxz///INatWrh1VdfRevWreHj44N79+7h0KFDWLduHerXr2/1vJHDhw9HdnY2Xn75ZRw/fhxPPfUUoqOjrRo1rdOtWzeMHz8en3/+Of7++2/UrFkTr7zyCtq0aYPAwECkpqbi4sWL2LNnD3755Rc8fPgQQ4cONZtf1apV0alTJ0RFRcl9Cfv16wdPT0+ry2SJsptEUlKSUf9jncjISKxevVoOjArS/1g5z/TEiRPx1ltvoWLFivK3HVWrVi3UVltrLF68GHPnzgWQG0x8+OGHuHz5ssVjnnjiCYtTbJnzzTffoHnz5sjKykK3bt3w2muvoUePHvD09MTRo0fxwQcfyO+XyZMnm/xqfPTo0fjyyy+RnZ2NHj16YNasWWjTpg2ysrKwb98+LFy4EBqNBjVq1DD7dauzszOaNm2Kffv2YcWKFWjYsCEaNGggv6aAgAC532pB3nN169bFG2+8gY8++gjHjh1Do0aNMHXqVDRs2BCpqanYsmULFi1ahJycHLi6uhZaf29bBQUFYfXq1Xj++edx+/ZtNGnSBMOGDcPTTz+NsLAwaDQa3LhxAwcPHsTGjRvx33//YcuWLfmeIrJBgwaYO3cu3n77bZw/fx5169bFuHHj0L59ewQGBiIxMRFxcXH45Zdf4OTkJA8YBXK/rfroo48wbtw43LhxA40bN8a0adPQqlUrZGdnY8eOHfjoo4/w6NEjSJKEr7/+2qZ7NS/R0dF499138eSTT6J79+6oV68eypcvj/T0dJw/fx5Lly6VA9ERI0YYva9btWqFmJgY/Pvvv/jggw/w9NNPy882d3d3hIaG5rtMXbp0QYUKFZCQkICZM2fiypUr6N27N8qVK4eLFy/KH0xbt25dpN1orGHtezAuLg7Dhw9H06ZN0aNHDzRq1AjBwcHQaDS4fPkyVq5cKY8TevbZZws8w5PdFe+6JGRvQ4cOlVfCsebHnJycHPHSSy9ZPHbEiBEiJyfHbB7KFYDy+nF2dhbPPfecuHLlis2vXblCVmxsrNXHZWVlCT8/PwFA1K1bV29fdna2ePXVV+Ulhs39GK7QZM3KRV999ZWcb+PGjY2WxtUdb24VM61WK+bOnSuvvGbpx9PT0+xy0jrKVe8A5LnSX35Vr15dzrt27dom06xcuVKvDJs3bzabX14r6QkhRL9+/cxek8uXL8vpLK0Ep5Sf1flMUa7YZ+2Pspz5tX37duHj42Mx/3Hjxll8Hy9cuNDssQEBAeLvv//O8/pt3brV7HvI8P629T0nRO5za+zYsRaP8/X1Nbn0thDW3VNC6N+nlurHmvvlt99+EwEBAXneByqVSuzcudNsPnmZN29ens8Kc/U3b948eblqUz9qtVqsXr3a5LH5WUnR3Ip3yuto6adnz54mn3M3btwwe42Vrzk/ZRVCiG3btgk3Nzez5YmMjBQnT560mKe1K+nltQqg7m+/uZVrrXkPGj5/zf20atXK5KqApQ27WJBNVCoVli9fjt9//x09e/ZESEgIXF1dERISgp49e+KPP/7At99+a9OgIg8PD1SsWBF169bFiy++iEWLFuHq1av4+eefUaVKFZvKm5GRgQ0bNgDIbZnLa7EAJRcXF/To0QMAcOLECb1+ek5OTli8eDEOHTqE0aNH4/HHH4enpydcXFwQHByMzp07Y+HChfj444/zXeZXXnkFixcvBpA7n2bnzp0tfqVsSJIkzJo1C+fPn8eUKVPQpEkTBAQEwMnJCd7e3qhVqxYGDRqE1atX4/bt23ku5tKrVy/5K9+IiAibV/ozx3DGClOU21UqFdq2bVugc/7www/48MMP0axZM/j6+hbaILjSonPnzrh48SJmzJghT5enVqtRuXJlDBo0CHv27MEXX3xh8bpMnDgR27ZtQ5cuXeDv7w+1Wo3w8HCMGzcOR48etaqPePfu3REdHS0/Syy1MhbkPadSqbBkyRL8/fffGDRoECpXrgy1Wg0fHx80aNAAM2bMwIULF9C5c2frLmAx6NGjBy5fvoyPP/4YHTp0QFBQEFxcXODu7o7w8HA888wzWLhwIa5cuYL27dvbfJ4ZM2bg9OnTmDBhAurUqQMfHx84OzujfPnyaNeuHd577z2Tq4rqjj169ChGjRqFxx57DO7u7vD09ETNmjUxfvx4nD17FkOGDLG5bHmZPHkyfv75Z4wZMwYtWrRA5cqV4ebmBjc3N1StWhX9+vXD1q1bsXnzZpPPudDQUBw8eBAjRoxA9erV9frkFkSXLl1w6NAhvPjii/J9rbueX3/9NaKjowvtW7iCsuY9OHDgQPzxxx+YOHEi2rRpg/DwcHh4eMDV1RVhYWF49tlnsWbNGuzZs0dvNqXSShKiFEy8R0R2d+HCBXm2iQULFhR4sQoiIqKSyrGaS4jIZrrZDpydnYu0NYiIiMjeGCATUZ4SExPx9ddfA8jtaqFbNpaIiKgs4iwWRGRSQkICkpOTcevWLcyZMwcPHjyAJEnydHtERERlVZlsQf7777/Ro0cPhISEQJIkbN68Oc9jdu3ahUaNGkGtVqN69epWT8lFVFZNmTIFNWrUQLt27eSpncaOHSuvSEVERFRWlckAOTU1FfXr18eSJUusSn/58mV0794d7du3R1xcHCZMmICRI0cW2epFRKWJq6sratasiYULF+Kzzz6zd3GIiIiKXJmfxUKSJGzatMni6j1Tp07F77//rrdYxYABA5CYmIht27YVQymJiIiIqKRgH2QAsbGxRuvDd+nSxeISyZmZmXprlmu1Wjx48ACBgYHySlxEREREVHIIIZCSkoKQkBCLc7wzQAYQHx+PoKAgvW1BQUFITk5Genq6yYnF58+fLy8JS0RERESlx/Xr1xEWFmZ2PwNkG02fPh2TJk2Sf09KSkLlypVx/vx5BAQE2LFkjkuj0SAmJgbt27e3uBIXFR3Wgf2xDuyL19/+WAf2V5LrICUlBeHh4fD29raYjgEygODgYNy5c0dv2507d+Dj42N2+V21Wg21Wm20PSAgoEwssVgaaTQaeHh4IDAwsMS9IR0F68D+WAf2xetvf6wD+yvJdaArT17dYcvkLBb51bJlS0RHR+tti4qKQsuWLe1UIiIiIiKylzIZID969AhxcXGIi4sDkDuNW1xcHK5duwYgt3uEcqncV155Bf/99x+mTJmCs2fP4ssvv8T69esxceJEexSfiIiIiOyoTAbIhw4dQsOGDdGwYUMAwKRJk9CwYUPMmjULAHD79m05WAaA8PBw/P7774iKikL9+vXxySef4Ntvv0WXLl3sUn4iIiIisp8y2Qc5MjISlqZ3NrVKXmRkJI4ePVqEpSIiIiKi0qBMtiATEREREdmKATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZFCmQ2QlyxZgqpVq8LNzQ3NmzfHwYMHzaZdtWoVJEnS+3FzcyvG0hIRERFRSVEmA+R169Zh0qRJmD17No4cOYL69eujS5cuSEhIMHuMj48Pbt++Lf9cvXq1GEtMRERERCVFmQyQFy5ciFGjRmH48OGoVasWli5dCg8PD6xYscLsMZIkITg4WP4JCgoqxhITERERUUnhbO8CFLasrCwcPnwY06dPl7epVCp06tQJsbGxZo979OgRqlSpAq1Wi0aNGuH9999H7dq1zabPzMxEZmam/HtycjIAQKPRQKPRFMIrofzSXXdef/thHdgf68C+eP3tj3VgfyW5DqwtkySEEEVclmJ169YthIaG4p9//kHLli3l7VOmTMHu3btx4MABo2NiY2Nx4cIF1KtXD0lJSfj444/x999/49SpUwgLCzN5njlz5mDu3LlG29euXQsPD4/Ce0FEREREVCjS0tLwwgsvICkpCT4+PmbTlbkWZFu0bNlSL5hu1aoVatasiWXLluHdd981ecz06dMxadIk+ffk5GRUqlQJ7du3R2BgYJGXmYxpNBpERUXhqaeegouLi72L45BYB/bHOrAvXn/7Yx3YX0muA903/nkpcwFyuXLl4OTkhDt37uhtv3PnDoKDg63Kw8XFBQ0bNsTFixfNplGr1VCr1SaPLWk3g6NhHdgf68D+WAf2xetvf6wD+yuJdWBtecrcID1XV1c0btwY0dHR8jatVovo6Gi9VmJLcnJycOLECVSsWLGoiklEREREJVSZa0EGgEmTJmHo0KFo0qQJmjVrhs8++wypqakYPnw4AGDIkCEIDQ3F/PnzAQDvvPMOWrRogerVqyMxMREfffQRrl69ipEjR9rzZRARERGRHZTJALl///64e/cuZs2ahfj4eDRo0ADbtm2Tp267du0aVKr/NZ4/fPgQo0aNQnx8PPz9/dG4cWP8888/qFWrlr1eAhERERHZSZkMkAHg1Vdfxauvvmpy365du/R+//TTT/Hpp58WQ6mIiIiIqKQrc32QiYiIiIgKggEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkYJzYWSi1Wpx+PBhXL16FWlpaRgyZEhhZEtEREREVOwK3IK8ePFiVKxYES1atED//v0xfPhwvf0PHz5EnTp1EBERgTt37hT0dERERERERapAAfK4ceMwYcIE3L17F97e3pAkySiNv78/GjVqhAsXLmDDhg0FOR0RERERUZGzOUDetm0bvvrqK3h5eWHTpk1ITExE+fLlTaZ94YUXIITAjh07bC4oEREREVFxsDlAXrp0KSRJwjvvvIOePXtaTNuyZUsAwIkTJ2w9HRERERFRsbA5QD5w4AAA4KWXXsozra+vL3x8fBAfH2/r6YiIiIiIioXNAfKDBw/g6+sLb29v606kUkGr1dp6OiIiIiKiYmFzgOzj44Pk5GRoNJo80z548ABJSUkoV66cracjIiIiIioWNgfIdevWhRBC7mphyY8//gghBJo0aWLr6fJtyZIlqFq1Ktzc3NC8eXMcPHjQYvoNGzYgIiICbm5uqFu3Lv74449iKikRERERlSQ2LxTy/PPPY9euXZgzZw7++usvqFSmY+1jx45h5syZkCQJAwcOtLmg+bFu3TpMmjQJS5cuRfPmzfHZZ5+hS5cuOHfuHCpUqGCU/p9//sHAgQMxf/58PPPMM1i7di169eqFI0eOoE6dOvk6971792zqSuLl5QV3d3ezeQoh8p0nAHh4eMDT09PkvgcPHiAnJ8emfN3c3Mx2r0lMTLTqmwVTXF1d4evra3JfUlISsrKyzB6r0WiQlJSEu3fvwsXFRW+fi4sL/Pz8TB6XkpKCjIwMm8rr5OSEgIAAk/tSU1ORlpZmU76SJJn9xiU9PR2PHj2yKV8AZmebyczMRHJyss35BgYGmtyelZWFpKQkm/P19/eHs7Pxoyo7OxsPHz60OV9fX1+4uroabddqtbh//77N+fr4+ECtVpvcd/fuXZvztfYZYel9YIojPSMsKaxnhOH15zPifwIDA03GCtY+I7JztMjI1gIGfw/9DJ4RGk02MrKBh48ykPrI9vecj4VnxIMCPCO8LTwj7hXgGeFp4RlxvwBxhHsezwitiWeEJjsbCfeTcPnGbbiYeH4DgDqPZ0S2jc8IFzPPCE9XZ6SkJFv9fJeEjVdMo9GgYcOGOHPmDCIjIzFx4kS89NJLuH//Ps6ePYsrV65gy5YtWL58OdLT09GyZUvs3bvX5FzJha158+Zo2rQpvvjiCwC5N3OlSpXw2muvYdq0aUbp+/fvj9TUVGzdulXe1qJFCzRo0ABLly616pzJyclmH9rW+OKLLzBu3DiT+8qXL4979+7ZlO/s2bMxZ84ck/tq166N06dP25Tv2LFjsWTJEpP7IiMjsXv3bpvyff75583Ol923b19s3LjRpnzbtWuHXbt2mdw3btw4fPnllzblW6tWLZw6dcrkvjlz5mDu3Lk25VuuXDmzwdSSJUvw6quv2pQvALMPyQ0bNqBfv34255uQkAA/Pz/88ccf6Natmxyc7dq1C+3bt7c535MnT6J27dpG20+dOpXvD7BKMTExiIyMNNp+9+5dkx+krfXuom/RoZvpmX1aVzcdeFij1vMTUPXJ50zui5rRA5pU2z6EVO8yDDWefgmm7oq9HwzBo/grNuVbqXUv1Hp+osl9B78Yj4eX4mzKN6h+OzQY9o7JfXGrZuHOMduePf6PNUDTcZ9DmLgSZ37+DDf2bbYpX8+gqmg5ZZXJfZe2r8Tlv1bblK+Lpy/avvOryX3X927C+U2f25QvAHT4OMbk/ZBwbBdOfW/bMw0AWs/eBFcvP6PtDy/FIW6p6XvFGhVfWgLX8lWMtmfdvYrbK0z/XbVG0MD34Va5ntH2nLQk3Fg8yOZ8y/WcBs+INib3XV3wjM35Bjz1CrwbmT7++qIXoE237cONb+uB8Gtj+vXe+nYsNPev2ZSvV8PuCOw8xuS++LXTkHn9pE35ejzRGuV7TTfavmtyJN4cM0yOI5KSkuDj42M2H5tbkF1cXPD777+ja9euiImJ0Qs+IiIi5P8LIVC3bl38/PPPxRIcZ2Vl4fDhw5g+/X8XR6VSoVOnToiNjTV5TGxsLCZNmqS3rUuXLti8ebPZ82RmZiIzM1P+vSCfqgEgJyfH5hYVW/O19dMkkPuho7jzLcggTyFEsedra8ubTnHnm52dXeB8dXkrz1HQfM/eSkSm5z1oBaDVCjzKzMbDNA1On7btoayzfM8lRD80btlLTXpQoHw/j76Ib2/m3fUsv+KTMpB6y/RzJkdr+3vu7qNMaG6bzjdTY/t742GaBmfjU0zuS8uy/Z5Iycg2m29Khu35pmVl49wd0/kmptn+bM7K1uJCgunW3MRU21q7gdw6v2gm35RHmSa3W+vS3VST21OTC5bv1QdpcMow/lYjIym9QPkSWSM7Oztff+9tDpABoEqVKjh8+DA++eQTrFixAlevXtXbHxoailGjRuGNN94w2zxf2O7du4ecnBwEBQXpbQ8KCsLZs2dNHhMfH28yvaVp6ebPn29z66App06dMtvv2davDAHgwoULZvMtyFdwV69eNZtvQb6ejo+PN5tvQaYJvH//vtl8De/b/Hj06JHZfC9cuGBzvllZWWbzPXnSdIu1terM2W5ye8qZYwXKt82HMXD28AXgBBzYKW9Pu1qw+c9f++kYXMsnGm3Pumt7vQFA1JkEuKVeN9qek2Z7dxAiIiqZdu3ala84okABMpDbf+3tt9/G22+/jVu3buHWrVvIyclBcHAwqlQx/sqjrJg+fbpeq3NycjIqVapkc361a9dGt27dTO4z1QfKWjVq1DCbr7KVPb+qVKliNt+FCxfanG9wcLDZfL/77jub8w0MDDSb77Zt22zO18vLy2y+hw4dsjlfV1dXdOvWDUII3ErKwImbyTh1KxknbiZj9828+5Vakplj+pscjbZg3/BotBK0JvLItr1xk4iIqFBERkbizw3WxxE290EuqbKysuDh4YGNGzeiV69e8vahQ4ciMTERv/5q3GercuXKmDRpEiZMmCBvmz17NjZv3oxjx6xrVdP1QT579qzZARmWcJBeroIO0tuxYwc6depU6IP09l68iw+3ncOtRIOvAiUVnNxNXwdtVgZEtu1fSTp5+ELAaDwKtJpMCI1tgwp1+ZoisjXQZtk2YAgAVO7ekCTjATgiRwNtZgHydfOCpHIyzlebA22G7d+CqNQekJyMP2wIoYU23fRX7TquTip4uRmXCQCc1B5QOZv+UKtJTYSTSkIFLzWCfNwR5KtGsI8b/DxcYWacs8zN3QNqN9PPiOSHD+S+szk5Wpw6eRK169SBs5PpMiqp3dzh5u5htF2SgJSkh9Dm2NbNwlWthrunl8l9j1KSkKPJX3cIXRc9ZxcXeHqb7jeYmpKc58Aecx8DnZyd4eVj+r2RlvoImkz997K5HoM5OTk4ceIE6tatCycnJ6icVPD29TeZNiMtDVk2DhCGBPj6/29grLI8mRnpyLB18B8A3wDTA241WZlIs/abRxPXx8fP3+QgPY0mC6kppt9zamcV3Fyc4O7qBHcXJ7i5qKAyuPi+fvqD9LKzs7Fz5060bdsWqY8sv5ct8fYxP0jv4QPbvyn18jY/SO/+PdsH6Xl4mo8jHtwvwCA9dw94mIkjHj40PUgvOzsHe/buwZNtnoSzs+nnkFrtBi8zcURSUgEG6bm4wsdEHFHeW420Rym4f/8+HnvssaLrg1xSubq6onHjxoiOjpYDZK1Wi+joaLMDm1q2bIno6Gi9ADkqKkpeIjs/ypUrZ3Y0v62Kav5oWwJ5a5gLRAsqr0GQGo0Gvr6+KF++vFWj93W8vb0tLnhzOykdM/44jHSNKyR3Ew9Lc88cFzdILm5Wl8PafFUuasDF9MO1ICRnFzg52z7Q1Gy+Ti5mg/I8j5UAJ0mCSpKgUuX+393VGf4eLvD3cIWvRwj83F3g4mzzjJVW8XFzQdVAD1QJ9ETVch4I8naDSlX0YyqsFyb/T6PR4I+0a+gWWSdf7wPTQgt4fHHnG1JE+VpPo9HA99FVdGtbuxCuP9lCo9HAXw1UCvSCS7DpDycFFRZQNN1GQ/wqF1G+tn/DbTlf0+9ljUaDa+d8UeexMJveByF+pgP9gvL19bV6PJzNAfK1a7YNkKlcuWgqX2nSpEkYOnQomjRpgmbNmuGzzz5Damoqhg8fDgAYMmQIQkNDMX/+fADA+PHj0a5dO3zyySfo3r07fvrpJxw6dAhff/11kZeVSr4Ff55FuqZgA+OKgr+HC+qE+qJemC/qhvqivLc1QXPeDwZrnh3mkmRnZ+Off/5B69atTU7NpjzezcVJbh1SO6vg6qyCk0qSg+Lc/5ekIJSIiByFzQFyeHh4vo+RJKnAI9qt0b9/f9y9exezZs1CfHw8GjRogG3btskD8a5du6b3NU+rVq2wdu1azJw5EzNmzECNGjWwefPmAk0hRWXDkWsPsTnulr2LAV93F9QL880NiENz/w3zdy9xAaRGo8Etb6B+mC9bz4iIqNSyOUC2pS9LcXZ3fvXVV812qTA1H27fvn3Rt2/fIi4VlSZarcDcLfrzRHu7OWP1S83g7pJ3387C4u3mjFC/khcMExERlVU2B8iXL1+2uD8pKQkHDhzAp59+irt37+L7779HzZo1bT0dUbHbHHcTx64n6m0b37EGGlUumj5tREREVDLYHCBbM4VbvXr1MHjwYHTs2BEjRozA0aNHbT0dUbFKzczGgm3682ZXK+eJIS2r2qdAREREVGyKdug3cqcDW7RoEW7fvo158+YV9emICsXS3Zdwx2DVqLe614RrEc+WQERERPZXLNO8NW7cGJ6entiyZQs+/9z29eHJcaVlZeNWYjpuPEzHzcR03E7MQIbBzBJarRb/XVHh2J/nTM61aS0B4If9+qu0PVmjHDpEVLA5TyIiIio9iiVA1mq1yMnJwe3bt4vjdFSGHL32EDM3n8SpW8lWHqHCrtsFW4LYkJNKwqxnanGQHBERkYMolu+LY2JikJGRUWQLSFDZJITAlI3H8xEcF40Xm1dGjSDzC4kQERFR2VKkAbJGo8H69esxdOhQSJKEDh06FOXpqIy58TAdFxJsX0q4MAT7uGFCp8ftWgYiIiIqXjZ3sahWrZrF/RkZGUhISIAQAkII+Pr6Yvbs2baejhzQ8RtJRtsCPF0R6ueOED83+Li56K36ptUKXL9+HZUqVSqUZYDLe6sxoGll+HsaLy9NREREZZfNAfKVK1esTtumTRssXrwYjz/Oljiy3rEbiXq/t3u8PFa/1Mxseo1Ggz/+uIpu3WpzFTciIiKymc0B8sqVKy1n7OwMf39/1K9fH6GhobaehhxYnMEiHfUr+dmlHERERORYbA6Qhw4dWpjlINKToxU4eVO/i0WDSr52Kg0RERE5Eq56QCXSxYRHSMvSn+e4XpiffQpDREREDoUBMpVIxwy6V4T6uaOcl9o+hSEiIiKHYlUXi2vXrhXaCStXrlxoeVHZFWcwQK8B+x8TERFRMbEqQA4PDy+Uk0mShOzs7ELJi8q24wYBcr0w9j8mIiKi4mFVgCyEKJSTFVY+VLZlaHJw9naK3jbOYEFERETFxaoA+fLly0VdDiLZqVvJyNb+78OUSgLqhrIFmYiIiIqHVQFylSpVirocRDLDAXrVK3jBU23zjIRERERE+cJZLKjEMex/XJ/TuxEREVExYoBMJc6xG/oLhLD/MRERERWnQvveOiEhATdu3EBqaqrFwXht27YtrFNSGZSUpsHle6l62zjFGxERERWnAgfIX3zxBRYtWoRLly7lmZbTvFFejt9M1Pvd1VmFJ4K97VMYIiIickgFCpAHDBiADRs2WD19G6d5o7wYDtCrHeIDFyf2BCIiIqLiY3Pk8dNPP2H9+vXw8fHBxo0bkZqa+7V4cHAwsrOzcePGDaxcuRLVq1dHuXLlEB0dDa1WW2gFp7Ip7rpB/2MO0CMiIqJiZnOAvGrVKkiShHfffRfPPfcc3N3d/5epSoWQkBAMHToUR44cQaVKldCrVy9cvHixUApNZZMQAscMZ7CoxPmPiYiIqHjZHCAfPXoUAPDiiy/qbTdsJfby8sIXX3yBlJQULFiwwNbTkQOIT87A3ZRMvW1sQSYiIqLiZnOAnJiYCG9vb/j5+cnbXFxc5K4WSi1btoSHhwd27Nhh6+nIARj2P/Zxc0bVQE/7FIaIiIgcls2D9AIDA5Genq63zc/PD/fu3UNiYqJe4KwTHx9v6+molDt1KwnbT8YjISUTD9Oy8DBVg4dpWUjJ+N+sJqmZ+jOc1Avzg0olFXdRiYiIyMHZHCCHhobiyJEjePToEby8vAAANWvWxJ49exATE4PevXvLaY8cOYK0tDT4+/sXvMRU6py6lYTeX/6DrOz8DdJk/2MiIiKyB5u7WDRq1AgA8O+//8rbunfvDiEEJk+ejH///RcajQaHDh3C0KFDIUkSWrduXfASU6mz/dSdfAfHANC4Cj9QERERUfGzOUDWBcMbNmyQt40ZMwahoaG4fPkyWrRoATc3NzRv3hynTp2Cs7Mz3nrrrUIpNJUuD1Iz805koFnVALStUb4ISkNERERkmc1dLLp164aYmBh4eHjI27y8vLBz504MGzYMsbGx8vbKlStjyZIlaN68ecFKS6WSsp8xADQPD8DTdYLh7+kKHzcXo37GPm7OqBPqC2cuEEJERER2YHOA7OzsjHbt2hltr1GjBvbt24cbN27g+vXr8PX1Rc2aNSFJHGzlqAwD5LaPl8ew1uF2Kg0RERGRZQVaatqSsLAwhIWFFVX2VIokp2v0fvdxK7LbjoiIiKjAbP4Oe82aNUbTvBGZYtiC7O3mYqeSEBEREeXN5gB58ODBCA4OxksvvYSYmJjCLBOVMSkZBi3I7mxBJiIiopLL5gDZ3d0dKSkpWL16NTp16oQqVargrbfewtmzZwuzfFQGsAWZiIiIShObA+SEhASsWrUK7du3hyRJuH79Oj744APUrl0bTZs2xRdffIH79+8XZlmpFNJqBR5lGQbIbEEmIiKiksvmANnT0xNDhgzBjh07cO3aNTk4FkLg8OHDGD9+PEJCQtCrVy9s3LgRWVlZhVluKiUeZWVDCP1tbEEmIiKikqxQJpoNCQnBlClTcPz4cRw9ehQTJ05EUFAQNBoNfvvtN/Tv3x8VK1bEmDFjCuN0VIoYdq8A2IJMREREJVuhr8RQv359fPLJJ7hx4wa2bduGQYMGwcPDAw8fPsTXX39d2KejEs5wijdJArxcGSATERFRyVVkS5WpVCo0aNAADRs2ROXKlYvqNFTCGbYge7k6G62cR0RERFSSFHpTXkZGBjZt2oTvv/8eO3bsQE5ODsT/d0Jt0KBBYZ+OSjjjKd7Y/5iIiIhKtkILkGNiYvD999/j559/xqNHj+SgOCQkBC+88AKGDBmCOnXqFNbpqJQwnuKN3SuIiIioZCtQtHLmzBl8//33WLNmDW7cuAEAEELAw8MDvXv3xpAhQ9CpUydIEr9Sd1SGLcgMkImIiKikszlaadKkCY4ePQogNyhWqVSIjIzEkCFD0KdPH3h6ehZaIan0SuYiIURERFTK2BwgHzlyBABQq1YtDB48GIMGDUJYWFihFYzKBnaxICIiotLG5mjl9ddfx+DBg9G4cePCLA+VMcnsYkFERESljM3RymeffVaIxaCyyrgFmV0siIiIqGQrsnmQiQAT07wxQCYiIqISjgEyFSn2QSYiIqLShgEyFSlO80ZERESlDQNkKlKGLcjsYkFEREQlHQNkKlLJ6WxBJiIiotKFATIVmRytQGpWjt42zmJBREREJR0DZCoyjwy6VwBsQSYiIqKSjwEyFRnDRUIABshERERU8jFApiJjOEBPJQGergyQiYiIqGQrlGjl1q1bOHHiBB48eACNxrjVUGnIkCGFcUoqBQynePNSO0OlkuxUGiIiIiLrFChAPnHiBF577TXs2bPHqvSSJDFAdiBcZpqIiIhKI5sD5HPnzuHJJ59ESkoKhBBwdXVF+fLl4ezMr9Apl2EfZPY/JiIiotLA5ohlzpw5SE5ORkhICJYuXYqnn34aTk5OhVk2KuW4SAgRERGVRjYHyDExMZAkCd999x06dOhQmGWiMoLLTBMREVFpZPMsFklJSVCr1YiMjCzE4lBZYtwHmQEyERERlXw2B8gVK1aEk5MTVCrOFEemJRt2sXBnFwsiIiIq+WyObnv06IG0tDQcPXq0MMtDZQi7WBAREVFpZHOA/NZbb6FcuXKYMGECMjMzC7NMBfLgwQMMGjQIPj4+8PPzw4gRI/Do0SOLx0RGRkKSJL2fV155pZhKXHZxmjciIiIqjWxu0svIyMDKlSsxePBgNGrUCJMnT0azZs3g7e1t8bjKlSvbekqrDBo0CLdv30ZUVBQ0Gg2GDx+O0aNHY+3atRaPGzVqFN555x35dw8PjyItpyPgNG9ERERUGtkcsYSHh8v/T0xMxMiRI/M8RpIkZGdn55nOVmfOnMG2bdvw77//okmTJgCAxYsXo1u3bvj4448REhJi9lgPDw8EBwcXWdkcEVuQiYiIqDSyOUAWQhTLMfkRGxsLPz8/OTgGgE6dOkGlUuHAgQPo3bu32WPXrFmDH374AcHBwejRowfefvtti63ImZmZel1LkpOTAQAajSbP5bYdRUq6/nXwcEaRXhtd3rz+9sM6sD/WgX3x+tsf68D+SnIdWFsmmwPky5cv23pokYmPj0eFChX0tjk7OyMgIADx8fFmj3vhhRdQpUoVhISE4Pjx45g6dSrOnTuHX375xewx8+fPx9y5c422x8TEsHvG/3uY6gRAkn8/efRfpF4s+vNGRUUV/UnIItaB/bEO7IvX3/5YB/ZXEusgLS3NqnQ2B8hVqlSx9dB8mzZtGhYsWGAxzZkzZ2zOf/To0fL/69ati4oVK6Jjx464dOkSHnvsMZPHTJ8+HZMmTZJ/T05ORqVKldC+fXsEBgbaXJayIjtHi/GxO/S2dY5sixpBXkV2To1Gg6ioKDz11FNwcWF3DntgHdgf68C+eP3tj3VgfyW5DnTf+OelVIyaeuONNzBs2DCLaapVq4bg4GAkJCTobc/OzsaDBw/y1b+4efPmAICLFy+aDZDVajXUarXRdhcXlxJ3M9hDqibLaFuAt3uxXBvWgf2xDuyPdWBfvP72xzqwv5JYB9aWp1AD5KtXr8oBaoUKFQqtlbl8+fIoX758nulatmyJxMREHD58GI0bNwYA7Ny5E1qtVg56rREXFwcgdzEUsk1yuvFgTM5iQURERKVBgZfBu337Nl5//XVUqFAB1apVQ4sWLdCiRQtUq1YNFSpUwIQJE3D79u3CKGueatasia5du2LUqFE4ePAg9u3bh1dffRUDBgyQZ7C4efMmIiIicPDgQQDApUuX8O677+Lw4cO4cuUKfvvtNwwZMgRt27ZFvXr1iqXcZZHhFG9OKgkerk52Kg0RERGR9QoUIO/btw/16tXDkiVLcO/ePQgh9H7u3buHxYsXo379+vjnn38Kq8wWrVmzBhEREejYsSO6deuGNm3a4Ouvv5b3azQanDt3Tu6k7erqih07dqBz586IiIjAG2+8gT59+mDLli3FUt6yynCKNy+1MyRJMpOaiIiIqOSw+TvvhIQEPPvss3j48CF8fHzwyiuv4KmnnkJYWBgA4MaNG9ixYweWLVuGe/fu4dlnn8Xp06eNZpkobAEBARYXBalataredHOVKlXC7t27i7RMjojLTBMREVFpZXPU8sknn+Dhw4eIiIhAVFQUQkND9fY/8cQT6NixI1577TV06tQJ586dw8KFC/HBBx8UuNBU8nGRECIiIiqtbO5i8fvvv0OSJHzzzTdGwbFSSEgIvvnmGwghsHXrVltPR6UMW5CJiIiotLI5QL5y5Qo8PT3RunXrPNO2bt0anp6euHr1qq2no1LGsAXZhy3IREREVEoUeBaL/Cjqpaap5DCcxcKHLchERERUStgcIFetWhWpqanYv39/nmljY2ORmpqKqlWr2no6KmWM+yAzQCYiIqLSweYA+emnn4YQAqNHj8bdu3fNpktISMDo0aMhSRK6detm6+molOEgPSIiIiqtbG7Wmzx5MpYvX45Tp06hZs2aGDNmDDp27CgP2Ltx4waio6OxbNky3L9/H35+fnjjjTcKreBUshl2sWALMhEREZUWNkctQUFB2LRpE3r37o0HDx7g/fffx/vvv2+UTggBPz8/bN68GUFBQQUqLJUebEEmIiKi0qpAg/TatWuH48eP4+WXX4a/v7/RSnr+/v4YM2YMTpw4gbZt2xZWmakU4DRvREREVFoVOGoJCwvDV199ha+++gqXL19GQkICAKBChQoIDw8vcAGpdDKa5s2dLchERERUOhRqs154eDiDYgLAPshERERUehXrPMjkGDQ5WmRotHrbOA8yERERlRYMkKnQGXavADhIj4iIiEoPqwJkJycnODk5oXbt2kbb8vPj7MxWREdgOEAPYBcLIiIiKj2silp0S0Qrl4rmstFkjmELspNKgruLk51KQ0RERJQ/VgXIMTExAAAPDw+jbUSGTA3QkyTJTqUhIiIiyh+rAuR27dpZtY0IMLVICLtXEBERUenBQXpU6JLT9VuQfThAj4iIiEoRmwPkDh06oG/fvlanHzhwIDp27Gjr6agUYQsyERERlWY2Ry67du1CcHCw1en379+Pa9eu2Xo6KkWMA2S2IBMREVHpUWxdLLRaLQdqOQjDad7YgkxERESlSbEEyDk5OUhISICnp2dxnI7szLAFmX2QiYiIqDSxumkvOTkZiYmJettycnJw/fp1s3MiCyGQmJiIlStXIjMzE/Xq1StQYal0SMlkCzIRERGVXlZHLp9++ineeecdvW337t1D1apVrTpekiQMHjw4X4Wj0ik5nYP0iIiIqPTKV+SibCmWJMnq1fRCQ0Pxyiuv4NVXX81f6ahUMuyDzC4WREREVJpYHSBPmDABw4YNA5AbKFerVg3ly5fHwYMHzR6jUqng4+MDX1/fAheUSg/OYkFERESlmdUBsq+vr16g27ZtW5QrVw5VqlQpkoJR6ZXMeZCJiIioFCvQPMhEpnCaNyIiIirNCjTNW3JyMh49epRnukePHiE5Obkgp6JSIitbi8xsrd42drEgIiKi0sTmAPmXX36Bv78/Ro8enWfaF198Ef7+/vjtt99sPR2VEoatxwDgwxZkIiIiKkVsDpA3bNgAABgxYkSeaUeNGgUhBNavX2/r6aiUMOx/DLAFmYiIiEoXmwPko0ePQqVSoXXr1nmm7dChA1QqFY4cOWLr6aiUMGxBdlZJcHMpthXNiYiIiArM5u++b968CT8/P7i5ueWZ1t3dHX5+frh586atp6MSIDEtC6v+uYJrD9LMprn3KEvvdx93F0iSVNRFIyIiIio0NgfIkiQhLc18oGQoPT2dgVIpN+aHI4j9736+juEMFkRERFTa2Pzdd6VKlZCRkYETJ07kmfbYsWNIT09HaGioracjO0vNzM53cAwAAZ6uRVAaIiIioqJjc4AcGRkJIQRmz56dZ9o5c+ZAkiS0b9/e1tORnaVmGQ++s8az9UMKuSRERERERcvm779fe+01LFu2DL/++itefPFFfPLJJwgKCtJLc+fOHUycOBG//vornJyc8Prrrxe4wGQfmRqt0bYRbcLh7GS624xKktCwkh+eqhVkcj8RERFRSWVzgBwREYF58+Zh+vTp+PHHH7Fx40Y0btxYXnr66tWrOHToELKzc1se33vvPdSqVatwSk3FLjM7x2jbjG414aRiv3IiIiIqWwo0gmrq1Knw8fHBtGnTkJKSgtjYWOzfvx8AIIQAAPj4+ODDDz+0akERKrkyDFqQXZwkBsdERERUJhV4ioExY8Zg4MCB2LhxI/755x/Ex8dDkiQEBwejVatW6Nu3L3x8fAqjrGRHhi3IamcnO5WEiIiIqGgVyhxcfn5+GDlyJEaOHFkY2VEJZNgHmYt/EBERUVnFKIesksEWZCIiInIQDJDJKoYtyGq2IBMREVEZVeAuFpcuXcL69etx/PhxPHjwABqNxmxaSZIQHR1d0FOSHbAFmYiIiBxFgQLkuXPn4r333oNWq5VnrbCES02XXuyDTERERI7C5gB5zZo1mDt3LgAgJCQEXbp0QUhICJydC2XcH5UwGRrDFmQGyERERFQ22RzNLlmyBADw7LPPYv369XB1dS20QlHJk5lt2ILMLhZERERUNtncDHjy5ElIkoQvv/ySwbEDMFwohC3IREREVFbZHOVIkgQfHx+EhIQUZnmohOJCIUREROQobA6QIyIikJaWhszMzMIsD5VQhi3IHKRHREREZZXNUc7IkSOh0WiwYcOGwiwPlVBsQSYiIiJHYXOAPGrUKDz77LN4/fXX8ffffxdmmagEMh6kxxZkIiIiKptsnsXinXfeQf369bFnzx60b98erVu3RvPmzeHt7W3xuFmzZtl6SrIj42ne2IJMREREZZPNAfKcOXPkhT+EENi7dy/27duX53EMkEsntiATERGRo7A5QG7bti1XxnMgbEEmIiIiR2FzgLxr165CLAaVdGxBJiIiIkfBKIeskskWZCIiInIQDJDJKoYtyGq2IBMREVEZxSiHrMI+yEREROQobO6D3KFDh3wfI0kSoqOjbT0l2RFbkImIiMhRFPkgPeVUcJz1ovQybEF2YwsyERERlVE2B8izZ8+2uD8pKQkHDhxAbGwsAgMDMWbMGDg5MagqrdiCTERERI6iyAJknZ07d+K5557D6dOnsXHjRltPR3ZmNM0bW5CJiIiojCryZsAOHTrg888/x6ZNm/Dtt98W9emoCGhytMjRCr1tbEEmIiKisqpYopz+/fvDycmJAXIpZdh6DABuLmxBJiIiorKpWAJkNzc3eHp64syZM8VxOipkhgP0AEDtzBZkIiIiKpuKJcq5efMmkpKSIITIO3EBzZs3D61atYKHhwf8/PysOkYIgVmzZqFixYpwd3dHp06dcOHChaItaCnCFmQiIiJyJEUeIKenp2Ps2LEAgLp16xb16ZCVlYW+fftizJgxVh/z4YcfYtGiRVi6dCkOHDgAT09PdOnSBRkZGUVY0tKDLchERETkSGyexeKdd96xuD8jIwPXr1/H9u3bcf/+fUiShHHjxtl6OqvNnTsXALBq1Sqr0gsh8Nlnn2HmzJno2bMnAOC7775DUFAQNm/ejAEDBhRVUUuNTI1+C7KTSoKLEwNkIiIiKptsDpDnzJlj1cIfQgioVCrMnDkTL7zwgq2nKzKXL19GfHw8OnXqJG/z9fVF8+bNERsbazZAzszMRGZmpvx7cnIyAECj0UCj0RRtoYvZo4xMvd/VzqoS+Rp1ZSqJZXMUrAP7Yx3YF6+//bEO7K8k14G1ZbI5QG7btq3FANnZ2Rn+/v6oX78++vXrhxo1ath6qiIVHx8PAAgKCtLbHhQUJO8zZf78+XJrtVJMTAw8PDwKt5B2diFJAvC/PseSNht//PGH/QqUh6ioKHsXweGxDuyPdWBfvP72xzqwv5JYB2lpaValK/KlpgvDtGnTsGDBAotpzpw5g4iIiGIqETB9+nRMmjRJ/j05ORmVKlVC+/btERgYWGzlKA67z98FTh+Vf/f2cEO3bu3sWCLTNBoNoqKi8NRTT8HFxcXexXFIrAP7Yx3YF6+//bEO7K8k14HuG/+82BwgF6c33ngDw4YNs5imWrVqNuUdHBwMALhz5w4qVqwob79z5w4aNGhg9ji1Wg21Wm203cXFpcTdDAWVLfT7G7u5OJXo11gW66C0YR3YH+vAvnj97Y91YH8lsQ6sLY/VAbJKpULFihVx8+ZNo31nzpyBRqNBvXr1rC9hPpQvXx7ly5cvkrzDw8MRHByM6OhoOSBOTk7GgQMH8jUTRlmWma0/iwWneCMiIqKyLF9TEZibx7hDhw5o1KhRoRSooK5du4a4uDhcu3YNOTk5iIuLQ1xcHB49eiSniYiIwKZNmwAAkiRhwoQJeO+99/Dbb7/hxIkTGDJkCEJCQtCrVy87vYqSxXAWC07xRkRERGVZoXWxKI5FQKwxa9YsrF69Wv69YcOGAHIHz0VGRgIAzp07h6SkJDnNlClTkJqaitGjRyMxMRFt2rTBtm3b4ObmVqxlL6kMW5DVbEEmIiKiMqxU9EHOj1WrVuU5B7JhMC9JEt55550853Z2VBlsQSYiIiIHwkiH8sQ+yERERORIGCBTntiCTERERI6EkQ7liS3IRERE5EgYIFOe2IJMREREjiRfg/Tu3LkDJyfzrYeW9gG5g+Gys7Pzc0oqAYxmsXBmCzIRERGVXfkKkEvKVG5UvAxbkN1c2IJMREREZZfVAfLs2bOLshxUgrEFmYiIiBwJA2TKU2Y2W5CJiIjIcTDSoTxlaAxbkHnbEBERUdnFSIfyZNyCzC4WREREVHYxQKY8GU3zxi4WREREVIYx0qE8GS0UwkF6REREVIYxQKY8ZbIFmYiIiBwIIx3KE1uQiYiIyJEwQKY8sQ8yERERORJGOpQnLhRCREREjoQBMlmUoxXQ5OgvMc6FQoiIiKgsY6RDFhm2HgNsQSYiIqKyjQEyWWQ4gwXAPshERERUtjHSIYsy2IJMREREDoYBMllkqgWZfZCJiIioLGOkQxYZtiBLEuDqxNuGiIiIyi5GOmSR0Sp6zipIkmSn0hAREREVPQbIZFGGhnMgExERkWNhgEwWZWbrtyCz/zERERGVdYx2yCK2IBMREZGjYYBMFhm2IKudecsQERFR2cZohywybEF2c2ELMhEREZVtDJDJIrYgExERkaNhtEMWGQ/SYwsyERERlW0MkMki40F6vGWIiIiobGO0QxaxBZmIiIgcDQNksiiTLchERETkYBjtkEVGg/TYgkxERERlHANksoh9kImIiMjRMNohi9gHmYiIiBwNA2SyiC3IRERE5GgY7ZBFxn2QecsQERFR2cZohywyWmramV0siIiIqGxjgEwWsQWZiIiIHA2jHbIoM5styERERORYGCCTRRkatiATERGRY2G0QxaxBZmIiIgcDQNksogtyERERORoGO2QRZmGs1hwoRAiIiIq4xggk0UZhrNYcKEQIiIiKuMY7ZBZQghkcalpIiIicjAMkMkswzmQAbYgExERUdnHaIfMytSYCpDZgkxERERlGwNkMivDYIo3AHDjLBZERERUxjHaIbPYgkxERESOiAEymWW4SAjAPshERERU9jHaIbMMFwlxdVJBpZLsVBoiIiKi4sEAmcwybEHmKnpERETkCBjxkFlGy0yz/zERERE5AAbIZJZhCzJnsCAiIiJHwIiHzDJuQebtQkRERGUfIx4yy7gFmV0siIiIqOxjgExmsQWZiIiIHBEjHjLLaBYLDtIjIiIiB8AAmcwybEHmID0iIiJyBIx4yCy2IBMREZEjYoBMZmVmswWZiIiIHA8jHjIrQ8MWZCIiInI8DJDJLLYgExERkSNixENmGbUgcx5kIiIicgBlLkCeN28eWrVqBQ8PD/j5+Vl1zLBhwyBJkt5P165di7agpYBRCzLnQSYiIiIH4GzvAhS2rKws9O3bFy1btsTy5cutPq5r165YuXKl/LtarS6K4pUqmWxBJiIiIgdU5gLkuXPnAgBWrVqVr+PUajWCg4OLoESll2ELMlfSIyIiIkdQ5gJkW+3atQsVKlSAv78/OnTogPfeew+BgYFm02dmZiIzM1P+PTk5GQCg0Wig0WiKvLzFIT0rW+93ZxVK9GvTla0kl7GsYx3YH+vAvnj97Y91YH8luQ6sLZMkhBBFXBa7WLVqFSZMmIDExMQ80/7000/w8PBAeHg4Ll26hBkzZsDLywuxsbFwcjLdrWDOnDlya7XS2rVr4eHhUdDilwgfH3fC9VRJ/n3gYzloUaFM3i5ERETkANLS0vDCCy8gKSkJPj4+ZtOVigB52rRpWLBggcU0Z86cQUREhPx7fgJkQ//99x8ee+wx7NixAx07djSZxlQLcqVKlXD79m2LLc+lSbfF+3AhIVX+/dO+dfFMvYp2LJFlGo0GUVFReOqpp+Di4mLv4jgk1oH9sQ7si9ff/lgH9leS6yA5ORnlypXLM0AuFV0s3njjDQwbNsximmrVqhXa+apVq4Zy5crh4sWLZgNktVptciCfi4tLibsZbJWVo//ZycPNtVS8trJUB6UV68D+WAf2xetvf6wD+yuJdWBteUpFgFy+fHmUL1++2M5348YN3L9/HxUrltzW0uKQqTFcKISzWBAREVHZV+amJbh27Rri4uJw7do15OTkIC4uDnFxcXj06JGcJiIiAps2bQIAPHr0CG+++Sb279+PK1euIDo6Gj179kT16tXRpUsXe72MEiEj23Cp6TJ3uxAREREZKRUtyPkxa9YsrF69Wv69YcOGAICYmBhERkYCAM6dO4ekpCQAgJOTE44fP47Vq1cjMTERISEh6Ny5M959912HnwuZLchERETkiMpcgLxq1ao850BWjkt0d3fH9u3bi7hUpY8Qgi3IRERE5JAY8ZBJmhwBw/lN2IJMREREjoABMplk2HoMsAWZiIiIHAMjHjLJsP8xwBZkIiIicgwMkMmkDA1bkImIiMgxMeIhkzKzjVuQGSATERGRI2DEQyYZtiA7qyQ4O/F2ISIiorKPEQ+ZZNiCzNZjIiIichSMesikTINZLDhAj4iIiBwFA2QyyXAWC7YgExERkaNg1EMmsQWZiIiIHBUDZDIpw6AF2ZUtyEREROQgGPWQSWxBJiIiIkfFAJlMMmxBZh9kIiIichSMesgktiATERGRo2KATCaxBZmIiIgcFaMeMsmwBVnNFmQiIiJyEAyQySTDFmQ3tiATERGRg2DUQyYZtyDzViEiIiLHwKiHTDJuQWYXCyIiInIMDJDJpMxsg0F6bEEmIiIiB8Goh0zK1BhM88YWZCIiInIQDJDJpAy2IBMREZGDYtRDJhm1IHOaNyIiInIQDJDJJKMWZE7zRkRERA6CUQ+ZxBZkIiIiclTO9i4A2d/g5Qdw4PIDvW1ZbEEmIiIiB8UAmZCdI4wCYkNqzmJBREREDoLNgmSVUH93exeBiIiIqFgwQKY8PdcwFDUqeNm7GERERETFgl0sCPN610FqZo7JfYFergjxY+sxEREROQ4GyIRq5dk6TERERKTDLhZERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUmCATERERESkwACZiIiIiEiBATIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhERERGRAgNkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESmUqQD5ypUrGDFiBMLDw+Hu7o7HHnsMs2fPRlZWlsXjMjIyMG7cOAQGBsLLywt9+vTBnTt3iqnURERERFSSlKkA+ezZs9BqtVi2bBlOnTqFTz/9FEuXLsWMGTMsHjdx4kRs2bIFGzZswO7du3Hr1i0899xzxVRqIiIiIipJnO1dgMLUtWtXdO3aVf69WrVqOHfuHL766it8/PHHJo9JSkrC8uXLsXbtWnTo0AEAsHLlStSsWRP79+9HixYtiqXsRERERFQylKkA2ZSkpCQEBASY3X/48GFoNBp06tRJ3hYREYHKlSsjNjbWbICcmZmJzMxMvfMAwIMHDwqp5JRfGo0GaWlpuH//PlxcXOxdHIfEOrA/1oF98frbH+vA/kpyHaSkpAAAhBAW05XpAPnixYtYvHix2dZjAIiPj4erqyv8/Pz0tgcFBSE+Pt7scfPnz8fcuXONtj/++OM2l5eIiIiIil5KSgp8fX3N7i8VAfK0adOwYMECi2nOnDmDiIgI+febN2+ia9eu6Nu3L0aNGlXoZZo+fTomTZok/56YmIgqVarg2rVrFi84FZ3k5GRUqlQJ169fh4+Pj72L45BYB/bHOrAvXn/7Yx3YX0muAyEEUlJSEBISYjFdqQiQ33jjDQwbNsximmrVqsn/v3XrFtq3b49WrVrh66+/tnhccHAwsrKykJiYqNeKfOfOHQQHB5s9Tq1WQ61WG2339fUtcTeDo/Hx8WEd2BnrwP5YB/bF629/rAP7K6l1YE1DZqkIkMuXL4/y5ctblfbmzZto3749GjdujJUrV0KlsjxRR+PGjeHi4oLo6Gj06dMHAHDu3Dlcu3YNLVu2LHDZiYiIiKh0KVPTvN28eRORkZGoXLkyPv74Y9y9exfx8fF6fYlv3ryJiIgIHDx4EEDup4gRI0Zg0qRJiImJweHDhzF8+HC0bNmSM1gQEREROaBS0YJsraioKFy8eBEXL15EWFiY3j7daEWNRoNz584hLS1N3vfpp59CpVKhT58+yMzMRJcuXfDll1/m69xqtRqzZ8822e2CigfrwP5YB/bHOrAvXn/7Yx3YX1moA0nkNc8FEREREZEDKVNdLIiIiIiICooBMhERERGRAgNkIiIiIiIFBshERERERAoMkAvBkiVLULVqVbi5uaF58+byFHJU+ObPn4+mTZvC29sbFSpUQK9evXDu3Dm9NBkZGRg3bhwCAwPh5eWFPn364M6dO3Yqcdn3wQcfQJIkTJgwQd7GOih6N2/exIsvvojAwEC4u7ujbt26OHTokLxfCIFZs2ahYsWKcHd3R6dOnXDhwgU7lrhsycnJwdtvv43w8HC4u7vjsccew7vvvgvluHfWQeH6+++/0aNHD4SEhECSJGzevFlvvzXX+8GDBxg0aBB8fHzg5+eHESNG4NGjR8X4Kko3S3Wg0WgwdepU1K1bF56enggJCcGQIUNw69YtvTxKSx0wQC6gdevWYdKkSZg9ezaOHDmC+vXro0uXLkhISLB30cqk3bt3Y9y4cdi/fz+ioqKg0WjQuXNnpKamymkmTpyILVu2YMOGDdi9ezdu3bqF5557zo6lLrv+/fdfLFu2DPXq1dPbzjooWg8fPkTr1q3h4uKCP//8E6dPn8Ynn3wCf39/Oc2HH36IRYsWYenSpThw4AA8PT3RpUsXZGRk2LHkZceCBQvw1Vdf4YsvvsCZM2ewYMECfPjhh1i8eLGchnVQuFJTU1G/fn0sWbLE5H5rrvegQYNw6tQpREVFYevWrfj7778xevTo4noJpZ6lOkhLS8ORI0fw9ttv48iRI/jll19w7tw5PPvss3rpSk0dCCqQZs2aiXHjxsm/5+TkiJCQEDF//nw7lspxJCQkCABi9+7dQgghEhMThYuLi9iwYYOc5syZMwKAiI2NtVcxy6SUlBRRo0YNERUVJdq1ayfGjx8vhGAdFIepU6eKNm3amN2v1WpFcHCw+Oijj+RtiYmJQq1Wix9//LE4iljmde/eXbz00kt625577jkxaNAgIQTroKgBEJs2bZJ/t+Z6nz59WgAQ//77r5zmzz//FJIkiZs3bxZb2csKwzow5eDBgwKAuHr1qhCidNUBW5ALICsrC4cPH0anTp3kbSqVCp06dUJsbKwdS+Y4kpKSAAABAQEAgMOHD0Oj0ejVSUREBCpXrsw6KWTjxo1D9+7d9a41wDooDr/99huaNGmCvn37okKFCmjYsCG++eYbef/ly5cRHx+vVwe+vr5o3rw566CQtGrVCtHR0Th//jwA4NixY9i7dy+efvppAKyD4mbN9Y6NjYWfnx+aNGkip+nUqRNUKhUOHDhQ7GV2BElJSZAkCX5+fgBKVx2UqZX0itu9e/eQk5ODoKAgve1BQUE4e/asnUrlOLRaLSZMmIDWrVujTp06AID4+Hi4urrKb0adoKAgvSXHqWB++uknHDlyBP/++6/RPtZB0fvvv//w1VdfYdKkSZgxYwb+/fdfvP7663B1dcXQoUPl62zq2cQ6KBzTpk1DcnIyIiIi4OTkhJycHMybNw+DBg0CANZBMbPmesfHx6NChQp6+52dnREQEMA6KQIZGRmYOnUqBg4cCB8fHwClqw4YIFOpNW7cOJw8eRJ79+61d1EcyvXr1zF+/HhERUXBzc3N3sVxSFqtFk2aNMH7778PAGjYsCFOnjyJpUuXYujQoXYunWNYv3491qxZg7Vr16J27dqIi4vDhAkTEBISwjogh6fRaNCvXz8IIfDVV1/Zuzg2YReLAihXrhycnJyMRuffuXMHwcHBdiqVY3j11VexdetWxMTEICwsTN4eHByMrKwsJCYm6qVnnRSew4cPIyEhAY0aNYKzszOcnZ2xe/duLFq0CM7OzggKCmIdFLGKFSuiVq1aettq1qyJa9euAYB8nflsKjpvvvkmpk2bhgEDBqBu3boYPHgwJk6ciPnz5wNgHRQ3a653cHCw0QD67OxsPHjwgHVSiHTB8dWrVxEVFSW3HgOlqw4YIBeAq6srGjdujOjoaHmbVqtFdHQ0WrZsaceSlV1CCLz66qvYtGkTdu7cifDwcL39jRs3houLi16dnDt3DteuXWOdFJKOHTvixIkTiIuLk3+aNGmCQYMGyf9nHRSt1q1bG01veP78eVSpUgUAEB4ejuDgYL06SE5OxoEDB1gHhSQtLQ0qlf6fUCcnJ2i1WgCsg+JmzfVu2bIlEhMTcfjwYTnNzp07odVq0bx582Ivc1mkC44vXLiAHTt2IDAwUG9/qaoDe48SLO1++uknoVarxapVq8Tp06fF6NGjhZ+fn4iPj7d30cqkMWPGCF9fX7Fr1y5x+/Zt+SctLU1O88orr4jKlSuLnTt3ikOHDomWLVuKli1b2rHUZZ9yFgshWAdF7eDBg8LZ2VnMmzdPXLhwQaxZs0Z4eHiIH374QU7zwQcfCD8/P/Hrr7+K48ePi549e4rw8HCRnp5ux5KXHUOHDhWhoaFi69at4vLly+KXX34R5cqVE1OmTJHTsA4KV0pKijh69Kg4evSoACAWLlwojh49Ks+QYM317tq1q2jYsKE4cOCA2Lt3r6hRo4YYOHCgvV5SqWOpDrKyssSzzz4rwsLCRFxcnN7f6MzMTDmP0lIHDJALweLFi0XlypWFq6uraNasmdi/f7+9i1RmATD5s3LlSjlNenq6GDt2rPD39xceHh6id+/e4vbt2/YrtAMwDJBZB0Vvy5Ytok6dOkKtVouIiAjx9ddf6+3XarXi7bffFkFBQUKtVouOHTuKc+fO2am0ZU9ycrIYP368qFy5snBzcxPVqlUTb731ll4gwDooXDExMSaf/0OHDhVCWHe979+/LwYOHCi8vLyEj4+PGD58uEhJSbHDqymdLNXB5cuXzf6NjomJkfMoLXUgCaFY9oeIiIiIyMGxDzIRERERkQIDZCIiIiIiBQbIREREREQKDJCJiIiIiBQYIBMRERERKTBAJiIiIiJSYIBMRERERKTAAJmIiIiISIEBMhGVSZGRkZAkCXPmzLF3UewqLS0Nb7/9NmrWrAl3d3dIkgRJkhAXF2fvohWZOXPmQJIkREZG2rsoNhk2bBgkScKwYcPsXRQih8UAmciB6AIHSZLg4eGBW7dumU175coVOe2uXbuKr5BUqPr374/33nsPZ8+ehSRJCAoKQlBQEFxcXOxdNIeza9cuzJkzB6tWrbJ3UYgoDwyQiRxUeno65s6da+9iUBE6e/Ystm7dCgBYt24d0tLSEB8fj/j4eNSuXdvOpXM8u3btwty5c/MMkCtWrIgnnngCFStWLJ6CEZERBshEDmzFihU4f/68vYtBReTEiRMAgMDAQPTr18/OpSFrzZ8/H2fPnsX8+fPtXRQih8UAmcgBVapUCfXq1UN2djZmzJhh7+JQEUlLSwMAeHl52bkkRESlCwNkIgekUqnk1qmff/4ZBw8ezNfxyv7JV65cMZuuatWqkCTJ6Ctlw+OvXr2KUaNGoXLlynBzc8Njjz2GmTNnIjU1VT7m5MmTePHFF1GpUiW4ubmhRo0aeO+996DRaPIsb1ZWFj744APUq1cPnp6e8Pf3x1NPPYU///wzz2NPnjyJ0aNHo0aNGvDw8ICXlxfq1auHt956C/fu3TN5jOEgsZ9//hmdO3dGhQoVoFKp8j1wMCMjA5999hlatWoFf39/uLm5oUqVKhgyZIjJwXa68+sGeV29elW+3rYO/tq3bx9efPFFVKlSBW5ubvD19UWzZs2wYMECPHr0SC+tRqNBuXLlIEkSFi1aZDHfFStWQJIk+Pj4yAE9AMTHx2Px4sXo2bMnatasCV9fX7i7u6N69eoYOXIkTp06le/XAFg3eNPSIL+HDx9i+fLl6NevH+rWrYuAgAC5Pl544QXs37/f6Bjd/a7r0rR79269+jB8j1gzSG/Xrl3o27cvQkNDoVarUa5cOXTs2BErV65ETk6OVa8rOjoa3bt3R/ny5eHm5oaaNWti7ty5yMjIMHve7du347nnnkNYWBhcXV3h4+ODatWqoXPnzvj444/x4MEDs8cSlSqCiBzG7NmzBQBRpUoVIYQQ7dq1EwBE+/btjdJevnxZABAARExMjNl9ly9fNnu+KlWqCABi5cqVZo//+eefhZ+fnwAgfHx8hJOTk7zvySefFFlZWWLr1q3Cw8NDABC+vr5CkiQ5Tf/+/U2eW/fapk+fLp588kkBQDg7O8vn0v3Mnj3bbPkXLFggVCqVnNbDw0O4urrKv1esWFEcOXLE7HVu166dmDRpkgAgJEkS/v7+wsnJyeI5Dd24cUPUqVNHPqeLi4vw9fWVf1epVGLRokV6x3z00UciKChI+Pj4yGmCgoLkn9dff93q8+fk5IjXX39d75p5eXnp1dMTTzwhrly5onfcuHHjBADRpEkTi/lHRkYKAGLYsGF624cOHSrn7+zsLAICAoSzs7O8Ta1Wi40bN5rMU3n9DenuC0t1YOl43T4AwsnJSfj7+wu1Wi1vkyRJfP7553rHXLt2TQQFBQlPT0+5DpX1ERQUJH766Sej1z506FCT5Zs4caLe+fz8/PTqo0OHDiI5Odni6/rwww+FJEny8cr3VPv27UV2drbR8XPnztW7Dzw8PISXl5feNsNnBVFpxQCZyIEYBsixsbHyH7Y///xTL21xBch+fn6iY8eO4tSpU0IIIdLS0sSiRYvkP/gzZ84Uvr6+on///nIQlpKSIt566y05j6ioKKNz6wIhX19foVarxdKlS0V6eroQIjdgef755+Xjf/31V6Pjv/32WzkYnDdvnrh9+7YQQojs7Gxx6NAh0aFDBwFAhIWFiZSUFJPXWRc8TJ06VSQkJAghhMjIyDAKJs3Jzs4WzZs3l1/HDz/8IDIzM4UQQly6dEk888wzcpD0xx9/GB2/cuVKvfq2xcyZMwUAUaFCBbFkyRJx//59IYQQWVlZIiYmRjRs2FAAEI0aNRI5OTnycQcOHJCv75kzZ0zmffXqVTkw27lzp96+d999V3z00UfixIkTQqPRCCFyg/WTJ0+KQYMGCQDC09NT3Lx50yjfogyQly1bJmbPni0OHTok14VWqxX//fefGD9+vJAkSTg5OeX5wckSSwHy4sWL5es6evRo+b589OiR+PTTT+UPEaY+OOrO7+fnJ1QqlZg+fbq4e/euEEKIpKQkMWvWLDnv5cuX6x175coV+cPipEmT9K57YmKi2LNnjxg7dqw4dOiQxddGVFowQCZyIIYBshBC9O7dWwAQDRo0EFqtVt5eXAFy7dq1RUZGhtGxgwcPltM89dRTemXT0bUMjxgxwmifLhAy9cdeiNxgq23btnIZlJKTk+WW5m3btpl8bRqNRjRu3FgAEJ9++qnePmUr46RJk0web42ffvpJzmf79u0my6ALoOvUqWO0v6AB8uXLl4WTk5Nwd3cXcXFxJtMkJyeLsLAwAUBs2rRJb98TTzwht+Kb8v777wsAonLlyibr15Lu3bsLAOLdd9812leUAXJedC3npu7JggbIaWlpIiAgQAAQAwcONHnsokWL5HvGMFhV3pfmXv9zzz0nAIhOnTrpbV+3bp0AIB5//HGLZScqK9gHmcjBvf/++3ByckJcXBx+/PHHYj//xIkToVarjbZ36dJF/v+0adMgSZLZNMePHzebf6VKlTB8+HCj7SqVCjNnzgQAnDp1Sp7xAcjtM5yYmIiGDRvqlUPJ2dkZAwcOBJDbL9MUlUqFqVOnmi1bXtatWwcAaNmyJTp37myyDLNnzwaQ21da+RoKw6pVq5CTk4OuXbuifv36JtN4e3ujV69eAIyvw+DBgwEAa9asgRDC6Njvv/8eADBo0CCT9WtJ9+7dAQB79+7N13FFrSjLFRUVJffxNdeHeuzYsfL0cGvXrjWZRq1WY/LkySb39ezZE4Dxe8rPzw8AkJKSojc2gKisYoBM5OAiIiLkAPLtt9+2atBbYWrWrJnJ7UFBQfL/mzZtajHNw4cPzeavG5RlypNPPglnZ2cAwKFDh+Tt+/btAwCcOXMGwcHBZn/eeecdALmD4EypXr06KlSoYLZsedGVqVOnTmbTtG/fHk5OTkavoTDorsNff/1l8TqsXLkSgPF1GDx4MCRJwrVr17B79269fYcPH8aZM2cAAEOGDDF5/mPHjmHs2LGoV68efHx8oFKp5EFtY8eOBQDcuHGjUF+zNf777z9MnjwZjRs3hp+fH5ycnORydevWrcjKpavfSpUq4fHHHzeZxsnJCR06dNBLb6h27dpmZzYJCQkBAKPBds2aNUO5cuVw+/ZtNG/eHF988QXOnj1r8oMPUVngbO8CEJH9zZkzB2vWrMF///2HpUuX4rXXXiu2c3t7e5vcrgtcrUljKagPDQ01u8/NzQ2BgYG4c+cOEhIS5O26FQYzMjIsjujXUc6+oFSQ4BiAXKa8XkO5cuWMXkNh0F2H1NRUq1oNDa9D5cqV0a5dO+zatQvff/+93qwQutbjpk2bIiIiwiivL774AuPHj4dWqwUASJIEX19f+duG9PR0JCcnF3tr5qZNmzBw4EBkZmbK23x8fODm5gZJkpCVlYWHDx8WSbmsuR8AICwsTC+9IXPvJ+B/76ns7Gy97X5+fvjxxx/xwgsv4NSpU/IzwtfXF23btkW/fv3Qv39/rtBIZQZbkIkIoaGh8h+89957z2jaLkejmyarf//+ELljNSz+mJvqTteyW1rprsPUqVOtug6mliTXtQ5v3LgR6enpAHKDL113Hl03DKUzZ85gwoQJ0Gq16Nu3Lw4ePIiMjAw8fPhQXglw4cKFAFCsLZj379/HsGHDkJmZiQ4dOmDXrl1IS0tDUlIS7ty5g/j4eGzYsKHYylPcOnXqhMuXL+O7777D0KFDUaNGDSQlJWHLli0YPHgwGjZsiJs3b9q7mESFggEyEQHI7efr7++PhIQEfPLJJxbTKlt3LbWwJiUlFVr5bGXpD3ZmZibu378PQL+1Nzg4GID5rhPFRVcmS1/XZ2RkmHwNhaEwrsPzzz8Pd3d3JCcn49dffwWQ22UjISEBLi4ucj9upY0bNyInJwc1a9bETz/9hKZNm8LV1VUvTXx8vE3l0d27tty3f/zxB5KTk+Hv748tW7agXbt2cHd3L5RyWcOa+0G5v7DvBwDw9PTE4MGDsWrVKpw/fx43btzAggUL4ObmpteyTFTaMUAmIgCAv78/pk2bBgD45JNPcPfuXYtpda5fv24yzfnz55GYmFioZbTF7t27zbYy7tmzR/4quUmTJvL21q1bA8jtJ3v79u2iL6QZujJFR0ebTbNr1y75NZjrq20r3XXYsWOHVV1NTFEO4tN1q9D9+/TTT6NcuXJGx+juqfr160OlMv1naseOHTaVR3fvmrtvAeDAgQMmt+uOeeKJJ+Dh4ZHvculei62t3rr74caNG2aXiM/JyUFMTAyAwr8fTAkNDcWUKVPwxhtvAMgdSEhUFjBAJiLZa6+9hrCwMKSkpODdd981m87T0xOPPfYYgNwZH0yZN29ekZQxv65du4bVq1cbbddqtXj//fcBALVq1ULdunXlfX379oWfnx80Gg0mTZpkMaDRarVF9kFgwIABAIDY2Fj89ddfRvuzs7PlgYJ16tRBnTp1CvX8L730EpydnXHv3j15tgxzsrKyzHbN0XWz+Ouvv3DhwgW5Jdnc4DxfX18AwIkTJ0xe+z///NNkdw5r6Gbj2L59u8l+wjt37kRsbKzFcp0/f97kB4a4uDizM0cAuX2VAdh8vzz11FMIDAwEYH4Wi2XLlsl9x021zttK2efaFF1LurkPNESlDe9kIpK5u7vLf3i3bNliMa3uj++KFSvw5Zdfyv1Lr1+/jpEjR2LdunVmW9mKk6+vL8aMGYNvvvlGDmquX7+OgQMHyi1t7733nt4xfn5++OyzzwAAP/30E7p3744DBw7IA8a0Wi3OnDmDTz75BLVr18bWrVuLpOx9+vRB8+bNAQD9+vXD2rVr5QGJly9fRp8+feRg7sMPPyz08z/22GN4++235fyHDBmCkydPyvuzs7MRFxeHd955B9WrVze57DWQG9gFBwcjOzsbL7zwAtLT0+Hv749nnnnGZPquXbsCyJ1+b9y4cfKMCqmpqVi2bBmef/55OVDMr379+kGlUuH+/fsYOHCg3B0hPT0dq1evRu/evREQEGDy2M6dO0OlUuHBgwcYNGiQ3H0nKysL69evR+fOnS0OgNN9gDl16hT++eeffJdd+f788ccf8corr+DOnTsAcgdILlq0CBMmTACQ23++cePG+T6HOQsWLMDTTz+N77//Xq+LR2ZmJtavX4+PPvoIwP+muSMq9YptxmUisjtTC4UYys7OFhEREXkuH5uSkiJq1aolp1GpVPLiGi4uLuLHH3+0aqEQcwuNxMTEyGnMsbQQhnKp6TZt2sjl8vf313ttM2fONJv/V199pbe0tFqtFoGBgcLFxUUvjx9++EHvuIIsNGHoxo0bonbt2vK5XF1d9ZbLVqlURksb6xTGSnparVa8/fbbeksRu7u7i8DAQL3ljQGIvXv3ms1Ht+S27ufll1+2eN4BAwbopVcup9y4cWN5RTlTry2v669cMQ7/v0qhbgW6Xr16yasHmjp+6tSpRsfq7ofw8HCxZs0as/etRqORF08BIPz9/UWVKlVElSpVxIYNG+R0+V1q2t/fX28Z7vbt2+e51LQ55t53ykVGdPdAQECA3n1Rs2ZNeWU/otKOLchEpMfJyUnuemCJl5cX9u7di0mTJiE8PBzOzs5wcXGRWzV13QPszdXVFdHR0Xj//ffxxBNPIDMzE76+vujYsSN+//13i11JXnnlFZw7dw6TJ09G/fr1oVarkZiYCC8vLzRp0gSvvfYaoqKiCvWrbEOhoaE4dOgQFi5ciBYtWsDd3R1paWmoVKkSBg8ejMOHD+P1118vsvNLkoR33nkHx48fx9ixY1GzZk04OTkhKSkJ/v7+aNWqFd588038888/cp9lUwy7U5jrXqGzZs0afPbZZ6hXrx7UajVycnJQt25dzJ8/H/v27TM7j6815s6di++//x4tWrSAp6cncnJy0KBBAyxduhS//PKLxdlHPvjgA3z33Xdo1qwZ3N3dodFoUL16dcyYMQNHjx6V5xE2xdnZGdHR0Rg5ciTCw8ORmpqKq1ev4urVq/maOWbhwoXYuXMn+vTpg6CgIDx69Aje3t5o3749VqxYgaioKIst2bYYPXo0vv76awwcOBB16tSBh4eHPGDxySefxGeffYYjR47IAzuJSjtJCM7yTURERESkwxZkIiIiIiIFBshERERERAoMkImIiIiIFBggExEREREpMEAmIiIiIlJggExEREREpMAAmYiIiIhIgQEyEREREZECA2QiIiIiIgUGyERERERECgyQiYiIiIgUGCATERERESkwQCYiIiIiUvg/iP+GCB/S2E0AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from matplotlib import rc\n", + "\n", + "def ensure_tensor(x, tkwargs, dim=0):\n", + " \"\"\"Helper function to concatenate or convert to tensor to handle potential multiple restarts.\"\"\"\n", + " if isinstance(x, (list, tuple)):\n", + " # flatten if there's only one element\n", + " if len(x) == 1:\n", + " x = x[0]\n", + " else:\n", + " x = torch.cat([xi if torch.is_tensor(xi) else torch.as_tensor(xi) for xi in x], dim=dim)\n", + " # ensure final type is tensor\n", + " if not torch.is_tensor(x):\n", + " x = torch.as_tensor(x)\n", + " return x.to(**tkwargs)\n", + "\n", + "Y_all = ensure_tensor(Y_all, tkwargs)\n", + "C_all = ensure_tensor(C_all, tkwargs)\n", + "\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 6))\n", + "\n", + "score = Y_all.clone()\n", + "# Set infeasible to -inf\n", + "score[~(C_all <= 0).all(dim=-1)] = float(\"-inf\")\n", + "fx = np.maximum.accumulate(score.cpu())\n", + "plt.plot(fx, marker=\"\", lw=3)\n", + "\n", + "plt.plot([0, len(Y_all)], [obj.fun.optimal_value, obj.fun.optimal_value], \"k--\", lw=3)\n", + "plt.ylabel(\"Function value\", fontsize=18)\n", + "plt.xlabel(\"Number of evaluations\", fontsize=18)\n", + "plt.title(\"10D Ackley with 2 outcome constraints\", fontsize=20)\n", + "plt.xlim([0, len(Y_all)])\n", + "plt.ylim([np.floor(min(fx[(C_all <= 0).all(dim=-1)])), 1]) \n", + "\n", + "plt.grid(True)\n", + "plt.show()" + ] + } + ], + "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.10.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}