From 7ea2c9fd2a61d17b5f650ec1a6b7cfc53ae57c44 Mon Sep 17 00:00:00 2001 From: dcluo Date: Sun, 20 Aug 2023 19:48:34 -0500 Subject: [PATCH] added chance constraint tutorial --- tutorials/3_Poisson_ChanceConstrained.ipynb | 516 ++++++++++++++++++++ 1 file changed, 516 insertions(+) create mode 100644 tutorials/3_Poisson_ChanceConstrained.ipynb diff --git a/tutorials/3_Poisson_ChanceConstrained.ipynb b/tutorials/3_Poisson_ChanceConstrained.ipynb new file mode 100644 index 0000000..0c67c9a --- /dev/null +++ b/tutorials/3_Poisson_ChanceConstrained.ipynb @@ -0,0 +1,516 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# Tutorial 2: Minimizing the superquantile/conditional-value-at-risk\n", + "In this tutorial, we will:\n", + "- Setup and solve a chance constrained optimization problem using SOUPy and `scipy.optimize`\n", + "\n", + "To run this with MPI for parallel sampling, export this notebook as a python script.\n", + "\n", + "## Problem definition\n", + "\n", + "We use again use the Poisson equation with a log-normal coefficient field with a distributed source as the control. Since CVaR may be sample-intensive, \n", + "to reduce the notebook's runtime, we will pose the problem in the 1D domain $\\Omega = (0,1)$. Recall the PDE is given by,\n", + "\n", + "\\begin{align*}\n", + "\\nabla \\cdot (e^{m} \\nabla u ) + z &= 0 \\qquad x \\in \\Omega, \\\\\n", + "u &= 0 \\qquad \\text{at } x = 0, \\\\\n", + "u &= 1 \\qquad \\text{at } x = 1. \\\\\n", + "\\end{align*}\n", + "\n", + "Here $u$ is the PDE solution, $m$, the uncertain parameter, is the log-coefficient field, and $z$, the optimization variable, is a distributed source.\n", + "We will again use $R(u,m,z)$ to denote the residual of the PDE. The corresponding weak form is \n", + "$$\n", + "\\text{Find } u \\in \\mathcal{U} \\text{ s.t. } \n", + "\\int_{\\Omega} e^m \\nabla u \\cdot \\nabla v dx - \\int_{\\Omega} z v dx = 0 \\qquad \\forall v \\in \\mathcal{V}\n", + "$$\n", + "with the trial and test spaces \n", + "\n", + "\\begin{align*}\n", + "\\mathcal{U} := \\{u \\in H^1(\\Omega) : u = x \\text{ on } \\Gamma_{D} \\}, \\\\\n", + "\\mathcal{V} := \\{v \\in H^1(\\Omega) : v = 0 \\text{ on } \\Gamma_{D} \\}.\n", + "\\end{align*}\n", + "\n", + "As with tutorial 1, we choose the log-coefficient $m$ to be distributed as a Matern Gaussian random field, $m \\sim \\mathcal{N}(\\bar{m}, \\mathcal{C})$, where the covariance operator $\\mathcal{C} = \\mathcal{A}^{-2}$ is given by the squared-inverse of an elliptic operator \n", + "$$ \\mathcal{A} = -\\gamma \\Delta + \\delta \\mathcal{I} $$\n", + "with Homogeneous Neumann boundary conditions. We will assume $(\\gamma, \\delta) = (0.5, 10)$ and $\\bar{m} = -3$ is a constant.\n", + "\n", + "The goal of our optimization problem is to reach a target source profile $z = 1$ while minimally disturbing the state from $u = x$. \n", + "We can formulate this optimization problem as a chance-constrained optimization problem, \n", + "\n", + "$$ \\min_{z} P(z) \\qquad \\text{s.t. PDE constraint and } \\mathbb{P}[Q > \\tau] < p_{\\text{tol}} $$ \n", + "\n", + "using the objective function\n", + "\n", + "$$ P(z) = \\int_{\\Omega}(z-1)^2 dx, $$ \n", + "\n", + "and the QoI\n", + "\n", + "$$ Q = \\int_{\\Omega} (u - x)^2 dx. $$ \n", + "\n", + "Here, $\\tau$ is some threshold value, $p_{\\text{tol}}$ is the maximum allows probability of $Q$ exceeding the threshold.\n", + "\n", + "\n", + "Recall that the probability can be formulated in terms of an expectation over an indicator function \n", + "\n", + "$$ \\mathbb{P}[Q > \\tau] = \\mathbb{E}[\\mathbb{1}_{(0, \\infty)}(Q - \\tau)], $$\n", + "\n", + "which, similar to the maximum function from the CVaR tutorial, is unfortunately not differentiable. We can also make a smooth approximation here, e.g. \n", + "\n", + "$$ \\mathbb{1}_{\\epsilon}(t) = \\frac{1}{1 + \\exp(-2x/\\epsilon)}. $$\n", + "\n", + "where $\\epsilon > 0$ controls the approximation error. We can then formulate the chance-constrained optimization using the SAA for the probability of exceedance, \n", + "\n", + "$$ \\min_{z} P(z) \\qquad \\text{s.t. PDE constraint and } \\frac{1}{N}\\sum_{i=1}^{N} \\mathbb{1}_{\\epsilon}(Q_i - \\tau) < p_{\\mathrm{tol}} $$ " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Import libraries \n", + "Note: hippylib and soupy paths need to be appended if cloning the repos instead of installing via pip" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import sys \n", + "import os \n", + "sys.path.append(os.environ.get('HIPPYLIB_PATH')) # Needed if using cloned repo\n", + "sys.path.append(os.environ.get('SOUPY_PATH')) # Needed if using cloned repo\n", + "import time\n", + "\n", + "import logging \n", + "logging.getLogger('FFC').setLevel(logging.WARNING)\n", + "logging.getLogger('UFL').setLevel(logging.WARNING)\n", + "\n", + "import scipy.optimize \n", + "import numpy as np \n", + "import matplotlib.pyplot as plt \n", + "import dolfin as dl \n", + "import hippylib as hp \n", + "from mpi4py import MPI \n", + "\n", + "import soupy " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Setup the function space\n", + "We will set up the mesh in 1D. Note that we will implement the code to be amenable for parallel sampling (see tutorial 1a). \n", + "To do so, we give the mesh the `MPI.COMM_SELF` communicator and save the `MPI.COMM_WORLD` communicator for sampling." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "N_ELEMENTS_X = 16\n", + "\n", + "comm_mesh = MPI.COMM_SELF\n", + "comm_sampler = MPI.COMM_WORLD\n", + "\n", + "mesh = dl.UnitIntervalMesh(comm_mesh, N_ELEMENTS_X) # Using MPI.COMM_SELF for the mesh so it does not partition\n", + "\n", + "Vh_STATE = dl.FunctionSpace(mesh, \"CG\", 1)\n", + "Vh_PARAMETER = dl.FunctionSpace(mesh, \"CG\", 1)\n", + "Vh_CONTROL = dl.FunctionSpace(mesh, \"CG\", 1)\n", + "Vh = [Vh_STATE, Vh_PARAMETER, Vh_STATE, Vh_CONTROL] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Define the PDE problem and Prior\n", + "This is the same as in tutorial 1, except for the change in the definition of the boundaries. " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Define PDE \n", + "\n", + "def residual(u,m,v,z):\n", + " return dl.exp(m)*dl.inner(dl.grad(u), dl.grad(v))*dl.dx - z * v *dl.dx \n", + "\n", + "def boundary(x, on_boundary):\n", + " return on_boundary and (dl.near(x[0], 0) or dl.near(x[0], 1))\n", + "\n", + "boundary_value = dl.Expression(\"x[0]\", degree=1, mpi_comm=comm_mesh) # Need to use the same mpi_comm as the mesh\n", + "\n", + "bc = dl.DirichletBC(Vh_STATE, boundary_value, boundary)\n", + "bc0 = dl.DirichletBC(Vh_STATE, dl.Constant(0.0), boundary)\n", + "pde = soupy.PDEVariationalControlProblem(Vh, residual, [bc], [bc0], is_fwd_linear=True)\n", + "\n", + "# Define prior\n", + "PRIOR_GAMMA = 1.0\n", + "PRIOR_DELTA = 10.0\n", + "PRIOR_MEAN = -3.0\n", + "\n", + "mean_vector = dl.interpolate(dl.Constant(PRIOR_MEAN), Vh_PARAMETER).vector()\n", + "prior = hp.BiLaplacianPrior(Vh_PARAMETER, PRIOR_GAMMA, PRIOR_DELTA, mean=mean_vector, robin_bc=False)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Define the QoI and control model\n", + "Here we prescribe the new QoI in its variational form. Note that in this example, the QoI will be in the constraint. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Define QoI\n", + "target_expression = dl.Expression(\"x[0]\", degree=1, mpi_comm=comm_mesh)\n", + "\n", + "def l2_qoi_form(u,m,z):\n", + " return (u - target_expression)**2*dl.dx \n", + "\n", + "qoi = soupy.VariationalControlQoI(Vh, l2_qoi_form)\n", + "\n", + "# Combine into control model\n", + "control_model = soupy.ControlModel(pde, qoi)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 5. Define the probability of exceedance as a risk measure.\n", + "We will define the probability of exceedance using the `soupy.TransformedMeanRiskMeasureSAASettings` class, which supports a risk of the form\n", + "\n", + "$$ g\\left( \\mathbb{E}[f(q)] \\right).$$\n", + "\n", + "In the case of the chance constraint, we have $g(x) = x$ and $f(x) = \\mathbb{1}_{\\epsilon}(x - \\tau)$. These two functions can be supplied to the risk measure class in the parameter list, `soupy.transformedMeanRiskMeasureSAASettings`, as `inner_function` for $f$ and `outer_function` for $g$. \n", + "\n", + "These are to be provided using the `soupy.FunctionWrapper` class, which will use finite differences to compute the derivatives of the function if not supplied. \n", + "\n", + "Note that by default, the `inner_function` and `outer_functions` will be set to identity functions $f(x) = x$ and $g(x) = x$." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+jElEQVR4nO3de1yVZb7///fisEAUUFFRPKKZhzFToRLTtDRMzWpydu5dqTXqLzuMqdWUOTMdZnaklTkzeait5q7M7GCznbSS7yNPpZka2kEnbTRRgQxFQDBO6/r9gWvlkoUCAjfrXq/n47Eexc211vpcXBzeXvd1X7fDGGMEAABgE0FWFwAAAFCbCDcAAMBWCDcAAMBWCDcAAMBWCDcAAMBWCDcAAMBWCDcAAMBWCDcAAMBWCDcAAMBWCDcIWNu2bdOvf/1rdejQQWFhYYqNjVVSUpIeeughq0s7rz179ujJJ5/UDz/8UOFzQ4YMUa9eveqlDofDoSeffPK8bX744Qc5HA49//zztfrenTp10l133eX5eMOGDXI4HNqwYUOtvo9bRkaGnnzySe3atatOXt8Xd598PX7zm9/UWx2+rF27ttKxP3dsACuEWF0AYIU1a9bopptu0pAhQzRnzhy1adNGmZmZ2rFjh9566y298MILVpdYqT179uipp57SkCFD1KlTJ6vLaRD69eunrVu3qmfPnnXy+hkZGXrqqafUqVMn9enTp07eozLPPPOMrr32Wq9jMTEx9VrDudauXav58+f7DDjvv/++oqKi6r8o4CyEGwSkOXPmKD4+Xh9//LFCQn75MfjP//xPzZkzx8LKUBNRUVHq37+/1WVUW2FhoSIiIs7bpmvXrn7Vt759+1pdAsBpKQSm48ePq0WLFl7Bxi0oyPvHolOnTrrxxhv1wQcfqG/fvmrUqJF69OihDz74QJK0bNky9ejRQ40bN9aVV16pHTt2VHjN1atXKykpSREREYqMjNT111+vrVu3Vmj36aefaujQoYqMjFRERIQGDBigNWvWeD6/bNky/cd//Ick6dprr/Wcpli2bJnX62zfvl2DBg1SRESEOnfurGeffVYul8urTV5enh5++GHFx8fL6XSqbdu2mjZtmgoKCiq0mzx5smJiYtSkSRPdcMMN2rdv33m+uue3bNkyORwOrV+/Xvfee69atGihmJgY3XrrrcrIyPBqW1JSot///vdq3bq1IiIiNHDgQH3xxRcVXrOy01Lbtm3T6NGjFRMTo/DwcHXp0kXTpk3zfP7777/X3Xffra5duyoiIkJt27bV6NGj9fXXX3u99hVXXCFJuvvuuz1f87NnLaoyvk8++aQcDoe+/PJL/eY3v1GzZs3UpUuXGn4Vy1V2CmjIkCEaMmSIVx8cDodWrFihWbNmKS4uTlFRURo2bJi+++67Cs//6KOPNHToUEVHRysiIkI9evRQSkqKJOmuu+7S/PnzJcnrVJn7NKmvmtLT03XnnXeqVatWCgsLU48ePfTCCy94fU+efQpz7ty5io+PV5MmTZSUlKTPP//8or5OCDyEGwSkpKQkbdu2TVOnTtW2bdtUUlJy3va7d+/WzJkz9eijj2rVqlWKjo7WrbfeqieeeEKLFy/WM888o+XLlys3N1c33nijTp8+7Xnum2++qZtvvllRUVFasWKFlixZopycHA0ZMkSffvqpp93GjRt13XXXKTc3V0uWLNGKFSsUGRmp0aNHa+XKlZKkUaNG6ZlnnpEkzZ8/X1u3btXWrVs1atQoz+tkZWXpjjvu0J133qnVq1drxIgRmjlzpt544w1Pm8LCQg0ePFj/+7//q6lTp+rDDz/Uo48+qmXLlummm26SMUaSZIzRLbfcotdff10PPfSQ3n//ffXv318jRoy46DGYNGmSQkND9eabb2rOnDnasGGD7rzzTq82kydP1vPPP6/x48fr//7v/zRmzBjdeuutysnJueDrf/zxxxo0aJDS09M1d+5cffjhh/rDH/6gH3/80dMmIyNDMTExevbZZ/XRRx9p/vz5CgkJ0VVXXeX5o9+vXz+9+uqrkqQ//OEPnq/5pEmTJFV9fN1uvfVWXXLJJXrnnXe0aNGiC/bD5XKptLTU61FTjz/+uA4dOqTFixfrlVde0f79+zV69GiVlZV52ixZskQjR46Uy+XSokWL9M9//lNTp07VkSNHJEl//OMfPWt+3F+LrVu3qk2bNj7f86efftKAAQO0bt06/fnPf9bq1as1bNgwPfzww3rggQcqtJ8/f75SU1M1b948LV++XAUFBRo5cqRyc3Nr3G8EIAMEoOzsbDNw4EAjyUgyoaGhZsCAASYlJcXk5+d7te3YsaNp1KiROXLkiOfYrl27jCTTpk0bU1BQ4Dn+j3/8w0gyq1evNsYYU1ZWZuLi4sxll11mysrKPO3y8/NNq1atzIABAzzH+vfvb1q1auX1/qWlpaZXr16mXbt2xuVyGWOMeeedd4wks379+gr9Gjx4sJFktm3b5nW8Z8+eZvjw4Z6PU1JSTFBQkNm+fbtXu3fffddIMmvXrjXGGPPhhx8aSeavf/2rV7v//u//NpLME088UfGLe5aDBw8aSea5557zHHv11VeNJHPfffd5tZ0zZ46RZDIzM40xxuzdu9dIMtOnT/dqt3z5ciPJTJgwwXNs/fr1Fb4mXbp0MV26dDGnT58+b41nKy0tNcXFxaZr165e77t9+3Yjybz66qte7aszvk888YSRZP70pz9VqRZ3n3w99u/fb4wp/948++vgNnjwYDN48OAKrzVy5Eivdm+//baRZLZu3eqpOyoqygwcONDz/ebL/fffbyr783FuTY899pjP78l7773XOBwO89133xljfvleueyyy0xpaamn3RdffGEkmRUrVlRaD3AuZm4QkGJiYrR582Zt375dzz77rG6++Wbt27dPM2fO1GWXXabs7Gyv9n369FHbtm09H/fo0UNS+fT/2Wsm3McPHTokSfruu++UkZGhcePGeZ3uatKkicaMGaPPP/9chYWFKigo0LZt2/Sb3/xGTZo08bQLDg7WuHHjdOTIEZ+nD3xp3bq1rrzySq9jvXv39tQkSR988IF69eqlPn36eM0IDB8+3Ov0zvr16yVJd9xxh9fr3X777VWq5XxuuummCjVKv3ztKnvv2267zefpxLPt27dP//73vzVx4kSFh4dX2q60tFTPPPOMevbsKafTqZCQEDmdTu3fv1979+69YB+qOr5nGzNmzAVf92yzZ8/W9u3bvR7t27ev1mu4XehrvmXLFuXl5em+++6Tw+Go0Xuc65NPPlHPnj0rfE/eddddMsbok08+8To+atQoBQcHV1ojUBUsKEZAS0xMVGJioqTy9R2PPvqoXnzxRc2ZM8drYXHz5s29nud0Os97/Oeff5ZUvrZHks8p+7i4OLlcLuXk5MgYI2NMpe3Ofq0L8XUlTVhYmNepsh9//FHff/+9QkNDfb6GO9wdP35cISEhFV6zdevWVaqlOnWGhYVJkqdOd3/PfS9f9Zzrp59+kiS1a9fuvO1mzJih+fPn69FHH9XgwYPVrFkzBQUFadKkSV5fr8pUdXzPDsCVnb6pTOfOnT3foxfrQl/zqn7dquP48eM+r+qr7Pv6QjUCVUG4Ac4IDQ3VE088oRdffFHffPNNrbym+xd1ZmZmhc9lZGQoKChIzZo1kzFGQUFBlbaTpBYtWtRKTe7XatSokZYuXVrp56Xy+ktLS3X8+HGvPzpZWVm1Vktl3O+XlZXlNWvmrud8WrZsKUmedSKVeeONNzR+/HjPOia37OxsNW3atMo1Xmh8z1ZbMyKSFB4erqKiogrHs7Oza/T9UtWvW3XExMTU2/c14MZpKQQkX79sJXlORbj/VXmxunXrprZt2+rNN9/0LNKVpIKCAr333nueK2waN26sq666SqtWrfL6F6rL5dIbb7yhdu3a6dJLL5VUO/+SvfHGG/Xvf/9bMTExntmrsx/uf2m791dZvny51/PffPPNGr93Vbmv9jn3vd9+++0LLqq99NJL1aVLFy1dutTnH383h8Ph+Xq6rVmzRkePHvU6VtnXvKrjW1c6deqkr776yuvYvn37qnwK81wDBgxQdHS0Fi1a5NWfc1Xne3Do0KHas2ePvvzyS6/jr732mhwOR4U9fIDawMwNAtLw4cPVrl07jR49Wt27d5fL5dKuXbv0wgsvqEmTJnrwwQdr5X2CgoI0Z84c3XHHHbrxxht1zz33qKioSM8995xOnjypZ5991tM2JSVF119/va699lo9/PDDcjqdWrBggb755hutWLHC8y9+9w7Er7zyiiIjIxUeHq74+Phqbew2bdo0vffee7rmmms0ffp09e7dWy6XS+np6Vq3bp0eeughXXXVVUpOTtY111yj3//+9yooKFBiYqI+++wzvf7667Xy9TmfHj166M4779S8efMUGhqqYcOG6ZtvvtHzzz9fpU3i5s+fr9GjR6t///6aPn26OnTooPT0dH388ceewHTjjTdq2bJl6t69u3r37q2dO3fqueeeq3BapkuXLmrUqJGWL1+uHj16qEmTJoqLi1NcXFyVx7cujBs3Tnfeeafuu+8+jRkzRocOHdKcOXM8MzDV1aRJE73wwguaNGmShg0bpsmTJys2Nlbff/+9du/erZdeekmSdNlll0kqXw80YsQIBQcHq3fv3p7TsmebPn26XnvtNY0aNUpPP/20OnbsqDVr1mjBggW69957PaEdqFUWLmYGLLNy5Upz++23m65du5omTZqY0NBQ06FDBzNu3DizZ88er7YdO3Y0o0aNqvAaksz999/vdczX1UHGlF9FddVVV5nw8HDTuHFjM3ToUPPZZ59VeM3Nmzeb6667zjRu3Ng0atTI9O/f3/zzn/+s0G7evHkmPj7eBAcHe13FM3jwYPOrX/2qQvsJEyaYjh07eh07deqU+cMf/mC6detmnE6niY6ONpdddpmZPn26ycrK8rQ7efKk+e1vf2uaNm1qIiIizPXXX2/+9a9/XfTVUudeqeXriqeioiLz0EMPmVatWpnw8HDTv39/s3Xr1gpX5Ph6rjHGbN261YwYMcJER0ebsLAw06VLF6+roHJycszEiRNNq1atTEREhBk4cKDZvHlzhauNjDFmxYoVpnv37iY0NLRC36syvu6rpX766afzfs3O7dM777xTaRuXy2XmzJljOnfubMLDw01iYqL55JNPKr1a6tzXco/PuVeBrV271gwePNg0btzYREREmJ49e5rZs2d7Pl9UVGQmTZpkWrZsaRwOh5FkDh48aIzxfQXXoUOHzO23325iYmJMaGio6datm3nuuee8rjCr7GfHGFOl7zXgbA5jzjP3CAAA4GdYcwMAAGyFcAMAAGyFcAMAAGyFcAMAAGyFcAMAAGyFcAMAAGwl4Dbxc7lcysjIUGRkZK1ugw4AAOqOMUb5+fmKi4vzulGtLwEXbjIyMmp8R10AAGCtw4cPX/DmrgEXbiIjIyWVf3GqsoU7AACwXl5entq3b+/5O34+ARdu3KeioqKiCDcAAPiZqiwpYUExAACwFcINAACwFcINAACwFcINAACwFcINAACwFcINAACwFcINAACwFcINAACwFcINAACwFcINAACwFUvDzaZNmzR69GjFxcXJ4XDoH//4xwWfs3HjRiUkJCg8PFydO3fWokWL6r5QAADgNywNNwUFBbr88sv10ksvVan9wYMHNXLkSA0aNEhpaWl6/PHHNXXqVL333nt1XCkAAPAXlt44c8SIERoxYkSV2y9atEgdOnTQvHnzJEk9evTQjh079Pzzz2vMmDF1VCUA1C1jjErKjMpcRqUul1xGcrmMXMbIZco/L0lGkjGS0ZmPjTzHa/q+QF0IDnKoTXQjy97fr+4KvnXrViUnJ3sdGz58uJYsWaKSkhKFhoZWeE5RUZGKioo8H+fl5dV5nQAClzFGh44Xak9mnrJyf9ZPp4p0LK9IP50q0snCYhUWl+l0cZkKi0t1uqTME2oAO2kVGaYvZg2z7P39KtxkZWUpNjbW61hsbKxKS0uVnZ2tNm3aVHhOSkqKnnrqqfoqEUCAMcboy/QcbdqXrV2HT2r3kZM6WVhSq+/hcEgOSQ6H48x/pfL/k875D9AghIVae72SX4UbqfyH+2zuadVzj7vNnDlTM2bM8Hycl5en9u3b112BAAKCy2X0yb+OacGG7/Vl+kmvzzmDg9SjTaTaNY9QyyZhahUVppZNwtQswqmIsGBFOEMU4QxWo9BgOUOCFBLkUEhw+X+DzzyCHA4FOSr/3Qagcn4Vblq3bq2srCyvY8eOHVNISIhiYmJ8PicsLExhYWH1UR6AAFDmMvrn7gwt3PBvffdjviTJGRKkG37VWomdmqlP+6bq3jpKzhB22gCs4lfhJikpSf/85z+9jq1bt06JiYk+19sAQG0qLXPpdyvS9OE35f/IahIWojv6d9DEq+PVKirc4uoAuFkabk6dOqXvv//e8/HBgwe1a9cuNW/eXB06dNDMmTN19OhRvfbaa5KkKVOm6KWXXtKMGTM0efJkbd26VUuWLNGKFSus6gKAAFHmMnrond368JssOYOD9LvrLtH4AZ0U3Yh/WAENjaXhZseOHbr22ms9H7vXxkyYMEHLli1TZmam0tPTPZ+Pj4/X2rVrNX36dM2fP19xcXH629/+xmXgAOqUy2X0+Kqv9X+7MhQS5NCCO/ppWM/YCz8RgCUcJsA2OsjLy1N0dLRyc3MVFRVldTkAGjhjjJ5Y/a1e23pIQQ7p7//VT6N6V7wyE0Ddqs7fb1a8AcB5PPvhv/Ta1kNyOKTn/+Nygg3gBwg3AFCJnYdy9PKmA5KkZ359mW7t187iigBUBeEGACrxP2eCzX8ktNN/XdnB4moAVBXhBgB8OHS8QB/vKb/ke/I1nS2uBkB1EG4AwIelnx6UMdKQbi11aWyk1eUAqAbCDQCc42Rhsd7ecUSSNHkQszaAvyHcAMA5lm9L1+mSMvVoE6UBXXzf2gVAw0W4AYCzFJe69L9bfpAkTR4Uz40rAT9EuAGAs6zenaFj+UWKjQrTjb3jrC4HQA0QbgDgDGOMFm8uv/z7rgHx3Nkb8FP85ALAGZv3Z+tfWfmKcAbrdva1AfwW4QYAzljxRfmNem9LbK/oCO72Dfgrwg0AqPyU1PYfciRJN3L/KMCvEW4AQNKRnNPKPlWk0GCHerWNtrocABeBcAMAkr5ML5+16RkXrfDQYIurAXAxCDcAICkt/aQkqW/7ppbWAeDiEW4AQL/M3PTr2MziSgBcLMINgID3c0mZ9mTkSZL6dWhqbTEALhrhBkDA+/porkpdRi0jw9S2aSOrywFwkQg3AALel4fOnJLq0JR7SQE2QLgBEPDci4n7dWC9DWAHhBsAAc0Y41lM3JdwA9gC4QZAQMvI/VnH8osUEuRQ73Zs3gfYAeEGQEBzr7fpGRfF5n2ATRBuAAQ0zykpNu8DbINwAyCgeRYTs3kfYBuEGwAB6+eSMn2bkSuJK6UAOyHcAAhY32bkqqTMqEUTp9o1Y/M+wC4INwACludmmR2asXkfYCOEGwABy3OzTE5JAbZCuAEQsH6ZuWlqaR0AahfhBkBAysw9rczcnxXM5n2A7RBuAASkvZl5kqSurZoowhlicTUAahPhBkBAOppzWpLUvnmExZUAqG2EGwAB6cjJ8nDTtimXgAN2Q7gBEJDcMzfsbwPYD+EGQEA6yswNYFuEGwAByT1zE0e4AWyHcAMg4BSVlulYfpEkqS2npQDbIdwACDhZuT9LksJDgxTT2GlxNQBqG+EGQMA5+5QU95QC7IdwAyDgcBk4YG+EGwABh8vAAXsj3AAIOFwGDtgb4QZAwHHP3HClFGBPhBsAAeeXmRvuKwXYEeEGQEBxuYwyc5m5AeyMcAMgoBzLL1JJmVFwkEOxkWFWlwOgDhBuAASUoycLJUmto8IVEsyvQMCO+MkGEFCOnizfnZhTUoB9EW4ABBTPHjdcBg7YFuEGQEBxn5bibuCAfRFuAAQU9rgB7I9wAyCgsDsxYH+EGwABwxjDzA0QAAg3AAJG7ukSFRSXSWLmBrAzy8PNggULFB8fr/DwcCUkJGjz5s3nbb98+XJdfvnlioiIUJs2bXT33Xfr+PHj9VQtAH925MysTYsmToWHBltcDYC6Ymm4WblypaZNm6ZZs2YpLS1NgwYN0ogRI5Senu6z/aeffqrx48dr4sSJ+vbbb/XOO+9o+/btmjRpUj1XDsAfsd4GCAyWhpu5c+dq4sSJmjRpknr06KF58+apffv2Wrhwoc/2n3/+uTp16qSpU6cqPj5eAwcO1D333KMdO3bUc+UA/BHrbYDAYFm4KS4u1s6dO5WcnOx1PDk5WVu2bPH5nAEDBujIkSNau3atjDH68ccf9e6772rUqFGVvk9RUZHy8vK8HgACUwYzN0BAsCzcZGdnq6ysTLGxsV7HY2NjlZWV5fM5AwYM0PLlyzV27Fg5nU61bt1aTZs21d///vdK3yclJUXR0dGeR/v27Wu1HwD8B6elgMBg+YJih8Ph9bExpsIxtz179mjq1Kn605/+pJ07d+qjjz7SwYMHNWXKlEpff+bMmcrNzfU8Dh8+XKv1A/AfnnDTLMLiSgDUpRCr3rhFixYKDg6uMEtz7NixCrM5bikpKbr66qv1yCOPSJJ69+6txo0ba9CgQfrLX/6iNm3aVHhOWFiYwsLCar8DAPyOZ80NMzeArVk2c+N0OpWQkKDU1FSv46mpqRowYIDP5xQWFiooyLvk4ODyyzmNMXVTKABbOF1cpuMFxZJYUAzYnaWnpWbMmKHFixdr6dKl2rt3r6ZPn6709HTPaaaZM2dq/PjxnvajR4/WqlWrtHDhQh04cECfffaZpk6dqiuvvFJxcXFWdQOAH3CfkooMC1F0o1CLqwFQlyw7LSVJY8eO1fHjx/X0008rMzNTvXr10tq1a9WxY0dJUmZmpteeN3fddZfy8/P10ksv6aGHHlLTpk113XXXafbs2VZ1AYCfcIcb7gYO2J/DBNj5nLy8PEVHRys3N1dRUVFWlwOgnry5LV2Pv/+1ruveSkvvusLqcgBUU3X+flt+tRQA1IejJwslsZgYCASEGwABgd2JgcBBuAEQEDJO/iyJmRsgEBBuAASEn04VSZJio8ItrgRAXSPcAAgIx8+Em+aNnRZXAqCuEW4A2F5JmUt5P5dKItwAgYBwA8D2ThaWSJKCHGIDPyAAEG4A2N6JM7ddaBrhVHCQ7xvzArAPwg0A23OHm2YRzNoAgYBwA8D2cgrLw01M4zCLKwFQHwg3AGzPfTfwZo2ZuQECAeEGgO3lnAk3XCkFBAbCDQDb+2XNDeEGCASEGwC2515zw8wNEBgINwBs7wSnpYCAQrgBYHue01KEGyAgEG4A2J57QXEM4QYICIQbALZmjPnlUnAWFAMBgXADwNZOl5SpqNQliTU3QKAg3ACwNfd6G2dIkCKcwRZXA6A+EG4A2NqJs9bbOBzcNBMIBIQbALbGBn5A4CHcALA1NvADAg/hBoCtnSgokUS4AQIJ4QaArZ0oKJJEuAECCeEGgK25Z25YcwMEDsINAFvL8dxXKtTiSgDUF8INAFv75aaZYRZXAqC+EG4A2NqJQvdNM5m5AQIF4QaArf1yWoo1N0CgINwAsC2Xy7DPDRCACDcAbCv3dIlcpvz/uVoKCByEGwC25V5vExkeotBgft0BgYKfdgC2xXobIDARbgDY1nHCDRCQCDcAbMszc8N6GyCgEG4A2NYve9wQboBAQrgBYFsnTpWHmxjCDRBQCDcAbIuZGyAwEW4A2BZrboDARLgBYFvum2YycwMEFsINANs6wa0XgIBEuAFgWzkFJZIIN0CgIdwAsKWi0jKdKiqVxJobINAQbgDYknvWJjjIoahGIRZXA6A+EW4A2JJnMXGEUw6Hw+JqANQnwg0AW8rxLCYOtbgSAPWNcAPAlo6fNXMDILAQbgDYknsDv5gmhBsg0BBuANjSCWZugIBFuAFgS+5wwx43QOAh3ACwJXYnBgIX4QaALeUwcwMELMINAFtizQ0QuAg3AGyJNTdA4CLcALAdY8xZm/gRboBAY3m4WbBggeLj4xUeHq6EhARt3rz5vO2Lioo0a9YsdezYUWFhYerSpYuWLl1aT9UC8AenikpVUmYkcVoKCESW3k1u5cqVmjZtmhYsWKCrr75aL7/8skaMGKE9e/aoQ4cOPp9z22236ccff9SSJUt0ySWX6NixYyotLa3nygE0ZO5TUo1Cg9XIGWxxNQDqm6XhZu7cuZo4caImTZokSZo3b54+/vhjLVy4UCkpKRXaf/TRR9q4caMOHDig5s2bS5I6depUnyUD8AOstwECm2WnpYqLi7Vz504lJyd7HU9OTtaWLVt8Pmf16tVKTEzUnDlz1LZtW1166aV6+OGHdfr06Urfp6ioSHl5eV4PAPbGehsgsFk2c5Odna2ysjLFxsZ6HY+NjVVWVpbP5xw4cECffvqpwsPD9f777ys7O1v33XefTpw4Uem6m5SUFD311FO1Xj+AhutEQYkkqWkEdwQHApHlC4odDofXx8aYCsfcXC6XHA6Hli9friuvvFIjR47U3LlztWzZskpnb2bOnKnc3FzP4/Dhw7XeBwANy8kzMzdNWUwMBCTLZm5atGih4ODgCrM0x44dqzCb49amTRu1bdtW0dHRnmM9evSQMUZHjhxR165dKzwnLCxMYWFhtVs8gAYt7/SZmZtGzNwAgciymRun06mEhASlpqZ6HU9NTdWAAQN8Pufqq69WRkaGTp065Tm2b98+BQUFqV27dnVaLwD/kXsm3EQTboCAZOlpqRkzZmjx4sVaunSp9u7dq+nTpys9PV1TpkyRVH5Kafz48Z72t99+u2JiYnT33Xdrz5492rRpkx555BH99re/VaNGjazqBoAG5uRp1twAgczSS8HHjh2r48eP6+mnn1ZmZqZ69eqltWvXqmPHjpKkzMxMpaene9o3adJEqamp+t3vfqfExETFxMTotttu01/+8herugCgAXLP3EQxcwMEJIcxxlhdRH3Ky8tTdHS0cnNzFRUVZXU5AOrALfM/067DJ/XyuAQN/1Vrq8sBUAuq8/fb8qulAKC2saAYCGyEGwC2415zE82aGyAgEW4A2IoxxrPmpmkj9rkBAhHhBoCtFBSXqcxVvpSQS8GBwES4AWAr7t2JnSFBCg/lVxwQiPjJB2ArZ2/gV9mtXADYG+EGgK3kFrI7MRDoCDcAbCWXy8CBgEe4AWAr3FcKAOEGgK2wxw0Awg0AW2HmBgDhBoCtnGRBMRDwCDcAbIX7SgEg3ACwlVzW3AABj3ADwFZOni7foZj7SgGBi3ADwFbcMzdRnJYCAhbhBoCtsKAYAOEGgG2UuYzyfy6VJDVlzQ0QsEKq+4TvvvtOK1as0ObNm/XDDz+osLBQLVu2VN++fTV8+HCNGTNGYWFhdVErAJxX/s8lnv9n5gYIXFWeuUlLS9P111+vyy+/XJs2bdIVV1yhadOm6c9//rPuvPNOGWM0a9YsxcXFafbs2SoqKqrLugGgAvcpqcbOYIUGMzENBKoqz9zccssteuSRR7Ry5Uo1b9680nZbt27Viy++qBdeeEGPP/54rRQJAFXB7sQApGqEm/3798vpvPCllUlJSUpKSlJxcfFFFQYA1XWSK6UAqBqnpaoSbCSpsLCwWu0BoLa4Z25YTAwEthqdlB4yZIiOHDlS4fi2bdvUp0+fi60JAGokt7B8xpjTUkBgq1G4iYqKUu/evfXWW29Jklwul5588kldc801uummm2q1QACoKs/MDbsTAwGt2peCS9Lq1au1aNEiTZo0SatXr9YPP/yg9PR0rVmzRsOGDavtGgGgSrivFACphuFGkqZMmaJDhw5p9uzZCgkJ0YYNGzRgwIDarA0AqoXdiQFINTwtlZOTozFjxmjhwoV6+eWXddtttyk5OVkLFiyo7foAoMq4FByAVMOZm169eik+Pl5paWmKj4/X5MmTtXLlSt13331as2aN1qxZU9t1AsAFnSTcAFANZ26mTJmiTZs2KT4+3nNs7Nix2r17N/vbALBMHpeCA1ANZ27++Mc/+jzerl07paamXlRBAFBTnJYCIFVj5iY9Pb1aL3z06NFqFwMAF8O9oJhLwYHAVuVwc8UVV2jy5Mn64osvKm2Tm5ur//mf/1GvXr20atWqWikQAKqiqLRMp0vKJDFzAwS6Kp+W2rt3r5555hndcMMNCg0NVWJiouLi4hQeHq6cnBzt2bNH3377rRITE/Xcc89pxIgRdVk3AHhxn5JyOKTI8BrvcgHABqo8c3PkyBHNnj1bGRkZWrRokS699FJlZ2dr//79kqQ77rhDO3fu1GeffUawAVDv3IuJo8JDFRTksLgaAFaq8j9v+vbtq6ysLLVs2VIPPfSQtm/frpiYmLqsDQCqjMXEANyqPHPTtGlTHThwQJL0ww8/yOVy1VlRAFBdnsXEXAYOBLwqz9yMGTNGgwcPVps2beRwOJSYmKjg4GCfbd0hCADqCzM3ANyqHG5eeeUV3Xrrrfr+++81depUTZ48WZGRkXVZGwBUmXvmJopwAwS8al1ScMMNN0iSdu7cqQcffJBwA6DBcM/cNCXcAAGvRtdLvvrqq7VdBwBcFE5LAXCr0b2lAKChyeW+UgDOINwAsAVmbgC4EW4A2MLJwmJJUjT3lQICHuEGgC0wcwPAjXADwBYINwDcCDcA/J4xhgXFADwINwD83umSMpWUGUnM3AAg3ACwAffuxKHBDkU4fd8WBkDgINwA8Htnr7dxOBwWVwPAaoQbAH6P+0oBOBvhBoDf475SAM5GuAHg9/K4DBzAWQg3APzeydPluxM3jWB3YgCEGwA2wAZ+AM5GuAHg91hQDOBsloebBQsWKD4+XuHh4UpISNDmzZur9LzPPvtMISEh6tOnT90WCKDBY0ExgLNZGm5WrlypadOmadasWUpLS9OgQYM0YsQIpaenn/d5ubm5Gj9+vIYOHVpPlQJoyDgtBeBsloabuXPnauLEiZo0aZJ69OihefPmqX379lq4cOF5n3fPPffo9ttvV1JSUj1VCqAh475SAM5mWbgpLi7Wzp07lZyc7HU8OTlZW7ZsqfR5r776qv7973/riSeeqOsSAfgJZm4AnC3EqjfOzs5WWVmZYmNjvY7HxsYqKyvL53P279+vxx57TJs3b1ZISNVKLyoqUlFRkefjvLy8mhcNoEFyLygm3ACQGsCC4nPvA2OM8XlvmLKyMt1+++166qmndOmll1b59VNSUhQdHe15tG/f/qJrBtBwlJa5PDM3zRuzzw0AC8NNixYtFBwcXGGW5tixYxVmcyQpPz9fO3bs0AMPPKCQkBCFhITo6aef1u7duxUSEqJPPvnE5/vMnDlTubm5nsfhw4frpD8ArJFzZtbG4WATPwDlLDst5XQ6lZCQoNTUVP3617/2HE9NTdXNN99coX1UVJS+/vprr2MLFizQJ598onfffVfx8fE+3ycsLExhYWG1WzyABiOn8MzuxI1CFRzEHcEBWBhuJGnGjBkaN26cEhMTlZSUpFdeeUXp6emaMmWKpPJZl6NHj+q1115TUFCQevXq5fX8Vq1aKTw8vMJxAIHjREF5uGnGKSkAZ1gabsaOHavjx4/r6aefVmZmpnr16qW1a9eqY8eOkqTMzMwL7nkDILDlnAk3zTklBeAMhzHGWF1EfcrLy1N0dLRyc3MVFRVldTkALtIbnx/SH/7xja7vGav/GZ9odTkA6kh1/n5bfrUUAFwM98xNDKelAJxBuAHg104UsuYGgDfCDQC/doI1NwDOQbgB4Ne4WgrAuQg3APyae58b1twAcCPcAPBrOQXlOxQzcwPAjXADwK8dLyi/MS5rbgC4EW4A+K3TxWX6ucQlSWrehHADoBzhBoDfcl8G7gwOUmNnsMXVAGgoCDcA/NaJU+4rpULlcHDTTADlCDcA/JZnAz/W2wA4C+EGgN/y3HqB9TYAzkK4AeC3PBv4MXMD4CyEGwB+y3PrBfa4AXAWwg0Av+Vec0O4AXA2wg0Av5XDzA0AHwg3APzWcdbcAPCBcAPAbzFzA8AXwg0Av5XDmhsAPhBuAPgll8sop7D8juCEGwBnI9wA8Et5P5eozGUkSU0jQi2uBkBDQrgB4Jfce9xEhoUoLISbZgL4BeEGgF9yr7dpxikpAOcg3ADwS8dPEW4A+Ea4AeCXPFdKsd4GwDkINwD80okC95VSYRZXAqChIdwA8Eu/7HHDzA0Ab4QbAH6JNTcAKkO4AeCX3DM3MYQbAOcg3ADwSye4aSaAShBuAPilE9w0E0AlCDcA/JL7juCsuQFwLsINAL9TXOpSflGpJNbcAKiIcAPA75w8s5g4yCFFhXMpOABvhBsAfuf4WYuJg4IcFlcDoKEh3ADwOzksJgZwHoQbAH7nBHcEB3AehBsAfsdzGTh73ADwgXADwO+c4DJwAOdBuAHgd9xrbrgMHIAvhBsAfudEYYkkZm4A+Ea4AeB3ThQUSZKaN2aPGwAVEW4A+J0TBeUzN80bh1lcCYCGiHADwO/kcLUUgPMg3ADwK8aYs66W4rQUgIoINwD8SkFxmYrLXJLYoRiAb4QbAH7FfUoqPDRIEc4Qi6sB0BARbgD4FXYnBnAhhBsAfoXdiQFcCOEGgF85wR3BAVwA4QaAX/HM3HBaCkAlCDcA/Epm7s+SpNgoNvAD4BvhBoBfOXqyUJLUtmkjiysB0FARbgD4laMnT0uS2jaLsLgSAA0V4QaAX8k4WX5aipkbAJUh3ADwG4XFpZ4FxW2bEW4A+GZ5uFmwYIHi4+MVHh6uhIQEbd68udK2q1at0vXXX6+WLVsqKipKSUlJ+vjjj+uxWgBWyjhzSioyLETRjbivFADfLA03K1eu1LRp0zRr1iylpaVp0KBBGjFihNLT032237Rpk66//nqtXbtWO3fu1LXXXqvRo0crLS2tnisHYIUjOe71NszaAKicwxhjrHrzq666Sv369dPChQs9x3r06KFbbrlFKSkpVXqNX/3qVxo7dqz+9Kc/Val9Xl6eoqOjlZubq6ioqBrVDcAay7cd0qz3v9HQ7q205K4rrC4HQD2qzt9vy2ZuiouLtXPnTiUnJ3sdT05O1pYtW6r0Gi6XS/n5+WrevHldlAiggTnKzA2AKrDslrrZ2dkqKytTbGys1/HY2FhlZWVV6TVeeOEFFRQU6Lbbbqu0TVFRkYqKijwf5+Xl1axgAJbzXAbOlVIAzsPyBcUOh8PrY2NMhWO+rFixQk8++aRWrlypVq1aVdouJSVF0dHRnkf79u0vumYA1mDmBkBVWBZuWrRooeDg4AqzNMeOHaswm3OulStXauLEiXr77bc1bNiw87adOXOmcnNzPY/Dhw9fdO0ArMHMDYCqsCzcOJ1OJSQkKDU11et4amqqBgwYUOnzVqxYobvuuktvvvmmRo0adcH3CQsLU1RUlNcDgP8pKXPpx7wzG/gxcwPgPCxbcyNJM2bM0Lhx45SYmKikpCS98sorSk9P15QpUySVz7ocPXpUr732mqTyYDN+/Hj99a9/Vf/+/T2zPo0aNVJ0dLRl/QBQ97Jyf5bLSM7gILVozE0zAVTO0nAzduxYHT9+XE8//bQyMzPVq1cvrV27Vh07dpQkZWZmeu158/LLL6u0tFT333+/7r//fs/xCRMmaNmyZfVdPoB65N7jJq5puIKCLrwuD0DgsnSfGyuwzw3gn97deUQPv7NbV18So+WT+ltdDoB65hf73ABAdXiulGIxMYALINwA8AsZniulIiyuBEBDR7gB4Bc8l4FzpRSACyDcAPAL7HEDoKoINwAaPJfLeMJNO2ZuAFwA4QZAg5ddUKTiUpeCHFLr6HCrywHQwBFuADR47iulYqPCFRrMry0A58dvCQANHuttAFQH4QZAg8fdwAFUB+EGQIPHzA2A6iDcAGjwmLkBUB2EGwANnnvmJo6ZGwBVQLgB0OC5Z27aEW4AVAHhBkCDlnu6RPlFpZI4LQWgagg3ABo09w0zm0WEKsIZYnE1APwB4QZAg8ZiYgDVRbgB0KBxGTiA6iLcAGjQfgk3ERZXAsBfEG4ANGiclgJQXYQbAA3aEU5LAagmwg2ABs2zxw0zNwCqiHADoME6VVSq7FNFkpi5AVB1hBsADdZXh09KKg82zRo7rS0GgN8g3ABosL5Mz5Ek9enQ1NpCAPgVwg2ABist/aQkqV+HZtYWAsCvEG4ANEjGGKWdOS3Vj5kbANVAuAHQIB06XqgTBcVyBgepZ1yU1eUA8COEGwANknu9Ta+2UQoLCba4GgD+hHADoEFyh5u+rLcBUE2EGwANEouJAdQU4QZAg1NYXKp/ZeVLkvp1bGptMQD8DuEGQIPz1ZFclbmMWkeFq000OxMDqB7CDYAGx73ehlkbADVBuAHQ4Hx56KQk1tsAqBnCDYAGxRijXYfdV0o1tbYYAH6JcAOgQTl84rSyTxUrNNihX8VFW10OAD9EuAHQoLjX2/SMi1Z4KJv3Aag+wg2ABiXNvZiYU1IAaohwA6BB+ZLN+wBcJMINgAbjdHGZ9mbmSWIxMYCaI9wAaDC+PpqrUpdRq8gwtW3K5n0AaoZwA6DB8Gze16GZHA6HxdUA8FeEGwANgjFGq748Iknq37m5xdUA8GeEGwANwsZ9P2nfj6fUJCxEtya0s7ocAH6McAOgQVi8+aAkaewV7RUVHmpxNQD8GeEGgOX2ZOTp0++zFRzk0N1Xd7K6HAB+jnADwHKLNx+QJI28rI3aNYuwuBoA/o5wA8BSWbk/a/XuDEnS5EHxFlcDwA4INwAstWzLDyp1GV0Z31y92zW1uhwANkC4AWCZgqJSvbntkCRp8qDOFlcDwC4INwAs8/aOw8r7uVSdWzTW0O6trC4HgE0QbgBY4ueSMi35tPzy798OjFdQEDsSA6gdhBsA9a6otEz/3+s7dSTntGIaOzWmH5v2Aag9hBsA9aqkzKUH3kzTpn0/qVFosBbemaBGzmCrywJgI4QbAPWmzGU0feUupe75Uc6QIC2ekKgr47mPFIDaRbgBUC9cLqNH3t2tD77KVGiwQy/fmaCrL2lhdVkAbCjE6gIA2JsxRlv+fVx//2S/Pj9wQsFBDv39v/rqWq6OAlBHLJ+5WbBggeLj4xUeHq6EhARt3rz5vO03btyohIQEhYeHq3Pnzlq0aFE9VQqgOlwuo4++ydItC7bojsXb9PmBEwoNdmjubZfrhl5trC4PgI1ZOnOzcuVKTZs2TQsWLNDVV1+tl19+WSNGjNCePXvUoUOHCu0PHjyokSNHavLkyXrjjTf02Wef6b777lPLli01ZswYC3oAwK2kzKXvsvK16/BJ7T58UtsOnlD6iUJJUlhIkP7zivaafE1n7h0FoM45jDHGqje/6qqr1K9fPy1cuNBzrEePHrrllluUkpJSof2jjz6q1atXa+/evZ5jU6ZM0e7du7V169YqvWdeXp6io6OVm5urqKioi+/EGWUuo8zc07X2ekB1ne8n2f05IyNjJKPy00WS5DLl/+8ykssYlbnKH6Uul0rKjErLjIrLylRYXP44XVymguJSnThVrGP5Rfopv0g/nSrS4ROFKip1eb1vZHiIxid11N1Xx6tFk7A66jmAQFCdv9+WzdwUFxdr586deuyxx7yOJycna8uWLT6fs3XrViUnJ3sdGz58uJYsWaKSkhKFhoZWeE5RUZGKioo8H+fl5dVC9RUdLyjSwNnr6+S1AX8RGR6iPu2bqk/7prq8XVNd1bm5IsMr/lwCQF2yLNxkZ2errKxMsbGxXsdjY2OVlZXl8zlZWVk+25eWlio7O1tt2lQ8j5+SkqKnnnqq9go/j7AQy5cwwcYcVdjA1yHvRmc/x+E5dqaVo/yYw+FQcJBDQY7y/w9ySCFBQQoJdigkyKHQ4CA5Q4LUKDRYjZzBinAGq1FoiGKaONUqMkwtzzziohupQ/MIdhoGYDnLr5ZynPMb2xhT4diF2vs67jZz5kzNmDHD83FeXp7at29f03Ir1SoyXN/9ZUStvy4AAKgey8JNixYtFBwcXGGW5tixYxVmZ9xat27ts31ISIhiYmJ8PicsLExhYZzrBwAgUFh2HsXpdCohIUGpqalex1NTUzVgwACfz0lKSqrQft26dUpMTPS53gYAAAQeSxeJzJgxQ4sXL9bSpUu1d+9eTZ8+Xenp6ZoyZYqk8lNK48eP97SfMmWKDh06pBkzZmjv3r1aunSplixZoocfftiqLgAAgAbG0jU3Y8eO1fHjx/X0008rMzNTvXr10tq1a9WxY0dJUmZmptLT0z3t4+PjtXbtWk2fPl3z589XXFyc/va3v7HHDQAA8LB0nxsr1NU+NwAAoO5U5+831y4DAABbIdwAAABbIdwAAABbIdwAAABbIdwAAABbIdwAAABbIdwAAABbIdwAAABbIdwAAABbsfT2C1Zwb8icl5dncSUAAKCq3H+3q3JjhYALN/n5+ZKk9u3bW1wJAACorvz8fEVHR5+3TcDdW8rlcikjI0ORkZFyOBy1+tp5eXlq3769Dh8+bMv7Vtm9f5L9+0j//J/d+0j//F9d9dEYo/z8fMXFxSko6PyragJu5iYoKEjt2rWr0/eIioqy7TetZP/+SfbvI/3zf3bvI/3zf3XRxwvN2LixoBgAANgK4QYAANgK4aYWhYWF6YknnlBYWJjVpdQJu/dPsn8f6Z//s3sf6Z//awh9DLgFxQAAwN6YuQEAALZCuAEAALZCuAEAALZCuAEAALZCuKmmBQsWKD4+XuHh4UpISNDmzZvP237jxo1KSEhQeHi4OnfurEWLFtVTpTVTnf5t2LBBDoejwuNf//pXPVZcdZs2bdLo0aMVFxcnh8Ohf/zjHxd8jr+NX3X76E9jmJKSoiuuuEKRkZFq1aqVbrnlFn333XcXfJ4/jWFN+uhPY7hw4UL17t3bs7lbUlKSPvzww/M+x5/Gr7r986ex8yUlJUUOh0PTpk07bzsrxpBwUw0rV67UtGnTNGvWLKWlpWnQoEEaMWKE0tPTfbY/ePCgRo4cqUGDBiktLU2PP/64pk6dqvfee6+eK6+a6vbP7bvvvlNmZqbn0bVr13qquHoKCgp0+eWX66WXXqpSe38bP6n6fXTzhzHcuHGj7r//fn3++edKTU1VaWmpkpOTVVBQUOlz/G0Ma9JHN38Yw3bt2unZZ5/Vjh07tGPHDl133XW6+eab9e233/ps72/jV93+ufnD2J1r+/bteuWVV9S7d+/ztrNsDA2q7MorrzRTpkzxOta9e3fz2GOP+Wz/+9//3nTv3t3r2D333GP69+9fZzVejOr2b/369UaSycnJqYfqapck8/7775+3jb+N37mq0kd/HsNjx44ZSWbjxo2VtvH3MaxKH/15DI0xplmzZmbx4sU+P+fv42fM+fvnr2OXn59vunbtalJTU83gwYPNgw8+WGlbq8aQmZsqKi4u1s6dO5WcnOx1PDk5WVu2bPH5nK1bt1ZoP3z4cO3YsUMlJSV1VmtN1KR/bn379lWbNm00dOhQrV+/vi7LrFf+NH4Xyx/HMDc3V5LUvHnzStv4+xhWpY9u/jaGZWVleuutt1RQUKCkpCSfbfx5/KrSPzd/G7v7779fo0aN0rBhwy7Y1qoxJNxUUXZ2tsrKyhQbG+t1PDY2VllZWT6fk5WV5bN9aWmpsrOz66zWmqhJ/9q0aaNXXnlF7733nlatWqVu3bpp6NCh2rRpU32UXOf8afxqyl/H0BijGTNmaODAgerVq1el7fx5DKvaR38bw6+//lpNmjRRWFiYpkyZovfff189e/b02dYfx686/fO3sZOkt956S19++aVSUlKq1N6qMQy4u4JfLIfD4fWxMabCsQu193W8oahO/7p166Zu3bp5Pk5KStLhw4f1/PPP65prrqnTOuuLv41fdfnrGD7wwAP66quv9Omnn16wrb+OYVX76G9j2K1bN+3atUsnT57Ue++9pwkTJmjjxo2VBgB/G7/q9M/fxu7w4cN68MEHtW7dOoWHh1f5eVaMITM3VdSiRQsFBwdXmMU4duxYhVTq1rp1a5/tQ0JCFBMTU2e11kRN+udL//79tX///touzxL+NH61qaGP4e9+9zutXr1a69evV7t27c7b1l/HsDp99KUhj6HT6dQll1yixMREpaSk6PLLL9df//pXn239cfyq0z9fGvLY7dy5U8eOHVNCQoJCQkIUEhKijRs36m9/+5tCQkJUVlZW4TlWjSHhpoqcTqcSEhKUmprqdTw1NVUDBgzw+ZykpKQK7detW6fExESFhobWWa01UZP++ZKWlqY2bdrUdnmW8Kfxq00NdQyNMXrggQe0atUqffLJJ4qPj7/gc/xtDGvSR18a6hj6YoxRUVGRz8/52/j5cr7++dKQx27o0KH6+uuvtWvXLs8jMTFRd9xxh3bt2qXg4OAKz7FsDOt0ubLNvPXWWyY0NNQsWbLE7Nmzx0ybNs00btzY/PDDD8YYYx577DEzbtw4T/sDBw6YiIgIM336dLNnzx6zZMkSExoaat59912runBe1e3fiy++aN5//32zb98+880335jHHnvMSDLvvfeeVV04r/z8fJOWlmbS0tKMJDN37lyTlpZmDh06ZIzx//Ezpvp99KcxvPfee010dLTZsGGDyczM9DwKCws9bfx9DGvSR38aw5kzZ5pNmzaZgwcPmq+++so8/vjjJigoyKxbt84Y4//jV93++dPYVebcq6UayhgSbqpp/vz5pmPHjsbpdJp+/fp5XaI5YcIEM3jwYK/2GzZsMH379jVOp9N06tTJLFy4sJ4rrp7q9G/27NmmS5cuJjw83DRr1swMHDjQrFmzxoKqq8Z92eW5jwkTJhhj7DF+1e2jP42hr35JMq+++qqnjb+PYU366E9j+Nvf/tbz+6Vly5Zm6NChnj/8xvj/+FW3f/40dpU5N9w0lDF0GHNmZQ8AAIANsOYGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGAADYCuEGgN/76aef1Lp1az3zzDOeY9u2bZPT6dS6dessrAyAFbi3FABbWLt2rW655RZt2bJF3bt3V9++fTVq1CjNmzfP6tIA1DPCDQDbuP/++/X//t//0xVXXKHdu3dr+/btCg8Pt7osAPWMcAPANk6fPq1evXrp8OHD2rFjh3r37m11SQAswJobALZx4MABZWRkyOVy6dChQ1aXA8AizNwAsIXi4mJdeeWV6tOnj7p37665c+fq66+/VmxsrNWlAahnhBsAtvDII4/o3Xff1e7du9WkSRNde+21ioyM1AcffGB1aQDqGaelAPi9DRs2aN68eXr99dcVFRWloKAgvf766/r000+1cOFCq8sDUM+YuQEAALbCzA0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALCV/x/ytNDNHAPLKwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "TAU = 2.0 # QoI threshold value \n", + "EPSILON = 1e-1\n", + "smoothed_indicator_form = lambda x : 1/ (1 + np.exp(-2 * (x - TAU)/EPSILON)) # Smoothed indicator function\n", + "smoothed_indicator = soupy.FunctionWrapper(smoothed_indicator_form) \n", + "\n", + "x_plot = np.linspace(0, 2*TAU, 100)\n", + "\n", + "if comm_sampler.rank == 0:\n", + " plt.figure()\n", + " plt.plot(x_plot, smoothed_indicator_form(x_plot))\n", + " plt.xlabel(\"x\")\n", + " plt.ylabel(\"f(x)\")\n", + " plt.title(\"Smoothed Indicator Function\")\n", + "\n", + "SAMPLE_SIZE = 400\n", + "RANDOM_SEED = 1\n", + "\n", + "risk_settings = soupy.transformedMeanRiskMeasureSAASettings()\n", + "risk_settings[\"sample_size\"] = SAMPLE_SIZE \n", + "risk_settings[\"seed\"] = RANDOM_SEED\n", + "risk_settings[\"inner_function\"] = smoothed_indicator # Outer function will be identity by default \n", + "\n", + "risk_measure = soupy.TransformedMeanRiskMeasureSAA(control_model, prior, risk_settings, comm_sampler=comm_sampler)\n", + "\n", + "chance_constraint = soupy.RiskMeasureControlCostFunctional(risk_measure, penalization=None)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 6. Define the penalization and cost functional\n", + "We now proceed to define the optimization objective as a penalization term since it does not depend on the PDE solution. \n", + "\n", + "This is then converted to a cost functional using the class `soupy.PenalizationControlCostFunctional`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "Z_TARGET = 1.0 \n", + "def l2_objective_form(z):\n", + " return (z- dl.Constant(Z_TARGET))**2*dl.dx \n", + "\n", + "penalty = soupy.VariationalPenalization(Vh, l2_objective_form)\n", + "\n", + "l2_objective = soupy.PenalizationControlCostFunctional(Vh, penalty)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 7. Optimization using SciPy\n", + "We will now use the SciPy wrapper to define the cost functional and chance constraints. We will use the SLSQP algorithm, which uses the first derivatives of the cost and the constraints.\n", + "\n", + "We will need to set the upper bound on the constraint function to be $p_{\\mathrm{tol}}$, which we take to be $0.05$ in this example. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Optimization terminated successfully (Exit mode 0)\n", + " Current function value: 0.09339550500373343\n", + " Iterations: 20\n", + " Function evaluations: 20\n", + " Gradient evaluations: 20\n" + ] + } + ], + "source": [ + "# Interface with scipy.optimize.minimize\n", + "scipy_cost = soupy.ScipyCostWrapper(l2_objective, verbose=False)\n", + "scipy_chance_constraint = soupy.ScipyCostWrapper(chance_constraint, verbose=False)\n", + "\n", + "P_TOL = 0.05\n", + "constraint = scipy.optimize.NonlinearConstraint(scipy_chance_constraint.function(), lb=-np.inf, ub=P_TOL, jac=scipy_chance_constraint.jac())\n", + "\n", + "z_initial_np = chance_constraint.generate_vector(soupy.CONTROL).get_local()\n", + "DISPLAY_ITERATIONS = True\n", + "MAX_ITER = 200\n", + "options = {'disp' : DISPLAY_ITERATIONS, 'maxiter' : MAX_ITER}\n", + "\n", + "# Use the SLSQP algorithm for constrained optimization\n", + "results = scipy.optimize.minimize(scipy_cost.function(), z_initial_np, method=\"SLSQP\", constraints=constraint,\n", + " jac=scipy_cost.jac(), options=options)\n", + "\n", + "z_optimal_np = results['x']\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 8. Postprocessing\n", + "We will now plot a sample of the parameter and solution at the optimal control.\n", + "\n", + "Since the result is in the form of a numpy array for the augmented vector, we need to set it to a dolfin vector for solving the PDE. \n", + "Remember to exclude the last component when setting it to a vector representing the control variable. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTLUlEQVR4nOzdd3hUVcIG8Hd6Mum9d0pC7x0pIh1ElwU/XSmiK2JDZFdZVMTGiooFBctSVsWGgq6CCiodpPfQk0B672Uy5Xx/3GSSIT2kTCbv73nyZObeM3fO5CaZd849RSaEECAiIiKiNk/e2hUgIiIioqbBYEdERERkIxjsiIiIiGwEgx0RERGRjWCwIyIiIrIRDHZERERENoLBjoiIiMhGMNgRERER2QgGOyIiIiIbwWBHRM3mzz//xF//+lf4+flBrVbD19cX06dPx6FDh27puGvWrMHGjRurbI+Li4NMJqt2X0sIDQ3FnDlzWuW565KUlIQXX3wRp06dapbjW/NrJ2pPGOyIqFmsXr0aQ4cORUJCAlauXInffvsNb775JhITEzFs2DC8//77jT52TcHOz88Phw4dwqRJk26h5rYpKSkJy5cvb7ZgR0TWQdnaFSAi23PgwAEsXLgQEydOxNatW6FUVvyrueeee3DXXXfhySefRO/evTF06NAme16NRoNBgwY12fHas6KiImi12tauBhE1EFvsiKjJrVixAjKZDGvXrrUIdQCgVCqxZs0ayGQy/Pvf/zZvf/HFFyGTyXDy5EncfffdcHZ2houLC/72t78hPT3dXC40NBTnz5/Hnj17IJPJIJPJEBoaCqD6S7Hlxz1z5gz++te/wsXFBe7u7li0aBEMBgMuXbqE8ePHw8nJCaGhoVi5cqVFfUtKSvD000+jV69e5scOHjwYP/zwQ6N/PiaTCatXr0avXr1gb28PV1dXDBo0CP/73/8syqxcuRKRkZHQaDTw9vbGrFmzkJCQYHGskSNHolu3bjh69CiGDx8OrVaL8PBw/Pvf/4bJZAIA7N69G/379wcAzJ071/xze/HFFwEAc+bMgaOjI86ePYuxY8fCyckJt99+OwAgKysLCxYsQEBAANRqNcLDw7F06VLodLpGv34iaj5ssSOiJmU0GrFr1y7069cPgYGB1ZYJCgpC37598ccff8BoNEKhUJj33XXXXZgxYwbmz5+P8+fP4/nnn0d0dDQOHz4MlUqFrVu3Yvr06XBxccGaNWsASC11dZkxYwb+9re/4eGHH8bOnTuxcuVK6PV6/Pbbb1iwYAEWL16ML774As888ww6dOiAu+++GwCg0+mQlZWFxYsXIyAgAKWlpfjtt99w9913Y8OGDZg1a1aDf0Zz5szB559/jnnz5uGll16CWq3GiRMnEBcXZy7zyCOP4OOPP8Zjjz2GyZMnIy4uDs8//zx2796NEydOwNPT01w2JSUF9913H55++mksW7YMW7duxZIlS+Dv749Zs2ahT58+2LBhA+bOnYvnnnvOfKm68vkpLS3F1KlT8fDDD+PZZ5+FwWBASUkJRo0ahWvXrmH58uXo0aMH9u3bhxUrVuDUqVPYtm1bg187ETUzQUTUhFJSUgQAcc8999RabubMmQKASE1NFUIIsWzZMgFAPPXUUxblNm3aJACIzz//3Lyta9euYsSIEVWOGRsbKwCIDRs2mLeVH/ett96yKNurVy8BQGzZssW8Ta/XCy8vL3H33XfXWG+DwSD0er2YN2+e6N27t8W+kJAQMXv27Fpf9969ewUAsXTp0hrLXLhwQQAQCxYssNh++PBhAUD861//Mm8bMWKEACAOHz5sUbZLly5i3Lhx5vtHjx6t8rMpN3v2bAFArF+/3mL7hx9+KACIb775xmL766+/LgCIHTt2mLfV57UTUfPjpVgiahVCCACATCaz2H7fffdZ3J8xYwaUSiV27dp1S883efJki/tRUVGQyWSYMGGCeZtSqUSHDh1w/fp1i7KbN2/G0KFD4ejoCKVSCZVKhXXr1uHChQsNrsfPP/8MAHj00UdrLFP+Wm8eZTpgwABERUXh999/t9ju6+uLAQMGWGzr0aNHlddRl7/85S8W9//44w84ODhg+vTpFtvL63VzPYio9THYEVGT8vT0hFarRWxsbK3l4uLioNVq4e7ubrHd19fX4r5SqYSHhwcyMzNvqV43P49arYZWq4WdnV2V7SUlJeb7W7ZswYwZMxAQEIDPP/8chw4dwtGjR/HAAw9YlKuv9PR0KBSKKq+zsvLX6ufnV2Wfv79/lZ+Fh4dHlXIajQbFxcX1rpdWq4Wzs3OVevj6+lYJ397e3lAqlbd8Toio6bGPHRE1KYVCgVGjRuGXX35BQkJCtf3sEhIScPz4cUyYMMGifx0g9RcLCAgw3zcYDMjMzKw2vLSEzz//HGFhYfj6668tAk5jBw94eXnBaDQiJSWl2uAGVAS15OTkKj+/pKQki/51TeXm8FZej8OHD0MIYbE/LS0NBoOhWepBRLeGLXZE1OSWLFkCIQQWLFgAo9Fosc9oNOKRRx6BEAJLliyp8thNmzZZ3P/mm29gMBgwcuRI87aGtkbdCplMBrVabRFsUlJSGj0qtvzS79q1a2ssM3r0aABSqKzs6NGjuHDhgnnEakOUDzBpyM/t9ttvR0FBAb7//nuL7Z9++ql5PxFZF7bYEVGTGzp0KN555x0sXLgQw4YNw2OPPYbg4GDcuHEDH3zwAQ4fPox33nkHQ4YMqfLYLVu2QKlU4o477jCPiu3ZsydmzJhhLtO9e3d89dVX+PrrrxEeHg47Ozt07969WV7L5MmTsWXLFixYsADTp09HfHw8Xn75Zfj5+eHKlSsNPt7w4cNx//3345VXXkFqaiomT54MjUaDkydPQqvV4vHHH0fnzp3x97//HatXr4ZcLseECRPMo2KDgoLw1FNPNfh5IyIiYG9vj02bNiEqKgqOjo7w9/eHv79/jY+ZNWsWPvjgA8yePRtxcXHo3r079u/fj9deew0TJ07EmDFjGlwPImpeDHZE1Cwef/xx9O/fH2+99RaefvppZGZmwt3dHcOGDcP+/fsxePDgah+3ZcsWvPjii1i7di1kMhmmTJmCd955B2q12lxm+fLlSE5OxkMPPYT8/HyEhIRYTBXSlObOnYu0tDR8+OGHWL9+PcLDw/Hss88iISEBy5cvb9QxN27ciD59+mDdunXYuHEj7O3t0aVLF/zrX/8yl1m7di0iIiKwbt06fPDBB3BxccH48eOxYsWKRl2W1mq1WL9+PZYvX46xY8dCr9dj2bJl5rnsqmNnZ4ddu3Zh6dKleOONN5Ceno6AgAAsXrwYy5Yta8xLJ6JmJhPlQ9OIiFrRiy++iOXLlyM9PZ19t4iIGol97IiIiIhsBIMdERERkY3gpVgiIiIiG8EWOyIiIiIbwWBHREREZCMY7IiIiIhsRLucx85kMiEpKQlOTk7VLqNDREREZC2EEMjPz4e/vz/k8jra5EQr27Nnj5g8ebLw8/MTAMTWrVvrfMzu3btFnz59hEajEWFhYWLt2rUNes74+HgBgF/84he/+MUvfvGrzXzFx8fXmXFavcWusLAQPXv2xNy5c/GXv/ylzvKxsbGYOHEiHnroIXz++ec4cOAAFixYAC8vr3o9HgCcnJwAAPHx8XB2dr6l+hMRERE1p7y8PAQFBZnzS21aPdhNmDDBvCh2fXz44YcIDg7GO++8AwCIiorCsWPH8Oabb9Y72JVffnV2dmawIyIiojahPt3HWj3YNdShQ4cwduxYi23jxo3DunXroNfroVKpWqlmlowmgXd/u4y/dfeEt7Nd9YUUCsCu0r7CwpoPKJcD9vaNK1tUBNQ0XaFMBmi1jStbXAyYTDXXw8GhcWVLSgCjsWnKarVSvQFApwMMhqYpa28v/ZwBoLQU0OubpqydnfR70dCyer1UviYaDaBUNryswSD9LGqiVgPlf3MNKWs0SueuJiqVVL6hZU0m6XetKcoqldLPApD+JoqKmqZsQ/7u+T+i+rL8H9HwsvwfId1u4v8ReqUKXx25gd7BbugW4FJz2ZbUoM5pzQyou49dx44dxauvvmqx7cCBAwKASEpKqvYxJSUlIjc31/xV3scuNze3qapexReHr4uQZ34SQvoXWP3XxImWD9Jqay47YoRlWU/Pmsv262dZNiSk5rJduliW7dKl5rIhIZZl+/Wruaynp2XZESNqLqvVWpadOLH2n1tl06fXXragoKLs7Nm1l01Lqyi7YEHtZWNjK8ouXlx72XPnKsouW1Z72SNHKsquXFl72V27Ksq+/37tZX/6qaLshg21l/3mm4qy33xTe9kNGyrK/lTH7/v771eU3bWr9rIrV1aUPXKk9rLLllWUPXeu9rKLF1eUjY2tveyCBRVl09JqLzt7dkXZgoLay06fLizUVpb/I6Qv/o+o+OL/COnLCv5HxM2cLUa+sUuEPPOTuPeTQ8JkMonmkpubK+qbW9rkdCc3N0UKIardXm7FihVwcXExfwUFBTV7HaP8nNE72LXZn4eIiIha3p7LGYjNKISHgxpju/iiLIq0OqtaUkwmk2Hr1q2YNm1ajWVuu+029O7dG++++65529atWzFjxgwUFRVVeylWp9NBV6n5t7wTYm5ubrP2sRNCYPvha1i14xKScqRm4u4BznhmQiT6hrjzMktNZXmZpeFleZlFus1LsY0ry/8R0m3+j2h42Xb0P+JKaj5W7byM3ZfSpaIaFeaM6oyHbguHk13zdgPLy8uDi4tLvXJLmwt2zzzzDH788UdER0ebtz3yyCM4deoUDh06VK/nacgPqCmU6I1YfyAWa3ZdQ4FO+kcwsbsvnh0fhWAPbR2PJiIiotaSnFuMt3dexrfHE2ASgEIuw8z+QVh4e8ea+9A3sYbkllYfPFFQUICrV6+a78fGxuLUqVNwd3dHcHAwlixZgsTERHz66acAgPnz5+P999/HokWL8NBDD+HQoUNYt24dvvzyy9Z6CXWyUymwYGQH/LVvEFbtvIyvj97A9rMp+C06DXOHhuLR0R3g3Mxpn4iIiOovr0SPtbuvYf3+WOgMUmvy+K6++Mf4zojwcmzl2tWs1Vvsdu/ejVGjRlXZPnv2bGzcuBFz5sxBXFwcdu/ebd63Z88ePPXUUzh//jz8/f3xzDPPYP78+fV+zpZusbvZxZQ8vLrtAvZdyQAAuDuo8dSYjvi/AcFQKtpkt0ciIiKboDMY8fmfN/D+H1eQXSRd6u4X4oYlE8u6UbWCNnsptqW0drADpP53uy+l45Vt0biWLvWF6eDtiKWTojCykxeXOiMiImpBJpPAj2eS8Mavl5CQLfWti/BywDPjI3FHF59WfV9msKuDNQS7cnqjCV8duYFVOy+bPxkM7+iJ5yZ1QWffumeYJiIioluz/0oGVvx8AeeT8gAA3k4aPHVHJ/y1b6BVXEljsKuDNQW7crnFenyw6yo2HIiF3igglwH3DAjGU2M6wctJ09rVIyIisjnnk3Lx758vmrtGOWqUmD8iHA8MC4NW3erDEMwY7OpgjcGu3PXMQvz754v4+VwKAOmXbMGoCDwwNAx2KkUr146IiKjti88qwqqdl/H9qUQIAagUMtw3MASPj+4AD0fra0xhsKuDNQe7ckdis/DKtmicScgFAAS42uOZCZGY0sOP/e+IiIgaIbuwFB/suopPD11HqVEa6Tqlpz/+MbazVU8/xmBXh7YQ7ACpI+f3pxKx8pdLSMmTJl/sHeyK5yd3QZ9gt1auHRERUdtQojdiw4E4rNl9Ffkl0nyyQyI88OyESPQIdG3dytUDg10d2kqwK1dcasQn+2Kwdvc1FOulGdSn9vTHP8d3RqCb9X7CICIiak1Gk8B3JxLw9s7LSM6VGkgifZ3w7IRIjGhDM1Aw2NWhrQW7cql5JXjz10v49kQChADUSjkeHBaGR0ZGNPtyJkRERG2FEAK7LqXh9Z8v4VJqPgDA38UOT4/tjGm9A6CQt41AV47Brg5tNdiVO5eYi1e2RePPmCwAgKejGovu6IyZ/YPa3C8rERFRUzoVn4MV2y/gcKz0Hulir8KjoyIwa3Bomx2EyGBXh7Ye7ADp08jO6FSs+PkiYjOkCY4jfZ2wdFIUhnf0auXaERERtay4jEK88eslbDubDEC6qjV3SCgWjOwAF23bvqrFYFcHWwh25UoNJnz+53W8+/sV5BZLExyP6uyFpZOi0MGbExwTEZFtyyjQYfXvV7Dp8A0YTAIyGXB370AsGtsJAa72rV29JsFgVwdbCnblcopK8e7vV/DZoeswmAQUchnuGxiMJ2/vaJVz8hAREd2K/BI9PtkXi3X7YlBYKg0sHNnZC8+Mj0SUn228t5djsKuDLQa7cjHpBXht+0X8diEVgPXOok1ERNQYJXojPjt0HWt2XzUvxdkj0AXPTojEkAjPVq5d82Cwq4MtB7tyB69m4LWfL+BcYsW6d4vu6ITpVrLuHRERUUMYjCZ8ezwB7/5+xTx1SbiXA/4xtjPGd/NtM1OXNAaDXR3aQ7ADpAmOfzyThDd+vYSE7GIAQEdvRzwzPhK3R3nb9B8BERHZBpNJYPu5ZKzacRkxZYMF/V3ssHBMJ9zdJ6BdNFYw2NWhvQS7cjqD1Gz9/q6ryClrth4Q5o4lEyLRmytYEBGRFRJCYO+VDKz85SLOJ0lXn9wd1Hh0VAfcNzC4zU5d0hgMdnVob8GuXG6xHmt3X8P6A7EoNUhr5E3q7od/jOuMUE+HVq4dERGR5Pj1bKz85aJ5LjpHjRIPDg/DvGFh7XJCfga7OrTXYFcuKacYq3ZexndlK1go5TL8bVAIHh/dgSNoiYio1VxKyccbv14yDwBUK+W4f1AIFoyMaNfvTwx2dWjvwa7cheQ8vP7LRey+lA6AI2iJiKh13Mgswtu/Xcb3pxIhBCCXAX/tG4Qnx3SEv43MRXcrGOzqwGBn6cDVDKzgCFoiImphaXklWP3HVXx19Ab0RimOTOzui0V3dEYHb8dWrp31YLCrA4NdVRxBS0RELSW3WI+P9lzDhgNxKNZLkwsP7+iJf46LRPdAl1aunfVhsKsDg13NOIKWiIiaS3GpERsOxuLD3deQV2IAAPQOdsU/x0VicIRHK9fOejHY1YHBrm7lI2g3HIiFjiNoiYjoFuiNJnx1NB6rf7+CtHwdAKCTjyMWj+2MO7r48KpQHRjs6sBgV3/VjaC9b2AwHr+9Izzb8QglIiKqm8kk8L/TSVi18zJuZBUBAALd7LHojk64s1cAFHIGuvpgsKsDg13DcQQtERHVlxACf1xMwxu/XsLFlHwAgKejBo+P7oD/GxAMtZID8xqCwa4ODHaNxzVoiYioNodjMrHy10s4fj0bAOBkp8T8ERGYOzSUDQGNxGBXBwa7W8MRtEREdLNzibl449dL2HNZurKjUcoxZ2goHhkRAVetupVr17Yx2NWBwa5p6AxGfP7nDaz+4wpH0BIRtVNX0wrwzm+X8dOZZABSX+yZ/YPwxO0d4eNs18q1sw0MdnVgsGtaNY2gXTyuM8I4gpaIyCYdi8vCh3tizMt/AcDUnv5YdEcnzp7QxBjs6sBg1zxuHkErlwFTevpj/ogIRPnx50xE1NaZTAI7olPx8d5rOHEjx7z9ji4+WDimI7r6c3Lh5sBgVwcGu+Z1ITkPK3+5iF1lI2gBYFRnLzwysgMGhLm3Ys2IiKgxSvRGbDmRiP/si0FMRiEAQK2Q4+4+AXhweDiX/2pmDcktVjGMcc2aNQgLC4OdnR369u2Lffv21Vp+06ZN6NmzJ7RaLfz8/DB37lxkZma2UG2pLlF+ztgwdwB+enwYJvfwg1wG7LqUjhkfHcL0tQfx+4VUmEzt7vMEEVGbk1NUivf/uIJhr+/Cv7aeRUxGIZztlFgwMgL7nxmFf/+lB0OdlWn1Fruvv/4a999/P9asWYOhQ4fio48+wn/+8x9ER0cjODi4Svn9+/djxIgRePvttzFlyhQkJiZi/vz56NixI7Zu3Vqv52SLXcuKyyjER3tj8N3xBJQapT54nX2cMH9kOCb38IeK06QQEVmVhOwirNsfi6+PxqOoVFrL1d/FDg8MC8M9A4LhqOG0JS2pTV2KHThwIPr06YO1a9eat0VFRWHatGlYsWJFlfJvvvkm1q5di2vXrpm3rV69GitXrkR8fHy9npPBrnWk5ZVg3YFYbPrzBgp00hqBgW72+Ptt4fhr3yDYqxWtXEMiovbtfFIuPt4bg5/OJMNYdmUl0tcJD4/gB/HW1GaCXWlpKbRaLTZv3oy77rrLvP3JJ5/EqVOnsGfPniqPOXjwIEaNGoWtW7diwoQJSEtLw4wZMxAVFYUPP/yw2ufR6XTQ6XTm+3l5eQgKCmKwayW5xXp8/ud1rN8fi8zCUgCAh4Mac4eG4v7BoXCxV7VyDYmI2g8hBA5czcRHe69h35UM8/ahHTzw8G0RGN7Rk/OTtrKGBLtWbUvNyMiA0WiEj4+PxXYfHx+kpKRU+5ghQ4Zg06ZNmDlzJkpKSmAwGDB16lSsXr26xudZsWIFli9f3qR1p8ZzsVfh0VEdMG9YGDYfi8dHe2OQkF2MN3dcxod7YnDfwGDMGxYGb85/RETUbAxGE7adTcZHe2IQnSytJiSXAZN6+OPh28LRLYAjXNuiVm2xS0pKQkBAAA4ePIjBgwebt7/66qv47LPPcPHixSqPiY6OxpgxY/DUU09h3LhxSE5Oxj/+8Q/0798f69atq/Z52GJn3QxGE346k4y1u6/hUqq0pqBaIcdf+gbi4dvCOR8SEVETKtQZ8M2xePxnXywSc6TVg+xVCszsH4R5w8IQ5K5t5RrSzWz6Uuz999+PkpISbN682bxt//79GD58OJKSkuDn51fn87KPnXUSQmDXpTSs2XUNx8rWGJTLgAnd/fDIiAh+eiQiugXp+Tr892AcPvvzOnKLpdWCPBzUmD0kFPcPCoGbA5f9slZt5lKsWq1G3759sXPnTotgt3PnTtx5553VPqaoqAhKpWW1FQqp0307nJLPpshkMoyO9MHoSB8cjcvC2t3X8MfFNGw7k4xtZ5JxWycvPDIiAoPC3dnfg4ionmLSC/DJvlh8dyIBpWWrA4V6aPHQbeH4S59A2Kk4cM2WtPp45UWLFuH+++9Hv379MHjwYHz88ce4ceMG5s+fDwBYsmQJEhMT8emnnwIApkyZgoceeghr1641X4pduHAhBgwYAH9//9Z8KdSE+oe6o/8cd1xIzsNHe67hxzPJ2Hs5HXsvp6NXkCseGRmBO6J8IJcz4BERVefEjWx8tOcadkSnorzdo1eQK+aPCMcdXXyh4P9Pm9Tq050A0gTFK1euRHJyMrp164a3334bt912GwBgzpw5iIuLw+7du83lV69ejQ8//BCxsbFwdXXF6NGj8frrryMgIKBez8dLsW3PjcwifLIvBt8cizevR9vB2xHzR0Tgzl4cgk9EBEhLfv1xMQ0f7b2Go3HZ5u23R3rj4RER6B/qxisebVCb6WPXWhjs2q70fB02HozFp4euI79EmgvP38UODw4Pxz0DgqBVt3ojNBFRi9MZjPjhZBI+2nsN19KlJb9UChmm9QrA328LR0cfp1auId0KBrs6MNi1fXklenxx+AbW7Y9Fer404tlNq8KcIWGYPSQErlp2AiYi25ddWIovjtzAfw/GIa3sf6GTRon7BoVg7tBQ+HDaKJvAYFcHBjvbUb4w9Ud7r+F6ZhEAQKtW4P8GSHPh+bvat3INiYiaXkx6AdYfiMW3xxNQope6p/g622HesDDcMyAITnac6N2WMNjVgcHO9hiMJvx8LgVrd18zT7SplMswqYcfHhrOiTaJqO0TQuBQTCbW7YvF7xfTzNu7+DnjweFhmNzDH2ol+xvbIga7OjDY2S4hBPZeycDa3VfxZ0yWefugcHc8OCwcoyO9OZKWiNqUUoMJP51Jwn/2xZo/uALAmChvzBsWzimg2gEGuzow2LUPZxNysW6/tJi1oWwx63AvB8wbFsa5m4jI6uUUlWLTYcv+c3YqOf7aNwhzh4Yi3MuxlWtILYXBrg4Mdu1LUk4x/nswDl8cuWEeSevuoMbfBgbj/sGh8HLStHINiYgqVNd/zttJg9lDQnHvgGCuENEOMdjVgcGufSrQGfDN0XisPxCLhGxpfUS1Uo67egVg3vAwdOJ0AETUSsr7z63fH4vfLrD/HFlisKsDg137ZjCa8Ov5VHyyLwan4nPM20d08sJDw8MxtIMH+6sQUYuorf/cA8PCMDic/4+Iwa5ODHZU7vj1LHyyNxa/RqeYl9yJ9HXCg8PDMbUnPyETUfOoqf/c9L6BmDs0DBHsP0eVMNjVgcGObnY9sxAbDsThm2PxKCo1Aqjo03LfwGBOeExETYL956gxGOzqwGBHNckt0mPTkev478E4pOZJn6LtVQr8tV8gHhgahlBPh1auIRG1NUII/BmThXX7Y9h/jhqFwa4ODHZUl/J+L5/si8WFsn4vMhkwtosPHhwejn4hXEibiGpX/n9k3f5YnE9i/zlqPAa7OjDYUX0JIXDoWiY+2ReDXZfSzdt7BrniwWFhmNDNF0oFP2kTUQX2n6OmxmBXBwY7aowrqflYtz8WW04motQg9Y0JcLXH3KGhmNmfazMStXcx6QXYcCAO3x5PQLHesq8u+8/RrWCwqwODHd2KjAIdPjt0HZ/9eR1ZhaUAACeNEvcMCMKcoWEIcLVv5RoSUUvJLNDhtwup2H42BXuvpJtH17P/HDUlBrs6MNhRUyjRG7H1ZCL+sy8G19ILAQAKuQyTuvvhweFh6B7gwv4zRDYoPqsIv55PwY7oVByLy4Kp0rvo7ZHemDec/eeoaTHY1YHBjpqSySSw+3IaPtkbi0Mxmebt4Z4OuKOLD8Z29UGvIDco5PwnT9QWCSFwMSVfCnPnUy0mEgaAbgHOGNfFF5N6+HH9VmoWDHZ1YLCj5nIuMRfr9sdi25lklBpN5u2ejmqMifLBHV18MLSDJ+xUilasJRHVxWgSOHEjG7+ek1rmbmQVmffJZcCAMHeM6+qLO7r4INBN24o1pfaAwa4ODHbU3Ap0Buy5lI4d0Sn442Ia8ksM5n1atQK3dfTC2K4+GB3pzcmPiayEzmDEwauZ2BGdgp3RqcgoKDXv0yjlGN7RC+O6+uD2KB+4cyAEtSAGuzow2FFL0htNOByThZ3R0if/5NwS8z6FXIYBoe64o4vUmhfkzk/+RC0pv0SPXZfSseN8CnZfSkeBruJDmLOdErdH+WBcVx/c1skLWrWyFWtK7RmDXR0Y7Ki1CCFwPikPO8o6Xl9MybfYH+XnjLFlIa+rvzM7XxM1g/R8aSTrr+dTcPBqpkW3CR9nDcZ28cXYrj4YFO4BFeepJCvAYFcHBjuyFjcyi7DzQip2nE/B0ZtG1wW42ptb8gaEufMNhugWXM8sxI7zUpg7fiMbld/5wr0cMK6rL8Z28UHPQFfIOdCJrAyDXR0Y7MgaZRWW4o+LadgZnYI9l9PNC4QD0iWh0ZHeGNvVF7d18oKjhpeEiGojhEB0ch5+PS99cLq5dbxnoAvGdvXFuK4+6ODt1Eq1JKofBrs6MNiRtSvRG7H/SgZ2RKfg9wtpyCys6MStVsgxtIMH7ujiizFdvOHtZNeKNSWyHkaTwLG4LCnMRacgIbvYvE8hl2FQeMVIVj8XTiRObQeDXR0Y7KgtKZ92obxf3vXMimkXZDKgV5ArxnaR3qw6eHMOLWo/hBCIyyzC4ZhM/BmTib1XMsyrwQDS+qwjOnlhXFdfjkCnNo3Brg4MdtRWCSFwJa0AO6NTsSM6Fafjcyz2h3tJkyLfHumDnkEu0Cg5Xx7ZDiEEYjMKcTg2C3+WhbnUPJ1FGVetCrdHShOD39bRC/Zq/g1Q28dgVwcGO7IVKbkl+O2CFPIOXcuA3ljx56xWytEr0BX9w9zQP9QdfUPc4GSnasXaEjVMeZD7M6YiyKXlWwY5tUKOXsGuGBTugcHhHugf6gYlBxqRjWGwqwODHdmivBJ92aTIqTh4NcOiXx4gzZYf5eeM/qHuGBDmjn6hbuyfR1ZFCIGYjMKyECeFufRqglzvsiA3KNwDvYNduZIL2TwGuzow2JGtK3+DPBqbhSNxWTgWl22xJFK5UA8t+oe6o3+YOwaEuiPEQ8u586jFCCFwLb3Q3Br3Z0wWMgpuCnJKOfoEu2JgGIMctV9tLtitWbMGb7zxBpKTk9G1a1e88847GD58eI3ldTodXnrpJXz++edISUlBYGAgli5digceeKBez8dgR+1RSm4JjsZl4WhcFo7EZuFSaj5u/uv3ctJgQKg7+oe6oX+YOyJ9naHgnF7URKQgV4BDZa1xh2sJcuUtcr2CGOSI2lSw+/rrr3H//fdjzZo1GDp0KD766CP85z//QXR0NIKDg6t9zJ133onU1FS88sor6NChA9LS0mAwGDBkyJB6PSeDHRGQW6zHievZOBKXhaOxWTiTkGsxAz8AOGmU6BPihgFh7ugf6o4egS58k6V6E0LgalqBuTXucGymxfqrgBTk+ga7lQU5d/RkkCOqok0Fu4EDB6JPnz5Yu3ateVtUVBSmTZuGFStWVCn/yy+/4J577kFMTAzc3d0b9ZwMdkRVleiNOB2fU9aql43j17Mt1s0EpP5NPYNczJdv+4a4wZkDMqhM+ajt8ta4P2Myq/T11Cjl6BsiBbmBYQxyRPXRZoJdaWkptFotNm/ejLvuusu8/cknn8SpU6ewZ8+eKo9ZsGABLl++jH79+uGzzz6Dg4MDpk6dipdffhn29vWbcJLBjqhuRpPAheS8Spdvs6tcNpPJgChfZ/Ol2wGh7vB25oAMW2c0CSTlFONqegGupRUgJqMQ19IKcDk1H9lFeouylYPcoHAPTsND1AgNyS2tui5RRkYGjEYjfHx8LLb7+PggJSWl2sfExMRg//79sLOzw9atW5GRkYEFCxYgKysL69evr/YxOp0OOl3FG1JeXl7TvQgiG6WQy9AtwAXdAlwwd2iYeTLYo2WXbo/GZSEuswjRyXmITs7Dfw9dBwCEeGjRxc8ZwR5ahHo4IMRdi2APLfxc7Nlfr40p0BkQk16Aa+kFiEkvNH+PyShEqcFU7WPsVGVBLswDgyI80COQQY6oJVnFgpM3j8ITQtQ4Ms9kMkEmk2HTpk1wcXEBAKxatQrTp0/HBx98UG2r3YoVK7B8+fKmrzhROyKTyRDm6YAwTwfM6BcEAEjLK8HRuGxzq96F5DxczyyyWB2jnFohR6C7PUI9HBDsrkVIWfAL9tAi0M2eb/6txGQSSMotxrX0wioh7ubJfytTK+QI83RAhLcDwj0dzd8j/Zx4LolaUasGO09PTygUiiqtc2lpaVVa8cr5+fkhICDAHOoAqU+eEAIJCQno2LFjlccsWbIEixYtMt/Py8tDUFBQE70KovbL29kOk3r4YVIPPwDSXHonb+TgWloBbmQV4XpmIa5nFiE+uwilRpPU2pNeWOU4Mhng72KPEA8tQjwcpO/uFbcdNFbxGbRNK9QZEJshBbaKEFeI2IwClOirb30DAE9HDSK8HBDu5YgILwdEeDkiwssRAW5sgSWyRq3631KtVqNv377YuXOnRR+7nTt34s4776z2MUOHDsXmzZtRUFAAR0dpXczLly9DLpcjMDCw2sdoNBpoNJqmfwFEZMHZToURnbwwopOXxXajSSA5t9jcmlce+K6Xhb+iUiMSc4qRmFOMg9cyqxzX01Ethbyyy7rlLX2hHg5w06o4914ZIQSSc0sqXTatCHFJuSU1Pk6lkCHUwwHhZcGtPMSFeznCxZ6DY4jaklYfFVs+3cmHH36IwYMH4+OPP8Ynn3yC8+fPIyQkBEuWLEFiYiI+/fRTAEBBQQGioqIwaNAgLF++HBkZGXjwwQcxYsQIfPLJJ/V6Tg6eILIeQghkFJRWCXvXM4twI6vIYlH36jhplBZhL8hNCyc7JexVCtirFbBTKcy3tZXuqxQyqwqEQgjoDCbkleiRX2JAfokBecXltyu+55XvK9FX2l5RxmCq+V+6h4O6UnirCHFBbvZchovIirWZwRMAMHPmTGRmZuKll15CcnIyunXrhu3btyMkJAQAkJycjBs3bpjLOzo6YufOnXj88cfRr18/eHh4YMaMGXjllVda6yUQ0S2QyWTwctLAy0mDfqFVpzDKK9HjRmYR4srDXtntG1lFSM4tQb7OgPNJeTif1LBBUQq5DPaqsqCnlpeFPyXsVfKqoVBVFgrVFfdvDo3l32WARdgqD2p5NwW0inBWUbbyWr+NpZTLEOyhtQhvEV5S/zc3B/UtH5+IrFurt9i1BrbYEdmGEr0R8VnS5d3ysJeYXYzCUgOK9SaUlBpRrDeiqNSIEr0RRaUG1NKgZRVkMsBRo4SznQpOdko42VW+ffN3aZ+zfcU2T0cNVGx9I7IpbarFjoiosexUCnT0cUJHH6d6lRdCQG8UKNZLQa+4LPgV640oKZUCoPl+Nful2yYUlxqk26XS/cplTSZhEb6c7StCWOVg5mwR2irKOKiVkHNQAhE1EoMdEbUbMpkMaqUMaqWcgwKIyCaxvZ6IiIjIRrTLFrvyboVcgYKIiIisXXleqc+wiHYZ7PLz8wGAkxQTERFRm5Gfn2+xQEN12uWoWJPJhKSkJDg5OTXrPFblK1zEx8dz9K2V4DmxPjwn1onnxfrwnFinljgvQgjk5+fD398fcnntvejaZYtdbatUNAdnZ2f+EVoZnhPrw3NinXherA/PiXVq7vNSV0tdOQ6eICIiIrIRDHZERERENoLBrhlpNBosW7YMGo2mtatCZXhOrA/PiXXiebE+PCfWydrOS7scPEFERERki9hiR0RERGQjGOyIiIiIbASDHREREZGNYLAjIiIishEMdrdozZo1CAsLg52dHfr27Yt9+/bVWn7Pnj3o27cv7OzsEB4ejg8//LCFatp+NOScbNmyBXfccQe8vLzg7OyMwYMH49dff23B2rYPDf07KXfgwAEolUr06tWreSvYTjX0vOh0OixduhQhISHQaDSIiIjA+vXrW6i27UNDz8mmTZvQs2dPaLVa+Pn5Ye7cucjMzGyh2tq+vXv3YsqUKfD394dMJsP3339f52Na/X1eUKN99dVXQqVSiU8++URER0eLJ598Ujg4OIjr169XWz4mJkZotVrx5JNPiujoaPHJJ58IlUolvv322xauue1q6Dl58sknxeuvvy6OHDkiLl++LJYsWSJUKpU4ceJEC9fcdjX0nJTLyckR4eHhYuzYsaJnz54tU9l2pDHnZerUqWLgwIFi586dIjY2Vhw+fFgcOHCgBWtt2xp6Tvbt2yfkcrl49913RUxMjNi3b5/o2rWrmDZtWgvX3HZt375dLF26VHz33XcCgNi6dWut5a3hfZ7B7hYMGDBAzJ8/32JbZGSkePbZZ6st/89//lNERkZabHv44YfFoEGDmq2O7U1Dz0l1unTpIpYvX97UVWu3GntOZs6cKZ577jmxbNkyBrtm0NDz8vPPPwsXFxeRmZnZEtVrlxp6Tt544w0RHh5use29994TgYGBzVbH9qw+wc4a3ud5KbaRSktLcfz4cYwdO9Zi+9ixY3Hw4MFqH3Po0KEq5ceNG4djx45Br9c3W13bi8ack5uZTCbk5+fD3d29OarY7jT2nGzYsAHXrl3DsmXLmruK7VJjzsv//vc/9OvXDytXrkRAQAA6deqExYsXo7i4uCWqbPMac06GDBmChIQEbN++HUIIpKam4ttvv8WkSZNaospUDWt4n1e2yLPYoIyMDBiNRvj4+Fhs9/HxQUpKSrWPSUlJqba8wWBARkYG/Pz8mq2+7UFjzsnN3nrrLRQWFmLGjBnNUcV2pzHn5MqVK3j22Wexb98+KJX8F9UcGnNeYmJisH//ftjZ2WHr1q3IyMjAggULkJWVxX52TaAx52TIkCHYtGkTZs6ciZKSEhgMBkydOhWrV69uiSpTNazhfZ4tdrdIJpNZ3BdCVNlWV/nqtlPjNfSclPvyyy/x4osv4uuvv4a3t3dzVa9dqu85MRqNuPfee7F8+XJ06tSpparXbjXkb8VkMkEmk2HTpk0YMGAAJk6ciFWrVmHjxo1stWtCDTkn0dHReOKJJ/DCCy/g+PHj+OWXXxAbG4v58+e3RFWpBq39Ps+Pw43k6ekJhUJR5ZNUWlpalbReztfXt9rySqUSHh4ezVbX9qIx56Tc119/jXnz5mHz5s0YM2ZMc1azXWnoOcnPz8exY8dw8uRJPPbYYwCkQCGEgFKpxI4dOzB69OgWqbsta8zfip+fHwICAuDi4mLeFhUVBSEEEhIS0LFjx2ats61rzDlZsWIFhg4din/84x8AgB49esDBwQHDhw/HK6+8wqtArcAa3ufZYtdIarUaffv2xc6dOy2279y5E0OGDKn2MYMHD65SfseOHejXrx9UKlWz1bW9aMw5AaSWujlz5uCLL75g35Qm1tBz4uzsjLNnz+LUqVPmr/nz56Nz5844deoUBg4c2FJVt2mN+VsZOnQokpKSUFBQYN52+fJlyOVyBAYGNmt924PGnJOioiLI5ZZv4wqFAkBFKxG1LKt4n2+xYRo2qHxo+rp160R0dLRYuHChcHBwEHFxcUIIIZ599llx//33m8uXD4N+6qmnRHR0tFi3bh2nO2liDT0nX3zxhVAqleKDDz4QycnJ5q+cnJzWegk2p6Hn5GYcFds8Gnpe8vPzRWBgoJg+fbo4f/682LNnj+jYsaN48MEHW+sl2JyGnpMNGzYIpVIp1qxZI65duyb2798v+vXrJwYMGNBaL8Hm5Ofni5MnT4qTJ08KAGLVqlXi5MmT5ilorPF9nsHuFn3wwQciJCREqNVq0adPH7Fnzx7zvtmzZ4sRI0ZYlN+9e7fo3bu3UKvVIjQ0VKxdu7aFa2z7GnJORowYIQBU+Zo9e3bLV9yGNfTvpDIGu+bT0PNy4cIFMWbMGGFvby8CAwPFokWLRFFRUQvX2rY19Jy89957okuXLsLe3l74+fmJ++67TyQkJLRwrW3Xrl27an2PsMb3eZkQbK8lIiIisgXsY0dERERkIxjsiIiIiGwEgx0RERGRjWCwIyIiIrIRDHZERERENoLBjoiIiMhGMNgRERER2QgGOyIiIiIbwWBHREREZCMY7IiIiIhsBIMdERERkY1gsCMiIiKyEQx2RERERDaCwY6IiIjIRjDYEREREdkIBjsiIiIiG8FgR0RERGQjGOyIqE3YuHEjZDKZ+UupVCIwMBBz585FYmJia1ev2SQlJeHFF1/EqVOnWrsqRNQGKFu7AkREDbFhwwZERkaiuLgYe/fuxYoVK7Bnzx6cPXsWDg4OrV29JpeUlITly5cjNDQUvXr1au3qEJGVY7AjojalW7du6NevHwBg1KhRMBqNePnll/H999/jvvvua/Rxi4uLYW9v31TVtHrFxcWws7ODTCZr7aoQURPipVgiatMGDRoEALh+/TqWL1+OgQMHwt3dHc7OzujTpw/WrVsHIYTFY0JDQzF58mRs2bIFvXv3hp2dHZYvXw4A+OCDD3DbbbfB29sbDg4O6N69O1auXAm9Xm9xjJEjR6Jbt244dOgQhgwZAnt7e4SGhmLDhg0AgG3btqFPnz7QarXo3r07fvnllyp1v3LlCu699154e3tDo9EgKioKH3zwgXn/7t270b9/fwDA3LlzzZehX3zxRXOZY8eOYerUqXB3d4ednR169+6Nb775xuJ5yi9j79ixAw888AC8vLyg1Wqh0+ka+VMnImvFFjsiatOuXr0KAPDy8sLBgwfx8MMPIzg4GADw559/4vHHH0diYiJeeOEFi8edOHECFy5cwHPPPYewsDDzZdxr167h3nvvRVhYGNRqNU6fPo1XX30VFy9exPr16y2OkZKSgrlz5+Kf//wnAgMDsXr1ajzwwAOIj4/Ht99+i3/9619wcXHBSy+9hGnTpiEmJgb+/v4AgOjoaAwZMgTBwcF466234Ovri19//RVPPPEEMjIysGzZMvTp0wcbNmzA3Llz8dxzz2HSpEkAgMDAQADArl27MH78eAwcOBAffvghXFxc8NVXX2HmzJkoKirCnDlzLOr7wAMPYNKkSfjss89QWFgIlUrVtCeDiFqfICJqAzZs2CAAiD///FPo9XqRn58vfvrpJ+Hl5SWcnJxESkqKRXmj0Sj0er146aWXhIeHhzCZTOZ9ISEhQqFQiEuXLtX6nOXH+PTTT4VCoRBZWVnmfSNGjBAAxLFjx8zbMjMzhUKhEPb29iIxMdG8/dSpUwKAeO+998zbxo0bJwIDA0Vubq7Fcz722GPCzs7O/FxHjx4VAMSGDRuq1C8yMlL07t1b6PV6i+2TJ08Wfn5+wmg0WvzsZs2aVevrJaK2j5diiahNGTRoEFQqFZycnDB58mT4+vri559/ho+PD/744w+MGTMGLi4uUCgUUKlUeOGFF5CZmYm0tDSL4/To0QOdOnWqcvyTJ09i6tSp8PDwMB9j1qxZMBqNuHz5skVZPz8/9O3b13zf3d0d3t7e6NWrl7llDgCioqIASJeLAaCkpAS///477rrrLmi1WhgMBvPXxIkTUVJSgj///LPWn8PVq1dx8eJFc7/Cm4+RnJyMS5cuWTzmL3/5S10/XiJq43gplojalE8//RRRUVFQKpXw8fGBn58fAODIkSMYO3YsRo4ciU8++QSBgYFQq9X4/vvv8eqrr6K4uNjiOOWPq+zGjRsYPnw4OnfujHfffRehoaGws7PDkSNH8Oijj1Y5hru7e5VjqNXqKtvVajUAKdABQGZmJgwGA1avXo3Vq1dX+zozMjJq/TmkpqYCABYvXozFixfX6xjVvWYisi0MdkTUpkRFRZlHxVb21VdfQaVS4aeffoKdnZ15+/fff1/tcaobDfr999+jsLAQW7ZsQUhIiHl7U88h5+bmBoVCgfvvvx+PPvpotWXCwsJqPYanpycAYMmSJbj77rurLdO5c2eL+xwBS2T7GOyIyCaUT1qsUCjM24qLi/HZZ5816BgAoNFozNuEEPjkk0+arqIAtFotRo0ahZMnT6JHjx7mFr3qlNfl5tbCzp07o2PHjjh9+jRee+21Jq0fEbVdDHZEZBMmTZqEVatW4d5778Xf//53ZGZm4s0337QIaXW54447oFar8X//93/45z//iZKSEqxduxbZ2dlNXt93330Xw4YNw/Dhw/HII48gNDQU+fn5uHr1Kn788Uf88ccfAICIiAjY29tj06ZNiIqKgqOjI/z9/eHv74+PPvoIEyZMwLhx4zBnzhwEBAQgKysLFy5cwIkTJ7B58+YmrzcRWTcOniAimzB69GisX78eZ8+exZQpU7B06VJMnz4dzz77bL2PERkZie+++w7Z2dm4++678fjjj6NXr1547733mry+Xbp0wYkTJ9CtWzc899xzGDt2LObNm4dvv/0Wt99+u7mcVqvF+vXrkZmZibFjx6J///74+OOPAUgTNB85cgSurq5YuHAhxowZg0ceeQS//fYbxowZ0+R1JiLrJxPippk7iYiIiKhNYosdERERkY1gsCMiIiKyEQx2RERERDaCwY6IiIjIRjDYEREREdkIBjsiIiIiG9EuJyg2mUxISkqCk5MTl9ghIiIiqyaEQH5+Pvz9/SGX194m1y6DXVJSEoKCglq7GkRERET1Fh8fj8DAwFrLtMtg5+TkBED6ATk7O7dybYiIiIhqlpeXh6CgIHN+qU27DHbll1+dnZ0Z7IiIiKhNqE/3sXYZ7IiIiNq6Er0RqXklSM3TIaNAB3uVAi5aFVztVXDVquFsp4RSwTGS7Q2DHRERkRUxmQSyikqRkluC1LwSpJSFt9Tc8tvS95wifZ3HctIopbCnVcHVXg0Xe1Wl8Cdtcy6/XXbfVauCnUrRAq+UmgODHRERUQsp0RuRUjmgVbqdmqdDSm4J0vJLoDeKeh3PTiWHr7MdPB010BlMyCkuRU6RHvklBgBAvs6AfJ0BCdnFDaqnWim3CH8uWhVc7CsCoYtWDVf7sm1lZdwd1XBQKzjbRCtjsCMiIrpFJpNARqEOaWXhrLrglpJbgryywFUXmQzwcNDA10UDX2c7+DjbSd9d7CzuO9srqw1SBqMJeSUG5BbrkVNUipxiPXKLKt0uv3/T/txiPQwmgVKDCWn5OqTl6xr0c1Ar5XDXquHuUL8vN60aCjmDYFNisCMiImoEo0ngcGwmfjydjF/OJSO7HpdGAcBepYCvix18nDUWYa3ybS8nDVS30D9OqZCbwxPgUO/HCSFQWGqUwl5Z0JPCoR45xaVl4bDsdtn23GI9sgpLoTOYUGowIaXsUnF9yGSAi71KqmtZIPRwlAJf5QDo4aCBm4MKHg4a2Kt5mbg2DHZERET1ZDIJnIzPxo+nk7HtbDLSK7VoyWSAl6MGvi528Hays2xtc6kIbk6a6lvZrIFMJoOjRglHjRKBbg17bHGpEZmFOmQVltb+VSR9zynSQwhIQbFIjxgU1ut57FRyeDhopBY/BzU8HNSwUykghIBJCJgEYBICQqBsW8X9yt9N5v0CArC4bzIBApXLVH/M8v2Dwj3w4tSuDf+BNwMGOyIioloIIXA+KQ8/nk7CT2eSkZhT0V/NxV6FCd18MaWnPwaEud9SK1tbZ69WIFCtRaCbtl7lDUYTcspa++r7VWo0oURvQmJOscV5aG2BbvatXQUzBjsiIqJqXEnNx4+nk/DjmWTEZlS0JjlqlBjbxQdTevpjaAdPqJXtN8zdCqVCDk9HDTwdNfUqX36ZOKugvNVPh6xCPbIKddDpTZDLZZDJALlMBhnKvpfdl8uk1siK7+X7Kt0HIJeXP66sLKp7bMWxy+9Ll7ytA4MdERFRmbiMQvx0RmqZu5iSb96uUcoxJsoHU3r6YWRnb04H0goqXyYO9qhfq2B7xGBHRETtWlJOMbadScaPZ5JwJiHXvF2lkGFEJy9M6emP26N84KjhWyZZP/6WEhFRu5Oer8PP55Lx4+kkHI3LNm9XyGUYEuGBKT39Ma6LL1y0qlasJVHDMdgREVG7kFNUil/Pp+DH08k4eC0DprI5gGUyoH+oO6b09MeEbr717vNFZI0Y7IiIyGYV6AzYGS2FuX1X0i1WdOgZ5IopPfwwqYcf/FysZ1Qj0a1gsCMiIptSojfij4tp+PF0Ev64mAadwWTeF+nrhCk9/TGlhz874JNNYrAjIqI2r9Rgwr4r6fjxdBJ2RqeisNRo3hfu6SCFuZ5+6ODt1Iq1JGp+DHZERNRmleiN+OZYPD7cfQ1JuRXLWAW42pvDXBc/Z6td6YGoqTXbrIpxcXGYN28ewsLCYG9vj4iICCxbtgylpaW1Pi41NRVz5syBv78/tFotxo8fjytXrliU0el0ePzxx+Hp6QkHBwdMnToVCQkJzfVSiIjIyhSXGrF+fyxuW7kLL/xwHkm5JfB01GDu0FBsWTAE+58ZhWcnRKKrvwtDHbUrzdZid/HiRZhMJnz00Ufo0KEDzp07h4ceegiFhYV48803q32MEALTpk2DSqXCDz/8AGdnZ6xatQpjxoxBdHQ0HBykhYwXLlyIH3/8EV999RU8PDzw9NNPY/LkyTh+/DgUCk4aSURkqwp1Bmw6fB0f741BRoHUUODnYocFIyPw135BnDiY2j2ZEELUXaxpvPHGG1i7di1iYmKq3X/58mV07twZ586dQ9eu0mK6RqMR3t7eeP311/Hggw8iNzcXXl5e+OyzzzBz5kwAQFJSEoKCgrB9+3aMGzeuznrk5eXBxcUFubm5cHZ2broXSEREzSK/RI9PD13Huv2xyCqUAl2gmz0WjOyAv/QNgEbJQEe2qyG5pUX72OXm5sLd3b3G/TqdDgBgZ2dn3qZQKKBWq7F//348+OCDOH78OPR6PcaOHWsu4+/vj27duuHgwYPVBjudTmc+NiD9gIiIyPrlFuux8UAc1h+IRW6xHgAQ6qHFglEdcFfvAKgUXKeVqLIWC3bXrl3D6tWr8dZbb9VYJjIyEiEhIViyZAk++ugjODg4YNWqVUhJSUFycjIAICUlBWq1Gm5ubhaP9fHxQUpKSrXHXbFiBZYvX950L4aIiJpVdmEp1h+IxcYDccjXGQAA4V4OeHx0B0zp4Q8lAx1RtRr8l/Hiiy9CJpPV+nXs2DGLxyQlJWH8+PH461//igcffLDGY6tUKnz33Xe4fPky3N3dodVqsXv3bkyYMKHOvnNCiBo7yC5ZsgS5ubnmr/j4+Ia+bCIiagGZBTq8/stFDHv9D6z+4yrydQZ08nHEe//XGzufGoG7egcy1BHVosEtdo899hjuueeeWsuEhoaabyclJWHUqFEYPHgwPv744zqP37dvX5w6dQq5ubkoLS2Fl5cXBg4ciH79+gEAfH19UVpaiuzsbItWu7S0NAwZMqTaY2o0Gmg0XCKGiMhapeWX4JO9Mfj8zxso1ktz0EX5OeOJ0R0wrqsv5HKObCWqjwYHO09PT3h6etarbGJiIkaNGoW+fftiw4YNkMvr/ynLxcUFAHDlyhUcO3YML7/8MgAp+KlUKuzcuRMzZswAACQnJ+PcuXNYuXJlA18NERG1ppTcEny45xq+PHLDvEJE9wAXPHF7R4yJ8uZUJUQN1Gx97JKSkjBy5EgEBwfjzTffRHp6unmfr6+v+XZkZCRWrFiBu+66CwCwefNmeHl5ITg4GGfPnsWTTz6JadOmmQdLuLi4YN68eXj66afh4eEBd3d3LF68GN27d8eYMWOa6+UQEVETSswpxtrdV/HN0QSUGqVA1zvYFU/c3hEjO3kx0BE1UrMFux07duDq1au4evUqAgMDLfZVnmHl0qVLyM3NNd9PTk7GokWLkJqaCj8/P8yaNQvPP/+8xePffvttKJVKzJgxA8XFxbj99tuxceNGzmFHRGTlbmQWYc3uq/juRAL0Rum9YECoO564vSOGdvBgoCO6RS06j5214Dx2REQtKzajEB/suoqtJxNhNElvO0MiPPDE7R0xKNyjlWtHZN2sdh47IiJqX66m5eP9P67if6eTUJbnMLyjJ564vSP6h9Y8rykRNQ6DHRERNbmLKXlY/cdVbD+bjPLrQqMjvfH46A7oHexW+4OJqNEY7IiIqMmcS8zF6j+u4NfzqeZtY7v44PHRHdE90KUVa0bUPjDYERHRLdEZjNh7OQNfHbmB3y+mAQBkMmBiNz88NroDovzYl5mopTDYERFRg5UaTNh/NR0/nU7GzuhU87Jfchkwpac/HhvVAR19nFq5lkTtD4MdERHVi95owoGrGfjpTDJ2nE9BXonBvM/X2Q4Tu/vhvkHBiPBybMVaErVvDHZERFQjg9GEQzGZ+Ol0Mn6NTkFOkd68z8tJg0nd/TC5hx/6BLtx2S8iK8BgR0REFgxGE47EZuHHM8n49XwKsgpLzfs8HdWY0E0Kc/1C3aFgmCOyKgx2REQEo0ngSGwWtp1Nwi/nUpBRUBHm3B3UmNDNF5N6+GFgmAfDHJEVY7AjImqnTCaBY9ezse1MErafS0F6vs68z1WrksJcd38MCneHUiFvxZoSUX0x2BERtSMmk8DJ+Gz8dCYZ288mIzWvIsw52ykxvpsvJvXwx5AID6gY5ojaHAY7IiIbJ4TAqfgcbCsLc0m5JeZ9TnZKjO3ii8k9/DC0gyfUSoY5oraMwY6IyAYJIXA2MRfbziTjpzPJSMwpNu9z1ChxRxcfTOruh+GdPKFRKlqxpkTUlBjsiIhshBAC55PysO1sMradScaNrCLzPq1agTFRPpjcww+3dfKCnYphjsgWMdgREbVBQgik5ulwPikX55PycD4pF+cS8yxa5uxVCoyO8saUHn4Y2dmbYY6oHWCwIyKyciaTQGxmoTnARSflITopD5mV5pcrp1HKMTrSG5N7+GNUpBe0av6bJ2pP+BdPRGRFdAYjLqcUSAEuOQ/nk/JwITkPRaXGKmUVchkivBzQ1d8FXf2d0cXfGT0DXeGg4b92ovaKf/1ERK0kr0SPC0l5ZS1xUmvc1bQCGEyiSlk7lRyRvs7o6u9sDnKdfZ14eZWILDDYERG1gLS8EnN4Kw9ylQc3VOaqVZkDXBc/KcyFeTpwkmAiqhODHRFREzKZBK5nFVkEuOikPGQU6Kot7+9ihy5lLXBd/Z3RNcAF/i52kMm4bBcRNRyDHRHRLcgqLMWJ69k4cSMbx69n43xSHgp0hirl5DIg3MvRHOC6+Elhzs1B3Qq1JiJbxWBHRFRPRpPApZR8nLghBbkT17MRl1n1cqpaKUeUrxO6+DubW+OifJ1hr2Z/OCJqXs0W7OLi4vDyyy/jjz/+QEpKCvz9/fG3v/0NS5cuhVpd8yfU1NRUPPPMM9ixYwdycnJw2223YfXq1ejYsaO5zMiRI7Fnzx6Lx82cORNfffVVc70cImqHcopKcfJGjjnInbqRg8JqRqdGeDmgT7Ab+oS4oXewKzp4ObI/HBG1imYLdhcvXoTJZMJHH32EDh064Ny5c3jooYdQWFiIN998s9rHCCEwbdo0qFQq/PDDD3B2dsaqVaswZswYREdHw8HBwVz2oYcewksvvWS+b29v31wvhYjaAZNJ4Epagbkl7viNbMSkF1Yp56BWoFewK/oGu6F3iBt6B7nCVcvLqURkHZot2I0fPx7jx4833w8PD8elS5ewdu3aGoPdlStX8Oeff+LcuXPo2rUrAGDNmjXw9vbGl19+iQcffNBcVqvVwtfXt7mqT9QoV1Lzse1sMk7cyIGjRgFXrRruWjVctSq4O6jhplXDzUENN60Kbg5qOGmU7CTfSnKL9TgVn2PuH3fqRg7yq+kbF+7pgN7BbugT4oo+wW7o5OMEhZznjIisU4v2scvNzYW7u3uN+3U6adSYnZ2deZtCoYBarcb+/fstgt2mTZvw+eefw8fHBxMmTMCyZcvg5ORU43HLjw0AeXl5t/pSiMwup+Zj25lkbD+bjCtpBQ16rFIuk8Kfg8ocAt0cVFIAvCkEupXtd7JTQs5g0SAmk0BMRgFOXM8xD3K4ml4AcdN0cVq1Aj0DXc0hrnewG9w5uIGI2pAWC3bXrl3D6tWr8dZbb9VYJjIyEiEhIViyZAk++ugjODg4YNWqVUhJSUFycrK53H333YewsDD4+vri3LlzWLJkCU6fPo2dO3dWe9wVK1Zg+fLlTf6aqP2qKcypFDIM7+iFUZHeMJkEsgpLkVNUiqwivfS9sBQ5RXpkFZaiWG+EwSSQUaCrcSqM6shlsAx91YRAV3sVXOxVcK703UGtaDetg/klepyOzzX3jTt5Iwe5xfoq5UI8tFLfuGBX9A52Q6SvE/vGEVGbJhPi5s+stXvxxRfrDElHjx5Fv379zPeTkpIwYsQIjBgxAv/5z39qfezx48cxb948nD59GgqFAmPGjIFcLv2j3b59e42P6devH44fP44+ffpU2V9di11QUBByc3Ph7Oxca32IAKn/5+XUAmw7K4W5q5XCnFohx/COnpjY3Q9juvjAxV5Vr2OW6I3ILipFdqFe+l5UiuzCUmRVvl+kR3Zhxb7qOu7Xl0Iug7OdEs72KjjblQc+pfTdTgp/0j6lOQxWLqdRNv+ITiEESvQmFJYaUKQzSt9LDSjUGS2/lxpRpCv7ftP+zEIdrqRVbY2zU8nRI9DVHOT6hLjB01HT7K+JiOhW5eXlwcXFpV65pcHBLiMjAxkZGbWWCQ0NNV9OTUpKwqhRozBw4EBs3LjRHNLqkpubi9LSUnh5eWHgwIHo168fPvjgg2rLCiGg0Wjw2WefYebMmXUeuyE/IGq/hBC4lJqP7WeSse1sMq5V6kivVshxW6eKMOdsV78wd6t0BiNyiqTgV7n1T2oNtAyE+cV65JXokVush97YoD/zammUcstWwLKQWBEMK0KiUYhKwcyIQt1N3y2Cm+X2hv1Hqlmgmz36BLuhb4gb+gS7IdLPCSq2xhFRG9SQ3NLgS7Genp7w9PSsV9nExESMGjUKffv2xYYNG+od6gDAxcUFgDSg4tixY3j55ZdrLHv+/Hno9Xr4+fnV+/hE1RFC4GJKPraflcJcjBWEuco0SgV8nBXwcbaru3CZ8law8pCXVynw5RUbat5WIm3P10lhS2cwIS1fh7T8+l82vhVatQJatRIOmrLvagW0mrLvNW3XKOFkp0RXf2d4O9X/Z0REZCsa3GJXX+WXX4ODg/Hpp59Coai4jFN5NGtkZCRWrFiBu+66CwCwefNmeHl5ITg4GGfPnsWTTz6Jvn374rvvvgMg9dXbtGkTJk6cCE9PT0RHR+Ppp5+Gvb09jh49avE8NWGLHVVWHubK+8zFZFQKc0o5buvohUk9fHF7VOuEudZmNAkU6AzIK7YMhpUDYMV2aZtCLqsziDlqlNCqFXAo/65WQquRvturFBwgQkRUpllb7Oprx44duHr1Kq5evYrAwECLfZWz5KVLl5Cbm2u+n5ycjEWLFiE1NRV+fn6YNWsWnn/+efN+tVqN33//He+++y4KCgoQFBSESZMmYdmyZfUKdUSA9Dt4IVlqmasuzI3o5IVJ3f1we5Q3nNphmKtMIZfBpeySa1BrV4aIiGrVbC121owtdu2TEALRyXllYS4FsTeFuZGdvDCphx9GRzLMERGR9bCKFjsiayCEwPmkPHPLXOV1PTVKOUZ29sLE7n64PcoHjhr+ORARUdvGdzKySeeTcvHTmWT8XE2YG9XZGxPLWuYY5oiIyJbwXY1sht5owi/nUrDhQCxO3Mgxb9co5Rgd6Y2J3aUw58AwR0RENorvcNTmZReW4sujN/DZoetIzi0BIK0AMSbKh2GOiIjaFb7bUZt1OTUfGw7EYuvJRJToTQAAT0c17hsYgvsGBXMeMyIiancY7KhNMZkEdl1Kw4YDcdh/tWIFlK7+zpg7NAxTevq1yNJXRERE1ojBjtqEAp0Bm4/F478H48yDIeQyYFxXX8wdGob+oW7tZoF7IiKimjDYkVW7kVmEjQfjsPlYPPJ1BgCAs50S9wwIxv2DQhDkrm3lGhIREVkPBjuyOkIIHLqWifUH4vD7xVTzovARXg6YMzQMf+kTAK2av7pEREQ347sjWY0SvRE/nErEhgNxuJiSb94+opMX5g4NxW0dvbh+KBERUS0Y7KjVpeSW4LM/4/DF4RvILtIDAOxVCkzvG4jZQ0LRwduxlWtIRETUNjDYUas5eSMbGw7EYfvZZBhM0vXWAFd7zB4Sgpn9guGi5XqtREREDcFgRy1KbzTh57LVIU5WWh1iQKg7HhgWijFRPlAq5K1XQSIiojaMwY5aRFZhKb48Iq0OkZInrQ6hVsgxpac/5g4NRbcAl1auIRERUdvHYEfN6lJKxeoQOkP56hAa/G1QMO4bGAIvJ00r15CIiMh2MNhRk8st1mPv5XR8eeQGDl7LNG/vFuCMB4aGYVIPrg5BRETUHBjs6JYJIXApNR+7LqZj16U0HL+eDWPZYAi5DBjfTVodol8IV4cgIiJqTgx21CiFOgMOXsvErktp2H0xDUm5JRb7I7wcMLarL+4bGIxAN64OQURE1BIY7KjeYtILsOtSOnZfSsPhmCyUGk3mfRqlHEMiPDAq0hsjO3kj2INhjoiIqKUx2FGNSvRGHI7Nwq6Ladh9KQ1xmUUW+4Pc7TG6szdGRnpjcLgH7FTsN0dERNSaGOzIQkJ2kdQqdzENB65loERf0SqnUsgwIMwdozp7Y2Rnb0R4ObDPHBERkRVhsGvn9EYTjsVlY9elNOy6mIYraQUW+32d7TAq0gsjO3tjaAdPOGr4K0NERGStmvVdeurUqTh16hTS0tLg5uaGMWPG4PXXX4e/v3+NjxFCYPny5fj444+RnZ2NgQMH4oMPPkDXrl3NZXQ6HRYvXowvv/wSxcXFuP3227FmzRoEBgY258uxGWl5Jdh9SRrBuu9KBgp0BvM+hVyGvsFuGBnphVGdvRHp68RWOSIiojZCJoQQzXXwt99+G4MHD4afnx8SExOxePFiAMDBgwdrfMzrr7+OV199FRs3bkSnTp3wyiuvYO/evbh06RKcnJwAAI888gh+/PFHbNy4ER4eHnj66aeRlZWF48ePQ6Gou59XXl4eXFxckJubC2dn56Z5sVbMaBI4FZ9tno7kfFKexX4PBzVGdJaC3G0dvbhGKxERkRVpSG5p1mB3s//973+YNm0adDodVKqq4UEIAX9/fyxcuBDPPPMMAKl1zsfHB6+//joefvhh5ObmwsvLC5999hlmzpwJAEhKSkJQUBC2b9+OcePG1VkPWw52RpNAZoEOafk6XEnLx+5L6dhzOR05RXpzGZkM6BHoilFlYa57gAvkcrbKERERWaOG5JYW6zCVlZWFTZs2YciQIdWGOgCIjY1FSkoKxo4da96m0WgwYsQIHDx4EA8//DCOHz8OvV5vUcbf3x/dunXDwYMHqw12Op0OOp3OfD8vL69KGWtnMJqQWViK1LwSpOXpkJovfU/Lt7yfUaCDqZqo7mynxG2dvDA60hu3dfKCpyOX8iIiIrI1zR7snnnmGbz//vsoKirCoEGD8NNPP9VYNiUlBQDg4+Njsd3HxwfXr183l1Gr1XBzc6tSpvzxN1uxYgWWL19+Ky+j2eiNJmQU6MpCmk4Kbvk6pJV9L7+fWUNgq45cJq3H6u9qb55brneQK5QKefO+GCIiImpVDQ52L774Yp0h6ejRo+jXrx8A4B//+AfmzZuH69evY/ny5Zg1axZ++umnWjvk37xPCFFnB/7ayixZsgSLFi0y38/Ly0NQUFCtx7tVeqMJ6fnVhLWbWtsyC0tR34vhCrkMXo4aeDtr4O1kV/ZdAx9nO4vvHo4aKHhplYiIqN1pcLB77LHHcM8999RaJjQ01Hzb09MTnp6e6NSpE6KiohAUFIQ///wTgwcPrvI4X19fAFKrnJ+fn3l7WlqauRXP19cXpaWlyM7Otmi1S0tLw5AhQ6qtj0ajgUbTspcevz2egCVbztarrFIug5eTBt5lwezmsOZV9t3dQc3ARkRERDVqcLArD2qNUT5Oo3J/t8rCwsLg6+uLnTt3onfv3gCA0tJS7NmzB6+//joAoG/fvlCpVNi5cydmzJgBAEhOTsa5c+ewcuXKRtWrOXg7aaBSyODtJAWzKmHNWQOfslY3d62agxeIiIjoljVbH7sjR47gyJEjGDZsGNzc3BATE4MXXngBERERFq11kZGRWLFiBe666y7IZDIsXLgQr732Gjp27IiOHTvitddeg1arxb333gsAcHFxwbx58/D000/Dw8MD7u7uWLx4Mbp3744xY8Y018tpsJGdvXHp5QkMbERERNRimi3Y2dvbY8uWLVi2bBkKCwvh5+eH8ePH46uvvrK4LHrp0iXk5uaa7//zn/9EcXExFixYYJ6geMeOHeY57ABpfjylUokZM2aYJyjeuHFjveawaym8ZEpEREQtrUXnsbMWubm5cHV1RXx8vM3NY0dERES2pXzQZ05ODlxcXGot2y4X/szPzweAZh8ZS0RERNRU8vPz6wx27bLFzmQyISkpCU5OzbsOannCZsug9eA5sT48J9aJ58X68JxYp5Y4L0II5Ofnw9/fH3J57XPStssWO7lcjsDAwBZ7PmdnZ/4RWhmeE+vDc2KdeF6sD8+JdWru81JXS105LkVAREREZCMY7IiIiIhsBINdM9JoNFi2bFmLr3pBNeM5sT48J9aJ58X68JxYJ2s7L+1y8AQRERGRLWKLHREREZGNYLAjIiIishEMdkREREQ2gsGOiIiIyEYw2N2iNWvWICwsDHZ2dujbty/27dtXa/k9e/agb9++sLOzQ3h4OD788MMWqmn70ZBzsmXLFtxxxx3w8vKCs7MzBg8ejF9//bUFa9s+NPTvpNyBAwegVCrRq1ev5q1gO9XQ86LT6bB06VKEhIRAo9EgIiIC69evb6Hatg8NPSebNm1Cz549odVq4efnh7lz5yIzM7OFamv79u7diylTpsDf3x8ymQzff/99nY9p9fd5QY321VdfCZVKJT755BMRHR0tnnzySeHg4CCuX79ebfmYmBih1WrFk08+KaKjo8Unn3wiVCqV+Pbbb1u45raroefkySefFK+//ro4cuSIuHz5sliyZIlQqVTixIkTLVxz29XQc1IuJydHhIeHi7Fjx4qePXu2TGXbkcacl6lTp4qBAweKnTt3itjYWHH48GFx4MCBFqy1bWvoOdm3b5+Qy+Xi3XffFTExMWLfvn2ia9euYtq0aS1cc9u1fft2sXTpUvHdd98JAGLr1q21lreG93kGu1swYMAAMX/+fIttkZGR4tlnn622/D//+U8RGRlpse3hhx8WgwYNarY6tjcNPSfV6dKli1i+fHlTV63dauw5mTlzpnjuuefEsmXLGOyaQUPPy88//yxcXFxEZmZmS1SvXWroOXnjjTdEeHi4xbb33ntPBAYGNlsd27P6BDtreJ/npdhGKi0txfHjxzF27FiL7WPHjsXBgwerfcyhQ4eqlB83bhyOHTsGvV7fbHVtLxpzTm5mMpmQn58Pd3f35qhiu9PYc7JhwwZcu3YNy5Yta+4qtkuNOS//+9//0K9fP6xcuRIBAQHo1KkTFi9ejOLi4paoss1rzDkZMmQIEhISsH37dgghkJqaim+//RaTJk1qiSpTNazhfV7ZIs9igzIyMmA0GuHj42Ox3cfHBykpKdU+JiUlpdryBoMBGRkZ8PPza7b6tgeNOSc3e+utt1BYWIgZM2Y0RxXbncackytXruDZZ5/Fvn37oFTyX1RzaMx5iYmJwf79+2FnZ4etW7ciIyMDCxYsQFZWFvvZNYHGnJMhQ4Zg06ZNmDlzJkpKSmAwGDB16lSsXr26JapM1bCG93m22N0imUxmcV8IUWVbXeWr206N19BzUu7LL7/Eiy++iK+//hre3t7NVb12qb7nxGg04t5778Xy5cvRqVOnlqpeu9WQvxWTyQSZTIZNmzZhwIABmDhxIlatWoWNGzey1a4JNeScREdH44knnsALL7yA48eP45dffkFsbCzmz5/fElWlGrT2+zw/DjeSp6cnFApFlU9SaWlpVdJ6OV9f32rLK5VKeHh4NFtd24vGnJNyX3/9NebNm4fNmzdjzJgxzVnNdqWh5yQ/Px/Hjh3DyZMn8dhjjwGQAoUQAkqlEjt27MDo0aNbpO62rDF/K35+fggICICLi4t5W1RUFIQQSEhIQMeOHZu1zrauMedkxYoVGDp0KP7xj38AAHr06AEHBwcMHz4cr7zyCq8CtQJreJ9ni10jqdVq9O3bFzt37rTYvnPnTgwZMqTaxwwePLhK+R07dqBfv35QqVTNVtf2ojHnBJBa6ubMmYMvvviCfVOaWEPPibOzM86ePYtTp06Zv+bPn4/OnTvj1KlTGDhwYEtV3aY15m9l6NChSEpKQkFBgXnb5cuXIZfLERgY2Kz1bQ8ac06Kioogl1u+jSsUCgAVrUTUsqzifb7FhmnYoPKh6evWrRPR0dFi4cKFwsHBQcTFxQkhhHj22WfF/fffby5fPgz6qaeeEtHR0WLdunWc7qSJNfScfPHFF0KpVIoPPvhAJCcnm79ycnJa6yXYnIaek5txVGzzaOh5yc/PF4GBgWL69Oni/PnzYs+ePaJjx47iwQcfbK2XYHMaek42bNgglEqlWLNmjbh27ZrYv3+/6NevnxgwYEBrvQSbk5+fL06ePClOnjwpAIhVq1aJkydPmqegscb3eQa7W/TBBx+IkJAQoVarRZ8+fcSePXvM+2bPni1GjBhhUX737t2id+/eQq1Wi9DQULF27doWrrHta8g5GTFihABQ5Wv27NktX3Eb1tC/k8oY7JpPQ8/LhQsXxJgxY4S9vb0IDAwUixYtEkVFRS1ca9vW0HPy3nvviS5dugh7e3vh5+cn7rvvPpGQkNDCtbZdu3btqvU9whrf52VCsL2WiIiIyBawjx0RERGRjWCwIyIiIrIRDHZERERENoLBjoiIiMhGMNgRERER2QgGOyIiIiIbwWBHREREZCMY7IiIiIhsBIMdERERkY1gsCMiIiKyEQx2RERERDaCwY6IiIjIRjDYEREREdkIBjsiIiIiG8FgR0RERGQjGOyIiIiIbASDHREREZGNYLAjonbt8OHDuOuuuxAcHAyNRgMfHx8MHjwYTz/9tLnMmjVrsHHjxlt6ntdeew3ff//9rVWWiKgOMiGEaO1KEBG1hm3btmHq1KkYOXIkHnroIfj5+SE5ORnHjh3DV199hYSEBABAt27d4Onpid27dzf6uRwdHTF9+vRbDohERLVRtnYFiIhay8qVKxEWFoZff/0VSmXFv8N77rkHK1eubMWaERE1Di/FElG7lZmZCU9PT4tQV04ul/49hoaG4vz589izZw9kMhlkMhlCQ0MBACUlJXj66afRq1cvuLi4wN3dHYMHD8YPP/xgcSyZTIbCwkL897//NR9j5MiR5v0pKSl4+OGHERgYCLVajbCwMCxfvhwGg6HZXjsR2Sa22BFRuzV48GD85z//wRNPPIH77rsPffr0gUqlsiizdetWTJ8+HS4uLlizZg0AQKPRAAB0Oh2ysrKwePFiBAQEoLS0FL/99hvuvvtubNiwAbNmzQIAHDp0CKNHj8aoUaPw/PPPAwCcnZ0BSKFuwIABkMvleOGFFxAREYFDhw7hlVdeQVxcHDZs2NBSPw4isgHsY0dE7VZmZiamTZuG/fv3AwBUKhX69++PKVOm4LHHHoOjoyOA+vexMxqNEEJg/vz5OHHiBE6cOGHeV1Mfu/nz52PTpk04f/48goODzdvfeustLF68GOfPn0eXLl2a5gUTkc3jpVgiarc8PDywb98+HD16FP/+979x55134vLly1iyZAm6d++OjIyMOo+xefNmDB06FI6OjlAqlVCpVFi3bh0uXLhQrzr89NNPGDVqFPz9/WEwGMxfEyZMAADs2bPnll4jEbUvDHZE1O7169cPzzzzDDZv3oykpCQ89dRTiIuLq3MAxZYtWzBjxgwEBATg888/x6FDh3D06FE88MADKCkpqddzp6am4scff4RKpbL46tq1KwDUK1wSEZVjHzsiokpUKhWWLVuGt99+G+fOnau17Oeff46wsDB8/fXXkMlk5u06na7ez+fp6YkePXrg1VdfrXa/v79/vY9FRMRgR0TtVnJyMvz8/KpsL7+MWh6qNBoNiouLq5STyWRQq9UWoS4lJaXKqNjajjF58mRs374dERERcHNza/RrISICOHiCiNqxHj16IDAwEFOmTEFkZCRMJhNOnTqFt956C/n5+Th48CC6d++OOXPm4KuvvsJ///tfhIeHw87ODt27d8eGDRvwwAMP4JFHHsH06dMRHx+Pl19+GXK5HFeuXEHlf68jR47EhQsX8J///Ad+fn5wcnJC586dkZycjMGDB8Pe3h5PPPEEOnfujJKSEsTFxWH79u348MMPERgY2Io/JSJqSxjsiKjd+uabb/DDDz/g6NGjSE5Ohk6ng5+fH0aMGIElS5YgKioKAHD9+nX8/e9/x6FDh5Cfn4+QkBDExcUBAF5//XV8+OGHSE5ORnh4OBYtWoSEhAQsX77cItidPn0ajz76KE6ePImioiKMGDHCPMo2IyMDL7/8Mn788UckJCTAyckJYWFhGD9+PJYsWQIHB4eW/tEQURvFYEdERERkIzgqloiIiMhGMNgRERER2QgGOyIiIiIbwWBHREREZCMY7IiIiIhsBIMdERERkY1olytPmEwmJCUlwcnJyWLGeCIiIiJrI4RAfn4+/P39IZfX3ibXLoNdUlISgoKCWrsaRERERPUWHx9f50o07TLYOTk5AZB+QM7Ozq1cGyIiIqKa5eXlISgoyJxfatMug1355VdnZ2cGOyIiImoT6tN9rE0GuxUrVmDLli24ePEi7O3tMWTIELz++uvo3Llza1eNiKjZCCFQajShuNSIolIjikoNZd8rbpsEoJDJoJCXfwHyyvfLbsvlMijlMot9cpm0rXy/QiaDXA4o5XLz7ZuPw37KRNalTQa7PXv24NFHH0X//v1hMBiwdOlSjB07FtHR0Vwsm4handEkUKwvC1s6KXgV6w0ovOl2cakRhaUGc1CrfLs8qJWXKb9tMFnX8t4yWUWQVMplcLZXwcVeBVdt2Xd7tXRbW3Hb1b7svlYNV3sVtGoFAyJRE5EJIazrv0QjpKenw9vbG3v27MFtt91WZ/m8vDy4uLggNzeXl2KJqF6EEMgt1iM1T4fUvBKk5pUgLb/idmqeDml5JcgqKkWJ3tTs9VEpZNCqldCqFbBXK+CgVsJerYBcBphMgFEIGE0VX6by+5W2m8z3AaPJVFYOVco1N6VcVhEEy8KeRRAs21d5v6tWBSc7FRRyBkKyfQ3JLW2yxe5mubm5AAB3d/cmPa7RaIRer2/SY1JVKpUKCoWitatB7ZQQAvk6A9LKwllqpe9p+WWBrex7qaFhgU0mA7QqBbSasgCmUsDhptv2aoVFmer2mW+rFeYwp1K03DSkppuCnlGUhcKbthuMAvklBuQUlyKnSI+cYj1yiyrdLtYjt0hfsb9Ij1KjCQaTQEZBKTIKSgEU1rteMhngbKcytwK6O6gR4GaPIDctgty1Zd/t4WKvYosgtRttPtgJIbBo0SIMGzYM3bp1q7aMTqeDTqcz38/Ly6vzmCkpKcjJyWnKqlItXF1d4evry3++1KSKSg2WLWzlt8ta2srDXLHeWO9jumlV8HG2g7ezHXycNPBxtoOPs0a672wHDwc1tGopjGmUcpv4nZbLZZBDBlUTf/4SQqBEbzIHvdxifdn3ijBYfr98X3m5Ap0BQkAKi8V6XK/leZw0SinwVQp75eEv0M0eDpo2/1ZIZNbmf5sfe+wxnDlzBvv376+xzIoVK7B8+fJ6H7M81Hl7e0Or1drEP2ZrJYRAUVER0tLSAAB+fn6tXCNqS4QQSMkrwdW0AvNXXGYhUnKlEJevM9T7WE52SnNI83EqC27OlYKbkx28nDSwa+p0047JZDLYqxWwV9vDz8W+QY/VG00WQTC3WI+0PB0SsosRn12E+KwixGcXIz1f+j24mJKPiyn51R7L3UGNIDd7BJYFv0BzCLRHgJs9NEqec2o72nQfu8cffxzff/899u7di7CwsBrLVddiFxQUVO21aqPRiMuXL8Pb2xseHh7NVneylJmZibS0NHTq1ImXZakKg9GE61lF5vB2La0A19ILcC29EAV1hDd7lQK+LnbwrtS6dnOLm7ezBlp1m/+cS9Uo0RvNYS+hLOxJoa8I8VnFyC2uvbuNTAb4ONkhyN0egW5aiwAY5C4FUvbzo+Zm833shBB4/PHHsXXrVuzevbvWUAcAGo0GGo2mXscu71On1WpvuZ5Uf+U/b71ez2DXjhWXGssCW0GVVji9sfrPoAq5DCHuWkR4O6KDtyMivBzh71LR4uaoUbLVvR2zUynQoex3ozp5JXokZFW08iXcFPyK9Uak5JUgJa8ER+OyqzxeKZfBz9UOQW5aRHg5oqu/M7r4O6OTjxNbd6lVtMlg9+ijj+KLL77ADz/8ACcnJ6SkpAAAXFxcYG/fsOb8mvCNoGXx592+ZBeW4mql8FYe5BJzilHTNQQ7lRwRXtIbdIfy796OCPFwgFrZcgMJyLY426nQxV+FLv5VW0GEEMgqLK3SypeQLQXAxOxilBpNiM8qRnxWMQ5eyzQ/VimXoYO3I7r4O6Orv4s58DnbqVry5VE71CYvxdYUAjZs2IA5c+bU+fjamjRLSkoQGxuLsLAw2NnZNUV1qR74c7c9Qggk51bq/5ZecRk1s7C0xse5aVXm0Bbh5Si1xHk5IsDVHnJe8iIrYjIJpOaXID6rGDeyinA5NR/nk3JxPikPOUXVX+INdteiq7+zOeh19XeBt5OGH26pVu3iUiwRWQ+D0YRLqfk4eSMHp+JzcDk1H9fSClBYWvNo0wBXe0R4OyLCy8GiFc7DsX7dJoham1wug5+L1M9uQFjFdFvlH2rOJ+WZg150Uh4Sc6QAeCOrCD+fSzGX93RUo0tZq17XsrAX4q7lBxlqlDYZ7Kiquj7tzZ49Gxs3bmyZytwkNDQUCxcuxMKFC1vl+anppeWX4NSNHJy4kYOTN7JxJiG32ilDlHIZQjy05hY4KcA5IdzLgVNMkM2SyWTwd7WHv6s97ujiY96eXViKC8l5FoHvWnoBMgpKsfdyOvZeTjeXdVArEOVXEfTK++2x2wHVhf9ZbURycrL59tdff40XXngBly5dMm9raN/D0tJSqNXqJqsftV2lBhOik/Nw4no2TsZLQS4hu7hKOSeNEr2CXdE7yBVd/F3K+r9pW3QiXSJr5uagxpAOnhjSwdO8rbjUiIspUtiLLgt9F5PzUFhqxLHr2Th2vWLAhkohQ0dvp4qWvQAXRPk5w5EfkqgS/jbYCF9fX/NtFxcXyGQy87bMzEzMnz8f+/btQ1ZWFiIiIvCvf/0L//d//2d+zMiRI9GtWzeo1Wp8+umn6Nq1K/bs2YP//e9/ePrpp5GQkIBBgwZhzpw5mDNnDrKzs+Hq6goAOHjwIJ599lkcPXoUnp6euOuuu7BixQo4ODhg5MiRuH79Op566ik89dRTAHgp3ZqVX0I6cSMbJ8ta484l5VVZcUEmAzp5O6F3sCv6BLuhd7ArIrwceemIqIHs1Qr0DnZD72A38zaD0YSYjEKpVS+xooUvr8SA6GQpAG4+XnGMUA8tuvq7oH+oG4Z28EQHb0f22WvHGOwaorCWpW4UCqByp//aysrlQOUWtJrKOjg0rH41KCkpQd++ffHMM8/A2dkZ27Ztw/3334/w8HAMHDjQXO6///0vHnnkERw4cABCCMTFxWH69Ol48skn8eCDD+LkyZNYvHixxbHPnj2LcePG4eWXX8a6deuQnp6Oxx57DI899hg2bNiALVu2oGfPnvj73/+Ohx56qEleDzWd4lIjzibm4mR5kIvPRmqerko5N60KvYPd0CfYFb2D3dAj0AVOHN1H1CyUCjk6+Tihk48T7uotbRNCIDGnuCzk5SG67FJucm4J4jKLEJdZhG1npSs3Xk4aDInwwNAITwzp4IFAN07f1Z60yVGxt6rRo2Jr+wQ0cSKwbVvFfQcHoKio+rIjRgC7d1fc9/ICMjKqlmvkqdm4cSMWLlxY65JokyZNQlRUFN58800AUotdbm4uTp48aS7z7LPPYtu2bTh79qx523PPPYdXX33V3GI3a9Ys2Nvb46OPPjKX2b9/P0aMGIHCwkLY2dnVq48dR8U2PyEErmcW4WR8eWtcDi4k58Fw0yLvCrkMUX5O5pa43kFuCPHgCixE1iirsBTnk3JxJiEXh65l4mhcFnQ3tbAHu2sxtIMHhkR4YnCEBzw5QKnNsflRsdQwRqMR//73v/H1118jMTHRvBKHw00tgv369bO4f+nSJfTv399i24ABAyzuHz9+HFevXsWmTZvM24QQMJlMiI2NRVRUVBO/GqqvAp0Bp8v6xEmtcTnIqmaaES8nDfqYL6m6oXuAC+zVnFiVqC1wd1BjeEcvDO/ohUdHdYDOYMSJ6zk4eC0DB69l4lR8jjQS90gRvjwSDwCI9HXCkAhPDO3ggQFh7mx9tzEMdg1RUFDzvptXSyhb+7Ra8ps6k8fFNbpK9fHWW2/h7bffxjvvvIPu3bvDwcEBCxcuRGmp5Zv8zUFPCFGllebmBl6TyYSHH34YTzzxRJXnDQ4ObqJXQPVRVGrAoWuZ2HM5HUdis3ApNb9Ko69aIUfXAGf0DnJDnxDpsqq/ix1b44hshEapwOAIDwyO8MDTkD7gHYnNxIGrmTh4LRMXkvPM6+auPxALhVyGHoEu0mXbCA/0CXHjihltHINdQzSkz1tzlW2Effv24c4778Tf/vY3AFIYu3LlSp2taZGRkdi+fbvFtmPHjlnc79OnD86fP48OHTrUeBy1Wg2jseb5zKhxhBC4klaAPZfSsftyGo7GZqPUaHkJJsDV3mKAQxd/Zy5oTtSOOGqUGB3pg9GR0rQrmQU6HIqRQt7BqxmIyywyd814f9dVaJRy9At1w5CyoNc9wAVKjmxvUxjs2oEOHTrgu+++w8GDB+Hm5oZVq1YhJSWlzmD38MMPY9WqVXjmmWcwb948nDp1yjwXXnkLzzPPPINBgwbh0UcfxUMPPQQHBwdcuHABO3fuxOrVqwFI89jt3bsX99xzDzQaDTw9PWt6SqpDXokeB65kYM/ldOy5nI7k3BKL/YFu9hjZ2QvDOniiT7AbvJ3ZX5GIKng4ajC5hz8m9/AHACTmFOPgVemy7YGrGUjL1+HAVamFD5CmMRoY7iENxujgiU4+HHFr7Rjs2oHnn38esbGxGDduHLRaLf7+979j2rRpyM3NrfVxYWFh+Pbbb/H000/j3XffxeDBg7F06VI88sgj0Gikzrc9evTAnj17sHTpUgwfPhxCCERERGDmzJnm47z00kt4+OGHERERAZ1Ox+lOGsBkEohOzpOC3KV0HL+RDWOlwQ4apRyDwj0wopMXRnT2QrinA//pElG9Bbja46/9gvDXfkEQQuBaeoE55B26lom8EgN+u5CK3y6kApBWyRgc4YmhZUEvyJ0jbq0NR8VyrdgGefXVV/Hhhx8iPj6+SY/Ln3uFrMJS7LsiBbm9V9KRUWDZFzLcywEjOnlhZGdvDAxzZ38YImoWRpNAdFIeDlzLwIGrGTgal4USvWV3j0A3ewyN8MSwjp4YFenNyZKbCUfFUpNZs2YN+vfvDw8PDxw4cABvvPEGHnvssdaulk0xmgROJ+Rg9yXp8uqZhByLQQ9atQJDIjwxsrMXRnTy4idkImoRCrkM3QNd0D3QBfNHREBnMOLUjRwcKOufdyo+BwnZxfj6WDy+PhYPtVKOUZ29MLmHP26P8oZWzYjRGvhTp1pduXIFr7zyCrKyshAcHIynn34aS5Ysae1qtXlpeSXmfnL7rmQgt1hvsT/S1wkjyoJcvxB3rg9JRK1Oo1RgYLgHBoZ7YNEdnVCoM+BIXBYOXs3A7xfSEJNRiF/Pp+LX86mwVykwOsobk7v7YVSkN68stCBeiuWlWKtg6z93vdGE49ezsedyOnZfSseF5DyL/c52SgzvKAW52zp5wdfF9n4GRGS7hBC4kJyPn84k4aczybiRVTFBv4NagTFdfDCpux9u6+TFkNcIvBRLZAUSsovMgx4OXstEgc5gsb9HoEtZXzkv9Ax05ZQCRNRmyWQydPF3Rhd/Z/xjXGecS8wzh7zEnGL8cCoJP5xKgpNGiTu6+GByTz8M6+DFqxHNgMGOqAldTcvHD6eS8PO5FFxNs5zQ2sNBjds6Sa1ywzp6clkfIrJJMllF37xnJ0TiVHwOfjqTjG1nkpGSV4ItJxOx5WQinO2UGNfVF5N7+mNIhAdU/HDbJBjsamAymeouRE2mLf+8k3OL8eNp6dPo+aSKS6xyGdAn2M08FUk3fxfI5ZyKhIjaD5lMht5lyxUunRiFEzeypZB3Nhnp+TpsPp6AzccT4KZVYXw3X0zu4Y+BYe68gnEL2MfupmvV5asyKBQKeHl5Qa1Wc16wZiSEQGlpKdLT02E0GtGxY0fIb15yzQrlFJVi+9kU/HAqEUfissyjWJVyGUZ08sLUXv4Y2ckbLlquwUhEdDOjSeBIbBa2nU3Cz2dTkFlpHWtPR7U55PUPdYeCH4gb1MeOwa6aH1BpaSmSk5NRVFRUzaOpOWi1Wvj5+UGtVrd2VWpUXGrEbxdS8cOpJOy5nAa9seJPZ0CoO+7s7Y+J3fzg5mC9r4GIyNoYjCb8GVMW8s6lIKeoYpYAbycNJnb3w+QefugT7NZur3ow2NWhPj8gIQQMBgPXOG0BCoUCSqXSKltG9UYTDlzNwA+nkrDjfAoKSyt+H6L8nHFnL39M6emPAFf7VqwlEZFt0BtNOHgtEz+dTsKv51OQV1Ix6MzPxc4c8noFuVrle0ZzYbCrQ0N+QNT+CCFw4kY2fjiVhG1nki0uEQS62ePOXv64s1cAOvk4tWItiYhsW6nBhP1X0/HT6WTsiE61mFkgwNUek3v4YXIPf3QLcLb5kMdgVwcGO6rO5dR8fH8yEf87nYSE7GLzdg8HNSb38MPUXgHoE9y+PiUSEVmDEr0Rey+n46czyfjtQiqKKl09CfHQ4s5eAbh3QLDNzgHKYFcHBjsql5BdhB9PJ+OHU4m4mJJv3u6gVmBcV1/c2TsAQyM8OEKLiMhKFJcasftSGn46k4zfL6aa169VymUY180XsweHon+om019CGewqwODXfuWVViK7WelMHc0Ltu8XaWQYWRnb9zZyx+3R/rAXs3Z0YmIrFmhzoDfLqRi0583cCQuy7w9ys8ZsweH4M5eATbxv5zBrg4Mdu1PUakBO6OlEa17L6fDYJJ+7WUyYFCYB+7s5Y8J3fw4PQkRURsVnZSHTw/F4ftTieZWPBd7FWb2D8LfBoYg2EPbyjVsPAa7OjDYtQ96own7rqSXjWhNRbG+ok9GtwBn3NkzAJN7+sHPhSNaiYhsRU5RKTYfS8Cnf8YhPkvqLy2TAbdHemPW4FAM6+DZ5qZNYbCrA4Odbbucmo9PD8Vh25lkZFeaDynEQ4s7e/pjaq8AdPB2bMUaEhFRczOaBHZfSsPGg3HYdyXDvD3c0wGzBofgL30D4WTXNq7SMNjVgcHO9phMArsvp2H9/jjsv1rxB+zpqMGUnn64s1cAega62FRnWiIiqp9r6QX47NB1fHs8wTxtioNagb/0DcSswSHo4G3d01cx2NWBwc52FOgM+O54AjYejENsRiEAaY3WcV19ce/AYAwO54hWIiKSFOgM2HIiAf89GIdr6YXm7UM7eGD24FDcHuVjlUuYMdjVgcGu7YvPKsJ/D8bh66PxyC/79OVkp8T/DQjG/YNCEOTedjvJEhFR8xJC4MDVTPz3UBx+v5CKsvF0CHC1x/2DQzCzX5BVLQ/JYFcHBru2SQiBw7FZ2HAgFjujK/4Qw70cMHdIKO7uEwgHjbJ1K0lERG1KfFYRPj98HV8fjTevU6tRyjG1pz9mDwlFtwCXVq5hOwh2e/fuxRtvvIHjx48jOTkZW7duxbRp0+r9eAa7tqVEb8SPp5Ow/kAcLiTnmbff1skLc4eGYkRHrzY3womIiKxLid6I/51KwsaDcYiu9F7TN8QNs4eEYnxXX6iVrdO1pyG5pU02bxQWFqJnz56YO3cu/vKXv7R2daiZpOWX4PM/b2DTn9fN67XaqeT4S59AzB0aavWdXYmIqO2wUykwo38Q/tovEMevZ+O/h67j57PJOH49G8evZ8PbSYN7Bwbj3gHB8Ha23qXL2mSLXWUymYwtdjbmbEIuNhyIxY9nkqA3Sr+efi52mDU4FP83IAiuWuvp90BERLYrNa8EXxy+gS+O3EB6vg6AtHTZxO5+mD0kBH2CW2bpMptvsSPbYzCasDM6FesPxFos89U3xA1zh4ZiXFdfqDi6lYiIWpCPsx2euqMTHh3VAT+fS8anh67j+PVs/O90Ev53Ogl9vDSY29kJd4zrBzuVdSxd1i6CnU6ng06nM9/Py8urpTS1pNwiPb46egOfHrqOxBxphnClXIbJPfwwd2gYega5tm4FiYiofcnLA/bsAdLSgPR0ID0d6vR03JmWhjvT05F05wy8HTUeP5xOQt6FK9Bt/A53XnkBvywcbhVzpbaLYLdixQosX768tatBlVxNK8DGg7H47niieakvdwc17hsYjL8NCoGPFfdfICIiKycEYDQCyrKYk5kJfP+9FNQqBTbz7QULgCVLpLKJicDUqTUe2n/IELzxQk8smRiFH347De3J7zGhu69VhDqgnQS7JUuWYNGiReb7eXl5CAoKasUatU9CCOy9koH1+2Ox53K6eXukrxMeGBqGqb38raYpm4iIrIgQQG4uIJcD5X3MUlOB9esrwlnlwJaeLgW1F1+UyqanAw8+WPPxExMrbvv4AP37A15e1X9FRgKQGiPm3tkfhsk/Y5TReoYrtItgp9FooNFoWrsa7VZRqQFbTiRiw4FY80zf0oLMPnhgWCgGh3tYzScdIiJqAeVBLS1NCmq+vtL2hATgrbeqD2t6PbBsWUVYy84G/vWvmp8jvaIBAb6+wKRJFeHM29vydnBwRVl3d+DIkXq/FKVCDqUVtUm0yWBXUFCAq1evmu/Hxsbi1KlTcHd3R3Dlk0OtKjGnGJ8eisOXh28gr0RaHcJRo8SMfkGYPSQEIR4OrVxDIiJqEiYTkJNTEcL8/YHwcGlfbCywdKllSMvIkIIaALzwAlDeXaqwEHjnnZqfJ7ticB38/IC5cy1b024ObOVcXYGffmrCF2y92mSwO3bsGEaNGmW+X36Zdfbs2di4cWMr1YrKxWYU4u2dl7HtbDKMZctDhHhoMWdIKKb3DYSTnaqVa0hERLUqD2qVw1inTkC3btL+K1eARx6paFXLyAAMhorHP/888NJL0m29Hvjyy+qfx8lJeq5y/v7AM8/UHNbs7SvKurhIl2LJQpsMdiNHjkQbn37PJqXkluDd36/gm2Px5kA3JMIDDwwNw6hIb6tcWJmIqF0wmaTWrpsvb/bpAwwYIJWJjgZmzqxoUTMaLY/x3HMVwQ4Afv+96vM4OUlBzNGxYltAgHR59eaQ5uUF2NlVffy//900r7mdapPBjqxLdmEpPtxzDRsPxkFnkD55jY70xqI7OlnFGntERDbHZAKysqoOGEhPB267DRg5Uip3+jQwdqw0KvTmoAZIl0jLg51aDZw7Z7nf2bkijPn7V2wPDAQ++6xqUKuuP7uDA1BpACM1LwY7arRCnQHr98fi470xyNdJTfD9Q93wz/GR6B/q3sq1IyJqQ0wmKXxVN7ozPR2YOFH6AoBjx4CBAy0vYVb2r39VBDtHR+l45Vxdqx3hCQAICgJ+/bWiZc3Ts/qgBkiXRP/2t1t91dQMGOyowXQGI748fAPv77qKjAJpDddIXyc8Mz4SIzt7cYQrEREAlJQA167VHNZmzACmT5fKHj0KDBpU87GcnSuCnZtbRahzda3aata/f8XjgoOBU6ek7Z6eUqtcTTQaqXWP2jQGO6o3o0ng+5OJePu3y0jIllaJCPHQYtEdnTClhz/k7ENHRLYuLw84frzqBLflXw88AMyaJZU9c0ZqWatJRERFsPPykr67uVUdNODtXdECBwAhIUBSkhTUVHUMRlOpgJ49G/1yqe1hsKM6CSGwMzoVb+64hMupBQAAbycNnri9I2b2D+IarkTUtmVmAr/9Vv2KBOnpwJNPAvPnS2UvXQJGj675WIMHV9z29pbmRLu5Ra08rFUOfaGhQGlp3UENkFZT8PNr1Esl28dgR7U6dC0TK3+9iJM3cgAALvYqPDIyArMHh8JebUUzMhJR+2YwSF/loyyTk4Fvv605rP3rX8DTT0tlY2OBe+6p+dhxcRW3fX2BqKjqg5qXF9CjR0XZ0FApNNaHXC59Ed0iBjuq1tmEXKz89SL2XckAANirFHhgWCj+flsEXOw5Dx0RNTO9Xgpg9vbS5UkAuHED+OST6sNadjawciXwj39IZRMTgSeeqPn4qakVt/39pZGkNbWsdexYUTYoSJoWhMhKMdiRhWvpBVi1Q5pcGABUChn+b0AwHhvdAd5OdnU8moioBqWl0txoaWlSWCqfOuPqVSmQ3RzUcnKk/a+/Dvzzn9Lt9HTglVdqfo7KS0gFBgJ//WvNLWsBARVl/f2BPXua9OUStRYGOwIAJOcW493frmDz8QQYTQIyGTCtVwCeGtMJwR7a1q4eEVmb0lLLMBYRIX0BUovWv/5l2bKWm1vx2H//W1pdAJAGI3zySfXPIZdLS0yVCw4GHn20+hUJvLyk/mzlfH2Bb75p2tdM1AYw2LVzWYWlWLv7Kv576DpKyyYXHhPljcXjOiPS17mVa0dELaakpOp0HL16Vaw0cOpUxRJS6elSIKtsxQrg2Wel2zod8MMPVZ9DLpdGclbuSxYaKq0TWlNQq1zWywt4//0mfNFEtofBrp0q0Bmwbl8sPtkXg4KyyYUHhLnjmfGd0TeEkwsTtXnlQe3m6TiGD6+Y5+zIEeD//k/anp9f9RivvWa5hNSff1ruVyikoObtLa3bWS48HFi7tuolUDe3qgME3N2lReCJqEkw2LUzOoMRm/68gQ92XUVmoTS5cFd/Z/xjXGeM6MTJhYmsVnGxZUCrHNgmTpQ6/wPA/v3AhAlAQUH1x3nttYpgp1YDMTEV+5TKiqDm5SX1UysXEQF8951lq5qra/UjOV1cKqYHIaIWxWDXThiMJmw5mYh3f7uCxBxpcuEwTwcsuqMTJnX34+TCRC1Nr5cmma1pVYKZM4Fx46Syu3bVPneas3NFsHNyqgh1SmXVy5uVl5Dq1AnYt6+ijKsrUNOHOycn4O67b/llE1HzYrCzcUII/Ho+BW/uuIyradI/ex9nDRaO6YTpfQM5uTBRUyooAC5cqDmsPfggMG2aVHbfPuD222s+VkRERbDz8JC+q1TVT8lReQmpyEhpEt3yy6O1tcJrtcCwYbf0konIujDY2bADVzOw8tdLOB2fAwBw1aqwYGQEZg0OhZ2KkwsT1UtmptS3rKaw9tRTUj81QFrvs7aWtUGDKoKdl5d0KbSmEZ6jRlU8rksXafoPZ+fagxogrffZqdOtvGIiasMY7GxQYk4xlv1wDr9dSAMAaNUKzBsWhoduC4ezHScXpnZICGnRdEXZB5rkZOCXX6oPa2lpUmf+efOksmfOAJMn13zsq1crbvv4SP3Sqgtr3t7AgAEVZbt1kwY41Kdfq1JpOTiBiKgGDHY2xGgS2HgwDm/tuISiUiNUChnuGxiCR0d1gJeTprWrR9R0hJAue8rlgIODtC0uDvjqq5oXZ1+1SpquA5AuVT7wQM3HT0ysuB0QAPTtW/0kt15elqNGu3QB4uPr9xo4UImImgGDnY04l5iLJVvO4myiNAlo/1A3vHZXd3T0cWrlmhHVgxDSdBvp6VIH/vI+ZZcuAR9/XP0l0JIS4IMPgAULpLLXrwNLltT8HJVXJQgOlkaO1rQqQXh4RdlOnYBjx5r8JRMRNQcGuzauUGfA2zsvY/2BWJgE4GynxJKJUZjZL4gjXan1CCFNYFsewoKDK5ZwOnPGcgmp8i+dTtr//vvS6gKAtJ7nqlU1P0/lBdZDQ4FZs6pe/iy/7etbUTY8HNi+vUlfMhGRNWCwa8P+uJiK578/b56+ZEpPfzw/OYprulLTE0JaEqryZc7u3SuWkDp2zHIJqYwMacmpcqtXA489Jt3OzQU2bar+ebRay8dFREiLutc0wKD8MiwAhIQA//1v075uIqI2hsGuDUrLK8HyH6Ox7WwyACDQzR4vT+uGUZ29W7lm1GYIIY2yvHnQwNChQNeuUpkDB6SWs/L9er3lMd57D3j8cel2SQmwc2fV53FwkAKYqtKgnU6dgDfeqD6oaW9alzggQGrdIyKiemGwa0NMJoEvjtzA679cRH6JAQq5DA8OC8OTYzpCq+apbNdMpoqgdnN/tAkTpM7/APDHH8B990ktagZD1eO8915FsAOA06ct9zs6VoSxyguuR0VJrWU3t6zZ21d9Dh8fYPHiW37JRERUFdNAG3EpJR//2noWx69nAwB6Brrgtbu7o6s/p0CwSeVBrXJIq3x75kypdQ0Afv4ZmDIFMBqrP5azc0Wws7MDUlIq9jk5WbaaBQVV7OvWTTp2eVDz9Kw+qAHSYIdZs275ZRMR0a1hsLNyJXojVv9xBR/tiYHBJOCgVuAf4zrj/sGhUHBwRNtSXCxNyVHdCM/0dOChhypWIti2DZg6teZjRURUBDtX14pQ5+xctdUsKqricT17AsePVwQ1u1r6Y7q4AOPH38orJiKiFsZgZ8UOXM3A0q1nEZdZBAAY28UHy+/sCj+XGlpNqOVlZ0ujPKub5DY9XbrkOGmSVPb336WWtZoMGlQR7Ly8pO8uLtVPyVF5Cak+faS507y8pFUHauPgIJUnIiKbxGBnhTILdHh12wVsOSlNkurrbIcXp3bF+G6+dTySmkRKirSO582XP8tvv/QS8Je/SGUPHKg9rE2dWhHsvL2l1rXqViTw8gJGjqx4XL9+0oCEuoIaIJUJDGzsqyUiIhvCYGdFhBD49ngCXtt+AdlFeshkwKxBIVg8rjOcuBRYwxgM0sjP8tGYcXHS5c2awto770iDCgDpUuWMGTUf+/r1itv+/tIoz+pGeHp5WS4hNWCA1MJXH0ql9EVERNQAfOewEjHpBVi69RwOxUgTrkb6OmHF3d3RO9itlWtmJQwGaSSno6P0BQAXLwJffll9WMvKAjZurOjQHx1dMY9adVJTK24HBwPDhtW83mflPmt9+kirIxAREVkBBrtWVmow4cM91/D+rqsoNZhgp5Jj4ZhOmDcsDCqFvLWr13zKg1pamtTq5ekpbT99Gvjww+qDGiCFtdmzpduxsdJl0ZpUXkIqPFy6fFrTep/BwRVlu3eXLsUSERG1MQx2rehoXBaWbDmLq2kFAIDbOnnhlTu7IdhDW8cjrZBeLwW18jAWFSUFNgA4fBh4/XXLlrXKlyQ3bADmzJFuJydLwa46Mpm0TFW5jh2lRd1rCmvl640CQGQk8O23TfqSiYiIrA2DXSvILdLj379cwJdH4gEAno5qPD+5C6b29IdMZiVTmJSWWga1tDRp1Gb54uh79gBLl1YEtZwcy8dXDmvZ2cDWrVWfQyaTwlfl+deiooBly6ofCerubtnvrEMHYM2apnzVREREbVqbDnZr1qzBG2+8geTkZHTt2hXvvPMOhg8f3trVqpEQAj+eScZLP0Yjo0Ba8Pye/kF4dkIkXLXq5n3y0tLqp+MYN66iz9iOHVI/tLQ0aT3Pm61fXxHsSkqkEaGVyeXSJdWbp93o0QP44IOqLWru7oBCYXmMkBDgxReb7GUTERG1J2022H399ddYuHAh1qxZg6FDh+Kjjz7ChAkTEB0djeDK/aWsRHxWEZ77/hz2XJb6fUV4OeC1u7pjYLhHHY+sgU5XdYLb8rA2YwbQq5dU7n//A+6/3/ISZmXr1lkOBrhypeJ25aB286XNPn2AzZstw5q7u/SYm/n7AwsWNO51EhERUb3JhBCitSvRGAMHDkSfPn2wdu1a87aoqChMmzYNK1asqPWxeXl5cHFxQW5uLpydnZu1ngajCesPxGLVzsso0ZugVsjx6KgOmD8yHBplpdaq0lJpZGZ1y0eVr0owcKBU9ptvpCWlavKf/wDz5km3d+4Exo6VbisUUlCr3Go2bx5wxx3S/uxs4Oz/t3e3MTX/fRzA39XpTpe6Io6TY5GVxIS6oszf/kZ2sZoHpo0ZxubMTDRcGZPM1jB2ydxMc/PgKprIPAidB6QbY5LNHFum3LSklfV3iFCf64H/KUeHOqfOr9PvvF/bedDX99f5nL13/D59f3ePu5u14GDbjRoREREpxp6+ZUiu2H358gVVVVXIyMiwGk9KSkJlZeUgVdXTY9Mr5J4twV+vGvDvtr8ww/8rknUaBF8oBHL+firBH398n3z5MrBixa9/2ezZ3Y2d5eHrGk33itqPzVpUVPd2CQnA06fdN8f9XaMWHNxdDxEREQ05Q7Kxa25uRkdHB7RardW4VqtF448POP9be3s72tvbu35+/6vDkgPofEUdHh48iZxrh349KTm5u5EaPfp7o/are6dZHuIOAHPnAi0tvTdqwPd7vv3Y6BEREZFqDcnGzuLnK0hFxOZVpdnZ2cjKylKqLADAvyaMwM3hI9EaPAr/GBcKjXZ0z6s8f1wd+/PP74dj+3JVrK9v3x41RURERG5lSDZ2ISEh8PLy6rE619TU1GMVDwB27tyJ9PT0rp/fv3+PcePGObXGKaFByP7vJvzzf//p2wY8l42IiIj6aUh2Ez4+PoiNjYXRaLQaNxqNSExM7DHf19cXgYGBVi8ljA8JUOR9iIiIiIAhumIHAOnp6Vi1ahXi4uKQkJCA06dP49WrVzAYDL1ua7kQWIlz7YiIiIj6w9Kv9OVGJkO2sUtNTUVLSwv27duHN2/eYOrUqSguLkZYWFiv25rNZgBw+uFYIiIiooFiNpsRFBT02zlD9j52/dHZ2YmGhgYMHz7cqY/wspzL9/r1a8UO/9LvMRPXw0xcE3NxPczENSmRi4jAbDYjNDQUnr2ckz9kV+z6w9PTE3q9XrH3U/K8PuobZuJ6mIlrYi6uh5m4Jmfn0ttKncWQvHiCiIiIiHpiY0dERESkEmzsnMjX1xeZmZnw5c2EXQYzcT3MxDUxF9fDTFyTq+XilhdPEBEREakRV+yIiIiIVIKNHREREZFKsLEjIiIiUgk2dkREREQqwcaun06cOIEJEybAz88PsbGxKCsr++380tJSxMbGws/PD+Hh4Th16pRClboPezK5cuUKFi5ciFGjRiEwMBAJCQm4efOmgtW6B3u/JxYVFRXQaDSYPn26cwt0U/bm0t7ejl27diEsLAy+vr6YOHEizp49q1C17sHeTPLy8hATE4Nhw4ZBp9Nh7dq1aGlpUaha9btz5w6Sk5MRGhoKDw8PXL16tddtBn0/L+Swixcvire3t+Tm5orJZJK0tDQJCAiQly9f2pxfW1srw4YNk7S0NDGZTJKbmyve3t5SWFiocOXqZW8maWlpcuDAAbl//77U1NTIzp07xdvbWx4+fKhw5eplbyYWra2tEh4eLklJSRITE6NMsW7EkVxSUlJk1qxZYjQapa6uTu7duycVFRUKVq1u9mZSVlYmnp6ecvToUamtrZWysjKZMmWKLF26VOHK1au4uFh27dolly9fFgBSVFT02/musJ9nY9cP8fHxYjAYrMaioqIkIyPD5vwdO3ZIVFSU1diGDRtk9uzZTqvR3dibiS3R0dGSlZU10KW5LUczSU1Nld27d0tmZiYbOyewN5fr169LUFCQtLS0KFGeW7I3k0OHDkl4eLjVWE5Ojuj1eqfV6M760ti5wn6eh2Id9OXLF1RVVSEpKclqPCkpCZWVlTa3uXv3bo/5ixYtwoMHD/D161en1eouHMnkZ52dnTCbzRgxYoQzSnQ7jmZy7tw5PH/+HJmZmc4u0S05ksu1a9cQFxeHgwcPYuzYsYiMjMS2bdvw6dMnJUpWPUcySUxMRH19PYqLiyEiePv2LQoLC7FkyRIlSiYbXGE/r1HkXVSoubkZHR0d0Gq1VuNarRaNjY02t2lsbLQ5/9u3b2huboZOp3Nave7AkUx+dvjwYXz8+BHLly93Rolux5FMnj17hoyMDJSVlUGj4X9RzuBILrW1tSgvL4efnx+KiorQ3NyMjRs34t27dzzPbgA4kkliYiLy8vKQmpqKz58/49u3b0hJScGxY8eUKJlscIX9PFfs+snDw8PqZxHpMdbbfFvj5Dh7M7G4cOEC9u7di4KCAowePdpZ5bmlvmbS0dGBFStWICsrC5GRkUqV57bs+a50dnbCw8MDeXl5iI+Px+LFi3HkyBGcP3+eq3YDyJ5MTCYTNm/ejD179qCqqgo3btxAXV0dDAaDEqXSLwz2fp5/DjsoJCQEXl5ePf6Sampq6tGtW4wZM8bmfI1Gg5EjRzqtVnfhSCYWBQUFWLduHS5duoQFCxY4s0y3Ym8mZrMZDx48QHV1NTZt2gTge0MhItBoNCgpKcH8+fMVqV3NHPmu6HQ6jB07FkFBQV1jkydPhoigvr4eERERTq1Z7RzJJDs7G3PmzMH27dsBANOmTUNAQADmzp2L/fv38yjQIHCF/TxX7Bzk4+OD2NhYGI1Gq3Gj0YjExESb2yQkJPSYX1JSgri4OHh7ezutVnfhSCbA95W6NWvWID8/n+emDDB7MwkMDMTjx4/x6NGjrpfBYMCkSZPw6NEjzJo1S6nSVc2R78qcOXPQ0NCADx8+dI3V1NTA09MTer3eqfW6A0cyaWtrg6en9W7cy8sLQPcqESnLJfbzil2moUKWS9PPnDkjJpNJtmzZIgEBAfLixQsREcnIyJBVq1Z1zbdcBr1161YxmUxy5swZ3u5kgNmbSX5+vmg0Gjl+/Li8efOm69Xa2jpYH0F17M3kZ7wq1jnszcVsNoter5dly5bJkydPpLS0VCIiImT9+vWD9RFUx95Mzp07JxqNRk6cOCHPnz+X8vJyiYuLk/j4+MH6CKpjNpulurpaqqurBYAcOXJEqquru25B44r7eTZ2/XT8+HEJCwsTHx8fmTlzppSWlnb92+rVq2XevHlW82/fvi0zZswQHx8fGT9+vJw8eVLhitXPnkzmzZsnAHq8Vq9erXzhKmbv9+RHbOycx95cnj59KgsWLBB/f3/R6/WSnp4ubW1tCletbvZmkpOTI9HR0eLv7y86nU5Wrlwp9fX1CletXrdu3frtPsIV9/MeIlyvJSIiIlIDnmNHREREpBJs7IiIiIhUgo0dERERkUqwsSMiIiJSCTZ2RERERCrBxo6IiIhIJdjYEREREakEGzsiIiIilWBjR0RERKQSbOyIiIiIVIKNHREREZFKsLEjIiIiUon/A7AgO6ddok11AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "# Initialize vectors for the state, parameter, adjoint (not used) and control variables\n", + "x = control_model.generate_vector()\n", + "\n", + "# Initialize the noise vector \n", + "noise = dl.Vector(comm_mesh)\n", + "prior.init_vector(noise, \"noise\")\n", + "\n", + "# Use hippylib's rng to sample Gaussian white noise \n", + "rng = hp.Random(seed=111)\n", + "\n", + "# This is sampling noise with 1.0 standard dev to the noise vector\n", + "rng.normal(1.0, noise) \n", + "\n", + "# The prior then turns the noise into a parameter sample\n", + "prior.sample(noise, x[soupy.PARAMETER])\n", + "\n", + "# Also set the CONTROL component of x to the optimal control z\n", + "x[soupy.CONTROL].set_local(z_optimal_np)\n", + "\n", + "# Solve the forward problem\n", + "control_model.solveFwd(x[soupy.STATE], x)\n", + "\n", + "# Convert to functions and plot\n", + "u_fun = dl.Function(Vh[soupy.STATE], x[soupy.STATE])\n", + "m_fun = dl.Function(Vh[soupy.PARAMETER], x[soupy.PARAMETER])\n", + "z_fun = dl.Function(Vh[soupy.CONTROL], x[soupy.CONTROL])\n", + "\n", + "x_obs = np.linspace(0, 1, 100)\n", + "if comm_sampler.rank == 0:\n", + " plt.figure()\n", + " plt.subplot(311)\n", + " dl.plot(z_fun, title=\"Optimal control\")\n", + " plt.plot(x_obs, Z_TARGET*np.ones_like(x_obs), '--r', label=\"Target\")\n", + " plt.legend()\n", + "\n", + " plt.subplot(312)\n", + " dl.plot(m_fun, title=\"Parameter\")\n", + "\n", + " plt.subplot(313)\n", + " dl.plot(u_fun, title=\"State\")\n", + " plt.plot(x_obs, x_obs, '--r', label=\"Target\")\n", + " plt.legend()\n", + "\n", + " plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Distribution of QoIs at optimal controls\n", + "We can now sample the QoI distribution using the optimal control. We will compare this against directly using the target value of $z = 1$. \n", + "As with the CVaR tutorial, we will make a new risk measure as the evaluation risk measure to do the sampling. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "SAMPLE_SIZE = 2000\n", + "RANDOM_SEED = 111\n", + "\n", + "risk_settings_evaluation = soupy.meanVarRiskMeasureSAASettings()\n", + "risk_settings_evaluation[\"sample_size\"] = SAMPLE_SIZE \n", + "risk_settings_evaluation[\"seed\"] = RANDOM_SEED\n", + "risk_evaluation = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings_evaluation, comm_sampler=comm_sampler)\n", + "\n", + "# Constrained \n", + "z = risk_evaluation.generate_vector(soupy.CONTROL)\n", + "z.set_local(z_optimal_np)\n", + "risk_evaluation.computeComponents(z, order=0)\n", + "qoi_samples_constrained = risk_evaluation.gather_samples()\n", + "\n", + "# Unconstrained, use z = 1\n", + "z = dl.interpolate(dl.Constant(Z_TARGET), Vh[soupy.CONTROL]).vector()\n", + "risk_evaluation.computeComponents(z, order=0)\n", + "qoi_samples_unconstrained = risk_evaluation.gather_samples()\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We plot the QoI distributions at the control from solving the chance-constrained optimization problem, $z^*$, and at the target control $z = 1$. \n", + "Observe that when using the target control, the probability of exceedance is much larger than using the chance-constraint, which attempts to respect the $p_{p_\\mathrm{tol}} = 0.05$ threshold. \n", + "\n", + "Note that the presence of errors in the smoothing approximation of the indicator and sample average approximation of the expectation means this constraint is not perfectly respected.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Probability of exceedance if using z = 1 (unconstrained): 0.996\n", + "Probability of exceedance if using chance constraint: 0.069\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGdCAYAAADJ6dNTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKd0lEQVR4nO3de1yO9/8H8Netw925lI466kDkkNqhjMrmEDMbWzZGOc4XMzVM83PcCBti5nyIMcM2ZvuaNCsiTOk2FDlE5CbHTlKq+/eHn/vn1sF9p1xX9Xo+Hvdj7uv4ui/Z/e5zfT6fS6JQKBQgIiIiEokmQgcgIiIiehqLEyIiIhIVFidEREQkKixOiIiISFRYnBAREZGosDghIiIiUWFxQkRERKLC4oSIiIhERVvoAOooLy/H9evXYWxsDIlEInQcIiIiUoNCoUB+fj7s7OzQpIn67SH1oji5fv06HBwchI5BRERENXD16lXY29urvX29KE6MjY0BPP5wJiYmAqehupCTk4Pt27cjJCQEVlZWQschIqJakJeXBwcHB+X3uLok9eHZOnl5eTA1NUVubi6LEyIionqipt/f7BBLonDv3j3s2LED9+7dEzoKEREJjMUJiUJmZiZCQkKQmZkpdBQiIhIYixMiIiISlXrRIZaISKFQoLS0FGVlZUJHIaL/o6WlBW1t7Vqf5oPFCRGJXklJCeRyOR48eCB0FCJ6hoGBAWxtbaGrq1trx2RxQqKgr68Pb29v6OvrCx2FRKa8vByZmZnQ0tKCnZ0ddHV1ORkjkQgoFAqUlJTg1q1byMzMhLu7u0YTrVWHxQmJgqenJ06cOCF0DBKhkpISlJeXw8HBAQYGBkLHIaKn6OvrQ0dHB1euXEFJSQn09PRq5bjsEEtE9UJt/UZGRLWrLv5t8l87iUJqaiqkUilSU1OFjkJERAJjcUKi8OTeZT2YsJiIiOoY+5wQUb20OC7jpZ4vvJtHnZ8jISEBQUFBuHfvHszMzKrcztnZGRMmTMCECRPqPNOLiomJwYQJE3D//n2hozR69ennhi0nRES1bOXKlTA2NkZpaalyWUFBAXR0dNC5c2eVbRMTEyGRSJCRkQF/f3/I5XKYmpoCePzFXl2RUtecnZ0RHR0t2Pkbu8DAwFotJI4fP45Ro0apvX1CQgIkEokghSWLEyKiWhYUFISCggIkJycrlyUmJsLGxgbHjx9Xma8lISEBdnZ28PDwgK6uLmxsbDhUmtT2ZHJCdVhaWtabEW8sTkgUPD09cfr0aXh6egodheiFtWzZEnZ2dkhISFAuS0hIQN++feHq6oqkpCSV5UFBQco/P/lNNSEhAUOHDkVubi4kEgkkEglmzpyp3O/BgwcYNmwYjI2N4ejoiNWrV6tkOHXqFLp27Qp9fX1YWFhg1KhRKCgoUK6v7Lfyd999F2FhYcr1V65cQXh4uPL8Vbl//z5GjRoFa2tr6OnpwcvLC3/88YfKNrGxsfD09ISRkRF69uwJuVyuXHf8+HF069YNzZo1g6mpKQICAipMLSCRSLB27Vq89957MDAwgLu7O3bv3q2yzZkzZ9C7d2+YmJjA2NgYnTt3xsWLF5XrN2zYAE9PT+jp6aFVq1ZYvnx5lZ8JeDzHzvz58+Hm5gapVApHR0fMmTNH7WscFhaGd999F99++y1sbW1hYWGBsWPH4tGjR8ptli9fDnd3d+jp6cHa2hrvv/++ct8DBw5gyZIlyut/+fJl5c9IbGwsfH19IZVKkZiYiIsXL6Jv376wtraGkZERXnnlFfz1118qn+fZlrDqrunly5eVP5dNmzaFRCJR/my8DCxOGrL4KNWXiOnr66NNmzachI0ajMDAQMTHxyvfx8fHIzAwEAEBAcrlJSUlOHLkiPJL4Gn+/v6Ijo6GiYkJ5HI55HI5Jk6cqFy/cOFC+Pr6IjU1FWPGjMF//vMfnD17FsDjwqVnz55o2rQpjh8/jh07duCvv/7CuHHj1M7/66+/wt7eHrNnz1aevzLl5eUIDg5GUlISNm/ejLS0NMybNw9aWlrKbR48eIBvv/0WP/zwAw4ePIisrCyVz5Kfn4/Q0FAkJibi6NGjcHd3R69evZCfn69yrlmzZiEkJAT//vsvevXqhUGDBuHu3bsAgOzsbHTp0gV6enr4+++/kZKSgmHDhilbFdasWYOpU6dizpw5SE9Px9y5czFt2jRs3LixymsQGRmJ+fPnY9q0aUhLS8OPP/4Ia2trja5xfHw8Ll68iPj4eGzcuBExMTGIiYkBACQnJ2P8+PGYPXs2zp07h71796JLly4AgCVLlsDPzw8jR45UXn8HBwflcSdPnoyoqCikp6ejXbt2KCgoQK9evfDXX38hNTUVPXr0QJ8+fZCVlVXt33NV19TBwQG//PILAODcuXOQy+VYsmRJtceqTewQS6Jw5coVfPXVV5g2bRqcnJyEjkP0wgIDAxEeHo7S0lIUFRUhNTUVXbp0QVlZGZYuXQoAOHr0KIqKiiotTnR1dWFqagqJRAIbG5sK63v16oUxY8YAAL744gssXrwYCQkJaNWqFbZs2YKioiJs2rQJhoaGAIBly5ahT58+mD9/vvILtjrm5ubQ0tKCsbFxped/4q+//sI///yD9PR0eHg87jTcokULlW0ePXqElStXwtXVFQAwbtw4zJ49W7m+a9euKtuvWrUKTZs2xYEDB/D2228rl4eFheGjjz4CAMydOxffffcd/vnnH/Ts2RPff/89TE1N8dNPP0FHRwcAlHkA4KuvvsLChQvRr18/AICLiwvS0tKwatUqhIaGVvhc+fn5WLJkCZYtW6Zc7+rqijfeeAMA1L7GTZs2xbJly6ClpYVWrVqhd+/e2L9/P0aOHImsrCwYGhri7bffhrGxMZycnODt7Q0AMDU1ha6uLgwMDCq9/rNnz0a3bt2U7y0sLNC+fXvl+6+//ho7d+7E7t27qy1Kq7um5ubmAAArK6uX3veJLSckCnfu3MG6detw584doaMQ1YqgoCAUFhbi+PHjSExMhIeHB6ysrBAQEIDjx4+jsLAQCQkJcHR0rPBlro527dop//ykgMnJyQEApKeno3379sovTQDo1KkTysvLce7cuRf/cE+RyWSwt7dXKQSeZWBgoCxMAMDW1laZFQBycnIwevRoeHh4wNTUFKampigoKKjwW//Tn9nQ0BDGxsbK48hkMnTu3FlZmDzt1q1buHr1KoYPHw4jIyPl6+uvv1a57fO09PR0FBcX480336xyvTrXuE2bNiqtSE9/9m7dusHJyQktWrTA4MGDsWXLFrWfH+Xr66vyvrCwEJMnT0br1q1hZmYGIyMjnD179rktJ9VdUyGx5YSIqA64ubnB3t4e8fHxuHfvHgICAgAANjY2cHFxweHDhxEfH1+h1UBdz34JSyQSlJeXA3jcSbKqPiJPljdp0qTCvEJP94VQlzq3YivL+vS5w8LCcOvWLURHR8PJyQlSqRR+fn4oKSl57nGefObqcjzZZs2aNXjttddU1j1dODzteZ9LnWv8vMzGxsY4ceIEEhISsG/fPkyfPh0zZ87E8ePHn9tS8XRRBACTJk1CbGwsvv32W7i5uUFfXx/vv/9+hWv4rOryCYktJ0REdSQoKAgJCQlISEhAYGCgcnlAQABiY2Nx9OjRSm/pPKGrq4uysjKNz9u6dWvIZDIUFhYqlx0+fBhNmjRRtnBYWlqq9CMpKyvD6dOnNT5/u3btcO3aNWRk1HzemcTERIwfPx69evVCmzZtIJVKcfv2bY2O0a5dOyQmJlZaYFlbW6N58+a4dOkS3NzcVF4uLi6VHs/d3R36+vrYv39/pevVucbq0NbWxltvvYUFCxbg33//xeXLl/H3338D0OzvPzExEWFhYXjvvffQtm1b2NjY4PLly2rnqMyTpwzX5GfwRWlcnBw8eBB9+vSBnZ0dJBIJdu3a9dx9iouLMXXqVGVF7OrqivXr19ckLxFRvREUFIRDhw5BJpMpW06Ax8XJmjVr8PDhw2qLE2dnZxQUFGD//v24ffu22k3+gwYNgp6eHkJDQ3H69GnEx8fj008/xeDBg5V9Ibp27Yr//ve/+O9//4uzZ89izJgxFeazcHZ2xsGDB5GdnV1lsRAQEIAuXbqgf//+iIuLQ2ZmJv7880/s3btXrazA41amH374Aenp6Th27BgGDRqkcef4cePGIS8vDx9++CGSk5Nx/vx5/PDDD8pbLDNnzkRUVBSWLFmCjIwMnDp1Chs2bMCiRYsqPZ6enh6++OILTJ48GZs2bcLFixdx9OhRrFu3DoB61/h5/vjjDyxduhQymQxXrlzBpk2bUF5ejpYtWwJ4fP2PHTuGy5cv4/bt29W2aLi5ueHXX3+FTCbDyZMnMXDgwBduAXFycoJEIsEff/yBW7duqYxEqmsa39YpLCxE+/btMXToUPTv31+tfUJCQnDz5k2sW7cObm5uyMnJUXtcNjUO1tbWmDJlitr/qIlexoytLyooKAhFRUVo1aqVys92QEAA8vPz4erqqjIC41n+/v4YPXo0BgwYgDt37mDGjBkqw4mrYmBggNjYWHz22Wd45ZVXYGBggP79+6t8EQ8bNgwnT57EkCFDoK2tjfDw8AqF0uzZs/HJJ5/A1dUVxcXFVT5e4pdffsHEiRPx0UcfobCwEG5ubpg3b95zcz6xfv16jBo1Ct7e3nB0dMTcuXNVRvOow8LCAn///TcmTZqEgIAAaGlpoUOHDujUqRMAYMSIETAwMMA333yDyZMnw9DQEG3btq12krNp06ZBW1sb06dPx/Xr12Fra4vRo0cDUO8aP4+ZmRl+/fVXzJw5Ew8fPoS7uzu2bt2KNm3aAAAmTpyI0NBQtG7dGkVFRcjMzKzyWIsXL8awYcPg7++PZs2a4YsvvkBeXp7aWSrTvHlzzJo1C1OmTMHQoUMxZMgQ5UijuiZRvMDDTCQSCXbu3Il33323ym327t2LDz/8EJcuXVL2/NVUXl4eTE1NkZubCxMTkxqmbYSeHT4cFClMDqIX8PDhQ2RmZsLFxaXWHsdORLWnun+jNf3+rvM+J7t374avry8WLFiA5s2bw8PDAxMnTkRRUVGV+xQXFyMvL0/lRQ1bfn4+EhISKsxrQEREjU+dFyeXLl3CoUOHcPr0aezcuRPR0dH4+eefMXbs2Cr3iYqKUg4nMzU1rbbZkxqG8+fPIygoCOfPnxc6ChERCazOi5Py8nJIJBJs2bIFr776Knr16oVFixYhJiamytaTyMhI5ObmKl9Xr16t65hEREQkEnU+z4mtrS2aN2+ufMom8Pg5KgqFAteuXYO7u3uFfaRSKaRSaV1HIyIiIhGq85aTTp064fr16ypDkDIyMtCkSRPY29vX9emJiIiontG4OCkoKIBMJoNMJgMAZGZmQiaTKafIjYyMxJAhQ5TbDxw4EBYWFhg6dCjS0tJw8OBBTJo0CcOGDeND3khJR0cHzZs3r3TqaSIialw0vq2TnJysMhY+IiICABAaGoqYmBjI5XKVufyNjIwQFxeHTz/9FL6+vrCwsEBISAi+/vrrWohPDUXbtm1x7do1oWMQEZEIaFycBAYGVjkRD4BKJ2hp1aoV4uLiND0VERERNUJ8tg6JwqlTp2Bvb49Tp04JHYVIMAkJCZBIJBWmkX+Ws7MzoqOja+28ly9fhkQiUd6uJ+GEhYVVO7FpY8GnEpMoPHr0CNnZ2TV6Kio1Us/OgFzXNJhheeXKlZg0aRLu3bsHbe3H/5stKChA06ZN8frrryMxMVG5bWJiIrp06YJz587B398fcrlcOboxJiYGEyZMeG6xQsKZOXMmdu3aVWuF3ZIlS6q9O1EZdWZrr2/YckJEVMuCgoJQUFCA5ORk5bLExETY2Njg+PHjKg/wS0hIgJ2dHTw8PKCrqwsbGxtIJBIhYlMdUvcXL1NTU5iZmdVtmHqAxQkRUS1r2bIl7OzskJCQoFyWkJCAvn37wtXVFUlJSSrLnwwyePq2TkJCAoYOHYrc3FxIJBJIJBKVh/49ePAAw4YNg7GxMRwdHbF69epqM5WXl2P+/Plwc3ODVCqFo6Mj5syZo7LNpUuXEBQUBAMDA7Rv3x5HjhxRrrtz5w4++ugj2Nvbw8DAAG3btsXWrVtV9g8MDMT48eMxefJkmJubw8bGpsKDCu/fv49Ro0bB2toaenp68PLywh9//KFcn5SUhC5dukBfXx8ODg4YP348CgsLq/1sTx6Toqenh2bNmqFfv37Kdffu3cOQIUPQtGlTGBgYIDg4WGUm6piYGJiZmSE2Nhaenp4wMjJCz549IZfLldskJCTg1VdfhaGhIczMzNCpUydcuXIFMTExmDVrFk6ePKn8O3rS71IikWDlypXo27cvDA0N8fXXX6OsrAzDhw+Hi4sL9PX10bJlSyxZskTlszx7W+d519TZ2RkA8N5770EikSjf13csToiI6kBgYCDi4+OV7+Pj4xEYGIiAgADl8pKSEhw5cqTC04CBx08kjo6OhomJCeRyOeRyucqTehcuXAhfX1+kpqZizJgx+M9//oOzZ89WmScyMhLz58/HtGnTkJaWhh9//LHCU8CnTp2KiRMnQiaTwcPDAx999JHyCfIPHz6Ej48P/vjjD5w+fRqjRo3C4MGDcezYMZVjbNy4EYaGhjh27BgWLFiA2bNnKwdElJeXIzg4GElJSdi8eTPS0tIwb948aGlpAXjc96xHjx7o168f/v33X2zbtg2HDh3CuHHjqvxc//3vf9GvXz/07t0bqamp2L9/P3x9fZXrw8LCkJycjN27d+PIkSNQKBTo1auXSkvGgwcP8O233+KHH37AwYMHkZWVpbzWpaWlePfddxEQEIB///0XR44cwahRoyCRSDBgwAB8/vnnaNOmjfLvaMCAAcrjzpgxA3379sWpU6cwbNgwlJeXw97eHtu3b0daWhqmT5+OL7/8Etu3b6/y8z3vmh4/fhwAsGHDBsjlcuX7+o59TkgU3N3dER8fX+mMwUT1UWBgIMLDw1FaWoqioiKkpqaiS5cuKCsrw9KlSwEAR48eRVFRUaXFia6uLkxNTSGRSGBjY1Nhfa9evTBmzBgAwBdffIHFixcjISEBrVq1qrBtfn4+lixZgmXLliE0NBQA4OrqijfeeENlu4kTJ6J3794AgFmzZqFNmza4cOECWrVqhebNm6sUR59++in27t2LHTt24LXXXlMub9euHWbMmAHg8b/rZcuWYf/+/ejWrRv++usv/PPPP0hPT4eHhwcAoEWLFsp9v/nmGwwcOBATJkxQ7r906VIEBARgxYoVlT6Ves6cOfjwww8xa9Ys5bL27dsDePzMrt27d+Pw4cPw9/cHAGzZsgUODg7YtWsXPvjgAwCPb7msXLkSrq6uAIBx48Zh9uzZAB4/VTc3Nxdvv/22cr2np6fyXEZGRtDW1q7072jgwIEYNmyYyrKnc7q4uCApKQnbt29HSEhIhf3VuaaWlpYAADMzs0oz1FcsTkgUjI2NERgYKHQMoloTFBSEwsJCHD9+HPfu3YOHhwesrKwQEBCAwYMHo7CwEAkJCXB0dFT5glZXu3btlH9+UsDk5ORUum16ejqKi4vx5ptvqn1MW1tbAEBOTg5atWqFsrIyzJs3D9u2bUN2djaKi4tRXFwMQ0PDKo/x5DhPcslkMtjb2ysLk2elpKTgwoUL2LJli3KZQqFAeXk5MjMzVYqCJ2QyGUaOHFnl59bW1lYpniwsLNCyZUukp6crlxkYGCgLj2czm5ubIywsDD169EC3bt3w1ltvISQkRHl9qvN0C84TK1euxNq1a3HlyhUUFRWhpKQEHTp0qPY41V3Thoq3dUgUsrOzERkZiezsbKGjENUKNzc32NvbIz4+HvHx8QgICAAA2NjYwMXFBYcPH0Z8fDy6du1ao+M/O5uyRCJBeXl5pduqOxv308d80in3yTEXLlyIxYsXY/Lkyfj7778hk8nQo0cPlJSUqJ3reTnKy8vxySefKGchl8lkOHnyJM6fP69SPKj72aoa9aJQKFQ6HVeW+el9N2zYgCNHjsDf3x/btm2Dh4cHjh49Wu1nAVChcNu+fTvCw8MxbNgw7Nu3DzKZDEOHDq1wDZ+lyd91Q8HihETh5s2bmDdvHm7evCl0FKJaExQUhISEBCQkJKi0DAYEBCA2NhZHjx6t9JbOE7q6uigrK3vhHO7u7tDX18f+/ftrfIzExET07dsXH3/8Mdq3b48WLVqodCxVR7t27XDt2jVkZGRUur5jx444c+YM3NzcKrx0dXWrPGZVn6t169YoLS1V6Rdz584dZGRkVNoKUx1vb29ERkYiKSkJXl5e+PHHHwFo9neUmJgIf39/jBkzBt7e3nBzc8PFixc1ylEZHR2dWvk5ERMWJ0REdSQoKAiHDh2CTCZTtpwAj4uTNWvW4OHDh9UWJ87OzigoKMD+/ftx+/ZtlSHImtDT08MXX3yByZMnY9OmTbh48SKOHj2KdevWqX0MNzc3xMXFISkpCenp6fjkk09w48YNjXIEBASgS5cu6N+/P+Li4pCZmYk///wTe/fuBfC478yRI0cwduxYyGQyZZ+RTz/9tMpjzpgxA1u3bsWMGTOQnp6OU6dOYcGCBQAeF2V9+/bFyJEjcejQIZw8eRIff/wxmjdvjr59+6qVOTMzE5GRkThy5AiuXLmCffv2qRQ3zs7OymfM3b59G8XFxVUey83NDcnJyYiNjUVGRgamTZtWKx1YnZ2dsX//fty4cQP37t174eOJAYsTIqI6EhQUhKKiIri5uamMjAkICEB+fj5cXV3h4OBQ5f7+/v4YPXo0BgwYAEtLS+WXbk1MmzYNn3/+OaZPnw5PT08MGDBAo34L06ZNQ8eOHdGjRw8EBgbCxsamRpN+/fLLL3jllVfw0UcfoXXr1pg8ebLyt/527drhwIEDOH/+PDp37gxvb29Mmzat2v4dgYGB2LFjB3bv3o0OHTqga9euKi0lGzZsgI+PD95++234+flBoVBgz549aj9k1MDAAGfPnkX//v3h4eGBUaNGYdy4cfjkk08AAP3790fPnj0RFBQES0vLCsOrnzZ69Gj069cPAwYMwGuvvYY7d+4oOzW/iIULFyIuLg4ODg7w9vZ+4eOJgUSh6VR0AsjLy4OpqSlyc3NhYmIidJz649kZNDWY4fJlO3HiBHx8fJCSkoKOHTsKHYdE5OHDh8jMzISLi0ulozWISFjV/Rut6fc3W05IFCwsLDB8+HBYWFgIHYWIiATGocQkCk5OTli7dq3QMYiISATYckKiUFRUhDNnzqCoqEjoKEREJDAWJyQK6enp8PLyUpkYiYiIGicWJ0RERCQqLE6IiIhIVFicEBERkaiwOCFRkEgk0NXVVXneBRERNU4cSkyi4O3tXe20z0RE1Hiw5YSIiIhEhcUJiUJ6ejo6duzIocRERMTihMShqKgIqampnISN6DkCAwMxYcKEenledY4h1OcjcWFxQkRUx2QyGT788EPY2NhAV1cXrq6umDlzJh49elTlPvySpsaMxQkRUR1av349Xn31VVhbW+OPP/5Aeno6pk2bhqVLlyIsLOylZCgpKXkp5yGqLSxOiIjqyIEDBzBy5EisWbMGS5Ysga+vL1xdXREWFoZvvvkGP/74I86ePVthv7CwMBw4cABLliyBRCKBRCLB5cuXlevLy8sxefJkmJubw8bGBjNnzlTZPzAwEOPGjUNERASaNWuGbt26QaFQYMGCBWjRogX09fXRvn17/Pzzzyr7/fzzz2jbti309fVhYWGBt956C4WFhWqft7i4GOPHj4eVlRX09PTwxhtv4Pjx41Ven8LCQgwZMgRGRkawtbXFwoUL1b+41KCxOCFRcHFxwfbt2+Hi4iJ0FKJaEx4ejuDgYISGhlZYFxQUBAA4efJkhXVLliyBn58fRo4cCblcDrlcDgcHB+X6jRs3wtDQEMeOHcOCBQswe/ZsxMXFqRxj48aN0NbWxuHDh7Fq1Sr8z//8DzZs2IAVK1bgzJkzCA8Px8cff4wDBw4AAORyOT766CMMGzYM6enpSEhIQL9+/aBQKNQ+7+TJk/HLL79g48aNOHHiBNzc3NCjRw/cvXu30uszadIkxMfHY+fOndi3bx8SEhKQkpKiwRWmhorznJAoNG3aFB988IHQMageefKl/bSmTZvCxcUFDx8+RFpaWoV9OnbsCAA4d+6cSosAADg7O8Pc3By3bt3C1atXVdYZGxvD3d1do3ynTp1Camoq5syZU+n6J52/tbUr/m/Y1NQUurq6MDAwgI2NTYX17dq1w4wZMwAA7u7uWLZsGfbv349u3bopt3Fzc8OCBQsAPG6hWLRoEf7++2/4+fkBAFq0aIFDhw5h1apVCAgIgFwuR2lpKfr16wcnJycAQNu2bdU+b2FhIVasWIGYmBgEBwcDANasWYO4uDisW7cOkyZNUjlWQUEB1q1bh02bNilzb9y4Efb29tVdVmokWJzUd/FRqu+DIoXJ8YJu3ryJLVu2YNCgQbC2thY6DtUDq1atwqxZs1SWDRo0CJs3b8a1a9fg4+NTYZ8nrQBhYWE4evSoyroffvgBH3/8MbZv345x48aprOvevTtiY2M1yieTyQAAHTp0qHT9iRMnAACnT59WaVn5888/0blz52qP3a5dO5X3tra2yMnJUVnm6+ur/HNaWhoePnyoUrwAj/uieHt7AwDat2+PN998E23btkWPHj3QvXt3vP/++2jatKla57148SIePXqETp06Kdfr6Ojg1VdfrXSKgIsXL6KkpERZLAGAubk5WrZsWe1np8aBxQmJQnZ2Nj7//HMEBgayOCG1fPLJJ3jnnXdUlj35IrW3t6/29kBMTEylLScAEBISovKFCTxuOdHUk06oenp6la7//vvv4e/vj4iICAwaNEi5vHnz5s89to6Ojsp7iUSC8vJylWWGhobKPz9Z99///rfC8aVSKQBAS0sLcXFxSEpKwr59+/Ddd99h6tSpOHbsmPJ2a3XnfVL4PfsICoVCUeljKZ6+XUT0LBYnRFQv2drawtbWttJ1enp6yls4lanut3NLS0tYWlq+cL4nrQwHDhzAu+++q7Ju4cKFOHXqFA4ePAhjY+NKix9dXV2UlZW9cA4AaN26NaRSKbKyshAQEFDldhKJBJ06dUKnTp0wffp0ODk5YefOnYiIiHjuOdzc3KCrq4tDhw5h4MCBAIBHjx4hOTm50iHRbm5u0NHRwdGjR+Ho6AgAuHfvHjIyMqrNSI2DxsXJwYMH8c033yAlJQVyuRw7d+6s8A+vKocPH0ZAQAC8vLyUTZ5ERA3RK6+8gp49e2Ls2LF49OgRfH19cfPmTaxduxZbt27Fr7/+WumtpyecnZ1x7NgxXL58GUZGRjA3N0eTJjUbw2BsbIyJEyciPDwc5eXleOONN5CXl4ekpCQYGRkhNDQUx44dw/79+9G9e3dYWVnh2LFjuHXrFjw9PdU6h6GhIf7zn/9g0qRJMDc3h6OjIxYsWIAHDx5g+PDhFbY3MjLC8OHDMWnSJFhYWMDa2hpTp06t8WekhkXj4qSwsBDt27fH0KFD0b9/f7X3y83NxZAhQ/Dmm2/i5s2bmp6WiKje+eWXXzBr1ixMmjQJ165dQ1lZGXr27ImMjIzn3r6ZOHEiQkND0bp1axQVFSEzM1N566kmvvrqK1hZWSEqKgqXLl2CmZkZOnbsiC+//BIAYGJigoMHDyI6Ohp5eXlwcnLCwoULlZ1b1TFv3jyUl5dj8ODByM/Ph6+vL2JjY1X6rTztm2++QUFBAd555x0YGxvj888/R25ubo0/IzUcEsUL3PiTSCRqt5x8+OGHcHd3h5aWFnbt2qVRy0leXh5MTU2Rm5sLExOTmsZtmKrrEFuPOstevHgR4eHhWLx4MVxdXYWOQyLy8OFDZGZmwsXFpcr+G/XFqFGj8NdffyElJaXKL2yi+qa6f6M1/f5+Ke1nGzZswMWLF5VD0J6nuLgYeXl5Ki9q2FxdXbF7924WJtSgfffddxg+fLhypA4RVa7OO8SeP38eU6ZMQWJiYqXj+SsTFRVVYYggNWyPHj3C/fv3YWZmVmFEAFFDIZVKMXXqVKFjEIlenbaclJWVYeDAgZg1axY8PDzU3i8yMhK5ubnK17MTIlHDc+rUKVhZWeHUqVNCRyEiIoHVactJfn4+kpOTkZqaqpzUqLy8HAqFAtra2ti3bx+6du1aYT+pVKoce09ERESNS50WJyYmJhV+E16+fDn+/vtv/Pzzz3yOChEREVWgcXFSUFCACxcuKN9nZmZCJpMpx7VHRkYiOzsbmzZtQpMmTeDl5aWy/5OnVT67nIioOpxRlEic6uLfpsbFSXJysvJpmgCUMweGhoYiJiYGcrkcWVlZtZeQiBq1Jx2kHzx4AH19fYHTENGzHjx4AKDi4w1exAvNc/KycJ6TajSQeU7KyspQWFgIQ0NDaGlpCR2HREYul+P+/fuwsrKCgYFBpc9qIaKXS6FQ4MGDB8jJyYGZmVmlj5Oo6fc3n61DoqClpcXCk6pkY2MDABWevEtEwjMzM1P+G60tLE5IFM6fP49x48Zh2bJlcHd3FzoOiYxEIoGtrS2srKzw6NEjoeMQ0f/R0dGpk9ZuFickCvn5+di3bx/y8/OFjkIipqWlxdt+RI0AH/9IREREosLihIiIiESFt3UammdH6BAREdUzbDkhUXBwcMCyZcvg4OAgdBQiIhIYW07qg3o0X0lNWVpaYuzYsULHICIiEWDLCYnC3bt3sXnzZty9e1foKEREJDAWJyQKly9fxuDBg3H58mWhoxARkcBYnBAREZGosDghIiIiUWFxQkRERKLC4oREwdDQEK+//joMDQ2FjkJERALjUGIShZYtW+LIkSNCxyAiIhFgywkRERGJCosTEoUTJ05AIpHgxIkTQkchIiKBsTghIiIiUWFxQkRERKLC4oSIiIhEhcUJERERiQqHEjdWTz/pWARPOW7dujXOnz8Pe3t7oaMQEZHAWJyQKOjp6cHNzU3oGEREJAK8rUOikJmZiY8//hiZmZlCRyEiIoGxOCFRuHfvHrZs2YJ79+4JHYWIiATG4oSIiIhEhcUJERERiQqLEyIiIhIVFickCra2tpgxYwZsbW2FjkJERALjUGISBVtbW8ycOVPoGEREJAIsTuqjpydQayDy8vJw5MgR+Pn5wcTEROg4REQkIN7WIVG4cOECevbsiQsXLggdhYiIBKZxcXLw4EH06dMHdnZ2kEgk2LVrV7Xb//rrr+jWrRssLS1hYmICPz8/xMbG1jQvERERNXAaFyeFhYVo3749li1bptb2Bw8eRLdu3bBnzx6kpKQgKCgIffr0QWpqqsZhiYiIqOHTuM9JcHAwgoOD1d4+Ojpa5f3cuXPx22+/4ffff4e3t7empyciIqIG7qV3iC0vL0d+fj7Mzc2r3Ka4uBjFxcXK93l5eS8jGglIKpXC1dUVUqlU6ChERCSwl94hduHChSgsLERISEiV20RFRcHU1FT5cnBweIkJSQht2rTBhQsX0KZNG6GjEBGRwF5qy8nWrVsxc+ZM/Pbbb7Cysqpyu8jISERERCjf5+XlsUCpDQ1wCDIRETU8L63lZNu2bRg+fDi2b9+Ot956q9ptpVIpTExMVF7UsP3777+wtLTEv//+K3QUIiIS2EspTrZu3YqwsDD8+OOP6N2798s4JdUzpaWluH37NkpLS4WOQkREAtP4tk5BQYHKRFmZmZmQyWQwNzeHo6MjIiMjkZ2djU2bNgF4XJgMGTIES5Ysweuvv44bN24AAPT19WFqalpLH4OIiIgaCo1bTpKTk+Ht7a0cBhwREQFvb29Mnz4dACCXy5GVlaXcftWqVSgtLcXYsWNha2urfH322We19BGIiIioIdG45SQwMBAKhaLK9TExMSrvExISND0FERERNWJ8tg6JgoeHB5KSkuDh4SF0FCIiEhifSkyiYGRkBD8/P6FjEBGRCLDlhETh2rVriIiIwLVr14SOQkREAmNxQqKQk5ODxYsXIycnR+goREQkMBYnREREJCosToiIiEhUWJwQERGRqLA4IVFo1qwZxowZg2bNmgkdhYiIBMahxCQKjo6O+P7774WOQUREIsCWExKFBw8e4MSJE3jw4IHQUYiISGAsTkgUzp49Cx8fH5w9e1boKEREJDAWJ0RERCQqLE6IiIhIVFicEBERkaiwOCFRaNKkCYyNjdGkCX8kiYgaOw4lJlHo0KED8vLyhI5BREQiwF9TiYiISFRYnJAopKWloU2bNkhLSxM6ChERCYy3daii+Kj//3NQ5Es55cOHD5GWloaHDx++lPMREZF4seWEiIiIRIXFCREREYkKixMiIiISFRYnJAotWrTAb7/9hhYtWggdhYiIBMYOsSQKZmZmeOedd4SOQUREIsCWExKFGzduICoqCjdu3BA6ChERCYzFCYnC9evX8eWXX+L69etCRyEiIoGxOCEiIiJRYZ8TsXp6IjQiIqJGhC0nREREJCosTkgUzMzM8P7778PMzEzoKEREJDDe1iFRaNGiBXbs2CF0DCIiEgG2nJAolJSU4Nq1aygpKRE6ChERCUzj4uTgwYPo06cP7OzsIJFIsGvXrufuc+DAAfj4+EBPTw8tWrTAypUra5KVGrDTp0/DwcEBp0+fFjoKEREJTOPipLCwEO3bt8eyZcvU2j4zMxO9evVC586dkZqaii+//BLjx4/HL7/8onFYIiIiavg07nMSHByM4OBgtbdfuXIlHB0dER0dDQDw9PREcnIyvv32W/Tv31/T0xMREVEDV+d9To4cOYLu3burLOvRoweSk5Px6NGjSvcpLi5GXl6eyouIiIgahzovTm7cuAFra2uVZdbW1igtLcXt27cr3ScqKgqmpqbKl4ODQ13HJCIiIpF4KaN1JBKJynuFQlHp8iciIyORm5urfF29erXOM5KwOnTogIcPH6JDhw5CRyEiIoHV+TwnNjY2FZ40m5OTA21tbVhYWFS6j1QqhVQqretoJCJNmjTh3zkREQF4CS0nfn5+iIuLU1m2b98++Pr6QkdHp65PT/VERkYGAgMDkZGRIXQUIiISmMbFSUFBAWQyGWQyGYDHQ4VlMhmysrIAPL4lM2TIEOX2o0ePxpUrVxAREYH09HSsX78e69atw8SJE2vnE1CDUFBQgAMHDqCgoEDoKEREJDCNb+skJycjKChI+T4iIgIAEBoaipiYGMjlcmWhAgAuLi7Ys2cPwsPD8f3338POzg5Lly7lMGIx4ROQiYhIRDQuTgIDA5UdWisTExNTYVlAQABOnDih6amIiIioEeKzdYiIiEhUWJyQKDg6OmLNmjVwdHQUOgoREQmszocSE6mjWbNmGDFihNAxiIhIBNhyQqJw+/ZtrF27tspZg4mIqPFgcUKikJWVhZEjR6qM9CIiosaJxQkRERGJCosTIiIiEhUWJ0RERCQqLE5IFIyMjBAQEAAjIyOhoxARkcA4lJhEwcPDAwkJCULHICIiEWDLCYlCeXk5iouLUV5eLnQUIiISGIsTEgWZTAY9PT3l066JiKjxYnFCREREosLihIiIiESFxQkRERGJCkfriEl8lNAJiIiIBMfihETBy8sLV69ehZWVldBRiIhIYCxOSBR0dXVhb28vdAwiIhIB9jkhUbh06RI++OADXLp0SegoREQkMBYnJAr379/Hzz//jPv37wsdhYiIBMbihIiIiESFxQkRERGJCjvEUvWeHd4cFClMDiIiajTYckKiYGdnh7lz58LOzk7oKEREJDC2nJAo2NjYIDKSrTJERMSWExKJ+/fvY/fu3RytQ0RELE5IHC5duoS+fftynhMiImJxQkREROLC4oSIiIhEhcUJERERiQqLExIFPT09tG7dGnp6ekJHISIigdWoOFm+fDlcXFygp6cHHx8fJCYmVrv9li1b0L59exgYGMDW1hZDhw7FnTt3ahSYGqbWrVvjzJkzaN26tdBRiIhIYBoXJ9u2bcOECRMwdepUpKamonPnzggODkZWVlal2x86dAhDhgzB8OHDcebMGezYsQPHjx/HiBEjXjg8ERERNTwaFyeLFi3C8OHDMWLECHh6eiI6OhoODg5YsWJFpdsfPXoUzs7OGD9+PFxcXPDGG2/gk08+QXJy8guHp4ZDJpPBxMQEMplM6ChERCQwjYqTkpISpKSkoHv37irLu3fvjqSkpEr38ff3x7Vr17Bnzx4oFArcvHkTP//8M3r37l3z1CQ+8VGqLw2Vl5cjPz8f5eXldRCOiIjqE42Kk9u3b6OsrAzW1tYqy62trXHjxo1K9/H398eWLVswYMAA6OrqwsbGBmZmZvjuu++qPE9xcTHy8vJUXkRERNQ41KhDrEQiUXmvUCgqLHsiLS0N48ePx/Tp05GSkoK9e/ciMzMTo0ePrvL4UVFRMDU1Vb4cHBxqEpOIiIjqIY2Kk2bNmkFLS6tCK0lOTk6F1pQnoqKi0KlTJ0yaNAnt2rVDjx49sHz5cqxfvx5yubzSfSIjI5Gbm6t8Xb16VZOYREREVI9pVJzo6urCx8cHcXFxKsvj4uLg7+9f6T4PHjxAkyaqp9HS0gLwuMWlMlKpFCYmJiovathatWqFlJQUtGrVSugoREQkMG1Nd4iIiMDgwYPh6+sLPz8/rF69GllZWcrbNJGRkcjOzsamTZsAAH369MHIkSOxYsUK9OjRA3K5HBMmTMCrr74KOzu72v00VG8ZGBigY8eOQscgIiIR0Lg4GTBgAO7cuYPZs2dDLpfDy8sLe/bsgZOTEwBALperzHkSFhaG/Px8LFu2DJ9//jnMzMzQtWtXzJ8/v/Y+BdV7WVlZmD9/Pr744gs4OjoKHYeIiAQkUVR1b0VE8vLyYGpqitzc3IZ9i6cGQ3BfuqDIypc/m72q7apw4sQJ+Pj4ICUlhS0oREQNRE2/v/lsHSIiIhIVFidEREQkKixOiIiISFQ07hBLtag+9DF51tOZNexXUh0rKyuEh4fDysqq1o5JRET1E4sTEgV7e3ssWrRI6BiNzuK4jCrXhXfzeIlJiIj+H2/rkCgUFBTgyJEjKCgoEDoKEREJjMUJiUJGRgb8/f2RkVH1b/JERNQ4sDghIiIiUWFxQkRERKLC4oSIiIhEhcUJiYK2tjaaNWsGbW0OICMiauz4TUCi0K5dO9y6dUvoGEREJAJsOSEiIiJRYXFConDmzBm4ubnhzJkzQkchIiKBsTghUSguLsbFixdRXFwsdBQiIhIYixMiIiISFRYnREREJCosToiIiEhUWJyQKLi5uWHv3r1wc3MTOgoREQmM85yQKJiYmKBHjx5CxyAiIhFgywmJglwux8yZMyGXy4WOQkREAmNxQqIgl8sxa9YsFidERMTihIiIiMSFxQkRERGJCosTIiIiEhUWJyQKTZs2xaBBg9C0aVOhoxARkcA4lJhEwcXFBZs3bxY6BhERiQBbTkgUHj58iAsXLuDhw4dCRyEiIoGxOCFRSEtLg7u7O9LS0oSOQkREAmNxQkRERKLC4oSIiIhEhcUJERERiUqNipPly5fDxcUFenp68PHxQWJiYrXbFxcXY+rUqXBycoJUKoWrqyvWr19fo8BERETUsGk8lHjbtm2YMGECli9fjk6dOmHVqlUIDg5GWloaHB0dK90nJCQEN2/exLp16+Dm5oacnByUlpa+cHhqODp27AiFQiF0DCIiEgGNi5NFixZh+PDhGDFiBAAgOjoasbGxWLFiBaKioipsv3fvXhw4cACXLl2Cubk5AMDZ2fnFUhMREVGDpVFxUlJSgpSUFEyZMkVleffu3ZGUlFTpPrt374avry8WLFiAH374AYaGhnjnnXfw1VdfQV9fv+bJqUE5d+4cwsLCEBMTg5YtWwodp95YHJdR7frwbh4vKQkRUe3RqDi5ffs2ysrKYG1trbLc2toaN27cqHSfS5cu4dChQ9DT08POnTtx+/ZtjBkzBnfv3q2y30lxcTGKi4uV7/Py8jSJSfVQYWEhjh49isLCQqGjEBGRwGrUIVYikai8VygUFZY9UV5eDolEgi1btuDVV19Fr169sGjRIsTExKCoqKjSfaKiomBqaqp8OTg41CQmERER1UMaFSfNmjWDlpZWhVaSnJycCq0pT9ja2qJ58+YwNTVVLvP09IRCocC1a9cq3ScyMhK5ubnK19WrVzWJSURERPWYRrd1dHV14ePjg7i4OLz33nvK5XFxcejbt2+l+3Tq1Ak7duxAQUEBjIyMAAAZGRlo0qQJ7O3tK91HKpVCKpVqEq3+iK/YaZioOs/rV0JE1NBofFsnIiICa9euxfr165Geno7w8HBkZWVh9OjRAB63egwZMkS5/cCBA2FhYYGhQ4ciLS0NBw8exKRJkzBs2DB2iCUlZ2dn/PDDDxzJRUREmg8lHjBgAO7cuYPZs2dDLpfDy8sLe/bsgZOTEwBALpcjKytLub2RkRHi4uLw6aefwtfXFxYWFggJCcHXX39de5+ChFGLrUDm5ub4+OOPa+149BhbXYioPpIo6sHMV3l5eTA1NUVubi5MTEyEjvNiGsttnaBIjTa/desWtm/fjpCQEFhaWtZRqPpJqAKDw5CJ6EXV9Pubz9YhUbh69SrGjRvHzs9ERMTihIiIiMRF4z4nRFS72C+EiEgVW06IiIhIVNhyQqJgbGyM7t27w9jYWOgopCY+14eI6gqLE6obT49KUmPkjru7O2JjY+swEGmKt5uISCi8rUOiUFZWhry8PJSVlQkdhYiIBMbihETh5MmTMDU1xcmTJ4WOQkREAmNxQkRERKLC4oSIiIhEhcUJERERiQqLEyIiIhIVDiUmUWjbti1ycnJgZmYmdBQiIhIYixMSBR0dHT6NmIiIALA4IZG4ePEiwsPDsXjxYri6ugodp1ZxMjMiIs2wzwmJQm5uLn7//Xfk5uYKHYWIiATG4oSIiIhEhbd1iOil40MDiag6bDkhIiIiUWFxQqLQvHlzLFy4EM2bNxc6ChERCYy3dUgUrK2tERERIXQMIiISAbackCjcu3cPO3bswL1794SOQkREAmNxQqKQmZmJkJAQZGZmCh2FiIgExuKEiIiIRIXFCREREYkKixMiIiISFRYnJAr6+vrw9vaGvr6+0FGIiEhgHEpMouDp6YkTJ04IHYOIiESALSdEREQkKixOSBRSU1MhlUqRmpoqdBQiIhIYb+uQKCgUCpSUlEChUAgdhWrJ8x7uR0RUFRYnJKz4qMf/zcgWNgcREYlGjW7rLF++HC4uLtDT04OPjw8SExPV2u/w4cPQ1tZGhw4danJaIiIiagQ0Lk62bduGCRMmYOrUqUhNTUXnzp0RHByMrKysavfLzc3FkCFD8Oabb9Y4LBERETV8GhcnixYtwvDhwzFixAh4enoiOjoaDg4OWLFiRbX7ffLJJxg4cCD8/PxqHJYaLk8nK5xe9xk8PT2FjkJERALTqDgpKSlBSkoKunfvrrK8e/fuSEpKqnK/DRs24OLFi5gxY4Za5ykuLkZeXp7Kixo2fakO2rhYcxI2IiLSrEPs7du3UVZWBmtra5Xl1tbWuHHjRqX7nD9/HlOmTEFiYiK0tdU7XVRUFGbNmqVJNKrnrty4h682x2Nai4FwcnISOo7GODKFiKj21KhDrEQiUXmvUCgqLAOAsrIyDBw4ELNmzYKHh4fax4+MjERubq7ydfXq1ZrEpHrkTt4DrNuTjDt37ggdhYiIBKZRy0mzZs2gpaVVoZUkJyenQmsKAOTn5yM5ORmpqakYN24cAKC8vBwKhQLa2trYt28funbtWmE/qVQKqVSqSTTxejJUloiIiNSiUcuJrq4ufHx8EBcXp7I8Li4O/v7+FbY3MTHBqVOnIJPJlK/Ro0ejZcuWkMlkeO21114sPRERETU4Gk/CFhERgcGDB8PX1xd+fn5YvXo1srKyMHr0aACPb8lkZ2dj06ZNaNKkCby8vFT2t7Kygp6eXoXlREREREANipMBAwbgzp07mD17NuRyOby8vLBnzx5lJ0a5XP7cOU+InmXd1AhTPgqo9PYgERE1LhJFPXiYSV5eHkxNTZGbmwsTExOh42iGfU6AoMiq1z17farbVsQ4Wqd2hXdTvwM9EYlXTb+/+WwdEoX8B8VIyciGj28+jI2NhY5DAquu2GPhQtTw1WgoMVFtO3/tNoIi1uL8+fNCRyEiIoGx5aQu8FYOERFRjbHlhIiIiESFLSdU9xpIp1ciIno52HJCoqCjrYXmzUygo6MjdBQiIhIYW05IFNq2sMG17VOAtm2FjkJERAJjcUIvHzsMExFRNXhbh0Th1KUbsA+Zh1OnTgkdhYiIBMbihEThUWkZsm/n4dGjR0JHISIigbE4ISIiIlFhnxMialBe5DlHnBqfSBzYckJERESiwuKERMHdvhniF42Au7u70FGIiEhgvK1DomBsIEVghxYAn0hMRNTosTghUci+lYtlu45inMcQNG/eXOg4FbxIPwYiItIMb+uQKNy8V4B5Ww/g5s2bQkchIiKBsTghIiIiUeFtHSKqV3iLjajhY3FC9H/4pUdEJA68rUOiYGFigOG9fGFhYSF0FCIiEhhbTkgUnGyaYu3EfoCTk9BRiIhIYGw5IVEoKn6EM5k3UVRUJHQUIiISGIsTEoX0KznwGr4E6enpQkchIiKBsTghIiIiUWGfExKv+Kj//3NQpHA5iIjopWLLCREREYkKixMSBYlEAl0dLUgkEqGjEBGRwHhbh0TB290OxbFfAd7eQkchIiKBseWEiIiIRIUtJyQK6VdyMGjOdmzZ/S48PT012/npjrMAO89SjT3vEQbh3TxeUhKixo0tJyQKRcWPkHrhOidhIyKimhUny5cvh4uLC/T09ODj44PExMQqt/3111/RrVs3WFpawsTEBH5+foiNja1xYCIiImrYNC5Otm3bhgkTJmDq1KlITU1F586dERwcjKysrEq3P3jwILp164Y9e/YgJSUFQUFB6NOnD1JTU184PBERETU8GhcnixYtwvDhwzFixAh4enoiOjoaDg4OWLFiRaXbR0dHY/LkyXjllVfg7u6OuXPnwt3dHb///vsLhyciIqKGR6MOsSUlJUhJScGUKVNUlnfv3h1JSUlqHaO8vBz5+fkwNzevcpvi4mIUFxcr3+fl5WkSk+ohF1tzbJ/+EVyu7QJyeduPiKgx06g4uX37NsrKymBtba2y3NraGjdu3FDrGAsXLkRhYSFCQkKq3CYqKgqzZs3SJBrVc02N9fFBYFuhYxBVi6N5iF6OGnWIfXYWT4VCodbMnlu3bsXMmTOxbds2WFlZVbldZGQkcnNzla+rV6/WJCbVIzfv5mPRjkO4eTdf6ChERCQwjVpOmjVrBi0trQqtJDk5ORVaU561bds2DB8+HDt27MBbb71V7bZSqRRSqVSTaFTPZd/Ow+cr9iCwvQuszY2FjkNERALSqDjR1dWFj48P4uLi8N577ymXx8XFoW/fvlXut3XrVgwbNgxbt25F7969a56W6AU8r0meiIjEQeMZYiMiIjB48GD4+vrCz88Pq1evRlZWFkaPHg3g8S2Z7OxsbNq0CcDjwmTIkCFYsmQJXn/9dWWri76+PkxNTWvxoxAREVFDoHFxMmDAANy5cwezZ8+GXC6Hl5cX9uzZAycnJwCAXC5XmfNk1apVKC0txdixYzF27Fjl8tDQUMTExLz4JyAiIqIGpUbP1hkzZgzGjBlT6bpnC46EhISanIIaGVNDPfTxawVTQz2hoxARkcD44D8SBdfmFtg9Z4jQMYiISAT44D8ShUelZbh1vwCPSsuEjkJERAJjcUKicOrSDVj1m4tTl9SbzI+IiBou3tahBuPIpTsAgKOlHDJM4sPZZYnUx+KEiKiWvMhcOixeiP4fb+sQERGRqLDlhOqH+CjV90GRVW76etZq5Z+POo6qq0RERFRHWJyQKLR3tUXu79NhqKcrdBQiIhIYixMSBS2tJjDhBGxERAT2OSGROH/tNnpM3oDz124LHYWIiATG4oREIf9BMfYln0f+g2KhoxARkcBYnBAREZGosDghIiIiUWGHWGrQnh5WDHBoMRFRfcCWExIFBytTLBvfBw5WpkJHISIigbHlpKaenRSMXoilmRHGvusndAwiegan1SchsOWEROFu3gNsjkvF3bwHQkchIiKBseWEROHyjXsYHLUDKSvHwtzEQOg4RKJTXQsGWy+ooWHLCREREYkKixMiIiISFd7WofqJHZKJiBosFickCoZ6uni9tQOfSkxUAxxRQw0NixNN8Lf1OtPS0RJHlv1H6BhERCQCLE6oXjly6Y7QEYgalOe1uhAJgR1iSRROZGRD0vVLnMjIFjoKEREJjC0n9FJV1fJxLjv3pZz/2WftPMFn7hDVDOdfobrAlhMiIiISFbacEIFPLyYiEhMWJ0REVCc4xJlqisUJiYKznQW2zxmG1s5WQkchanA4IofqGxYnJApSHW3YWzdF6rU8oaMQEZHA2CGWROH6rVzMXLMH12+9nFE7REQkXjVqOVm+fDm++eYbyOVytGnTBtHR0ejcuXOV2x84cAARERE4c+YM7OzsMHnyZIwePbrGoanhyX/wEPuOpeOj7j4ATIWOo4KdZYmEwWHKjZfGxcm2bdswYcIELF++HJ06dcKqVasQHByMtLQ0ODo6Vtg+MzMTvXr1wsiRI7F582YcPnwYY8aMgaWlJfr3718rH4KIiOofIfvCsPARN42Lk0WLFmH48OEYMWIEACA6OhqxsbFYsWIFoqIqPntm5cqVcHR0RHR0NADA09MTycnJ+Pbbb8VRnDz7vJygSGFykKhUNVnbs+uqa0Wpbrvqjv/0tuqei6ixaYwjgRrTZ9aoOCkpKUFKSgqmTJmisrx79+5ISkqqdJ8jR46ge/fuKst69OiBdevW4dGjR9DR0amwT3FxMYqLi5Xvc3Mf90PIy6uDzpKFD1XfV3eOZ7dtoP65fLfa9a86m9f42IVFxZUuLyouUf63qm3E6GFhQZXrnv4cz25X3Wd8etvqjkFEVYvadaLG+9bJd83/+f7vC9WuH9vVrcp1z/t/wPNyV3fu6s77Ip5kUigUmu2o0EB2drYCgOLw4cMqy+fMmaPw8PCodB93d3fFnDlzVJYdPnxYAUBx/fr1SveZMWOGAgBffPHFF1988dUAXlevXtWk3FDUqEOsRCJRea9QKCose972lS1/IjIyEhEREcr35eXluHv3LiwsLKo9j9jk5eXBwcEBV69ehYmJidBxRI3XSn28VurjtVIPr5P6eK3U9+RapaWlwc7OTqN9NSpOmjVrBi0tLdy4cUNleU5ODqytrSvdx8bGptLttbW1YWFhUek+UqkUUqlUZZmZmZkmUUXFxMSEP8Rq4rVSH6+V+nit1MPrpD5eK/U1b94cTZpoNnOJRlvr6urCx8cHcXFxKsvj4uLg7+9f6T5+fn4Vtt+3bx98fX0r7W9CREREjZvGk7BFRERg7dq1WL9+PdLT0xEeHo6srCzlvCWRkZEYMmSIcvvRo0fjypUriIiIQHp6OtavX49169Zh4sSJtfcpiIiIqMHQuM/JgAEDcOfOHcyePRtyuRxeXl7Ys2cPnJycAAByuRxZWVnK7V1cXLBnzx6Eh4fj+++/h52dHZYuXSqOYcR1TCqVYsaMGRVuUVFFvFbq47VSH6+Venid1Mdrpb4XuVYShULT8T1EREREdYfP1iEiIiJRYXFCREREosLihIiIiESFxQkRERGJCouTOnDw4EH06dMHdnZ2kEgk2LVrl9CRRCsqKgqvvPIKjI2NYWVlhXfffRfnzp0TOpborFixAu3atVNO/OTn54c///xT6Fj1QlRUFCQSCSZMmCB0FNGZOXMmJBKJysvGxkboWKKVnZ2Njz/+GBYWFjAwMECHDh2QkpIidCzRcXZ2rvBzJZFIMHbsWLWPweKkDhQWFqJ9+/ZYtmyZ0FFE78CBAxg7diyOHj2KuLg4lJaWonv37igsLBQ6mqjY29tj3rx5SE5ORnJyMrp27Yq+ffvizJkzQkcTtePHj2P16tVo166d0FFEq02bNpDL5crXqVOnhI4kSvfu3UOnTp2go6ODP//8E2lpaVi4cGG9nr28rhw/flzlZ+rJRKwffPCB2seo0bN1qHrBwcEIDg4WOka9sHfvXpX3GzZsgJWVFVJSUtClSxeBUolPnz59VN7PmTMHK1aswNGjR9GmTRuBUolbQUEBBg0ahDVr1uDrr78WOo5oaWtrs7VEDfPnz4eDgwM2bNigXObs7CxcIBGztLRUeT9v3jy4uroiICBA7WOw5YREJTc3FwBgbm4ucBLxKisrw08//YTCwkL4+fkJHUe0xo4di969e+Ott94SOoqonT9/HnZ2dnBxccGHH36IS5cuCR1JlHbv3g1fX1988MEHsLKygre3N9asWSN0LNErKSnB5s2bMWzYMI0e3MvihERDoVAgIiICb7zxBry8vISOIzqnTp2CkZERpFIpRo8ejZ07d6J169ZCxxKln376CSdOnEBUVJTQUUTttddew6ZNmxAbG4s1a9bgxo0b8Pf3x507d4SOJjqXLl3CihUr4O7ujtjYWIwePRrjx4/Hpk2bhI4mart27cL9+/cRFham0X68rUOiMW7cOPz77784dOiQ0FFEqWXLlpDJZLh//z5++eUXhIaG4sCBAyxQnnH16lV89tln2LdvH/T09ISOI2pP335u27Yt/Pz84Orqio0bNyIiIkLAZOJTXl4OX19fzJ07FwDg7e2NM2fOYMWKFSrPkyNV69atQ3BwMOzs7DTajy0nJAqffvopdu/ejfj4eNjb2wsdR5R0dXXh5uYGX19fREVFoX379liyZInQsUQnJSUFOTk58PHxgba2NrS1tXHgwAEsXboU2traKCsrEzqiaBkaGqJt27Y4f/680FFEx9bWtsIvAp6enirPkiNVV65cwV9//YURI0ZovC9bTkhQCoUCn376KXbu3ImEhAS4uLgIHaneUCgUKC4uFjqG6Lz55psVRpwMHToUrVq1whdffAEtLS2BkolfcXEx0tPT0blzZ6GjiE6nTp0qTHOQkZGhfOgtVfRkgEPv3r013pfFSR0oKCjAhQsXlO8zMzMhk8lgbm4OR0dHAZOJz9ixY/Hjjz/it99+g7GxMW7cuAEAMDU1hb6+vsDpxOPLL79EcHAwHBwckJ+fj59++gkJCQkVRjsRYGxsXKHPkqGhISwsLNiX6RkTJ05Enz594OjoiJycHHz99dfIy8tDaGio0NFEJzw8HP7+/pg7dy5CQkLwzz//YPXq1Vi9erXQ0USpvLwcGzZsQGhoKLS1a1BqKKjWxcfHKwBUeIWGhgodTXQqu04AFBs2bBA6mqgMGzZM4eTkpNDV1VVYWloq3nzzTcW+ffuEjlVvBAQEKD777DOhY4jOgAEDFLa2tgodHR2FnZ2dol+/foozZ84IHUu0fv/9d4WXl5dCKpUqWrVqpVi9erXQkUQrNjZWAUBx7ty5Gu0vUSgUitqpk4iIiIheHDvEEhERkaiwOCEiIiJRYXFCREREosLihIiIiESFxQkRERGJCosTIiIiEhUWJ0RERCQqLE6IiIhIVFicEBERkaiwOCEiIiJRYXFCREREosLihIiIiETlfwHLUVRNOFNDQAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if comm_sampler.rank == 0:\n", + " print(\"Probability of exceedance if using z = 1 (unconstrained): %g\" %(np.mean(qoi_samples_unconstrained > TAU)))\n", + " print(\"Probability of exceedance if using chance constraint: %g\" %(np.mean(qoi_samples_constrained > TAU)))\n", + "\n", + " plt.figure()\n", + " bar = plt.hist(qoi_samples_unconstrained, bins=50, density=True, alpha=0.5, label='Without chance constraint')\n", + " bar = plt.hist(qoi_samples_constrained, bins=50, density=True, alpha=0.5, label='With chance constraint')\n", + " plt.axvline(TAU, color='k', linestyle='dashed', linewidth=1, label=\"$Q$-threshold\")\n", + " plt.legend()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "soupy", + "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.6" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +}