{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Tutorial 3: Chance constrained optimization\n", "In this tutorial, we will:\n", "- Setup and solve a chance constrained optimization problem using SOUPy and `scipy.optimize`.\n", "- Show how to use cost functional objects as constraints.\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 strong residual of the PDE. The corresponding weak form, is \n", "\n", "$$\n", "\\text{Find } u \\in \\mathcal{U} \\text{ s.t. } \n", "r(u,m,v,z):= \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", "\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", "\n", "$$ \\mathcal{A} = -\\gamma \\Delta + \\delta \\mathcal{I} $$\n", "\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", "where $f$ and $g$ are functions of scalar variables. \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_function` will be set to identity functions $f(x) = x$ and $g(x) = x$." ] }, { "cell_type": "code", "execution_count": 5, "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/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACYoklEQVR4nOzdd3xT5eIG8Ce7TfduugeUlr03MkSQKXq54HUwRAURFYF7AVERF1dU1IsC6o/hQEUUXKCCCmXK3oUCHdC998g8vz/Spg3dpSNNn+/nk0+Tkzcnb3La5sl73iESBEEAEREREbV54tauABERERE1DQY7IiIiIivBYEdERERkJRjsiIiIiKwEgx0RERGRlWCwIyIiIrISDHZEREREVoLBjoiIiMhKMNgRERERWQkGOyJqNn///Tf++c9/QqVSQS6Xw9vbG1OnTsWxY8fuaL/r16/H1q1bq2yPj4+HSCSq9r6WEBQUhFmzZrXKc9clOTkZr7zyCs6dO9cs+7fk107UnjDYEVGzWLduHYYMGYLExESsWbMGf/zxB9555x0kJSVh6NCh+PDDDxu975qCnUqlwrFjxzBhwoQ7qLl1Sk5OxqpVq5ot2BGRZZC2dgWIyPocOXIECxcuxPjx47Fr1y5IpRX/ah588EHcf//9eO6559CrVy8MGTKkyZ5XoVBg4MCBTba/9qy4uBhKpbK1q0FEDcQWOyJqcqtXr4ZIJMKGDRvMQh0ASKVSrF+/HiKRCP/9739N21955RWIRCKcPXsWDzzwABwdHeHk5IRHHnkEGRkZpnJBQUG4fPkyIiMjIRKJIBKJEBQUBKD6U7Hl+71w4QL++c9/wsnJCa6urli0aBF0Oh2io6Nx7733wsHBAUFBQVizZo1ZfUtLS7F48WL07NnT9NhBgwbhxx9/bPT7YzAYsG7dOvTs2RO2trZwdnbGwIED8dNPP5mVWbNmDcLDw6FQKODp6YkZM2YgMTHRbF8jRoxA165dcfLkSQwbNgxKpRIhISH473//C4PBAAA4cOAA+vXrBwCYPXu26X175ZVXAACzZs2Cvb09Ll68iDFjxsDBwQF33303ACA7Oxvz58+Hr68v5HI5QkJCsGLFCqjV6ka/fiJqPmyxI6ImpdfrsX//fvTt2xd+fn7VlvH390efPn3w119/Qa/XQyKRmO67//77MW3aNMybNw+XL1/GSy+9hKioKBw/fhwymQy7du3C1KlT4eTkhPXr1wMwttTVZdq0aXjkkUcwd+5c7Nu3D2vWrIFWq8Uff/yB+fPnY8mSJfjqq6+wdOlSdOjQAQ888AAAQK1WIzs7G0uWLIGvry80Gg3++OMPPPDAA9iyZQtmzJjR4Pdo1qxZ+PLLLzFnzhy8+uqrkMvlOHPmDOLj401lnnrqKXzyySdYsGABJk6ciPj4eLz00ks4cOAAzpw5A3d3d1PZ1NRUPPzww1i8eDFWrlyJXbt2Yfny5fDx8cGMGTPQu3dvbNmyBbNnz8aLL75oOlVd+fhoNBpMnjwZc+fOxbJly6DT6VBaWoqRI0ciJiYGq1atQvfu3XHo0CGsXr0a586dw+7duxv82omomQlERE0oNTVVACA8+OCDtZabPn26AEBIS0sTBEEQVq5cKQAQnn/+ebNy27ZtEwAIX375pWlbly5dhOHDh1fZZ1xcnABA2LJli2lb+X7fffdds7I9e/YUAAg7d+40bdNqtYKHh4fwwAMP1FhvnU4naLVaYc6cOUKvXr3M7gsMDBRmzpxZ6+s+ePCgAEBYsWJFjWWuXLkiABDmz59vtv348eMCAOGFF14wbRs+fLgAQDh+/LhZ2c6dOwtjx4413T558mSV96bczJkzBQDC5s2bzbZv3LhRACB8++23ZtvfeustAYCwd+9e07b6vHYian48FUtErUIQBACASCQy2/7www+b3Z42bRqkUin2799/R883ceJEs9sREREQiUQYN26caZtUKkWHDh1w8+ZNs7I7duzAkCFDYG9vD6lUCplMhk2bNuHKlSsNrsevv/4KAHj66adrLFP+Wm8fZdq/f39ERETgzz//NNvu7e2N/v37m23r3r17lddRl3/84x9mt//66y/Y2dlh6tSpZtvL63V7PYio9THYEVGTcnd3h1KpRFxcXK3l4uPjoVQq4erqarbd29vb7LZUKoWbmxuysrLuqF63P49cLodSqYSNjU2V7aWlpabbO3fuxLRp0+Dr64svv/wSx44dw8mTJ/HYY4+ZlauvjIwMSCSSKq+zsvLXqlKpqtzn4+NT5b1wc3OrUk6hUKCkpKTe9VIqlXB0dKxSD29v7yrh29PTE1Kp9I6PCRE1PfaxI6ImJZFIMHLkSPz2229ITEystp9dYmIiTp8+jXHjxpn1rwOM/cV8fX1Nt3U6HbKysqoNLy3hyy+/RHBwMLZv324WcBo7eMDDwwN6vR6pqanVBjegIqilpKRUef+Sk5PN+tc1ldvDW3k9jh8/DkEQzO5PT0+HTqdrlnoQ0Z1hix0RNbnly5dDEATMnz8fer3e7D69Xo+nnnoKgiBg+fLlVR67bds2s9vffvstdDodRowYYdrW0NaoOyESiSCXy82CTWpqaqNHxZaf+t2wYUONZUaNGgXAGCorO3nyJK5cuWIasdoQ5QNMGvK+3X333SgsLMQPP/xgtv3zzz833U9EloUtdkTU5IYMGYL3338fCxcuxNChQ7FgwQIEBATg1q1b+Oijj3D8+HG8//77GDx4cJXH7ty5E1KpFPfcc49pVGyPHj0wbdo0U5lu3brhm2++wfbt2xESEgIbGxt069atWV7LxIkTsXPnTsyfPx9Tp05FQkICXnvtNahUKly/fr3B+xs2bBgeffRRvP7660hLS8PEiROhUChw9uxZKJVKPPPMM+jUqROefPJJrFu3DmKxGOPGjTONivX398fzzz/f4OcNDQ2Fra0ttm3bhoiICNjb28PHxwc+Pj41PmbGjBn46KOPMHPmTMTHx6Nbt244fPgw3nzzTYwfPx6jR49ucD2IqHkx2BFRs3jmmWfQr18/vPvuu1i8eDGysrLg6uqKoUOH4vDhwxg0aFC1j9u5cydeeeUVbNiwASKRCJMmTcL7778PuVxuKrNq1SqkpKTgiSeeQEFBAQIDA82mCmlKs2fPRnp6OjZu3IjNmzcjJCQEy5YtQ2JiIlatWtWofW7duhW9e/fGpk2bsHXrVtja2qJz58544YUXTGU2bNiA0NBQbNq0CR999BGcnJxw7733YvXq1Y06La1UKrF582asWrUKY8aMgVarxcqVK01z2VXHxsYG+/fvx4oVK/D2228jIyMDvr6+WLJkCVauXNmYl05EzUwklA9NIyJqRa+88gpWrVqFjIwM9t0iImok9rEjIiIishIMdkRERERWgqdiiYiIiKwEW+yIiIiIrASDHREREZGVYLAjIiIishLtch47g8GA5ORkODg4VLuMDhEREZGlEAQBBQUF8PHxgVhcR5uc0MoiIyOFiRMnCiqVSgAg7Nq1q87HHDhwQOjdu7egUCiE4OBgYcOGDQ16zoSEBAEAL7zwwgsvvPDCS5u5JCQk1JlxWr3FrqioCD169MDs2bPxj3/8o87ycXFxGD9+PJ544gl8+eWXOHLkCObPnw8PD496PR4AHBwcAAAJCQlwdHS8o/oTERERNaf8/Hz4+/ub8kttWj3YjRs3zrQodn1s3LgRAQEBeP/99wEAEREROHXqFN555516B7vy06+Ojo4MdkRERNQm1Kf7WKsHu4Y6duwYxowZY7Zt7Nix2LRpE7RaLWQyWSvVzJzeIOCDP67hkW7u8HS0qb6QRALYVLqvqKjmHYrFgK1t48oWFwM1TVcoEgFKZePKlpQABkPN9bCza1zZ0lJAr2+askqlsd4AoFYDOl3TlLW1Nb7PAKDRAFpt05S1sTH+XjS0rFZrLF8ThQKQShteVqczvhc1kcuB8r+5hpTV643HriYymbF8Q8saDMbftaYoK5Ua3wvA+DdRXNw0ZRvyd8//EdWX5f+Ihpfl/wjj9Sb+H6GVyvDNiVvoFeCCrr5ONZdtSQ3qnNbMgLr72HXs2FF44403zLYdOXJEACAkJydX+5jS0lIhLy/PdCnvY5eXl9dUVa/iq+M3hcClvwiC8V9g9Zfx480fpFTWXHb4cPOy7u41l+3b17xsYGDNZTt3Ni/buXPNZQMDzcv27VtzWXd387LDh9dcVqk0Lzt+fO3vW2VTp9ZetrCwouzMmbWXTU+vKDt/fu1l4+Iqyi5ZUnvZS5cqyq5cWXvZEycqyq5ZU3vZ/fsryn74Ye1lf/mlouyWLbWX/fbbirLfflt72S1bKsr+Usfv+4cfVpTdv7/2smvWVJQ9caL2sitXVpS9dKn2skuWVJSNi6u97Pz5FWXT02svO3NmRdnCwtrLTp0qmKmtLP9HGC/8H1Fx4f8I48UC/kfET58pjHh7vxC49BfhoU+PCQaDQWgueXl5Qn1zS5uc7uT2pkhBEKrdXm716tVwcnIyXfz9/Zu9jhEqR/QKcG725yEiIqKWF3ktE3GZRXCzk2NMZ2+URZFWZ1FLiolEIuzatQtTpkypscxdd92FXr164YMPPjBt27VrF6ZNm4bi4uJqT8Wq1WqoKzX/lndCzMvLa9Y+doIgYM/xGKzdG43kXGMzcTdfRywdF44+ga48zVJTWZ5maXhZnmYxXuep2MaV5f8I43X+j2h42Xb0P+J6WgHW7ruGA9EZxqIKGWaN7IQn7gqBg03zdgPLz8+Hk5NTvXJLmwt2S5cuxc8//4yoqCjTtqeeegrnzp3DsWPH6vU8DXmDmkKpVo/NR+Kwfn8MCtXGfwTju3lj2b0RCHBT1vFoIiIiai0peSV4b981fHc6EQYBkIhFmN7PHwvv7lhzH/om1pDc0uqDJwoLC3Hjxg3T7bi4OJw7dw6urq4ICAjA8uXLkZSUhM8//xwAMG/ePHz44YdYtGgRnnjiCRw7dgybNm3C119/3VovoU42Mgnmj+iAf/bxx9p917D95C3suZiKP6LSMXtIEJ4e1QGOzZz2iYiIqP7yS7XYcCAGmw/HQa0ztibf28Ub/763E0I97Fu5djVr9Ra7AwcOYOTIkVW2z5w5E1u3bsWsWbMQHx+PAwcOmO6LjIzE888/j8uXL8PHxwdLly7FvHnz6v2cLd1id7urqfl4Y/cVHLqeCQBwtZPj+dEd8a/+AZBK2mS3RyIiIqug1unx5d+38OFf15FTbDzV3TfQBcvHl3WjagVt9lRsS2ntYAcY+98diM7A67ujEJNh7AvTwdMeKyZEYESYB5c6IyIiakEGg4CfLyTj7d+jkZhj7FsX6mGHpfeG457OXq36ucxgVwdLCHbltHoDvjlxC2v3XTN9MxjW0R0vTuiMTt51zzBNREREd+bw9Uys/vUKLifnAwA8HRR4/p4w/LOPn0WcSWOwq4MlBbtyeSVafLT/BrYciYNWL0AsAh7sH4DnR4fBw0HR2tUjIiKyOpeT8/DfX6+aukbZK6SYNzwEjw0NhlLe6sMQTBjs6mCJwa7czawi/PfXq/j1UioA4y/Z/JGheGxIMGxkklauHRERUduXkF2Mtfuu4YdzSRAEQCYR4eEBgXhmVAe42VteYwqDXR0sOdiVOxGXjdd3R+FCYh4AwNfZFkvHhWNSdxX73xERETVCTpEGH+2/gc+P3YRGbxzpOqmHD/49ppNFTz/GYFeHthDsAGNHzh/OJWHNb9FIzTdOvtgrwBkvTeyM3gEurVw7IiKitqFUq8eWI/FYf+AGCkqN88kODnXDsnHh6O7n3LqVqwcGuzq0lWBXrkSjx6eHYrHhQAxKtMYZ1Cf38MF/7u0EPxfL/YZBRETUmvQGAd+fScR7+64hJc/YQBLu7YBl48IxvA3NQMFgV4e2FuzKpeWX4p3fo/HdmUQIAiCXivH40GA8NSK02ZczISIiaisEQcD+6HS89Ws0otMKAAA+TjZYPKYTpvTyhUTcNgJdOQa7OrTVYFfuUlIeXt8dhb9jswEA7vZyLLqnE6b3829zv6xERERN6VxCLlbvuYLjccbPSCdbGZ4eGYoZg4La7CBEBrs6tPVgBxi/jeyLSsPqX68iLtM4wXG4twNWTIjAsI4erVw7IiKilhWfWYS3f4/G7ospAIxntWYPDsL8ER3gpGzbZ7UY7OpgDcGunEZnwJd/38QHf15HXolxguORnTywYkIEOnhygmMiIrJumYVqrPvzOrYdvwWdQYBIBDzQyw+LxoTB19m2tavXJBjs6mBNwa5cbrEGH/x5HV8cuwmdQYBELMLDAwLw3N0dLXJOHiIiojtRUKrFp4fisOlQLIo0xoGFIzp5YOm94YhQWcdnezkGuzpYY7ArF5tRiDf3XMUfV9IAWO4s2kRERI1RqtXji2M3sf7ADdNSnN39nLBsXDgGh7q3cu2aB4NdHaw52JU7eiMTb/56BZeSKta9W3RPGKZayLp3REREDaHTG/Dd6UR88Od109QlIR52+PeYTri3q3ebmbqkMRjs6tAegh1gnOD45wvJePv3aCTmlAAAOnraY+m94bg7wtOq/wiIiMg6GAwC9lxKwdq91xBbNljQx8kGC0eH4YHevu2isYLBrg7tJdiVU+uMzdYf7r+B3LJm6/7Brlg+Lhy9uIIFERFZIEEQcPB6Jtb8dhWXk41nn1zt5Hh6ZAc8PCCgzU5d0hgMdnVob8GuXF6JFhsOxGDzkThodMY18iZ0U+HfYzshyN2ulWtHRERkdPpmDtb8dtU0F529QorHhwVjztDgdjkhP4NdHdprsCuXnFuCtfuu4fuyFSykYhEeGRiIZ0Z14AhaIiJqNdGpBXj792jTAEC5VIxHBwZi/ojQdv35xGBXh/Ye7MpdScnHW79dxYHoDAAcQUtERK3jVlYx3vvjGn44lwRBAMQi4J99/PHc6I7wsZK56O4Eg10dGOzMHbmRidUcQUtERC0sPb8U6/66gW9O3oJWb4wj47t5Y9E9ndDB076Va2c5GOzqwGBXFUfQEhFRS8kr0eLjyBhsORKPEq1xcuFhHd3xn7Hh6Obn1Mq1szwMdnVgsKsZR9ASEVFzKdHoseVoHDYeiEF+qQ4A0CvAGf8ZG45BoW6tXDvLxWBXBwa7upWPoN1yJA5qjqAlIqI7oNUb8M3JBKz78zrSC9QAgDAveywZ0wn3dPbiWaE6MNjVgcGu/qobQfvwgAA8c3dHuLfjEUpERFQ3g0HAT+eTsXbfNdzKLgYA+LnYYtE9Ybivpy8kYga6+mCwqwODXcNxBC0REdWXIAj462o63v49GldTCwAA7vYKPDOqA/7VPwByKQfmNQSDXR0Y7BqPa9ASEVFtjsdmYc3v0Th9MwcA4GAjxbzhoZg9JIgNAY3EYFcHBrs7wxG0RER0u0tJeXj792hEXjOe2VFIxZg1JAhPDQ+Fs1LeyrVr2xjs6sBg1zTUOj2+/PsW1v11nSNoiYjaqRvphXj/j2v45UIKAGNf7On9/PHs3R3h5WjTyrWzDgx2dWCwa1o1jaBdMrYTgjmClojIKp2Kz8bGyFjT8l8AMLmHDxbdE8bZE5oYg10dGOyax+0jaMUiYFIPH8wbHooIFd9nIqK2zmAQsDcqDZ8cjMGZW7mm7fd09sLC0R3RxYeTCzcHBrs6MNg1rysp+Vjz21XsLxtBCwAjO3ngqREd0D/YtRVrRkREjVGq1WPnmST836FYxGYWAQDkEjEe6O2Lx4eFcPmvZtaQ3GIRwxjXr1+P4OBg2NjYoE+fPjh06FCt5bdt24YePXpAqVRCpVJh9uzZyMrKaqHaUl0iVI7YMrs/fnlmKCZ2V0EsAvZHZ2Dax8cwdcNR/HklDQZDu/s+QUTU5uQWa/DhX9cx9K39eGHXRcRmFsHRRor5I0JxeOlI/Pcf3RnqLEyrt9ht374djz76KNavX48hQ4bg448/xv/93/8hKioKAQEBVcofPnwYw4cPx3vvvYdJkyYhKSkJ8+bNQ8eOHbFr1656PSdb7FpWfGYRPj4Yi+9PJ0KjN/bB6+TlgHkjQjCxuw9knCaFiMiiJOYUY9PhOGw/mYBijXEtVx8nGzw2NBgP9g+AvYLTlrSkNnUqdsCAAejduzc2bNhg2hYREYEpU6Zg9erVVcq/88472LBhA2JiYkzb1q1bhzVr1iAhIaFez8lg1zrS80ux6Ugctv19C4Vq4xqBfi62ePKuEPyzjz9s5ZJWriERUft2OTkPnxyMxS8XUqAvO7MS7u2AucP5Rbw1tZlgp9FooFQqsWPHDtx///2m7c899xzOnTuHyMjIKo85evQoRo4ciV27dmHcuHFIT0/HtGnTEBERgY0bN1b7PGq1Gmq12nQ7Pz8f/v7+DHatJK9Eiy//vonNh+OQVaQBALjZyTF7SBAeHRQEJ1tZK9eQiKj9EAQBR25k4eODMTh0PdO0fUgHN8y9KxTDOrpzftJW1pBg16ptqZmZmdDr9fDy8jLb7uXlhdTU1GofM3jwYGzbtg3Tp09HaWkpdDodJk+ejHXr1tX4PKtXr8aqVauatO7UeE62Mjw9sgPmDA3GjlMJ+PhgLBJzSvDO3mvYGBmLhwcEYM7QYHhy/iMiomaj0xuw+2IKPo6MRVSKcTUhsQiY0N0Hc+8KQVdfjnBti1q1xS45ORm+vr44evQoBg0aZNr+xhtv4IsvvsDVq1erPCYqKgqjR4/G888/j7FjxyIlJQX//ve/0a9fP2zatKna52GLnWXT6Q345UIKNhyIQXSacU1BuUSMf/Txw9y7QjgfEhFREypS6/DtqQT836E4JOUaVw+ylUkwvZ8/5gwNhr+rspVrSLez6lOxjz76KEpLS7Fjxw7TtsOHD2PYsGFITk6GSqWq83nZx84yCYKA/dHpWL8/BqfK1hgUi4Bx3VR4angovz0SEd2BjAI1Pjsajy/+vom8EuNqQW52cswcHIRHBwbCxY7LflmqNnMqVi6Xo0+fPti3b59ZsNu3bx/uu+++ah9TXFwMqdS82hKJsdN9O5ySz6qIRCKMCvfCqHAvnIzPxoYDMfjrajp2X0jB7gspuCvMA08ND8XAEFf29yAiqqfYjEJ8eigO359JhKZsdaAgNyWeuCsE/+jtBxsZB65Zk1Yfr7xo0SI8+uij6Nu3LwYNGoRPPvkEt27dwrx58wAAy5cvR1JSEj7//HMAwKRJk/DEE09gw4YNplOxCxcuRP/+/eHj49OaL4WaUL8gV/Sb5YorKfn4ODIGP19IwcFrGTh4LQM9/Z3x1IhQ3BPhBbGYAY+IqDpnbuXg48gY7I1KQ3m7R09/Z8wbHoJ7OntDwv+fVqnVpzsBjBMUr1mzBikpKejatSvee+893HXXXQCAWbNmIT4+HgcOHDCVX7duHTZu3Ii4uDg4Oztj1KhReOutt+Dr61uv5+Op2LbnVlYxPj0Ui29PJZjWo+3gaY95w0NxX08OwSciAoxLfv11NR0fH4zByfgc0/a7wz0xd3go+gW58IxHG9Rm+ti1Fga7tiujQI2tR+Pw+bGbKCg1zoXn42SDx4eF4MH+/lDKW70Rmoioxal1evx4NhkfH4xBTIZxyS+ZRIQpPX3x5F0h6Ojl0Mo1pDvBYFcHBru2L79Ui6+O38Kmw3HIKDCOeHZRyjBrcDBmDg6Es5KdgInI+uUUafDViVv47Gg80sv+FzoopHh4YCBmDwmCF6eNsgoMdnVgsLMe5QtTf3wwBjezigEASrkE/+pvnAvPx9m2lWtIRNT0YjMKsflIHL47nYhSrbF7irejDeYMDcaD/f3hYMOJ3q0Jg10dGOysj05vwK+XUrHhQIxpok2pWIQJ3VV4Yhgn2iSitk8QBByLzcKmQ3H482q6aXtnlSMeHxaMid19IJeyv7E1YrCrA4Od9RIEAQevZ2LDgRv4OzbbtH1giCseHxqCUeGeHElLRG2KRmfALxeS8X+H4kxfXAFgdIQn5gwN4RRQ7QCDXR0Y7NqHi4l52HTYuJi1rmwx6xAPO8wZGsy5m4jI4uUWa7DtuHn/ORuZGP/s44/ZQ4IQ4mHfyjWklsJgVwcGu/YlObcEnx2Nx1cnbplG0rrayfHIgAA8OigIHg6KVq4hEVGF6vrPeTooMHNwEB7qH8AVItohBrs6MNi1T4VqHb49mYDNR+KQmGNcH1EuFeP+nr6YMywYYZwOgIhaSXn/uc2H4/DHFfafI3MMdnVgsGvfdHoDfr+chk8PxeJcQq5p+/AwDzwxLARDOrixvwoRtYja+s89NjQYg0L4/4gY7OrEYEflTt/MxqcH4/B7VKppyZ1wbwc8PiwEk3vwGzIRNY+a+s9N7eOH2UOCEcr+c1QJg10dGOzodjezirDlSDy+PZWAYo0eQEWflocHBHDCYyJqEuw/R43BYFcHBjuqSV6xFttO3MRnR+ORlm/8Fm0rk+Cfff3w2JBgBLnbtXINiaitEQQBf8dmY9PhWPafo0ZhsKtDfd8gvV4PrVbbgjVrn2QyGSQSy5p6pLzfy6eH4nClrN+LSASM6eyFx4eFoG8gF9ImotqV/x/ZdDgOl5PZf44aj8GuDnW9QYIgIDU1Fbm5uS1fuXbK2dkZ3t7eFvdPThAEHIvJwqeHYrE/OsO0vYe/Mx4fGoxxXb0hlfCbNhFVYP85amoMdnWo6w1KSUlBbm4uPD09oVQqLS5sWBNBEFBcXIz09HQ4OztDpVK1dpVqdD2tAJsOx2Hn2SRodMa+Mb7Otpg9JAjT+3FtRqL2LjajEFuOxOO704ko0Zr31WX/OboTDHZ1qO0N0uv1uHbtGjw9PeHm5tZKNWx/srKykJ6ejrCwMIs7LXu7zEI1vjh2E1/8fRPZRRoAgINCigf7+2PWkGD4Otu2cg2JqKVkFarxx5U07LmYioPXM0yj69l/jpoSg10danuDSktLERcXh6CgINja8gO6pZSUlCA+Ph7BwcGwsbFp7erUS6lWj11nk/B/h2IRk1EEAJCIRZjQTYXHhwWjm68TW3uJrFBCdjF+v5yKvVFpOBWfDUOlT9G7wz0xZxj7z1HTakiwk7ZQndoc/kG2rLb4ftvIJPhX/wBM7+uPA9fS8enBOByLzcJP55Px0/lkhLjb4Z7OXhjTxQs9/V0gEbe910hExi4jV1MLjGHucprZRMIA0NXXEWM7e2NCdxXXb6VWx2BHdIfEYhFGhXthVLgXLiXlYdPhOOy+kILYzCJ8fDAWHx+Mhbu9HKMjvHBPZy8M6eAOG5lln24mau/0BgFnbuXg90vGlrlb2cWm+8QioH+wK8Z28cY9nb3g56JsxZoSmeOp2BpOxbalU4LWwNre90K1DpHRGdgblYq/rqajoFRnuk8pl+Cujh4Y08ULo8I9OfkxkYVQ6/Q4eiMLe6NSsS8qDZmFGtN9CqkYwzp6YGwXL9wd4QVXDoSgFsRTse1QXacyZ86cia1bt7ZMZW4TFBSEhQsXYuHCha3y/K3BXiHFhO4qTOiuglZvwPHYbOyLMn7zT8krxW+XU/Hb5VRIxCL0D3LFPZ2NrXn+rvzmT9SSCkq12B+dgb2XU3EgOgOF6oovYY42Utwd4YWxXbxwV5gHlHJ+ZJLl42+plUhJSTFd3759O15++WVER0ebtjV0IIhGo4Fczm+kTUEmEWNoR3cM7eiOVyZ3weXkfOwt63h9NbUAx2KzcCw2C6/+EoUIlSPGlIW8Lj6ObbLvIZGlyygwjmT9/XIqjt7IgkZvMN3n5ajAmM7eGNPFCwND3CDjPJXUxvBUrBWeit26dSsWLlxommA5KysLCxYswKFDh5CdnY3Q0FC88MIL+Ne//mV6zIgRI9C1a1fI5XJ8/vnn6NKlCyIjI/HTTz9h8eLFSExMxMCBAzFr1izMmjULOTk5cHZ2BgAcPXoUy5Ytw8mTJ+Hu7o77778fq1evhp2dHUaMGIHIyEiz+lX3K2cN73tj3Moqxr4radh7ORUnbxtd5+tsa2rJ6x/syg8YojtwM6sIey8bw9zpWzmo/G8oxMMOY7t4Y0xnL/Twc4aYA53IwvBUbHMpKqr5PokEqBxIaisrFgOVW9BqKmvXNOuSlpaWok+fPli6dCkcHR2xe/duPProowgJCcGAAQNM5T777DM89dRTOHLkCARBQHx8PKZOnYrnnnsOjz/+OM6ePYslS5aY7fvixYsYO3YsXnvtNWzatAkZGRlYsGABFixYgC1btmDnzp3o0aMHnnzySTzxxBNN8nqsSYCbEnOGBmPO0GBkF2nw19V07ItKReS1DCTllmDr0XhsPRoPRxspRoV7YkwXb9wV5gF7Bf90iWojCAKiUvLx+2XjF6erqQVm9/fwc8KYLt4Y28ULHTwdWqmWRE2PLXYNabGr7bTY+PHA7t0Vt+3sgOLi6ssOHw4cOFBx28MDyMysWq6Rh+b2FrvqTJgwAREREXjnnXcAGFvs8vLycPbsWVOZZcuWYffu3bh48aJp24svvog33njD1GI3Y8YM2Nra4uOPPzaVOXz4MIYPH46ioiLY2NjUq49de22xq0mpVo/D1zOxNyoVf15JR1ZRRSduuUSMIR3ccE9nb4zu7AlPB75fRIBxJOup+GxjmItKRWJOiek+iViEgSEVI1lVTpynlNoOttiRGb1ej//+97/Yvn07kpKSoFaroVarYXdbi2Dfvn3NbkdHR6Nfv35m2/r37292+/Tp07hx4wa2bdtm2iYIAgwGA+Li4hAREdHEr6Z9sJFJMLqzF0Z39jJNu1DeL+9mVjH2R2dgf3QGVvwA9PR3xpjOxg+rDp6cQ4vaD0EQEJ9VjOOxWfg7NgsHr2eaVoMBjOuzDg/zwNgu3hyBTu0Gg11DFBbWfN/ty2Clp9dcVnxbX6n4+EZXqT7effddvPfee3j//ffRrVs32NnZYeHChdBoNGblbg96giBU6bx/ewOvwWDA3Llz8eyzz1Z53oCAgCZ6Be2bRCxCvyBX9AtyxQvjI3A9vRD7otKwNyoN5xNycfaW8fLWb1cR4mGcFPnucC/08HeCQsr58sh6CIKAuMwiHI/Lxt9lYS4tX21Wxlkpw93hxonB7+roAVs5/waofWGwa4iG9HlrrrKNcOjQIdx333145JFHABjD2PXr1+tsTQsPD8eePXvMtp06dcrsdu/evXH58mV06NChxv3I5XLo9fpG1p4qE4lECPNyQJiXA54e2QGpeaX444ox5B2LyURsRhE+jozFx5GxkEvF6OnnjH7BLugX5Io+gS5wsJG19ksgqrfyIPd3bEWQSy8wD3JyiRg9A5wxMMQNg0Lc0C/IBVIONKJ2jMGuHejQoQO+//57HD16FC4uLli7di1SU1PrDHZz587F2rVrsXTpUsyZMwfnzp0zzYVX3pK3dOlSDBw4EE8//TSeeOIJ2NnZ4cqVK9i3bx/WrVsHwDiP3cGDB/Hggw9CoVDA3d29WV9ve+LtZINHBgbikYGByC/Vlk2KnIajNzKRVaTBifhsnIjPBhADsQiIUDmiX5Ar+ge7om+QC/vnkUURBAGxmUVlIc4Y5jKqCXK9yoLcwBA39Apw5kouRJUw2LUDL730EuLi4jB27FgolUo8+eSTmDJlCvLy8mp9XHBwML777jssXrwYH3zwAQYNGoQVK1bgqaeegkKhAAB0794dkZGRWLFiBYYNGwZBEBAaGorp06eb9vPqq69i7ty5CA0NhVqtrna6E7pzjjYyTOrhg0k9fEwfkCfjjMHuVHwObmUX43JyPi4n52Pr0XgAQJCb0niaN9gV/YNcEeim5Nx51GIEQUBMRpGpNe7v2GxkFt4W5KRi9A5wxoBgBjmi+rCIUbHr16/H22+/jZSUFHTp0gXvv/8+hg0bVmN5tVqNV199FV9++SVSU1Ph5+eHFStW4LHHHqvX81n7PHbN6Y033sDGjRuRkJDQpPvl+978UvNKcTI+Gyfjs3EiLhvRaQVVBl57OCjQP8gV/YJc0C/YFeHejpBwTi9qIsYgV4hjZa1xx2sJcuUtcj39GeSI2tSo2O3bt2PhwoVYv349hgwZgo8//hjjxo1DVFRUjZ3vp02bhrS0NGzatAkdOnRAeno6dDpdtWXpzqxfvx79+vWDm5sbjhw5grfffhsLFixo7WpRI3g72Zha9AAgr0SLMzdzcCI+GyfjsnEhMQ8ZBWrsvpiC3ReNK5k4KKToHeiC/sHGwRvd/Zz4IUv1JggCbqQXmlrjjsdlma2/ChiDXJ8Al7Ig54oeDHJEd6TVW+wGDBiA3r17Y8OGDaZtERERmDJlClavXl2l/G+//YYHH3wQsbGxcHV1bdRzssWu/p5//nls374d2dnZCAgIwKOPPorly5dDKm3a7wR831tfqVaP8wm5Za16OTh9M8ds3UzA2L+ph7+T6fRtn0AXOHJABpURBAHXy4Lc8bJWucpzMAKAQipGn0BjkBsQzCBHVB8NabFr1WCn0WigVCqxY8cO3H///abtzz33HM6dO1dlKSoAmD9/Pq5du4a+ffviiy++gJ2dHSZPnozXXnut3uuhMthZHr7vlkdvEHAlJb/S6ducKqfNRCIgwtvRdOq2f5ArPB15/Kyd3iAgObcENzIKEZNeiNjMIsSkF+JaWgFyirVmZSsHuYEhbpyGh6gR2syp2MzMTOj1enh5eZlt9/LyQmpqarWPiY2NxeHDh2FjY4Ndu3YhMzMT8+fPR3Z2NjZv3lztY8on5C2Xn5/fdC+CyEpJxCJ09XVCV18nzB4SbJoM9mTZqduT8dmIzypGVEo+olLy8dmxmwCAQDclOqscEeCmRJCbHQJdlQhwU0LlZMv+em1MoVqH2IxCxGQUIjajyPQzNrMIGp2h2sfYyMqCXLAbBoa6obsfgxxRS2r1PnYAqp0Et6aReQaDASKRCNu2bYOTkxMAYO3atZg6dSo++uijalvtVq9ejVWrVjV9xYnaEZFIhGB3OwS722FaX38AQHp+KU7G55ha9a6k5ONmVjFuZlVdTk8uEcPP1RZBbnYIcFUisCz4Bbgp4ediyw//VmIwCEjOK0FMRlGVEHf75L+VySViBLvbIdTTDiHu9qaf4SoHHkuiVtSqwc7d3R0SiaRK61x6enqVVrxyKpUKvr6+plAHGPvkCYKAxMREdOzYscpjli9fjkWLFplu5+fnw9/fv9a6GQzVfxul5sH3u23ydLTBhO4qTOiuAgDkl2px9lYuYtILcSu7GDezinAzqxgJOcXQ6A3G1p6Moir7EYkAHydbBLopEehmZ/zpWnHdTmER30HbtCK1DnGZxsBWEeKKEJdZiFJtzX9/7vYKhHrYIcTDHqEedgj1sEeohz18XdgCS2SJWvW/pVwuR58+fbBv3z6zPnb79u3DfffdV+1jhgwZgh07dqCwsBD29sZ1Ma9duwaxWAw/P79qH6NQKEzzrtWnTmKxGMnJyfDw8IBcLue8Xs1IEARoNBpkZGRALBZDLudajm2Zo40Mw8M8MDzMw2y73iAgJa/E1JpXHvhuloW/Yo0eSbklSMotwdGYrCr7dbeXG0Ne2Wnd8pa+IDc7uChl/BstIwgCUvJKK502rQhxyXmlNT5OJhEhyM0OIWXBrTzEhXjYw8mWg2OI2pJWHxW7fft2PProo9i4cSMGDRqETz75BJ9++ikuX76MwMBALF++HElJSfj8888BAIWFhYiIiMDAgQOxatUqZGZm4vHHH8fw4cPx6aef1us56+qEqNFokJKSguLiqqeTqHkolUqoVCoGu3ZIEARkFmqqhL2bWcW4lV1stqh7dRwUUrOw5++ihIONFLYyCWzlEtjIJKbrykq3ZRKRRQVCQRCg1hmQX6pFQakOBaU65JeUX6/4mV9+X6m20vaKMjpDzf/S3ezklcJbRYjzd7HlMlxEFqzNDJ4AgOnTpyMrKwuvvvoqUlJS0LVrV+zZsweBgYEAgJSUFNy6dctU3t7eHvv27cMzzzyDvn37ws3NDdOmTcPrr7/eZHWSy+UICAiATqfjGqctQCKRQCqVWtSHLLUckUgEDwcFPBwU6BtUdQqj/FItbmUVI7487JVdv5VdjJS8UhSodaYVNRpCIhbBVlYW9OTisvAnha1MXDUUyspCobzi9u2hsfynCDALW+VBLf+2gFYRzirKavV3/j1bKhYhwE1pFt5CPYz931zs+MWJyNq1eotda2hI8iUiy1Wq1SMh23h6tzzsJeWUoEijQ4nWgFKNHiVaPYo1epRq9SjW6FBLg5ZFEIkAe4UUjjYyONhI4WBT+frtP433OdpWbHO3V0DG1jciq9KmWuyIiBrLRiZBRy8HdPRyqFd5QRCg1Qso0RqDXklZ8CvR6lGqMQZA0+1q7jdeN6BEozNe1xhvVy5rMAhm4cvRtiKEVQ5mjmahraKMnVwKMQclEFEjMdgRUbshEokgl4ogl4o5KICIrBLb64mIiIisRLtssSvvVsgVKIiIiMjSleeV+gyLaJfBrqCgAADqnKSYiIiIyFIUFBSYLdBQnXY5KtZgMCA5ORkODg7NOsVG+QoXCQkJHH1rIXhMLA+PiWXicbE8PCaWqSWOiyAIKCgogI+PD8Ti2nvRtcsWu9pWqWgOjo6O/CO0MDwmlofHxDLxuFgeHhPL1NzHpa6WunIcPEFERERkJRjsiIiIiKwEg10zUigUWLlyJRQKRWtXhcrwmFgeHhPLxONieXhMLJOlHZd2OXiCiIiIyBqxxY6IiIjISjDYEREREVkJBjsiIiIiK8FgR0RERGQlGOzu0Pr16xEcHAwbGxv06dMHhw4dqrV8ZGQk+vTpAxsbG4SEhGDjxo0tVNP2oyHHZOfOnbjnnnvg4eEBR0dHDBo0CL///nsL1rZ9aOjfSbkjR45AKpWiZ8+ezVvBdqqhx0WtVmPFihUIDAyEQqFAaGgoNm/e3EK1bR8aeky2bduGHj16QKlUQqVSYfbs2cjKymqh2lq/gwcPYtKkSfDx8YFIJMIPP/xQ52Na/XNeoEb75ptvBJlMJnz66adCVFSU8Nxzzwl2dnbCzZs3qy0fGxsrKJVK4bnnnhOioqKETz/9VJDJZMJ3333XwjW3Xg09Js8995zw1ltvCSdOnBCuXbsmLF++XJDJZMKZM2dauObWq6HHpFxubq4QEhIijBkzRujRo0fLVLYdacxxmTx5sjBgwABh3759QlxcnHD8+HHhyJEjLVhr69bQY3Lo0CFBLBYLH3zwgRAbGyscOnRI6NKlizBlypQWrrn12rNnj7BixQrh+++/FwAIu3btqrW8JXzOM9jdgf79+wvz5s0z2xYeHi4sW7as2vL/+c9/hPDwcLNtc+fOFQYOHNhsdWxvGnpMqtO5c2dh1apVTV21dquxx2T69OnCiy++KKxcuZLBrhk09Lj8+uuvgpOTk5CVldUS1WuXGnpM3n77bSEkJMRs2//+9z/Bz8+v2erYntUn2FnC5zxPxTaSRqPB6dOnMWbMGLPtY8aMwdGjR6t9zLFjx6qUHzt2LE6dOgWtVttsdW0vGnNMbmcwGFBQUABXV9fmqGK709hjsmXLFsTExGDlypXNXcV2qTHH5aeffkLfvn2xZs0a+Pr6IiwsDEuWLEFJSUlLVNnqNeaYDB48GImJidizZw8EQUBaWhq+++47TJgwoSWqTNWwhM95aYs8ixXKzMyEXq+Hl5eX2XYvLy+kpqZW+5jU1NRqy+t0OmRmZkKlUjVbfduDxhyT27377rsoKirCtGnTmqOK7U5jjsn169exbNkyHDp0CFIp/0U1h8Ycl9jYWBw+fBg2NjbYtWsXMjMzMX/+fGRnZ7OfXRNozDEZPHgwtm3bhunTp6O0tBQ6nQ6TJ0/GunXrWqLKVA1L+Jxni90dEolEZrcFQaiyra7y1W2nxmvoMSn39ddf45VXXsH27dvh6enZXNVrl+p7TPR6PR566CGsWrUKYWFhLVW9dqshfysGgwEikQjbtm1D//79MX78eKxduxZbt25lq10TasgxiYqKwrPPPouXX34Zp0+fxm+//Ya4uDjMmzevJapKNWjtz3l+HW4kd3d3SCSSKt+k0tPTq6T1ct7e3tWWl0qlcHNza7a6theNOSbltm/fjjlz5mDHjh0YPXp0c1azXWnoMSkoKMCpU6dw9uxZLFiwAIAxUAiCAKlUir1792LUqFEtUndr1pi/FZVKBV9fXzg5OZm2RUREQBAEJCYmomPHjs1aZ2vXmGOyevVqDBkyBP/+978BAN27d4ednR2GDRuG119/nWeBWoElfM6zxa6R5HI5+vTpg3379plt37dvHwYPHlztYwYNGlSl/N69e9G3b1/IZLJmq2t70ZhjAhhb6mbNmoWvvvqKfVOaWEOPiaOjIy5evIhz586ZLvPmzUOnTp1w7tw5DBgwoKWqbtUa87cyZMgQJCcno7Cw0LTt2rVrEIvF8PPza9b6tgeNOSbFxcUQi80/xiUSCYCKViJqWRbxOd9iwzSsUPnQ9E2bNglRUVHCwoULBTs7OyE+Pl4QBEFYtmyZ8Oijj5rKlw+Dfv7554WoqChh06ZNnO6kiTX0mHz11VeCVCoVPvroIyElJcV0yc3Nba2XYHUaekxux1GxzaOhx6WgoEDw8/MTpk6dKly+fFmIjIwUOnbsKDz++OOt9RKsTkOPyZYtWwSpVCqsX79eiImJEQ4fPiz07dtX6N+/f2u9BKtTUFAgnD17Vjh79qwAQFi7dq1w9uxZ0xQ0lvg5z2B3hz766CMhMDBQkMvlQu/evYXIyEjTfTNnzhSGDx9uVv7AgQNCr169BLlcLgQFBQkbNmxo4Rpbv4Yck+HDhwsAqlxmzpzZ8hW3Yg39O6mMwa75NPS4XLlyRRg9erRga2sr+Pn5CYsWLRKKi4tbuNbWraHH5H//+5/QuXNnwdbWVlCpVMLDDz8sJCYmtnCtrdf+/ftr/YywxM95kSCwvZaIiIjIGrCPHREREZGVYLAjIiIishIMdkRERERWgsGOiIiIyEow2BERERFZCQY7IiIiIivBYEdERERkJRjsiIiIiKwEgx0RERGRlWCwIyIiIrISDHZEREREVoLBjoiIiMhKMNgRERERWQkGOyIiIiIrwWBHREREZCUY7IiIiIisBIMdERERkZVgsCOiNmHr1q0QiUSmi1QqhZ+fH2bPno2kpKTWrl6zSU5OxiuvvIJz5861dlWIqA2QtnYFiIgaYsuWLQgPD0dJSQkOHjyI1atXIzIyEhcvXoSdnV1rV6/JJScnY9WqVQgKCkLPnj1buzpEZOEY7IioTenatSv69u0LABg5ciT0ej1ee+01/PDDD3j44Ycbvd+SkhLY2to2VTUtXklJCWxsbCASiVq7KkTUhHgqlojatIEDBwIAbt68iVWrVmHAgAFwdXWFo6MjevfujU2bNkEQBLPHBAUFYeLEidi5cyd69eoFGxsbrFq1CgDw0Ucf4a677oKnpyfs7OzQrVs3rFmzBlqt1mwfI0aMQNeuXXHs2DEMHjwYtra2CAoKwpYtWwAAu3fvRu/evaFUKtGtWzf89ttvVep+/fp1PPTQQ/D09IRCoUBERAQ++ugj0/0HDhxAv379AACzZ882nYZ+5ZVXTGVOnTqFyZMnw9XVFTY2NujVqxe+/fZbs+cpP429d+9ePPbYY/Dw8IBSqYRarW7ku05ElootdkTUpt24cQMA4OHhgaNHj2Lu3LkICAgAAPz999945plnkJSUhJdfftnscWfOnMGVK1fw4osvIjg42HQaNyYmBg899BCCg4Mhl8tx/vx5vPHGG7h69So2b95sto/U1FTMnj0b//nPf+Dn54d169bhscceQ0JCAr777ju88MILcHJywquvvoopU6YgNjYWPj4+AICoqCgMHjwYAQEBePfdd+Ht7Y3ff/8dzz77LDIzM7Fy5Ur07t0bW7ZswezZs/Hiiy9iwoQJAAA/Pz8AwP79+3HvvfdiwIAB2LhxI5ycnPDNN99g+vTpKC4uxqxZs8zq+9hjj2HChAn44osvUFRUBJlM1rQHg4han0BE1AZs2bJFACD8/fffglarFQoKCoRffvlF8PDwEBwcHITU1FSz8nq9XtBqtcKrr74quLm5CQaDwXRfYGCgIJFIhOjo6Fqfs3wfn3/+uSCRSITs7GzTfcOHDxcACKdOnTJty8rKEiQSiWBrayskJSWZtp87d04AIPzvf/8zbRs7dqzg5+cn5OXlmT3nggULBBsbG9NznTx5UgAgbNmypUr9wsPDhV69eglardZs+8SJEwWVSiXo9Xqz927GjBm1vl4iavt4KpaI2pSBAwdCJpPBwcEBEydOhLe3N3799Vd4eXnhr7/+wujRo+Hk5ASJRAKZTIaXX34ZWVlZSE9PN9tP9+7dERYWVmX/Z8+exeTJk+Hm5mbax4wZM6DX63Ht2jWzsiqVCn369DHddnV1haenJ3r27GlqmQOAiIgIAMbTxQBQWlqKP//8E/fffz+USiV0Op3pMn78eJSWluLvv/+u9X24ceMGrl69aupXePs+UlJSEB0dbfaYf/zjH3W9vUTUxvFULBG1KZ9//jkiIiIglUrh5eUFlUoFADhx4gTGjBmDESNG4NNPP4Wfnx/kcjl++OEHvPHGGygpKTHbT/njKrt16xaGDRuGTp064YMPPkBQUBBsbGxw4sQJPP3001X24erqWmUfcrm8yna5XA7AGOgAICsrCzqdDuvWrcO6deuqfZ2ZmZm1vg9paWkAgCVLlmDJkiX12kd1r5mIrAuDHRG1KREREaZRsZV98803kMlk+OWXX2BjY2Pa/sMPP1S7n+pGg/7www8oKirCzp07ERgYaNre1HPIubi4QCKR4NFHH8XTTz9dbZng4OBa9+Hu7g4AWL58OR544IFqy3Tq1MnsNkfAElk/BjsisgrlkxZLJBLTtpKSEnzxxRcN2gcAKBQK0zZBEPDpp582XUUBKJVKjBw5EmfPnkX37t1NLXrVKa/L7a2FnTp1QseOHXH+/Hm8+eabTVo/Imq7GOyIyCpMmDABa9euxUMPPYQnn3wSWVlZeOedd8xCWl3uueceyOVy/Otf/8J//vMflJaWYsOGDcjJyWny+n7wwQcYOnQohg0bhqeeegpBQUEoKCjAjRs38PPPP+Ovv/4CAISGhsLW1hbbtm1DREQE7O3t4ePjAx8fH3z88ccYN24cxo4di1mzZsHX1xfZ2dm4cuUKzpw5gx07djR5vYnIsnHwBBFZhVGjRmHz5s24ePEiJk2ahBUrVmDq1KlYtmxZvfcRHh6O77//Hjk5OXjggQfwzDPPoGfPnvjf//7X5PXt3Lkzzpw5g65du+LFF1/EmDFjMGfOHHz33Xe4++67TeWUSiU2b96MrKwsjBkzBv369cMnn3wCwDhB84kTJ+Ds7IyFCxdi9OjReOqpp/DHH39g9OjRTV5nIrJ8IkG4beZOIiIiImqT2GJHREREZCUY7IiIiIisBIMdERERkZVgsCMiIiKyEgx2RERERFaCwY6IiIjISrTLCYoNBgOSk5Ph4ODAJXaIiIjIogmCgIKCAvj4+EAsrr1Nrl0Gu+TkZPj7+7d2NYiIiIjqLSEhAX5+frWWaZfBzsHBAYDxDXJ0dGzl2hARERHVLD8/H/7+/qb8Upt2GezKT786Ojoy2BEREVGbUJ/uY+0y2BEREbV1pVo90vJLkZavRmahGrYyCZyUMjjbyuCslMPRRgqphGMk2xsGOyIiIgtiMAjILtYgNa8UafmlSC0Lb2l55deNP3OLtXXuy0EhNYY9pQzOtnI42coqhT/jNsfy62W3nZUy2MgkLfBKqTkw2BEREbWQUq0eqZUDWqXraflqpOaVIr2gFFq9UK/92cjE8Ha0gbu9AmqdAbklGuQWa1FQqgMAFKh1KFDrkJhT0qB6yqVis/DnpJTBybYiEDop5XC2LdtWVsbVXg47uYSzTbQyBjsiIqI7ZDAIyCxSI70snFUX3FLzSpFfFrjqIhIBbnYKeDsp4O1oAy9HG+NPJxuz24620mqDlE5vQH6pDnklWuQWa5BbokVecaXr5bdvuz+vRAudQYBGZ0B6gRrpBeoGvQ9yqRiuSjlc7ep3cVHKIREzCDYlBjsiIqJG0BsEHI/Lws/nU/DbpRTk1OPUKADYyiTwdrKBl6PCLKxVvu7hoIDsDvrHSSViU3gC7Or9OEEQUKTRG8NeWdAzhkMtcks0ZeGw7HrZ9rwSLbKLNFDrDNDoDEgtO1VcHyIR4GQrM9a1LBC62RsDX+UA6GangIudDG52CtjKeZq4Ngx2RERE9WQwCDibkIOfz6dg98UUZFRq0RKJAA97BbydbODpYGPe2uZUEdwcFNW3slkCkUgEe4UU9gop/Fwa9tgSjR5ZRWpkF2lqvxQbf+YWayEIMAbFYi1iUVSv57GRieFmpzC2+NnJ4WYnh41MAkEQYBAEGATAIAgQBJRtq7hd+afBdL8AATC7bTAAAiqXqX6f5fcPDHHDK5O7NPwNbwYMdkRERLUQBAGXk/Px8/lk/HIhBUm5Ff3VnGxlGNfVG5N6+KB/sOsdtbK1dbZyCfzkSvi5KOtVXqc3ILesta++F43egFKtAUm5JWbHobX5udi2dhVMGOyIiIiqcT2tAD+fT8bPF1IQl1nRmmSvkGJMZy9M6uGDIR3cIZe23zB3J6QSMdztFXC3V9SrfPlp4uzC8lY/NbKLtMguUkOtNUAsFkEkAsQiEUQo+1l2WywytkZW/Cy/r9JtAGJx+ePKyqK6x1bsu/y28ZS3ZWCwIyIiKhOfWYRfLhhb5q6mFpi2K6RijI7wwqQeKozo5MnpQFpB5dPEAW71axVsjxjsiIioXUvOLcHuCyn4+UIyLiTmmbbLJCIMD/PApB4+uDvCC/YKfmSS5eNvKRERtTsZBWr8eikFP59Pxsn4HNN2iViEwaFumNTDB2M7e8NJKWvFWhI1HIMdERG1C7nFGvx+ORU/n0/B0ZhMGMrmABaJgH5BrpjUwwfjunrXu88XkSVisCMiIqtVqNZhX5QxzB26nmG2okMPf2dM6q7ChO4qqJwsZ1Qj0Z1gsCMiIqtSqtXjr6vp+Pl8Mv66mg61zmC6L9zbAZN6+GBSdx92wCerxGBHRERtnkZnwKHrGfj5fDL2RaWhSKM33RfibmcMcz1U6ODp0Iq1JGp+DHZERNRmlWr1+PZUAjYeiEFyXsUyVr7OtqYw11nlaLErPRA1tWabVTE+Ph5z5sxBcHAwbG1tERoaipUrV0Kj0dT6uLS0NMyaNQs+Pj5QKpW49957cf36dbMyarUazzzzDNzd3WFnZ4fJkycjMTGxuV4KERFZmBKNHpsPx+GuNfvx8o+XkZxXCnd7BWYPCcLO+YNxeOlILBsXji4+Tgx11K40W4vd1atXYTAY8PHHH6NDhw64dOkSnnjiCRQVFeGdd96p9jGCIGDKlCmQyWT48ccf4ejoiLVr12L06NGIioqCnZ1xIeOFCxfi559/xjfffAM3NzcsXrwYEydOxOnTpyGRcNJIIiJrVaTWYdvxm/jkYCwyC40NBSonG8wfEYp/9vXnxMHU7okEQRDqLtY03n77bWzYsAGxsbHV3n/t2jV06tQJly5dQpcuxsV09Xo9PD098dZbb+Hxxx9HXl4ePDw88MUXX2D69OkAgOTkZPj7+2PPnj0YO3ZsnfXIz8+Hk5MT8vLy4Ojo2HQvkIiImkVBqRafH7uJTYfjkF1kDHR+LraYP6ID/tHHFwopAx1Zr4bklhbtY5eXlwdXV9ca71er1QAAGxsb0zaJRAK5XI7Dhw/j8ccfx+nTp6HVajFmzBhTGR8fH3Tt2hVHjx6tNtip1WrTvgHjG0RERJYvr0SLrUfisflIHPJKtACAIDcl5o/sgPt7+UIm4TqtRJW1WLCLiYnBunXr8O6779ZYJjw8HIGBgVi+fDk+/vhj2NnZYe3atUhNTUVKSgoAIDU1FXK5HC4uLmaP9fLyQmpqarX7Xb16NVatWtV0L4aIiJpVTpEGm4/EYeuReBSodQCAEA87PDOqAyZ194GUgY6oWg3+y3jllVcgEolqvZw6dcrsMcnJybj33nvxz3/+E48//niN+5bJZPj+++9x7do1uLq6QqlU4sCBAxg3blydfecEQaixg+zy5cuRl5dnuiQkJDT0ZRMRUQvIKlTjrd+uYuhbf2HdXzdQoNYhzMse//tXL+x7fjju7+XHUEdUiwa32C1YsAAPPvhgrWWCgoJM15OTkzFy5EgMGjQIn3zySZ3779OnD86dO4e8vDxoNBp4eHhgwIAB6Nu3LwDA29sbGo0GOTk5Zq126enpGDx4cLX7VCgUUCi4RAwRkaVKLyjFpwdj8eXft1CiNc5BF6FyxLOjOmBsF2+IxRzZSlQfDQ527u7ucHd3r1fZpKQkjBw5En369MGWLVsgFtf/W5aTkxMA4Pr16zh16hRee+01AMbgJ5PJsG/fPkybNg0AkJKSgkuXLmHNmjUNfDVERNSaUvNKsTEyBl+fuGVaIaKbrxOevbsjRkd4cqoSogZqtj52ycnJGDFiBAICAvDOO+8gIyPDdJ+3t7fpenh4OFavXo37778fALBjxw54eHggICAAFy9exHPPPYcpU6aYBks4OTlhzpw5WLx4Mdzc3ODq6oolS5agW7duGD16dHO9HCIiakJJuSXYcOAGvj2ZCI3eGOh6BTjj2bs7YkSYBwMdUSM1W7Dbu3cvbty4gRs3bsDPz8/svsozrERHRyMvL890OyUlBYsWLUJaWhpUKhVmzJiBl156yezx7733HqRSKaZNm4aSkhLcfffd2Lp1K+ewIyKycLeyirH+wA18fyYRWr3xs6B/kCuevbsjhnRwY6AjukMtOo+dpeA8dkRELSsuswgf7b+BXWeToDcYP3YGh7rh2bs7YmCIWyvXjsiyWew8dkRE1L7cSC/Ah3/dwE/nk1GW5zCsozuevbsj+gXVPK8pETUOgx0RETW5q6n5WPfXDey5mILy80Kjwj3xzKgO6BXgUvuDiajRGOyIiKjJXErKw7q/ruP3y2mmbWM6e+GZUR3Rzc+pFWtG1D4w2BER0R1R6/Q4eC0T35y4hT+vpgMARCJgfFcVFozqgAgV+zITtRQGOyIiajCNzoDDNzLwy/kU7ItKMy37JRYBk3r4YMHIDujo5dDKtSRqfxjsiIioXrR6A47cyMQvF1Kw93Iq8kt1pvu8HW0wvpsKDw8MQKiHfSvWkqh9Y7AjIqIa6fQGHIvNwi/nU/B7VCpyi7Wm+zwcFJjQTYWJ3VXoHeDCZb+ILACDHRERmdHpDTgRl42fL6Tg98upyC7SmO5zt5djXFdjmOsb5AoJwxyRRWGwIyIi6A0CTsRlY/fFZPx2KRWZhRVhztVOjnFdvTGhuwoDgt0Y5ogsGIMdEVE7ZTAIOHUzB7svJGPPpVRkFKhN9zkrZcYw180HA0NcIZWIW7GmRFRfDHZERO2IwSDgbEIOfrmQgj0XU5CWXxHmHG2kuLerNyZ098HgUDfIGOaI2hwGOyIiKycIAs4l5GJ3WZhLzis13edgI8WYzt6Y2F2FIR3cIZcyzBG1ZQx2RERWSBAEXEzKw+4LKfjlQgqScktM99krpLinsxcmdFNhWJg7FFJJK9aUiJoSgx0RkZUQBAGXk/Ox+2IKdl9Iwa3sYtN9SrkEoyO8MLG7CneFecBGxjBHZI0Y7IiI2iBBEJCWr8bl5DxcTs7H5eQ8XErKN2uZs5VJMCrCE5O6qzCikyfDHFE7wGBHRGThDAYBcVlFpgAXlZyPqOR8ZFWaX66cQirGqHBPTOzug5HhHlDK+W+eqD3hXzwRkQVR6/S4llpoDHAp+bicnI8rKfko1uirlJWIRQj1sEMXHyd08XFEZx9H9PBzhp2C/9qJ2iv+9RMRtZL8Ui2uJOeXtcQZW+NupBdCZxCqlLWRiRHu7YguPo6mINfJ24GnV4nIDIMdEVELSM8vNYW38iBXeXBDZc5KmSnAdVYZw1ywux0nCSaiOjHYERE1IYNBwM3sYrMAF5Wcj8xCdbXlfZxs0LmsBa6LjyO6+DrBx8kGIhGX7SKihmOwIyK6A9lFGpy5mYMzt3Jw+mYOLifno1Ctq1JOLAJCPOxNAa6zyhjmXOzkrVBrIrJWDHZERPWkNwiITi3AmVvGIHfmZg7is6qeTpVLxYjwdkBnH0dTa1yEtyNs5ewPR0TNq9mCXXx8PF577TX89ddfSE1NhY+PDx555BGsWLECcnnN31DT0tKwdOlS7N27F7m5ubjrrruwbt06dOzY0VRmxIgRiIyMNHvc9OnT8c033zTXyyGidii3WIOzt3JNQe7crVwUVTM6NdTDDr0DXNA70AW9ApzRwcOe/eGIqFU0W7C7evUqDAYDPv74Y3To0AGXLl3CE088gaKiIrzzzjvVPkYQBEyZMgUymQw//vgjHB0dsXbtWowePRpRUVGws7MzlX3iiSfw6quvmm7b2to210shonbAYBBwPb3Q1BJ3+lYOYjOKqpSzk0vQM8AZfQJc0CvQBb38neGs5OlUIrIMzRbs7r33Xtx7772m2yEhIYiOjsaGDRtqDHbXr1/H33//jUuXLqFLly4AgPXr18PT0xNff/01Hn/8cVNZpVIJb2/v5qo+UaNcTyvA7ospOHMrF/YKCZyVcrgq5XBWyuBqJ4eLUg4XOzlclDK42MnhoJCyk3wrySvR4lxCrql/3LlbuSiopm9ciLsdegW4oHegM3oHuCDMywESMY8ZEVmmFu1jl5eXB1dX1xrvV6uNo8ZsbGxM2yQSCeRyOQ4fPmwW7LZt24Yvv/wSXl5eGDduHFauXAkHB4ca91u+bwDIz8+/05dCZHItrQC7L6Rgz8UUXE8vbNBjpWKRMfzZyUwh0MVOZgyAt4VAl7L7HWykEDNYNIjBICA2sxBnbuaaBjncyCiEcNt0cUq5BD38nE0hrleAC1w5uIGI2pAWC3YxMTFYt24d3n333RrLhIeHIzAwEMuXL8fHH38MOzs7rF27FqmpqUhJSTGVe/jhhxEcHAxvb29cunQJy5cvx/nz57Fv375q97t69WqsWrWqyV8TtV81hTmZRIRhHT0wMtwTBoOA7CINcos1yC7WGn8WaZBbrEV2kQYlWj10BgGZheoap8KojlgE89BXTQh0tpXByVYGx0o/7eSSdtM6WFCqxfmEPFPfuLO3cpFXoq1SLtBNaewbF+CMXgEuCPd2YN84ImrTRIJw+3fW2r3yyit1hqSTJ0+ib9++ptvJyckYPnw4hg8fjv/7v/+r9bGnT5/GnDlzcP78eUgkEowePRpisfEf7Z49e2p8TN++fXH69Gn07t27yv3Vtdj5+/sjLy8Pjo6OtdaHCDD2/7yWVojdF41h7kalMCeXiDGsozvGd1NhdGcvONnK6rXPUq0eOcUa5BRpjT+LNcgp0iC78u1iLXKKKu6rruN+fUnEIjjaSOFoK4OjTXngkxp/2hjDn/E+qSkMVi6nkDb/iE5BEFCqNaBIo0OxWm/8qdGhSK03/6nRo1hd9vO2+7OK1LieXrU1zkYmRnc/Z1OQ6x3oAnd7RbO/JiKiO5Wfnw8nJ6d65ZYGB7vMzExkZmbWWiYoKMh0OjU5ORkjR47EgAEDsHXrVlNIq0teXh40Gg08PDwwYMAA9O3bFx999FG1ZQVBgEKhwBdffIHp06fXue+GvEHUfgmCgOi0Auy5kILdF1MQU6kjvVwixl1hFWHO0aZ+Ye5OqXV65BYbg1/l1j9ja6B5ICwo0SK/VIu8Ei20+gb9mVdLIRWbtwKWhcSKYFgREvWCUCmY6VGkvu2nWXAz396w/0g183OxRe8AF/QJdEHvABeEqxwgY2scEbVBDcktDT4V6+7uDnd393qVTUpKwsiRI9GnTx9s2bKl3qEOAJycnAAYB1ScOnUKr732Wo1lL1++DK1WC5VKVe/9E1VHEARcTS3AnovGMBdrAWGuMoVUAi9HCbwcbeouXKa8Faw85OVXCnz5Jbqat5UatxeojWFLrTMgvUCN9IL6nza+E0q5BEq5FHaKsp9yCZSKsp81bVdI4WAjRRcfR3g61P89IiKyFg1usauv8tOvAQEB+PzzzyGRVJzGqTyaNTw8HKtXr8b9998PANixYwc8PDwQEBCAixcv4rnnnkOfPn3w/fffAzD21du2bRvGjx8Pd3d3REVFYfHixbC1tcXJkyfNnqcmbLGjysrDXHmfudjMSmFOKsZdHT0wobs37o5onTDX2vQGAYVqHfJLzINh5QBYsd24TSIW1RnE7BVSKOUS2JX/lEuhVBh/2sokHCBCRFSmWVvs6mvv3r24ceMGbty4AT8/P7P7KmfJ6Oho5OXlmW6npKRg0aJFSEtLg0qlwowZM/DSSy+Z7pfL5fjzzz/xwQcfoLCwEP7+/pgwYQJWrlxZr1BHBBh/B6+kGFvmqgtzw8M8MKGbCndHeMKhHYa5yiRiEZzKTrn6t3ZliIioVs3WYmfJ2GLXPgmCgKiU/LIwl4q428LciDAPTOiuwqhwhjkiIrIcFtFiR2QJBEHA5eR8U8tc5XU9FVIxRnTywPhuKtwd4QV7Bf8ciIiobeMnGVmly8l5+OVCCn6tJsyN7OSJ8WUtcwxzRERkTfipRlZDqzfgt0up2HIkDmdu5Zq2K6RijAr3xPhuxjBnxzBHRERWip9w1OblFGnw9clb+OLYTaTklQIwrgAxOsKLYY6IiNoVftpRm3UtrQBbjsRh19kklGoNAAB3ezkeHhCIhwcGcB4zIiJqdxjsqE0xGATsj07HliPxOHyjYgWULj6OmD0kGJN6qFpk6SsiIiJLxGBHbUKhWocdpxLw2dF402AIsQgY28Ubs4cEo1+QS7tZ4J6IiKgmDHZk0W5lFWPr0XjsOJWAArUOAOBoI8WD/QPw6MBA+LsqW7mGREREloPBjiyOIAg4FpOFzUfi8efVNNOi8KEedpg1JBj/6O0LpZy/ukRERLfjpyNZjFKtHj+eS8KWI/G4mlpg2j48zAOzhwThro4eXD+UiIioFgx21OpS80rxxd/x+Or4LeQUawEAtjIJpvbxw8zBQejgad/KNSQiImobGOyo1Zy9lYMtR+Kx52IKdAbj+VZfZ1vMHByI6X0D4KTkeq1EREQNwWBHLUqrN+DXstUhzlZaHaJ/kCseGxqE0RFekErErVdBIiKiNozBjlpEdpEGX58wrg6Rmm9cHUIuEWNSDx/MHhKErr5OrVxDIiKito/BjppVdGrF6hBqXfnqEAo8MjAADw8IhIeDopVrSEREZD0Y7KjJ5ZVocfBaBr4+cQtHY7JM27v6OuKxIcGY0J2rQxARETUHBju6Y4IgIDqtAPuvZmB/dDpO38yBvmwwhFgE3NvVuDpE30CuDkFERNScGOyoUYrUOhyNycL+6HQcuJqO5LxSs/tDPewwpos3Hh4QAD8Xrg5BRETUEhjsqN5iMwqxPzoDB6LTcTw2Gxq9wXSfQirG4FA3jAz3xIgwTwS4McwRERG1NAY7qlGpVo/jcdnYfzUdB6LTEZ9VbHa/v6stRnXyxIhwTwwKcYONjP3miIiIWhODHZlJzCk2tspdTceRmEyUaita5WQSEfoHu2JkJ0+M6OSJUA879pkjIiKyIAx27ZxWb8Cp+Bzsj07H/qvpuJ5eaHa/t6MNRoZ7YEQnTwzp4A57BX9liIiILFWzfkpPnjwZ586dQ3p6OlxcXDB69Gi89dZb8PHxqfExgiBg1apV+OSTT5CTk4MBAwbgo48+QpcuXUxl1Go1lixZgq+//holJSW4++67sX79evj5+TXny7Ea6fmlOBBtHMF66HomCtU6030SsQh9AlwwItwDIzt5Itzbga1yREREbYRIEAShuXb+3nvvYdCgQVCpVEhKSsKSJUsAAEePHq3xMW+99RbeeOMNbN26FWFhYXj99ddx8OBBREdHw8HBAQDw1FNP4eeff8bWrVvh5uaGxYsXIzs7G6dPn4ZEUnc/r/z8fDg5OSEvLw+Ojo5N82ItmN4g4FxCjmk6ksvJ+Wb3u9nJMbyTMcjd1dGDa7QSERFZkIbklmYNdrf76aefMGXKFKjVashkVcODIAjw8fHBwoULsXTpUgDG1jkvLy+89dZbmDt3LvLy8uDh4YEvvvgC06dPBwAkJyfD398fe/bswdixY+ushzUHO71BQFahGukFalxPL8CB6AxEXstAbrHWVEYkArr7OWNkWZjr5usEsZitckRERJaoIbmlxTpMZWdnY9u2bRg8eHC1oQ4A4uLikJqaijFjxpi2KRQKDB8+HEePHsXcuXNx+vRpaLVaszI+Pj7o2rUrjh49Wm2wU6vVUKvVptv5+flVylg6nd6ArCIN0vJLkZ6vRlqB8Wd6gfntzEI1DNVEdUcbKe4K88CocE/cFeYBd3su5UVERGRtmj3YLV26FB9++CGKi4sxcOBA/PLLLzWWTU1NBQB4eXmZbffy8sLNmzdNZeRyOVxcXKqUKX/87VavXo1Vq1bdyctoNlq9AZmF6rKQpjYGtwI10st+lt/OqiGwVUcsMq7H6uNsa5pbrpe/M6QScfO+GCIiImpVDQ52r7zySp0h6eTJk+jbty8A4N///jfmzJmDmzdvYtWqVZgxYwZ++eWXWjvk336fIAh1duCvrczy5cuxaNEi0+38/Hz4+/vXur87pdUbkFFQTVi7rbUtq0iD+p4Ml4hF8LBXwNNRAU8Hm7KfCng52pj9dLNXQMJTq0RERO1Og4PdggUL8OCDD9ZaJigoyHTd3d0d7u7uCAsLQ0REBPz9/fH3339j0KBBVR7n7e0NwNgqp1KpTNvT09NNrXje3t7QaDTIyckxa7VLT0/H4MGDq62PQqGAQtGypx6/O52I5Tsv1qusVCyCh4MCnmXB7Paw5lH209VOzsBGRERENWpwsCsPao1RPk6jcn+3yoKDg+Ht7Y19+/ahV69eAACNRoPIyEi89dZbAIA+ffpAJpNh3759mDZtGgAgJSUFly5dwpo1axpVr+bg6aCATCKCp4MxmFUJa44KeJW1urkq5Ry8QERERHes2frYnThxAidOnMDQoUPh4uKC2NhYvPzyywgNDTVrrQsPD8fq1atx//33QyQSYeHChXjzzTfRsWNHdOzYEW+++SaUSiUeeughAICTkxPmzJmDxYsXw83NDa6urliyZAm6deuG0aNHN9fLabARnTwR/do4BjYiIiJqMc0W7GxtbbFz506sXLkSRUVFUKlUuPfee/HNN9+YnRaNjo5GXl6e6fZ//vMflJSUYP78+aYJivfu3Wuaww4wzo8nlUoxbdo00wTFW7durdccdi2Fp0yJiIiopbXoPHaWIi8vD87OzkhISLC6eeyIiIjIupQP+szNzYWTk1OtZdvlwp8FBQUA0OwjY4mIiIiaSkFBQZ3Brl222BkMBiQnJ8PBoXnXQS1P2GwZtBw8JpaHx8Qy8bhYHh4Ty9QSx0UQBBQUFMDHxwdice1z0rbLFjuxWAw/P78Wez5HR0f+EVoYHhPLw2NimXhcLA+PiWVq7uNSV0tdOS5FQERERGQlGOyIiIiIrASDXTNSKBRYuXJli696QTXjMbE8PCaWicfF8vCYWCZLOy7tcvAEERERkTViix0RERGRlWCwIyIiIrISDHZEREREVoLBjoiIiMhKMNjdofXr1yM4OBg2Njbo06cPDh06VGv5yMhI9OnTBzY2NggJCcHGjRtbqKbtR0OOyc6dO3HPPffAw8MDjo6OGDRoEH7//fcWrG370NC/k3JHjhyBVCpFz549m7eC7VRDj4tarcaKFSsQGBgIhUKB0NBQbN68uYVq2z409Jhs27YNPXr0gFKphEqlwuzZs5GVldVCtbV+Bw8exKRJk+Dj4wORSIQffvihzse0+ue8QI32zTffCDKZTPj000+FqKgo4bnnnhPs7OyEmzdvVls+NjZWUCqVwnPPPSdERUUJn376qSCTyYTvvvuuhWtuvRp6TJ577jnhrbfeEk6cOCFcu3ZNWL58uSCTyYQzZ860cM2tV0OPSbnc3FwhJCREGDNmjNCjR4+WqWw70pjjMnnyZGHAgAHCvn37hLi4OOH48ePCkSNHWrDW1q2hx+TQoUOCWCwWPvjgAyE2NlY4dOiQ0KVLF2HKlCktXHPrtWfPHmHFihXC999/LwAQdu3aVWt5S/icZ7C7A/379xfmzZtnti08PFxYtmxZteX/85//COHh4Wbb5s6dKwwcOLDZ6tjeNPSYVKdz587CqlWrmrpq7VZjj8n06dOFF198UVi5ciWDXTNo6HH59ddfBScnJyErK6slqtcuNfSYvP3220JISIjZtv/973+Cn59fs9WxPatPsLOEz3meim0kjUaD06dPY8yYMWbbx4wZg6NHj1b7mGPHjlUpP3bsWJw6dQparbbZ6tpeNOaY3M5gMKCgoACurq7NUcV2p7HHZMuWLYiJicHKlSubu4rtUmOOy08//YS+fftizZo18PX1RVhYGJYsWYKSkpKWqLLVa8wxGTx4MBITE7Fnzx4IgoC0tDR89913mDBhQktUmaphCZ/z0hZ5FiuUmZkJvV4PLy8vs+1eXl5ITU2t9jGpqanVltfpdMjMzIRKpWq2+rYHjTkmt3v33XdRVFSEadOmNUcV253GHJPr169j2bJlOHToEKRS/otqDo05LrGxsTh8+DBsbGywa9cuZGZmYv78+cjOzmY/uybQmGMyePBgbNu2DdOnT0dpaSl0Oh0mT56MdevWtUSVqRqW8DnPFrs7JBKJzG4LglBlW13lq9tOjdfQY1Lu66+/xiuvvILt27fD09OzuarXLtX3mOj1ejz00ENYtWoVwsLCWqp67VZD/lYMBgNEIhG2bduG/v37Y/z48Vi7di22bt3KVrsm1JBjEhUVhWeffRYvv/wyTp8+jd9++w1xcXGYN29eS1SVatDan/P8OtxI7u7ukEgkVb5JpaenV0nr5by9vastL5VK4ebm1mx1bS8ac0zKbd++HXPmzMGOHTswevTo5qxmu9LQY1JQUIBTp07h7NmzWLBgAQBjoBAEAVKpFHv37sWoUaNapO7WrDF/KyqVCr6+vnBycjJti4iIgCAISExMRMeOHZu1ztauMcdk9erVGDJkCP79738DALp37w47OzsMGzYMr7/+Os8CtQJL+Jxni10jyeVy9OnTB/v27TPbvm/fPgwePLjaxwwaNKhK+b1796Jv376QyWTNVtf2ojHHBDC21M2aNQtfffUV+6Y0sYYeE0dHR1y8eBHnzp0zXebNm4dOnTrh3LlzGDBgQEtV3ao15m9lyJAhSE5ORmFhoWnbtWvXIBaL4efn16z1bQ8ac0yKi4shFpt/jEskEgAVrUTUsizic77FhmlYofKh6Zs2bRKioqKEhQsXCnZ2dkJ8fLwgCIKwbNky4dFHHzWVLx8G/fzzzwtRUVHCpk2bON1JE2voMfnqq68EqVQqfPTRR0JKSorpkpub21ovweo09JjcjqNim0dDj0tBQYHg5+cnTJ06Vbh8+bIQGRkpdOzYUXj88cdb6yVYnYYeky1btghSqVRYv369EBMTIxw+fFjo27ev0L9//9Z6CVanoKBAOHv2rHD27FkBgLB27Vrh7NmzpiloLPFznsHuDn300UdCYGCgIJfLhd69ewuRkZGm+2bOnCkMHz7crPyBAweEXr16CXK5XAgKChI2bNjQwjW2fg05JsOHDxcAVLnMnDmz5StuxRr6d1IZg13zaehxuXLlijB69GjB1tZW8PPzExYtWiQUFxe3cK2tW0OPyf/+9z+hc+fOgq2traBSqYSHH35YSExMbOFaW6/9+/fX+hlhiZ/zIkFgey0RERGRNWAfOyIiIiIrwWBHREREZCUY7IiIiIisBIMdERERkZVgsCMiIiKyEgx2RERERFaCwY6IiIjISjDYEREREVkJBjsiIiIiK8FgR0RERGQlGOyIiIiIrASDHREREZGVYLAjIiIishIMdkRERERWgsGOiIiIyEow2BERERFZCQY7IiIiIivBYEdE7drx48dx//33IyAgAAqFAl5eXhg0aBAWL15sKrN+/Xps3br1jp7nzTffxA8//HBnlSUiqoNIEAShtStBRNQadu/ejcmTJ2PEiBF44oknoFKpkJKSglOnTuGbb75BYmIiAKBr165wd3fHgQMHGv1c9vb2mDp16h0HRCKi2khbuwJERK1lzZo1CA4Oxu+//w6ptOLf4YMPPog1a9a0Ys2IiBqHp2KJqN3KysqCu7u7WagrJxYb/z0GBQXh8uXLiIyMhEgkgkgkQlBQEACgtLQUixcvRs+ePeHk5ARXV1cMGjQIP/74o9m+RCIRioqK8Nlnn5n2MWLECNP9qampmDt3Lvz8/CCXyxEcHIxVq1ZBp9M122snIuvEFjsiarcGDRqE//u//8Ozzz6Lhx9+GL1794ZMJjMrs2vXLkydOhVOTk5Yv349AEChUAAA1Go1srOzsWTJEvj6+kKj0eCPP/7AAw88gC1btmDGjBkAgGPHjmHUqFEYOXIkXnrpJQCAo6MjAGOo69+/P8RiMV5++WWEhobi2LFjeP311xEfH48tW7a01NtBRFaAfeyIqN3KysrClClTcPjwYQCATCZDv379MGnSJCxYsAD29vYA6t/HTq/XQxAEzJs3D2fOnMGZM2dM99XUx27evHnYtm0bLl++jICAANP2d999F0uWLMHly5fRuXPnpnnBRGT1eCqWiNotNzc3HDp0CCdPnsR///tf3Hfffbh27RqWL1+Obt26ITMzs8597NixA0OGDIG9vT2kUilkMhk2bdqEK1eu1KsOv/zyC0aOHAkfHx/odDrTZdy4cQCAyMjIO3qNRNS+MNgRUbvXt29fLF26FDt27EBycjKef/55xMfH1zmAYufOnZg2bRp8fX3x5Zdf4tixYzh58iQee+wxlJaW1uu509LS8PPPP0Mmk5ldunTpAgD1CpdEROXYx46IqBKZTIaVK1fivffew6VLl2ot++WXXyI4OBjbt2+HSCQybVer1fV+Pnd3d3Tv3h1vvPFGtff7+PjUe19ERAx2RNRupaSkQKVSVdlefhq1PFQpFAqUlJRUKScSiSCXy81CXWpqapVRsbXtY+LEidizZw9CQ0Ph4uLS6NdCRARw8AQRtWPdu3eHn58fJk2ahPDwcBgMBpw7dw7vvvsuCgoKcPToUXTr1g2zZs3CN998g88++wwhISGwsbFBt27dsGXLFjz22GN46qmnMHXqVCQkJOC1116DWCzG9evXUfnf64gRI3DlyhX83//9H1QqFRwcHNCpUyekpKRg0KBBsLW1xbPPPotOnTqhtLQU8fHx2LNnDzZu3Ag/P79WfJeIqC1hsCOiduvbb7/Fjz/+iJMnTyIlJQVqtRoqlQrDhw/H8uXLERERAQC4efMmnnzySRw7dgwFBQUIDAxEfHw8AOCtt97Cxo0bkZKSgpCQECxatAiJiYlYtWqVWbA7f/48nn76aZw9exbFxcUYPny4aZRtZmYmXnvtNfz8889ITEyEg4MDgoODce+992L58uWws7Nr6beGiNooBjsiIiIiK8FRsURERERWgsGOiIiIyEow2BERERFZCQY7IiIiIivBYEdERERkJRjsiIiIiKxEu1x5wmAwIDk5GQ4ODmYzxhMRERFZGkEQUFBQAB8fH4jFtbfJtctgl5ycDH9//9auBhEREVG9JSQk1LkSTbsMdg4ODgCMb5Cjo2Mr14aIiIioZvn5+fD39zfll9q0y2BXfvrV0dGRwY6IiIjahPp0H2uTwW716tXYuXMnrl69CltbWwwePBhvvfUWOnXq1NpVIyJqNoIgQKM3oESjR7FGj2KNruxnxXWDAEhEIkjE5RdAXPl22XWxWASpWGR2n1hk3FZ+v0QkglgMSMVi0/Xb98N+ykSWpU0Gu8jISDz99NPo168fdDodVqxYgTFjxiAqKoqLZRNRq9MbBJRoy8KW2hi8SrQ6FN12vUSjR5FGZwpqla+XB7XyMuXXdQbLWt5bJKoIklKxCI62MjjZyuCsLPtpKzdeV1Zcd7Ytu62Uw9lWBqVcwoBI1EREgiBY1n+JRsjIyICnpyciIyNx11131Vk+Pz8fTk5OyMvL46lYIqoXQRCQV6JFWr4aafmlSMsvRXpBxfW0fDXS80uRXaxBqdbQ7PWRSURQyqVQyiWwlUtgJ5fCVi6BWAQYDIBeEKA3VFwM5bcrbTeYbgN6g6GsHKqUa25SsagiCJaFPbMgWHZf5fudlTI42MggETMQkvVrSG5pky12t8vLywMAuLq6Nul+9Xo9tFptk+6TqpLJZJBIJK1dDWqnBEFAgVqH9LJwllbpZ3pBWWAr+6nRNSywiUSAUiaBUlEWwGQS2N123VYuMStT3X2m63KJKczJJC03DanhtqCnF8pC4W3bdXoBBaU65JZokFusRW6JFnnFla6XaJFXrK24v1gLjd4AnUFAZqEGmYUaAEX1rpdIBDjayEytgK52cvi62MLfRQl/V2XZT1s42crYIkjtRpsPdoIgYNGiRRg6dCi6du1abRm1Wg21Wm26nZ+fX+c+U1NTkZub25RVpVo4OzvD29ub/3ypSRVrdOYtbOXXy1raysNciVZf7326KGXwcrSBp6MNvBwU8HK0gZejwnjb0QZudnIo5cYwppCKreJ3WiwWQQwRZE38/UsQBJRqDaagl1eiLftZEQbLb5ffV16uUK2DIMAYFku0uFnL8zgopMbAVynslYc/Pxdb2Cna/EchkUmb/21esGABLly4gMOHD9dYZvXq1Vi1alW991ke6jw9PaFUKq3iH7OlEgQBxcXFSE9PBwCoVKpWrhG1JYIgIDW/FDfSC02X+KwipOYZQ1yBWlfvfTnYSE0hzcuhLLg5VgpuDjbwcFDApqnTTTsmEolgK5fAVm4LlZNtgx6r1RvMgmBeiRbp+Wok5pQgIacYCdnFSMgpQUaB8ffgamoBrqYWVLsvVzs5/F1s4VcW/PxMIdAWvi62UEh5zKntaNN97J555hn88MMPOHjwIIKDg2ssV12Lnb+/f7XnqvV6Pa5duwZPT0+4ubk1W93JXFZWFtLT0xEWFsbTslSFTm/AzexiU3iLSS9ETEYhYjKKUFhHeLOVSeDtZAPPSq1rt7e4eToqoJS3+e+5VI1Srd4U9hLLwp4x9BUjIbsEeSW1d7cRiQAvBxv4u9rCz0VpFgD9XY2BlP38qLlZfR87QRDwzDPPYNeuXThw4ECtoQ4AFAoFFApFvfZd3qdOqVTecT2p/srfb61Wy2DXjpVo9GWBrbBKK5xWX/13UIlYhEBXJUI97dHB0x6hHvbwcapocbNXSNnq3o7ZyCToUPa7UZ38Ui0Ssyta+RJvC34lWj1S80uRml+Kk/E5VR4vFYugcraBv4sSoR726OLjiM4+jgjzcmDrLrWKNhnsnn76aXz11Vf48ccf4eDggNTUVACAk5MTbG0b1pxfE34QtCy+3+1LTpEGNyqFt/Igl5RbgprOIdjIxAj1MH5Adyj/6WmPQDc7yKUtN5CArIujjQydfWTo7FO1FUQQBGQXaaq08iXmGANgUk4JNHoDErJLkJBdgqMxWabHSsUidPC0R2cfR3TxcTIFPkcbWUu+PGqH2uSp2JpCwJYtWzBr1qw6H19bk2ZpaSni4uIQHBwMGxubpqgu1QPfd+sjCAJS8ir1f8uoOI2aVaSp8XEuSpkptIV62Btb4jzs4etsCzFPeZEFMRgEpBWUIiG7BLeyi3EtrQCXk/NwOTkfucXVn+INcFWii4+jKeh18XGCp4OCX26pVu3iVCwRWQ6d3oDotAKcvZWLcwm5uJZWgJj0QhRpah5t6utsi1BPe4R62Jm1wrnZ16/bBFFrE4tFUDkZ+9n1D66Ybqv8S83l5HxT0ItKzkdSrjEA3souxq+XUk3l3e3l6FzWqtelLOwFuir5RYYapU0GO6qqrm97M2fOxNatW1umMrcJCgrCwoULsXDhwlZ5fmp66QWlOHcrF2du5eLsrRxcSMyrdsoQqViEQDelqQXOGOAcEOJhxykmyGqJRCL4ONvCx9kW93T2Mm3PKdLgSkq+WeCLyShEZqEGB69l4OC1DFNZO7kEEaqKoFfeb4/dDqgu/M9qJVJSUkzXt2/fjpdffhnR0dGmbQ3te6jRaCCXy5usftR2aXQGRKXk48zNHJxNMAa5xJySKuUcFFL0DHBGL39ndPZxKuv/pmzRiXSJLJmLnRyDO7hjcAd307YSjR5XU41hL6os9F1NyUeRRo9TN3Nw6mbFgA2ZRISOng4VLXu+TohQOcKeX5KoEv42WAlvb2/TdScnJ4hEItO2rKwszJs3D4cOHUJ2djZCQ0Pxwgsv4F//+pfpMSNGjEDXrl0hl8vx+eefo0uXLoiMjMRPP/2ExYsXIzExEQMHDsSsWbMwa9Ys5OTkwNnZGQBw9OhRLFu2DCdPnoS7uzvuv/9+rF69GnZ2dhgxYgRu3ryJ559/Hs8//zwAnkq3ZOWnkM7cysHZsta4S8n5VVZcEImAME8H9ApwRu8AF/QKcEaohz1PHRE1kK1cgl4BLugV4GLaptMbEJtZZGzVS6po4csv1SEqxRgAd5yu2EeQmxJdfJzQL8gFQzq4o4OnPfvstWMMdg1RVMtSNxIJULnTf21lxWKgcgtaTWXt7BpWvxqUlpaiT58+WLp0KRwdHbF79248+uijCAkJwYABA0zlPvvsMzz11FM4cuQIBEFAfHw8pk6diueeew6PP/44zp49iyVLlpjt++LFixg7dixee+01bNq0CRkZGViwYAEWLFiALVu2YOfOnejRoweefPJJPPHEE03yeqjplGj0uJiUh7PlQS4hB2n56irlXJQy9ApwQe8AZ/QKcEF3Pyc4cHQfUbOQSsQI83JAmJcD7u9l3CYIApJyS8pCXj6iyk7lpuSVIj6rGPFZxdh90XjmxsNBgcGhbhgS6o7BHdzg58Lpu9qTNjkq9k41elRsbd+Axo8Hdu+uuG1nBxQXV192+HDgwIGK2x4eQGZm1XKNPDRbt27FwoULa10SbcKECYiIiMA777wDwNhil5eXh7Nnz5rKLFu2DLt378bFixdN21588UW88cYbpha7GTNmwNbWFh9//LGpzOHDhzF8+HAUFRXBxsamXn3sOCq2+QmCgJtZxTibUN4al4srKfnQ3bbIu0QsQoTKwdQS18vfBYFuXIGFyBJlF2lwOTkPFxLzcCwmCyfjs6G+rYU9wFWJIR3cMDjUHYNC3eDOAUptjtWPiqWG0ev1+O9//4vt27cjKSnJtBKH3W0tgn379jW7HR0djX79+plt69+/v9nt06dP48aNG9i2bZtpmyAIMBgMiIuLQ0RERBO/GqqvQrUO58v6xBlb43KRXc00Ix4OCvQ2nVJ1QTdfJ9jKObEqUVvgaifHsI4eGNbRA0+P7AC1To8zN3NxNCYTR2OycC4h1zgS90Qxvj6RAAAI93bA4FB3DOnghv7Brmx9tzIMdg1RWFjzfbevllC29mm1xLd1Jo+Pb3SV6uPdd9/Fe++9h/fffx/dunWDnZ0dFi5cCI3G/EP+9qAnCEKVVprbG3gNBgPmzp2LZ599tsrzBgQENNEroPoo1uhwLCYLkdcycCIuG9FpBVUafeUSMbr4OqKXvwt6BxpPq/o42bA1jshKKKQSDAp1w6BQNyyG8QveibgsHLmRhaMxWbiSkm9aN3fzkThIxCJ093MynrYNdUPvQBeumNHGMdg1REP6vDVX2UY4dOgQ7rvvPjzyyCMAjGHs+vXrdbamhYeHY8+ePWbbTp06ZXa7d+/euHz5Mjp06FDjfuRyOfT6muczo8YRBAHX0wsRGZ2BA9fScTIuBxq9+SkYX2dbswEOnX0cuaA5UTtir5BiVLgXRoUbp13JKlTjWKwx5B29kYn4rGJT14wP99+AQipG3yAXDC4Let18nSDlyPY2hcGuHejQoQO+//57HD16FC4uLli7di1SU1PrDHZz587F2rVrsXTpUsyZMwfnzp0zzYVX3sKzdOlSDBw4EE8//TSeeOIJ2NnZ4cqVK9i3bx/WrVsHwDiP3cGDB/Hggw9CoVDA3d29pqekOuSXanHkeiYir2Ug8loGUvJKze73c7HFiE4eGNrBHb0DXODpyP6KRFTBzV6Bid19MLG7DwAgKbcER28YT9seuZGJ9AI1jtwwtvABxmmMBoS4GQdjdHBHmBdH3Fo6Brt24KWXXkJcXBzGjh0LpVKJJ598ElOmTEFeXl6tjwsODsZ3332HxYsX44MPPsCgQYOwYsUKPPXUU1AojJ1vu3fvjsjISKxYsQLDhg2DIAgIDQ3F9OnTTft59dVXMXfuXISGhkKtVnO6kwYwGAREpeQbg1x0Bk7fyoG+0mAHhVSMgSFuGB7mgeGdPBDibsd/ukRUb77OtvhnX3/8s68/BEFATEahKeQdi8lCfqkOf1xJwx9X0gAYV8kYFOqOIWVBz9+VI24tDUfFcq3YBnnjjTewceNGJCQkNOl++b5XyC7S4NB1Y5A7eD0DmYXmfSFDPOwwPMwDIzp5YkCwK/vDEFGz0BsERCXn40hMJo7cyMTJ+GyUas27e/i52GJIqDuGdnTHyHBPTpbcTDgqlprM+vXr0a9fP7i5ueHIkSN4++23sWDBgtaullXRGwScT8zFgWjj6dULiblmgx6UcgkGh7pjRCcPDA/z4DdkImoRErEI3fyc0M3PCfOGh0Kt0+PcrVwcKeufdy4hF4k5Jdh+KgHbTyVALhVjZCcPTOzug7sjPKGUM2K0Br7rVKvr16/j9ddfR3Z2NgICArB48WIsX768tavV5qXnl5r6yR26nom8Eq3Z/eHeDhheFuT6BrpyfUgianUKqQQDQtwwIMQNi+4JQ5FahxPx2Th6IxN/XklHbGYRfr+cht8vp8FWJsGoCE9M7KbCyHBPnlloQTwVy1OxFsHa33et3oDTN3MQeS0DB6IzcCUl3+x+RxsphnU0Brm7wjzg7WR97wERWS9BEHAlpQC/XEjGLxdScCu7YoJ+O7kEozt7YUI3Fe4K82DIawSeiiWyAIk5xaZBD0djslCo1pnd393PqayvnAd6+DlzSgEiarNEIhE6+ziis48j/j22Ey4l5ZtCXlJuCX48l4wfzyXDQSHFPZ29MLGHCkM7ePBsRDNgsCNqQjfSC/DjuWT8eikVN9LNJ7R2s5PjrjBjq9zQju5c1oeIrJJIVNE3b9m4cJxLyMUvF1Kw+0IKUvNLsfNsEnaeTYKjjRRju3hjYg8fDA51g4xfbpsEg10NDAZD3YWoybTl9zslrwQ/nzd+G72cXHGKVSwCege4mKYi6erjBLGYU5EQUfshEonQq2y5whXjI3DmVo4x5F1MQUaBGjtOJ2LH6US4KGW4t6s3Jnb3wYBgV57BuAPsY3fbueryVRkkEgk8PDwgl8s5L1gzEgQBGo0GGRkZ0Ov16NixI8S3L7lmgXKLNdhzMRU/nkvCifhs0yhWqViE4WEemNzTByPCPOGk5BqMRES30xsEnIjLxu6Lyfj1YiqyKq1j7W4vN4W8fkGukPALcYP62DHYVfMGaTQapKSkoLi4uJpHU3NQKpVQqVSQy+WtXZUalWj0+ONKGn48l4zIa+nQ6iv+dPoHueK+Xj4Y31UFFzvLfQ1ERJZGpzfg79iykHcpFbnFFbMEeDooML6bChO7q9A7wKXdnvVgsKtDfd4gQRCg0+m4xmkLkEgkkEqlFtkyqtUbcORGJn48l4y9l1NRpKn4fYhQOeK+nj6Y1MMHvs62rVhLIiLroNUbcDQmC7+cT8bvl1ORX1ox6EzlZGMKeT39nS3yM6O5MNjVoSFvELU/giDgzK0c/HguGbsvpJidIvBzscV9PX1wX09fhHk5tGItiYism0ZnwOEbGfjlfAr2RqWZzSzg62yLid1VmNjdB119Ha0+5DHY1YHBjqpzLa0AP5xNwk/nk5GYU2La7mYnx8TuKkzu6YveAe3rWyIRkSUo1epx8FoGfrmQgj+upKG40tmTQDcl7uvpi4f6B1jtHKAMdnVgsKNyiTnF+Pl8Cn48l4SrqQWm7XZyCcZ28cZ9vXwxJNSNI7SIiCxEiUaPA9Hp+OVCCv68mmZav1YqFmFsV2/MHBSEfkEuVvUlnMGuDgx27Vt2kQZ7LhrD3Mn4HNN2mUSEEZ08cV9PH9wd7gVbOWdHJyKyZEVqHf64koZtf9/Cifhs0/YIlSNmDgrEfT19reJ/OYNdHRjs2p9ijQ77oowjWg9ey4DOYPy1F4mAgcFuuK+nD8Z1VXF6EiKiNioqOR+fH4vHD+eSTK14TrYyTO/nj0cGBCLATdnKNWw8Brs6MNi1D1q9AYeuZ5SNaE1DibaiT0ZXX0fc18MXE3uooHLiiFYiImuRW6zBjlOJ+PzveCRkG/tLi0TA3eGemDEoCEM7uLe5aVMY7OrAYGfdrqUV4PNj8dh9IQU5leZDCnRT4r4ePpjc0xcdPO1bsYZERNTc9AYBB6LTsfVoPA5dzzRtD3G3w4xBgfhHHz842LSNszQMdnVgsLM+BoOAA9fSsflwPA7fqPgDdrdXYFIPFe7r6Ysefk5W1ZmWiIjqJyajEF8cu4nvTieapk2xk0vwjz5+mDEoEB08LXv6Kga7OjDYWY9CtQ7fn07E1qPxiMssAmBco3VsF288NCAAg0I4opWIiIwK1TrsPJOIz47GIyajyLR9SAc3zBwUhLsjvCxyCTMGuzow2LV9CdnF+OxoPLafTEBB2bcvBxsp/tU/AI8ODIS/a9vtJEtERM1LEAQcuZGFz47F488raSgbTwdfZ1s8OigQ0/v6W9TykAx2dWCwa5sEQcDxuGxsORKHfVEVf4ghHnaYPTgID/T2g51C2rqVJCKiNiUhuxhfHr+J7ScTTOvUKqRiTO7hg5mDg9DV16mVa9gOgt3Bgwfx9ttv4/Tp00hJScGuXbswZcqUej+ewa5tKdXq8fP5ZGw+Eo8rKfmm7XeFeWD2kCAM7+jR5kY4ERGRZSnV6vHTuWRsPRqPqEqfNX0CXTBzcBDu7eINubR1uvY0JLe0yeaNoqIi9OjRA7Nnz8Y//vGP1q4ONZP0glJ8+fctbPv7pmm9VhuZGP/o7YfZQ4IsvrMrERG1HTYyCab188c/+/rh9M0cfHbsJn69mILTN3Nw+mYOPB0UeGhAAB7qHwBPR8tduqxNtthVJhKJ2GJnZS4m5mHLkTj8fCEZWr3x11PlZIMZg4Lwr/7+cFZaTr8HIiKyXmn5pfjq+C18deIWMgrUAIxLl43vpsLMwYHoHdAyS5dZfYsdWR+d3oB9UWnYfCTObJmvPoEumD0kCGO7eEPG0a1ERNSCvBxt8Pw9YXh6ZAf8eikFnx+7idM3c/DT+WT8dD4ZvT0UmN3JAfeM7QsbmWUsXdYugp1arYZarTbdzs/Pr6U0taS8Yi2+OXkLnx+7iaRc4wzhUrEIE7urMHtIMHr4O7duBYmIqH3JzwciI4H0dCAjA8jIgDwjA/elp+O+jAwk3zcN70Xcix/PJyP/ynWot36P+66/jN8WDrOIuVLbRbBbvXo1Vq1a1drVoEpupBdi69E4fH86ybTUl6udHA8PCMAjAwPhZcH9F4iIyMIJAqDXA9KymJOVBfzwgzGoVQpspuvz5wPLlxvLJiUBkyfXuGufwYPx9ss9sHx8BH784zyUZ3/AuG7eFhHqgHYS7JYvX45FixaZbufn58Pf378Va9Q+CYKAg9czsflwHCKvZZi2h3s74LEhwZjc08dimrKJiMiCCAKQlweIxUB5H7O0NGDz5opwVjmwZWQYg9orrxjLZmQAjz9e8/6Tkique3kB/foBHh7VX8LDARgbI2bf1w+6ib9ipN5yhiu0i2CnUCigUChauxrtVrFGh51nkrDlSJxppm/jgsxeeGxoEAaFuFnMNx0iImoB5UEtPd0Y1Ly9jdsTE4F3360+rGm1wMqVFWEtJwd44YWanyOjogEB3t7AhAkV4czT0/x6QEBFWVdX4MSJer8UqUQMqQW1SbTJYFdYWIgbN26YbsfFxeHcuXNwdXVFQOWDQ60qKbcEnx+Lx9fHbyG/1Lg6hL1Ciml9/TFzcCAC3exauYZERNQkDAYgN7cihPn4ACEhxvvi4oAVK8xDWmamMagBwMsvA+XdpYqKgPffr/l5cioG10GlAmbPNm9Nuz2wlXN2Bn75pQlfsOVqk8Hu1KlTGDlypOl2+WnWmTNnYuvWra1UKyoXl1mE9/Zdw+6LKdCXLQ8R6KbErMFBmNrHDw42slauIRER1ao8qFUOY2FhQNeuxvuvXweeeqqiVS0zE9DpKh7/0kvAq68ar2u1wNdfV/88Dg7G5yrn4wMsXVpzWLO1rSjr5GQ8FUtm2mSwGzFiBNr49HtWKTWvFB/8eR3fnkowBbrBoW54bEgwRoZ7WuTCykRE7YLBYGztuv30Zu/eQP/+xjJRUcD06RUtanq9+T5efLEi2AHAn39WfR4HB2MQs7ev2Obrazy9entI8/AAbGyqPv6//22a19xOtclgR5Ylp0iDjZEx2Ho0Hmqd8ZvXqHBPLLonzCLW2CMisjoGA5CdXXXAQEYGcNddwIgRxnLnzwNjxhhHhd4e1ADjKdLyYCeXA5cumd/v6FgRxnx8Krb7+QFffFE1qFXXn93ODqg0gJGaF4MdNVqRWofNh+PwycFYFKiNTfD9glzwn3vD0S/ItZVrR0TUhhgMxvBV3ejOjAxg/HjjBQBOnQIGDDA/hVnZCy9UBDt7e+P+yjk7VzvCEwDg7w/8/ntFy5q7e/VBDTCeEn3kkTt91dQMGOyowdQ6Pb4+fgsf7r+BzELjGq7h3g5Yem84RnTy4AhXIiIAKC0FYmJqDmvTpgFTpxrLnjwJDBxY874cHSuCnYtLRahzdq7aatavX8XjAgKAc+eM293dja1yNVEojK171KYx2FG96Q0CfjibhPf+uIbEHOMqEYFuSiy6JwyTuvtAzD50RGTt8vOB06erTnBbfnnsMWDGDGPZCxeMLWs1CQ2tCHYeHsafLi5VBw14ela0wAFAYCCQnGwMarI6BqPJZECPHo1+udT2MNhRnQRBwL6oNLyzNxrX0goBAJ4OCjx7d0dM7+fPNVyJqG3LygL++KP6FQkyMoDnngPmzTOWjY4GRo2qeV+DBlVc9/Q0zol2e4taeVirHPqCggCNpu6gBhhXU1CpGvVSyfox2FGtjsVkYc3vV3H2Vi4AwMlWhqdGhGLmoCDYyi1oRkYiat90OuOlfJRlSgrw3Xc1h7UXXgAWLzaWjYsDHnyw5n3Hx1dc9/YGIiKqD2oeHkD37hVlg4KMobE+xGLjhegOMdhRtS4m5mHN71dx6HomAMBWJsFjQ4Pw5F2hcLLlPHRE1My0WmMAs7U1np4EgFu3gE8/rT6s5eQAa9YA//63sWxSEvDsszXvPy2t4rqPj3EkaU0tax07VpT19zdOC0JkoRjsyExMRiHW7jVOLgwAMokI/+ofgAWjOsDTwaaORxMR1UCjMc6Nlp5uDEvlU2fcuGEMZLcHtdxc4/1vvQX85z/G6xkZwOuv1/wclZeQ8vMD/vnPmlvWfH0ryvr4AJGRTfpyiVoLgx0BAFLySvDBH9ex43Qi9AYBIhEwpacvnh8dhgA3ZWtXj4gsjUZjHsZCQ40XwNii9cIL5i1reXkVj/3vf42rCwDGwQifflr9c4jFxiWmygUEAE8/Xf2KBB4exv5s5by9gW+/bdrXTNQGMNi1c9lFGmw4cAOfHbsJTdnkwqMjPLFkbCeEezu2cu2IqMWUlladjqNnz4qVBs6dq1hCKiPDGMgqW70aWLbMeF2tBn78sepziMXGkZyV+5IFBRnXCa0pqFUu6+EBfPhhE75oIuvDYNdOFap12HQoDp8eikVh2eTC/YNdsfTeTugTyMmFidq88qB2+3Qcw4ZVzHN24gTwr38ZtxcUVN3Hm2+aLyH199/m90skxqDm6Wlct7NcSAiwYUPVU6AuLlUHCLi6GheBJ6ImwWDXzqh1emz7+xY+2n8DWUXGyYW7+Dji32M7YXgYJxcmslglJeYBrXJgGz/e2PkfAA4fBsaNAwoLq9/Pm29WBDu5HIiNrbhPKq0Iah4exn5q5UJDge+/N29Vc3aufiSnk1PF9CBE1KIY7NoJnd6AnWeT8MEf15GUa5xcONjdDovuCcOEbipOLkzU0rRa4ySzNa1KMH06MHassez+/bXPneboWBHsHBwqQp1UWvX0ZuUlpMLCgEOHKso4OwM1fblzcAAeeOCOXzYRNS8GOysnCAJ+v5yKd/Zew4104z97L0cFFo4Ow9Q+fpxcmKgpFRYCV67UHNYefxyYMsVY9tAh4O67a95XaGhFsHNzM/6UyaqfkqPyElLh4cZJdMtPj9bWCq9UAkOH3tFLJiLLwmBnxY7cyMSa36NxPiEXAOCslGH+iFDMGBQEGxknFyaql6wsY9+ymsLa888b+6kBxvU+a2tZGziwIth5eBhPhdY0wnPkyIrHde5snP7D0bH2oAYY1/sMC7uTV0xEbRiDnRVKyi3Byh8v4Y8r6QAApVyCOUOD8cRdIXC04eTC1A4JgnHRdEnZF5qUFOC336oPa+npxs78c+YYy164AEycWPO+b9youO7lZeyXVl1Y8/QE+vevKNu1q3GAQ336tUql5oMTiIhqwGBnRfQGAVuPxuPdvdEo1ughk4jw8IBAPD2yAzwcFK1dPaKmIwjG055iMWBnZ9wWHw98803Ni7OvXWucrgMwnqp87LGa95+UVHHd1xfo06f6SW49PMxHjXbuDCQk1O81cKASETUDBjsrcSkpD8t3XsTFJOMkoP2CXPDm/d3Q0cuhlWtGVA+CYJxuIyPD2IG/vE9ZdDTwySfVnwItLQU++giYP99Y9uZNYPnymp+j8qoEAQHGkaM1rUoQElJRNiwMOHWqyV8yEVFzYLBr44rUOry37xo2H4mDQQAcbaRYPj4C0/v6c6QrtR5BME5gWx7CAgIqlnC6cMF8Canyi1ptvP/DD42rCwDG9TzXrq35eSovsB4UBMyYUfX0Z/l1b++KsiEhwJ49TfqSiYgsAYNdG/bX1TS89MNl0/Qlk3r4/H979x4U1ZXnAfzb0LwVIoIIoihGBFFRYEAgxiSlGOPIuDMpmdIymjJToVwnGMtkcckEcbNhDRXL4Iomjo+pWjSOKClTQWPXJCGIKR0JZpm0WY3ggxJE8NUGeZ/9407T3XTz6IZuum9/P1Vdudw+t/t361ed+/Pce87Bn34dxTVdafgJIS0JpX+bc9Ys3RJSFy8aLiHV1CQtOaW1axewYYO0/fAhUFRk+nu8vQ2PmzpVWtS9rwEG2tuwABAWBvzlL8N73kREDoaFnQNqfNSK3M/V+KK6HgAQOsYL/7F8Jp6fPm6EIyOHIYQ0yrL3oIGUFCA6WmpTUSH1nGnf7+gw/IyCAuCPf5S2W1sBlcr4e3x8pALMTW/QTkQEkJ9vulDz7rUu8YQJUu8eERENCgs7B9LdLXD4wk1sP/0TNK2dcHVR4LVnpiBz4TR4uzOVTq27W1eo9X4ebckS6eF/APjqK2DVKqlHrbPT+HMKCnSFHQD88IPh+6NG6Yox/QXXo6Kk3rLePWteXsbfERQEbN485FMmIiJjrAYcxP81aPDvJdWovHEfABAT6of3fzsL0SGcAkGWtIWafpGmv52eLvWuAcCpU8CyZUBXl+nP8vXVFXaenkBDg+690aMNe80mTtS9N3Om9NnaQi0gwHShBkiDHV55ZcinTUREQ8PCzs61dnRh11dX8XFZDTq7BXzcXfHW4ulYnTQZrhwc4ViePJGm5DA1wvPuXeAPf9CtRPDFF0BaWt+fNXWqrrB76ildUefra9xrFhWlOy4mBqis1BVqnv08j+nnB7z44lDOmIiIbIyFnR2r+LkJ2SXVuN7cAgBInRGE3N9EI9ivj14Tsr3796VRnqYmub17V7rluHSp1PZvf5N61voyb56usAsMlP7r52d6Sg79JaRiY6W50wIDpVUH+uPjI7UnIiJZYmFnh5oft+E/v7iME1XSJKnjfT2xNS0aL84cP8CRNCwaGqR1PHvf/tRub9sG/O53UtuKiv6LtbQ0XWE3bpzUu2ZqRYLAQOC553THxcdLAxIGKtQAqU1oqKVnS0REMsLCzo4IIVBcWYf3Sy/jfksHFArglXlh2Lx4OkZzKTDzdHZKIz+1ozGvX5dub/ZVrO3cKQ0qAKRblStW9P3ZN27otkNCpFGepkZ4BgYaLiGVkCD18A2GUim9iIiIzMArh52oufsY2SX/wHc10oSrkeNHI++3szB30pgRjsxOdHZKIzlHjZJeAPDTT8CRI6aLtXv3gEOHdA/0q9W6edRMuXNHtz1pEvDMM32v96n/zFpsrLQ6AhERkR1gYTfC2ju7sbfsGv7765/R3tkNTzcXbFwYgXXPTIGbq8tIh2c92kKtsVHq9QoIkPb/8AOwd6/pQg2QirU1a6Tt2lrptmhf9JeQCg+Xbp/2td7npEm6trNmSbdiiYiIHAwLuxH09+v3sOVENX5ufAwAeDYiEO/9ZiYmjfUe4Eg71NEhFWraYiwqSirYAOD8eWD7dsOeNf1bkgcPAmvXStv19VJhZ4pCIS1TpTVtmrSoe1/Fmna9UQCIjASKi4f1lImIiOwNC7sR8LClA/91+jKOXLgFAAgY5Y4//XoG0mJCoFDYyRQm7e2GhVpjozRqU7s4elkZkJ2tK9QePDA8Xr9Yu38fKCkx/g6FQiq+9Odfi4oCcnJMjwT19zd87uzpp4HCwuE8ayIiIofm0IVdYWEh8vPzUV9fj+joaOzcuRPz588f6bD6JITA5/9bj22fq9H0WFrw/Pe/moisJZF4ytvdul/e3m56Oo7Fi3XPjJ05Iz2H1tgorefZ24EDusKutVUaEarPxUW6pdp72o3Zs4Hdu4171Pz9AVdXw88ICwO2bh220yYiInImDlvYHT16FBs3bkRhYSFSUlLw8ccfY8mSJVCr1Zik/7yUnbh1rwXvfPYPlF2RnvuaGuiD9/9lFhLDxw5wZB/a2ownuNUWaytWAHPmSO1OngRWrza8halv/37DwQBXr+q29Qu13rc2Y2OBY8cMizV/f+mY3kJCgPXrLTtPIiIiGjSFEEKMdBCWSExMRGxsLPbs2dOzLyoqCsuXL0deXl6/xz569Ah+fn54+PAhfH19rRpnZ1c3DlTUYofqClo7uuHu6oJ/ff5pZDwXDg+lXm9Ve7s0MtPU8lHaVQkSE6W2f/2rtKRUX/78Z2DdOmlbpQJSU6VtV1epUNPvNVu3Dli0SHr//n2gulpXrI0ZY7pQIyIiIpsxp25xyB679vZ2VFZWIisry2B/amoqzp07N0JRGatW38S+A2fw8OZtLGl5iLleHVgWrMSYI8VAwT9XJXj2Wanx8ePAypV9f9i8ebrCTrv4ulKp61HTL9YiI3XHJSUBly/rJsftr1AbM0YXDxERETkchyzsmpqa0NXVhaCgIIP9QUFBaNBf4Pyf2tra0NbW1vP3o75uSw6jQxW1+P6DPSg4md93o2XLdIXUuHFSodbX3GnaRdwBYP58oLl54EINkOZ80y/0iIiISLYcsrDT6j2CVAhhclRpXl4ecnNzbRUWAOBXU/zx5eixeDAmEKMmhkAZNM54lKd+79jzz0u3YwczKtbDY3BLTREREZFTccjCLiAgAK6urka9c42NjUa9eACwZcsWbNq0qefvR48eYeLEiVaNMTrED3k7N+Cp//m3wR3AZ9mIiIhoiByymnB3d0dcXBxUKpXBfpVKheTkZKP2Hh4e8PX1NXjZwuQAH5t8DxERERHgoD12ALBp0yasXr0a8fHxSEpKwieffIKbN28iIyNjwGO1A4Ft8awdERER0VBo65XBTGTisIVdeno6mpubsW3bNtTX12PmzJkoLS1FWFjYgMdqNBoAsPrtWCIiIqLhotFo4Ofn128bh53Hbii6u7tx+/ZtjB492qpLeGmf5bt165bNbv9S/5gT+8Oc2Cfmxf4wJ/bJFnkRQkCj0SAkJAQuAzyT77A9dkPh4uKC0NBQm32fLZ/ro8FhTuwPc2KfmBf7w5zYJ2vnZaCeOi2HHDxBRERERMZY2BERERHJBAs7K/Lw8EBOTg48OJmw3WBO7A9zYp+YF/vDnNgne8uLUw6eICIiIpIj9tgRERERyQQLOyIiIiKZYGFHREREJBMs7IiIiIhkgoXdEBUWFmLKlCnw9PREXFwcysvL+21fVlaGuLg4eHp6Ijw8HHv37rVRpM7DnJycOHECixYtQmBgIHx9fZGUlIQvv/zShtE6B3N/J1oVFRVQKpWYM2eOdQN0Uubmpa2tDdnZ2QgLC4OHhwemTp2KAwcO2Cha52BuToqKihATEwNvb28EBwfj1VdfRXNzs42ilb9vv/0Wy5YtQ0hICBQKBT777LMBjxnx67wgi3366afCzc1N7Nu3T6jVapGZmSl8fHzEjRs3TLavqakR3t7eIjMzU6jVarFv3z7h5uYmiouLbRy5fJmbk8zMTLF9+3Zx4cIFceXKFbFlyxbh5uYmvv/+extHLl/m5kTrwYMHIjw8XKSmpoqYmBjbBOtELMlLWlqaSExMFCqVStTW1orz58+LiooKG0Ytb+bmpLy8XLi4uIiPPvpI1NTUiPLychEdHS2WL19u48jlq7S0VGRnZ4vjx48LAKKkpKTf9vZwnWdhNwQJCQkiIyPDYF9kZKTIysoy2f7tt98WkZGRBvtef/11MW/ePKvF6GzMzYkpM2bMELm5ucMdmtOyNCfp6eninXfeETk5OSzsrMDcvJw6dUr4+fmJ5uZmW4TnlMzNSX5+vggPDzfYV1BQIEJDQ60WozMbTGFnD9d53oq1UHt7OyorK5GammqwPzU1FefOnTN5zHfffWfUfvHixbh48SI6OjqsFquzsCQnvXV3d0Oj0cDf398aITodS3Ny8OBBXLt2DTk5OdYO0SlZkpeTJ08iPj4eH3zwASZMmICIiAhs3rwZT548sUXIsmdJTpKTk1FXV4fS0lIIIXDnzh0UFxdj6dKltgiZTLCH67zSJt8iQ01NTejq6kJQUJDB/qCgIDQ0NJg8pqGhwWT7zs5ONDU1ITg42GrxOgNLctLbhx9+iF9++QUrVqywRohOx5KcXL16FVlZWSgvL4dSyf9FWYMleampqcHZs2fh6emJkpISNDU1Yf369bh37x6fsxsGluQkOTkZRUVFSE9PR2trKzo7O5GWloZdu3bZImQywR6u8+yxGyKFQmHwtxDCaN9A7U3tJ8uZmxOtI0eOYOvWrTh69CjGjRtnrfCc0mBz0tXVhZUrVyI3NxcRERG2Cs9pmfNb6e7uhkKhQFFRERISEvDSSy9hx44dOHToEHvthpE5OVGr1XjjjTfw7rvvorKyEqdPn0ZtbS0yMjJsESr1YaSv8/znsIUCAgLg6upq9C+pxsZGo2pda/z48SbbK5VKjB071mqxOgtLcqJ19OhRrFu3DseOHcPChQutGaZTMTcnGo0GFy9eRFVVFTZs2ABAKiiEEFAqlThz5gxeeOEFm8QuZ5b8VoKDgzFhwgT4+fn17IuKioIQAnV1dZg2bZpVY5Y7S3KSl5eHlJQUvPXWWwCA2bNnw8fHB/Pnz8d7773Hu0AjwB6u8+yxs5C7uzvi4uKgUqkM9qtUKiQnJ5s8Jikpyaj9mTNnEB8fDzc3N6vF6iwsyQkg9dStXbsWhw8f5rMpw8zcnPj6+qK6uhqXLl3qeWVkZGD69Om4dOkSEhMTbRW6rFnyW0lJScHt27fx+PHjnn1XrlyBi4sLQkNDrRqvM7AkJy0tLXBxMbyMu7q6AtD1EpFt2cV13mbDNGRIOzR9//79Qq1Wi40bNwofHx9x/fp1IYQQWVlZYvXq1T3ttcOg33zzTaFWq8X+/fs53ckwMzcnhw8fFkqlUuzevVvU19f3vB48eDBSpyA75uakN46KtQ5z86LRaERoaKh4+eWXxY8//ijKysrEtGnTxGuvvTZSpyA75ubk4MGDQqlUisLCQnHt2jVx9uxZER8fLxISEkbqFGRHo9GIqqoqUVVVJQCIHTt2iKqqqp4paOzxOs/Cboh2794twsLChLu7u4iNjRVlZWU9761Zs0YsWLDAoP0333wj5s6dK9zd3cXkyZPFnj17bByx/JmTkwULFggARq81a9bYPnAZM/d3oo+FnfWYm5fLly+LhQsXCi8vLxEaGio2bdokWlpabBy1vJmbk4KCAjFjxgzh5eUlgoODxapVq0RdXZ2No5avr7/+ut9rhD1e5xVCsL+WiIiISA74jB0RERGRTLCwIyIiIpIJFnZEREREMsHCjoiIiEgmWNgRERERyQQLOyIiIiKZYGFHREREJBMs7IiIiIhkgoUdERERkUywsCMiIiKSCRZ2RERERDLBwo6IiIhIJv4fSxbMx9pOnCoAAAAASUVORK5CYII=", "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 }