tutorial/SMT_ExpandedLHS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/raul-rufato/smt/blob/newsmt/tutorial/SMT_MixedInteger_application.ipynb\" target=\"_blank\" rel=\"nofollow\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Augmented LHS Documentation for the Toolbox SMT\n",
    "\n",
    "<div class=\"jumbotron text-left\"><b>\n",
    "The French Aerospace Lab ONERA <br/>\n",
    "Information Processing and Systems Department (DTIS) <br/>\n",
    "Multidisciplinary Methods, Integrated Concepts (M2CI) Research Unit <br/>\n",
    "</div>\n",
    "\n",
    "Rubén CONDE ARENZANA ONERA/DTIS/M2CI\n",
    "**Latest update:** April 2023 - version 2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we present some examples in order to use the augmented LHS functionality within the Toolbox SMT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info fade in\" id=\"d110\">\n",
    "<p>Latest updates</p>\n",
    "<ol> - add possibility to extand a LHS for mixed integer or categorical variables  </ol>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p class=\"alert alert-success\" style=\"padding:1em\">\n",
    "To use SMT models, please follow this link: https://github.com/SMTorg/SMT/blob/master/README.md. The documentation is available here: http://smt.readthedocs.io/en/latest/\n",
    "</p>\n",
    "\n",
    "\n",
    "**Reference paper:** https://www.semanticscholar.org/paper/Adaptive-Response-Surface-Method-Using-Inherited-Wang/cb3505730d406b2baae25b34fa3fa74058461592"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cite us:\n",
    "\n",
    "M.-A. Bouhlel, J. T. Hwang, N. Bartoli, R. Lafage, J. Morlier, J .R.R.A Martins (2019), A Python surrogate modeling framework with derivatives, Advances in Engineering Software, 102662\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "@article{SMT2019,\n",
    "  title={A Python surrogate modeling framework with derivatives},\n",
    "  author={Mohamed Amine Bouhlel and John T. Hwang and Nathalie Bartoli and R{\\'e}mi Lafage and Joseph Morlier and Joaquim R. R. A. Martins},\n",
    "  journal={Advances in Engineering Software},\n",
    "  pages={102662},\n",
    "  year={2019},\n",
    "  publisher={Elsevier}\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Statement "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generation of optimized design of experiments (DoE) is an essential tool when dealing with the analysis and simulation processes that are computationally expensive. In this regard, Latin Hypercube Samplings (LHS) appears to be a sampling method that can cope with many input variables ensuring that each of them has all portions of its range represented. Additionally, this method is computionally cheap to generate. \n",
    "\n",
    "In this document, an expanded LHS is presented. To understand this functionality, we can start with a 2D LHS that counts on 10 points. For some purposes, for instance adaptive design or nested designs, it is necessary to add new points to the DoE but keeping the ensemble LHS. For the moment, we can just add a multiple of the initial number of points to ensure an LHS (in this example case we can add 10,20,30... points).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial LHS Sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we create a LHS sampling in 4D that will be the initial LHS to be augmented"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x14b40f3ffa0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAGwCAYAAAC5LOMAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA62ElEQVR4nO3de3RU9aH+/2fnQhIuGS6BZJBggsg9QSVQI3JRKgqtkcLpodojUKXroChqpCih3i1pqxXECxRF0VoVVxNszhFZUCFBfmIhGJpQERUTw9GJaaxMbMAAyf79YTNfhwxhT5iZPZO8X2vNOsyez955PnzOMQ/77L3HME3TFAAAAIA2RdkdAAAAAIgEFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABbE2B0g1Jqbm/X555+rR48eMgzD7jgAAMAC0zT19ddfq3///oqK4rwf7NHpivPnn3+u1NRUu2MAAIB2OHz4sAYMGGB3DHRSna449+jRQ9K3/4eXmJhocxoAAGBFfX29UlNTPb/HATt0uuLccnlGYmIixRkAgAjDZZawExcJAQAAABZQnAEAAAALKM4AAACABRRnAAAAwAKKMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACygOAMAAAAWUJwBAAAAC2wtzqtXr1ZmZqbn66+zs7P15ptvtrlPSUmJxowZo/j4eA0aNEhr1qwJUVoAAAB0ZrYW5wEDBujXv/61SktLVVpaqssvv1zXXHON/v73v/scX1lZqenTp2vChAkqKytTXl6eFi1apIKCghAnBwAAQGdjmKZp2h3iu3r37q1HHnlEN954Y6vP7rrrLhUVFenAgQOebQsWLNDf/vY37dq1y9Lx6+vr5XA45Ha7lZiYGLDcAAAgePj9jXAQNtc4NzU16dVXX1VDQ4Oys7N9jtm1a5emTp3qte3KK69UaWmpTpw44XOfxsZG1dfXe70AAAAAf8XYHaCiokLZ2dn65ptv1L17d23cuFEjRozwObampkbJycle25KTk3Xy5EnV1dXJ6XS22ic/P18PPPBAULL7MnLkSJ133nkh+3mBUF5erszMTLtj+IXMoUHm0IjEzFJk5iZzaBw6dOi0l10Ckcz24jx06FDt27dPR44cUUFBgebOnauSkpLTlmfDMLzet1xpcur2FkuXLlVubq7nfX19vVJTUwOUvrXzzjtPRUVFQTt+MOTk5JA5BMgcGmQOnUjMTebQyMnJsTsCEBS2F+cuXbpo8ODBkqSsrCzt2bNHjz/+uH7/+9+3GpuSkqKamhqvbbW1tYqJiVGfPn18Hj8uLk5xcXGBDw4AFrjcx1RZ16D0pG5yOhLsjgMAOAu2F+dTmaapxsZGn59lZ2frf/7nf7y2bdmyRVlZWYqNjQ1FPACwbMOeai0trFCzKUUZUv7MDM0eO9DuWACAdrL15sC8vDy9/fbbqqqqUkVFhZYtW6bi4mL99Kc/lfTtZRZz5szxjF+wYIE+/fRT5ebm6sCBA3ruuee0bt06LV682K4pAIBPLvcxT2mWpGZTyivcL5f7mL3BAADtZusZ5y+++ELXX3+9XC6XHA6HMjMztXnzZl1xxRWSJJfLperqas/49PR0bdq0SXfccYeeeuop9e/fX6tWrdKsWbPsmgIA+FRZ1+ApzS2aTFNVdUe5ZAMAIpStxXndunVtfr5+/fpW2yZNmqT33nsvSIkAIDDSk7opypBXeY42DKUldbUvFADgrITNc5wBoCNxOhKUPzND0f9+4k+0YWj5zFGcbQaACBZ2NwcCQEcxe+xATRzSV1V1R5WW1JXSDAARjuIMAEHkdCRQmAGgg+BSDQAAAMACijMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAAAsoDgDAAAAFlCcAQAAAAsozgAAAIAFFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABZQnAEAAAALKM4AAACABRRnAAAAwAKKMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACygOAMAAAAWUJwBAAAACyjOAAAAgAWGaZqm3SFCqb6+Xg6HQ263W4mJiQE/flpamjIzMwN+3GDavXu3xo0bZ3cMv5A5NMgcGpGYWYrM3GQOjfLyclVVVQX0mMH+/Q1YEWN3gI4mMzNTRUVFdsfwS05ODplDgMyhQebQicTcZA6NnJwcuyMAQcGlGgAAAIAFFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABZQnAEAAAALKM4AAACABRRnAAAAwAKKMzotl/uY3jlUJ5f7mN1RAABABIixOwBghw17qrW0sELNphRlSPkzMzR77EC7YwEAgDDGGWd0Oi73MU9plqRmU8or3M+ZZwAA0CaKMzqdyroGT2lu0WSaqqo7ak8gAAAQESjO6HTSk7opyvDeFm0YSkvqak8gAAAQESjO6HScjgTlz8xQtPFte442DC2fOUpOR4LNyQAAQDiztTjn5+dr7Nix6tGjh/r166cZM2bo4MGDbe5TXFwswzBavT744IMQpUZHMHvsQO28+zK98vOLtfPuy7gxEAAAnJGtT9UoKSnRwoULNXbsWJ08eVLLli3T1KlT9f7776tbt25t7nvw4EElJiZ63vft2zfYcdHBOB0JnGUGAACW2VqcN2/e7PX++eefV79+/bR3715NnDixzX379eunnj17BjEdAAAA8P+E1TXObrdbktS7d+8zjr3wwgvldDo1ZcoUbd++/bTjGhsbVV9f7/UCAAAA/BU2xdk0TeXm5urSSy/VqFGjTjvO6XRq7dq1KigoUGFhoYYOHaopU6Zox44dPsfn5+fL4XB4XqmpqcGaAgAAADqwsPnmwFtuuUXl5eXauXNnm+OGDh2qoUOHet5nZ2fr8OHDevTRR31e3rF06VLl5uZ63tfX11OeAQAA4LewOON86623qqioSNu3b9eAAQP83v/iiy/WRx995POzuLg4JSYmer0AAAAAf9l6xtk0Td16663auHGjiouLlZ6e3q7jlJWVyel0BjgdAAAA8P/YWpwXLlyol19+WX/+85/Vo0cP1dTUSJIcDocSEr59TNjSpUv12Wef6cUXX5QkrVy5UmlpaRo5cqSOHz+ul156SQUFBSooKLBtHgAAAOj4bC3Oq1evliRNnjzZa/vzzz+vefPmSZJcLpeqq6s9nx0/flyLFy/WZ599poSEBI0cOVJvvPGGpk+fHqrYAAAA6IRsv1TjTNavX+/1fsmSJVqyZEmQEgEAAAC+hcXNgQAAAEC4ozgDAAAAFlCcAQAAAAsozgAAAIAFhmnlDr0OpL6+Xg6HQ263OyhfhpKWlqbMzMyAHzeYdu/erXHjxtkdwy9kDg0yh0YkZpYiMzeZQ6O8vFxVVVUBPWawf38DVoTNV253FJmZmSoqKrI7hl9ycnLIHAJkDg0yh87pcrvcx1RZ16D0pG5yOhJsSHZ6kfh3HamZgY6I4gwACJgNe6q1tLBCzaYUZUj5MzM0e+xAu2MBQEBwjTMAICBc7mOe0ixJzaaUV7hfLvcxe4MBQIBQnAEAAVFZ1+ApzS2aTFNVdUftCQQAAUZxBgAERHpSN0UZ3tuiDUNpSV3tCQQAAUZxBgAEhNORoPyZGYo2vm3P0Yah5TNHhd0NggDQXtwcCAAImNljB2rikL6qqjuqtKSulGYAHQrFGQAQUE5HAoUZQIfEpRoAAACABRRnAAAAwAKKMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACygOAMAAAAWUJwBAAAACyjOAAAAgAUUZwAAAMACijMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAAAsoDgDAAAAFlCcAQAAAAsozgAAAIAFFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABZQnAEAAAALKM4AAACABYZpmqbdIUKpvr5eDodDbrdbiYmJAT9+WlqaMjMzA37cYNq9e7fGjRtndwy/kDk0yBwakZhZiszckZL5mxNNOnq8SV27RKu8bG9EZP6u8vJyVVVVBfSYwf79DVgRY3eAjiYzM1NFRUV2x/BLTk4OmUOAzKFB5tCJxNyRkHnDnmotLaxQsylFGdKAhpVhn/lUOTk5dkcAgoJLNQAACBMu9zFPaZakZlM64KqXy33M3mAAJFGcAQAIG5V1DZ7S3MKUVFV31JY8ALxRnAEACBPpSd0UZXhvMySlJXW1JQ8AbxRnAADChNORoPyZGYo2vm3P0Yah4c5EOR0JNicDIHFzIAAAYWX22IGaOKSvquqOKi2pq/77+jV2RwLwb7aecc7Pz9fYsWPVo0cP9evXTzNmzNDBgwfPuF9JSYnGjBmj+Ph4DRo0SGvW8B8VAEDH4XQkKPu8PpxpBsKMrcW5pKRECxcu1LvvvqutW7fq5MmTmjp1qhoaGk67T2VlpaZPn64JEyaorKxMeXl5WrRokQoKCkKYHAAAAJ2NrZdqbN682ev9888/r379+mnv3r2aOHGiz33WrFmjgQMHauXKlZKk4cOHq7S0VI8++qhmzZoV7MgAAADopMLq5kC32y1J6t2792nH7Nq1S1OnTvXaduWVV6q0tFQnTpxoNb6xsVH19fVeLwAAAMBfYVOcTdNUbm6uLr30Uo0aNeq042pqapScnOy1LTk5WSdPnlRdXV2r8fn5+XI4HJ5XampqwLMDAACg4wub4nzLLbeovLxcr7zyyhnHGob3Qy5N0/S5XZKWLl0qt9vteR0+fDgwgQEAANCphMXj6G699VYVFRVpx44dGjBgQJtjU1JSVFNT47WttrZWMTEx6tOnT6vxcXFxiouLC2heAAAAdD62nnE2TVO33HKLCgsLtW3bNqWnp59xn+zsbG3dutVr25YtW5SVlaXY2NhgRQUAAEAnZ2txXrhwoV566SW9/PLL6tGjh2pqalRTU6Njx455xixdulRz5szxvF+wYIE+/fRT5ebm6sCBA3ruuee0bt06LV682I4pAAAAoJOwtTivXr1abrdbkydPltPp9Lw2bNjgGeNyuVRdXe15n56erk2bNqm4uFgXXHCBHnroIa1atYpH0QEAACCobL3GueWmvrasX7++1bZJkybpvffeC0IiAAAAwLeweaoGAAAAEM4ozgAAAIAFFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABb4VZyPHTumnTt36v3332/12TfffKMXX3wxYMEAAACAcGK5OH/44YcaPny4Jk6cqIyMDE2ePFkul8vzudvt1s9+9rOghAQAAADsZrk433XXXcrIyFBtba0OHjyoxMREjR8/XtXV1cHMBwAAAIQFy8X5nXfe0fLly5WUlKTBgwerqKhI06ZN04QJE/TJJ58EMyMAAABguxirA48dO6aYGO/hTz31lKKiojRp0iS9/PLLAQ8HAAAAhAvLxXnYsGEqLS3V8OHDvbY/8cQTMk1TOTk5AQ8He7ncx1RZ16D0pG5yOhLsjgMAAGAry5dq/OhHP9Irr7zi87Mnn3xS1157rUzTDFgw2GvDnmqN//U2XffMXzX+19u0YQ/XsgMAgM7NMDtZ262vr5fD4ZDb7VZiYmLAj5+WlqbMzMyAHzeYdu/erXHjxnnef3OiSf/fx3X67v9iGJLGD05SfGx0yPP5cmrmSEDm0CBz6ERibjKHRnl5uaqqqgJ6zGD//gassHypBqzJzMxUUVGR3TH8kpOT45X5nUN1uu6Zv7Ya98DPL1b2eX1CGe20Ts0cCcgcGmQOnUjMTebQ4PJNdFR8cyBaSU/qpijDe1u0YSgtqas9gQAAAMIAZ5zRitORoPyZGcor3K8m01S0YWj5zFHcIAgAgKSmpiadOHHC7hgIkNjYWEVHW7sUleIMn2aPHaiJQ/qqqu6o0pK6UpoBAJ2eaZqqqanRkSNH7I6CAOvZs6dSUlJkGEab4yjOOC2nI4HCDADAv7WU5n79+qlr165nLFkIf6Zp6ujRo6qtrZUkOZ3ONse3qzh/+OGHKi4uVm1trZqbm70+u/fee9tzSAAAgLDV1NTkKc19+oTHjfIIjISEb08S1tbWql+/fm1etuF3cX7mmWd00003KSkpqdUpbcMwKM4AAKDDabmmuWtXbpTviFrW9cSJE4Etzg8//LB+9atf6a677mp/OgAAgAjE5Rkdk9V19ftxdF999ZV+/OMf+x0IAAAAiGR+F+cf//jH2rJlSzCyAAAAIAwYhqHXX3+9zTHz5s3TjBkzLB+zqqpKhmFo3759Z5XNTn5fqjF48GDdc889evfdd5WRkaHY2FivzxctWhSwcAAAADg78+bN05EjR85YhL/L5XKpV69ekr4tvOnp6SorK9MFF1zgGfP444/LNM2AZp08ebIuuOACrVy50ufnhmFo48aNrQr7qXOsra3VPffcozfffFNffPGFevXqpdGjR+v+++9XdnZ2u/P5XZzXrl2r7t27q6SkRCUlJa0mQ3EGAACIbCkpKWcc43A4QpCkfWbNmqUTJ07ohRde0KBBg/TFF1/orbfe0j//+c+zOq7fxbmysvKsfiAAAEBn53IfU2Vdg9KTuoX8OxMmT56szMxMxcfH69lnn1WXLl20YMEC3X///Z4x3z2zm56eLkm68MILJUmTJk1ScXFxq7O8mzdv1sMPP6z9+/crOjpa2dnZevzxx3XeeeeFdH5HjhzRzp07VVxcrEmTJkmSzj33XI0bN+6sj+33Nc7fZZpmwE/RAwAAdGQb9lRr/K+36bpn/qrxv96mDXuqQ57hhRdeULdu3fTXv/5Vv/3tb/Xggw9q69atPsfu3r1bkvSXv/xFLpdLhYWFPsc1NDQoNzdXe/bs0VtvvaWoqCj96Ec/avWdH8HWvXt3de/eXa+//roaGxsDeux2FecXX3xRGRkZSkhIUEJCgjIzM/WHP/whoMEAAAA6Gpf7mJYWVqj53+cdm00pr3C/XO5jIc2RmZmp++67T+eff77mzJmjrKwsvfXWWz7H9u3bV5LUp08fpaSkqHfv3j7HzZo1SzNnztT555+vCy64QOvWrVNFRYXef//9gGa/9tprPeW45fXHP/7R83lMTIzWr1+vF154QT179tT48eOVl5en8vLys/7Zfhfnxx57TDfddJOmT5+u1157TRs2bNBVV12lBQsWaMWKFWcdCAAAoKOqrGvwlOYWTaapqrqjIc2RmZnp9d7pdHq+drq9Dh06pOuuu06DBg1SYmKi5xKP6urAnlFfsWKF9u3b5/XKycnxGjNr1ix9/vnnKioq0pVXXqni4mJddNFFWr9+/Vn9bL+vcX7iiSe0evVqzZkzx7Ptmmuu0ciRI3X//ffrjjvuOKtAAAAAHVV6UjdFGfIqz9GGobSk0H4j4alPRTMM46wvqbj66quVmpqqZ555Rv3791dzc7NGjRql48ePn9VxT5WSkqLBgwd7bevRo4eOHDnitS0+Pl5XXHGFrrjiCt17772aP3++7rvvPs2bN6/dP9vvM84ul0uXXHJJq+2XXHKJXC5Xu4MAAAB0dE5HgvJnZij6399UF20YWj5zVMhvEPRHly5dJElNTU2nHfPll1/qwIED+uUvf6kpU6Zo+PDh+uqrr0IV0ZIRI0aooaHhrI7Rruc4v/baa8rLy/PavmHDBp1//vlnFQYAAKCjmz12oCYO6auquqNKS+oa1qVZkvr166eEhARt3rxZAwYMUHx8fKtH0fXq1Ut9+vTR2rVr5XQ6VV1drbvvvrtdP+8f//hHqy9JSUlJsfSIPOnbEv/jH/9YN9xwgzIzM9WjRw+Vlpbqt7/9ra655pp2ZWrhd3F+4IEHNHv2bO3YsUPjx4+XYRjauXOn3nrrLb322mtnFQYAAKAzcDoSwr4wt4iJidGqVav04IMP6t5779WECRNUXFzsNSYqKkqvvvqqFi1apFGjRmno0KFatWqVJk+e7PfPe/nll/Xyyy97bbvvvvu8HpfXlu7du+t73/ueVqxYoUOHDunEiRNKTU3Vz3/+81Ynfv3ld3GeNWuW/vrXv2rFihV6/fXXZZqmRowYod27d3ue7wcAAIDwcOoNcaeWXkmtvlXw1McNz58/X/Pnz2/zuN///vdbPUHju8dJS0s742OMfWVrK5evLHFxccrPz1d+fn6bx2oPv4uzJI0ZM0YvvfRSoLMAAAAAYctSca6vr1diYqLnz21pGQcAAAB0JJaKc69eveRyudSvXz/17NlTxr/vBP0u0zRlGEabd1wCAAAAkcpScd62bZvnW2K2b98e1EAAAABAOLJUnCdNmuTzzwAAAJ3JmW5uQ2Syuq5+fwHK5s2btXPnTs/7p556ShdccIGuu+66sHvQNQAAQCC0fNPe0aOh/WpshEbLup76jYqn8vupGr/4xS/0m9/8RpJUUVGh3Nxc3Xnnndq2bZtyc3P1/PPPtyMuAABA+IqOjlbPnj1VW1srSeratavPe74QWUzT1NGjR1VbW6uePXsqOjq6zfF+F+fKykqNGDFCklRQUKCrr75ay5cv13vvvafp06e3LzUAAECYa/nmupbyjI6jZ8+elr6Z0O/i3KVLF8/p7L/85S+aM2eOJKl3795nfFTdqXbs2KFHHnlEe/fulcvl0saNGzVjxozTji8uLtZll13WavuBAwc0bNgwv342AACAPwzDkNPpVL9+/XTixAm74yBAYmNjz3imuYXfxfnSSy9Vbm6uxo8fr927d2vDhg2SpA8//FADBgzw61gNDQ0aPXq0fvazn2nWrFmW9zt48KDX86L79u3r188FAABor+joaMtFCx2L38X5ySef1M0336w//elPWr16tc455xxJ0ptvvqmrrrrKr2NNmzZN06ZN8zeC53nSAAAAQKj4XZwHDhyo//3f/221fcWKFQEJZMWFF16ob775RiNGjNAvf/lLn5dvtGhsbFRjY6Pnvb+XkwAAAABSO4qzJDU3N+vjjz9WbW2tmpubvT6bOHFiQIL54nQ6tXbtWo0ZM0aNjY36wx/+oClTpqi4uPi0Pzc/P18PPPBA0DKdqry8XDk5OSH7eYGwe/duMocAmUODzKETibnJHBrl5eV2RwCCw/TTrl27zPT0dDMqKso0DMPrFRUV5e/hPCSZGzdu9Hu/H/7wh+bVV1992s+/+eYb0+12e16HDx82JZlut7vdWdvSVpZwRebQIHNokDl0IjE3mUMjGJndbndQf38DVvh9xnnBggXKysrSG2+8IafTafszDC+++GK99NJLp/08Li5OcXFxIUwEAACAjsjv4vzRRx/pT3/6kwYPHhyMPH4rKyuT0+m0OwYAAAA6OL+L8/e+9z19/PHHASnO//rXv/Txxx973ldWVmrfvn3q3bu3Bg4cqKVLl+qzzz7Tiy++KElauXKl0tLSNHLkSB0/flwvvfSSCgoKVFBQcNZZAAAAgLb4XZxvvfVW3XnnnaqpqVFGRkar7/TOzMy0fKzS0lKvJ2Lk5uZKkubOnav169fL5XKpurra8/nx48e1ePFiffbZZ0pISNDIkSP1xhtv8I2FAAAACDq/i3PLF5XccMMNnm2GYcg0TRmGoaamJsvHmjx5skzTPO3n69ev93q/ZMkSLVmyxL/AAAAAQAD4XZwrKyuDkQMAAAAIa34X53PPPTcYOQAAAICwFtWenf7whz9o/Pjx6t+/vz799FNJ39649+c//zmg4QAAAIBw4XdxXr16tXJzczV9+nQdOXLEc01zz549tXLlykDnAwAAAMKC38X5iSee0DPPPKNly5YpOjrasz0rK0sVFRUBDQcAAACEC7+Lc2VlpS688MJW2+Pi4tTQ0BCQUJBc7mN651CdXO5jdkcBAACA2nFzYHp6uvbt29fqJsE333xTI0aMCFiwzmzDnmotLaxQsylFGVL+zAzNHjvQ7lgAAACdmt/F+Re/+IUWLlyob775RqZpavfu3XrllVeUn5+vZ599NhgZOxWX+5inNEtSsynlFe7XxCF95XQk2BsOAACgE/O7OP/sZz/TyZMntWTJEh09elTXXXedzjnnHD3++OP6yU9+EoyMnUplXYOnNLdoMk1V1R2lOAMAANjI7+IsST//+c/185//XHV1dWpubla/fv0CnavTSk/qpihDXuU52jCUltTVvlAAAABo33OcWyQlJVGaA8zpSFD+zAxFG4akb0vz8pmjONsMAABgM7/POH/55Ze69957tX37dtXW1qq5udnr83/+858BC9dZzR47UBOH9FVV3VGlJXWlNAMAAIQBv4vzf/3Xf+nQoUO68cYblZycLOPfZ0YRWE5HAoUZAAAgjPhdnHfu3KmdO3dq9OjRwcgDAAAAhCW/r3EeNmyYjh3jSzkAAADQufhdnJ9++mktW7ZMJSUl+vLLL1VfX+/1AgAAADoivy/V6Nmzp9xuty6//HKv7aZpyjAMNTU1BSwcAAAAEC78Ls4//elP1aVLF7388svcHAgAAIBOw+/ivH//fpWVlWno0KHByAMAAACEJb+vcc7KytLhw4eDkQUAAAAIW36fcb711lt122236Re/+IUyMjIUGxvr9XlmZmbAwgEAAADhwu/iPHv2bEnSDTfc4NlmGAY3BwIAAKBD87s4V1ZWBiMHAAAAENb8Ls7nnntuMHIAAAAAYc1ScS4qKtK0adMUGxuroqKiNsfm5OQEJBgAAAAQTiwV5xkzZqimpkb9+vXTjBkzTjuOa5wBAADQUVkqzs3NzT7/DAAAAHQWhmmapt0hQqm+vl4Oh0Nut1uJiYkBP35aWlrEPZJv9+7dGjdunN0x/ELm0CBzaERiZikyc5M5NMrLy1VVVRXQYwb79zdghV83BzY3N2v9+vUqLCxUVVWVDMNQenq6/uM//kPXX389X7+tb59jfabrwMNNTk4OmUOAzKFB5sByuY+psq5B6Und5HQkeH0WzrlPh8yhwf1O6KgsF2fTNJWTk6NNmzZp9OjRysjIkGmaOnDggObNm6fCwkK9/vrrQYwKAAilDXuqtbSwQs2mFGVI+TMzNHvsQLtjAYBtLBfn9evXa8eOHXrrrbd02WWXeX22bds2zZgxQy+++KLmzJkT8JAAgNByuY95SrMkNZtSXuF+TRzSt9WZZwDoLKKsDnzllVeUl5fXqjRL0uWXX667775bf/zjHwMaDgBgj8q6Bk9pbtFkmqqqO2pPIAAIA5aLc3l5ua666qrTfj5t2jT97W9/C0goAIC90pO6KeqU21aiDUNpSV3tCQQAYcBycf7nP/+p5OTk036enJysr776KiChAAD2cjoSlD8zQ9H/vuk72jC0fOYoLtMA0KlZvsa5qalJMTGnHx4dHa2TJ08GJBQAwH6zxw7UxCF9VVV3VGlJXSnNADo9v56qMW/ePMXFxfn8vLGxMWChAADhwelIoDADwL9ZLs5z58494xieqAEAAICOynJxfv7554OZAwAAAAhrlm8OBAAAADozijMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAAAsoDgDAAAAFlCcAQAAAAtsLc47duzQ1Vdfrf79+8swDL3++utn3KekpERjxoxRfHy8Bg0apDVr1gQ/KAAAADo9W4tzQ0ODRo8erSeffNLS+MrKSk2fPl0TJkxQWVmZ8vLytGjRIhUUFAQ5KQAAADo7y98cGAzTpk3TtGnTLI9fs2aNBg4cqJUrV0qShg8frtLSUj366KOaNWtWkFICAAAAEXaN865duzR16lSvbVdeeaVKS0t14sQJn/s0Njaqvr7e6wUAAAD4K6KKc01NjZKTk722JScn6+TJk6qrq/O5T35+vhwOh+eVmpoaiqgAAADoYCKqOEuSYRhe703T9Lm9xdKlS+V2uz2vw4cPBz0jAAAAOh5br3H2V0pKimpqary21dbWKiYmRn369PG5T1xcnOLi4kIRDwAAAB1YRJ1xzs7O1tatW722bdmyRVlZWYqNjbUpFQAAADoDW4vzv/71L+3bt0/79u2T9O3j5vbt26fq6mpJ315mMWfOHM/4BQsW6NNPP1Vubq4OHDig5557TuvWrdPixYvtiA8AAIBOxNZLNUpLS3XZZZd53ufm5kqS5s6dq/Xr18vlcnlKtCSlp6dr06ZNuuOOO/TUU0+pf//+WrVqFY+iAwAAQNDZWpwnT57subnPl/Xr17faNmnSJL333ntBTAUAAAC0FlHXOAMAAAB2oTgDAAAAFlCcAQAAAAsozgAAAIAFhtnW3XkdUH19vRwOh9xutxITEwN+/LS0NGVmZgb8uMG0e/dujRs3zu4YfiFzaJA5NCIxsxSZuckcGuXl5aqqqgroMYP9+xuwIqK+OTASZGZmqqioyO4YfsnJySFzCJA5NMgcOpGYm8yhkZOTY3cEICi4VAMAAACwgOIMAAAAWEBxBgAAACygOAMAAAAWUJwBAAAACyjOAAAAgAUUZwAAAMACijMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAAAsoDgDCDqX+5jeOVQnl/uY3VEAAGi3GLsDAOjYNuyp1tLCCjWbUpQh5c/M0OyxA+2OBQCA3zjjDCBoXO5jntIsSc2mlFe4nzPPAICIRHEGEDSVdQ2e0tyiyTRVVXfUnkAAAJwFijOAoElP6qYow3tbtGEoLamrPYEAADgLFGcAQeN0JCh/ZoaijW/bc7RhaPnMUXI6EmxOBgCA/7g5EEBQzR47UBOH9FVV3VGlJXWlNAMAIhbFGUDQOR0JFGYAQMTjUg0AAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABZQnAEAAAALKM4AAACABRRnAAAAwAKKMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACwwTNM07Q4RSvX19XI4HHK73UpMTAz48dPS0pSZmRnw4wbT7t27NW7cOLtj+IXMoUHm0IjEzFJk5iZzaJSXl6uqqiqgxwz272/Aihi7A3Q0mZmZKioqsjuGX3JycsgcAmQODTKHTiTmJnNo5OTk2B0BCAou1QAAAAAsoDgDAAAAFlCcAQAAAAsozgAAAIAFFGcAAADAAtuL89NPP6309HTFx8drzJgxevvtt087tri4WIZhtHp98MEHIUwMAACAzsjW4rxhwwbdfvvtWrZsmcrKyjRhwgRNmzZN1dXVbe538OBBuVwuz+v8888PUWIAAAB0VrYW58cee0w33nij5s+fr+HDh2vlypVKTU3V6tWr29yvX79+SklJ8byio6NDlBgAAACdlW3F+fjx49q7d6+mTp3qtX3q1Kl655132tz3wgsvlNPp1JQpU7R9+/Y2xzY2Nqq+vt7rBQAAAPjLtuJcV1enpqYmJScne21PTk5WTU2Nz32cTqfWrl2rgoICFRYWaujQoZoyZYp27Nhx2p+Tn58vh8PheaWmpgZ0HgAAAOgcbP/KbcMwvN6bptlqW4uhQ4dq6NChnvfZ2dk6fPiwHn30UU2cONHnPkuXLlVubq7nfX19PeUZAAAAfrPtjHNSUpKio6NbnV2ura1tdRa6LRdffLE++uij034eFxenxMRErxcAAADgL9uKc5cuXTRmzBht3brVa/vWrVt1ySWXWD5OWVmZnE5noOMBAAAAXmy9VCM3N1fXX3+9srKylJ2drbVr16q6uloLFiyQ9O1lFp999plefPFFSdLKlSuVlpamkSNH6vjx43rppZdUUFCggoICO6cBAACATsDW4jx79mx9+eWXevDBB+VyuTRq1Cht2rRJ5557riTJ5XJ5PdP5+PHjWrx4sT777DMlJCRo5MiReuONNzR9+nS7pgAAAIBOwvabA2+++WbdfPPNPj9bv3691/slS5ZoyZIlIUgFAAAAeLP9K7cBAACASEBxBgAAACygOAMAAAAWUJwBAAAACyjOAAAAgAUUZwAAAMACijMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAK243Mf0zqE6udzH7I4ChA3bv3IbAACElw17qrW0sELNphRlSPkzMzR77EC7YwG244wzAADwcLmPeUqzJDWbUl7hfs48A6I4AwCA76isa/CU5hZNpqmquqP2BALCCMUZAAB4pCd1U5ThvS3aMJSW1NWeQEAYoTgDAAAPpyNB+TMzFG18256jDUPLZ46S05FgczLAfoZpmuaZh3Uc9fX1cjgccrvdSkxMDPjx09LSlJmZGfDjBtPu3bs1btw4u2P4hcyhQebQiMTMUmTmJrN135xo0rHjTUroEq342Gi/9i0vL1dVVVVA8wT79zdgBU/VCLDMzEwVFRXZHcMvOTk5ZA4BMocGmUMnEnOTOTRycnLsjgAEBZdqAAAAABZQnAEAAAALKM4AAACABRRnAAAAwAKKMwAAAGABxRkIEy73Mb1zqI6vtQUAIEzxODogDGzYU62lhRVqNqUoQ8qfmaHZYwfaHQsAAHwHZ5wBm7ncxzylWZKaTSmvcD9nngEACDMUZ8BmlXUNntLcosk0VVV31J5AAADAJ4ozYLP0pG6KMry3RRuG0pK62hMIAAD4RHEGbOZ0JCh/ZoaijW/bc7RhaPnMUXI6EmxOBgAAvoubA4EwMHvsQE0c0ldVdUeVltSV0gwAQBiiOANhwulIoDADABDGuFQDAAAAsIDiDAAAAFhAcQYAAAAsoDgDAAAAFlCcAQAAAAsozgAAAIAFFGcAAADAAoozAAAAYAHFGQAAALCA4gwAAABYQHEGAAAALKA4AwAAABZQnAEAAAALKM4AAACABRRnAAAAwALbi/PTTz+t9PR0xcfHa8yYMXr77bfbHF9SUqIxY8YoPj5egwYN0po1a0KUFAAAAJ2ZrcV5w4YNuv3227Vs2TKVlZVpwoQJmjZtmqqrq32Or6ys1PTp0zVhwgSVlZUpLy9PixYtUkFBQYiTAwAAoLOxtTg/9thjuvHGGzV//nwNHz5cK1euVGpqqlavXu1z/Jo1azRw4ECtXLlSw4cP1/z583XDDTfo0UcfDXFyAAAAdDa2Fefjx49r7969mjp1qtf2qVOn6p133vG5z65du1qNv/LKK1VaWqoTJ0743KexsVH19fVeLwAAAMBfMXb94Lq6OjU1NSk5Odlre3JysmpqanzuU1NT43P8yZMnVVdXJ6fT2Wqf/Px8PfDAA4ELfgaHDh1STk5OyH5eIJSXl5M5BMgcGmQOnUjMTebQOHTokN0RgKCwrTi3MAzD671pmq22nWm8r+0tli5dqtzcXM/7+vp6paamtjfuGf39738P2rEBAABgH9uKc1JSkqKjo1udXa6trW11VrlFSkqKz/ExMTHq06ePz33i4uIUFxcXmNAAAADotGy7xrlLly4aM2aMtm7d6rV969atuuSSS3zuk52d3Wr8li1blJWVpdjY2KBlBQAAAGx9qkZubq6effZZPffcczpw4IDuuOMOVVdXa8GCBZK+vcxizpw5nvELFizQp59+qtzcXB04cEDPPfec1q1bp8WLF9s1BQAAAHQStl7jPHv2bH355Zd68MEH5XK5NGrUKG3atEnnnnuuJMnlcnk90zk9PV2bNm3SHXfcoaeeekr9+/fXqlWrNGvWLLumAAAAgE7CMFvurusk6uvr5XA45Ha7lZiYaHccAABgAb+/EQ5s/8ptAAAAIBJQnAEAAAALKM4AAACABRRnAAAAwAKKMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACygOAMAAAAW2PqV23Zo+aLE+vp6m5MAAACrWn5vd7IvPEaY6XTF+euvv5Ykpaam2pwEAAD46+uvv5bD4bA7Bjopw+xk/3Rrbm7W559/rh49esgwjIAeu76+XqmpqTp8+LASExMDeuxw0NHnJ3X8OTK/yNfR58j8Il+w5miapr7++mv1799fUVFcaQp7dLozzlFRURowYEBQf0ZiYmKH/Q+i1PHnJ3X8OTK/yNfR58j8Il8w5siZZtiNf7IBAAAAFlCcAQAAAAsozgEUFxen++67T3FxcXZHCYqOPj+p48+R+UW+jj5H5hf5OsMc0Xl1upsDAQAAgPbgjDMAAABgAcUZAAAAsIDiDAAAAFhAcQYAAAAsoDj76emnn1Z6erri4+M1ZswYvf32222OLykp0ZgxYxQfH69BgwZpzZo1IUraPv7Mr7i4WIZhtHp98MEHIUxs3Y4dO3T11Verf//+MgxDr7/++hn3ibT183eOkbSG+fn5Gjt2rHr06KF+/fppxowZOnjw4Bn3i6Q1bM8cI2kNV69erczMTM8XY2RnZ+vNN99sc59IWj9/5xdJa+dLfn6+DMPQ7bff3ua4SFpD4Ewozn7YsGGDbr/9di1btkxlZWWaMGGCpk2bpurqap/jKysrNX36dE2YMEFlZWXKy8vTokWLVFBQEOLk1vg7vxYHDx6Uy+XyvM4///wQJfZPQ0ODRo8erSeffNLS+EhbP8n/ObaIhDUsKSnRwoUL9e6772rr1q06efKkpk6dqoaGhtPuE2lr2J45toiENRwwYIB+/etfq7S0VKWlpbr88st1zTXX6O9//7vP8ZG2fv7Or0UkrN2p9uzZo7Vr1yozM7PNcZG2hsAZmbBs3Lhx5oIFC7y2DRs2zLz77rt9jl+yZIk5bNgwr23//d//bV588cVBy3g2/J3f9u3bTUnmV199FYJ0gSXJ3LhxY5tjIm39TmVljpG8hrW1taYks6Sk5LRjIn0NrcwxktfQNE2zV69e5rPPPuvzs0hfP9Nse36RunZff/21ef7555tbt241J02aZN52222nHdsR1hD4Ls44W3T8+HHt3btXU6dO9do+depUvfPOOz732bVrV6vxV155pUpLS3XixImgZW2P9syvxYUXXiin06kpU6Zo+/btwYwZUpG0fmcrEtfQ7XZLknr37n3aMZG+hlbm2CLS1rCpqUmvvvqqGhoalJ2d7XNMJK+flfm1iLS1W7hwoX7wgx/o+9///hnHRvIaAr5QnC2qq6tTU1OTkpOTvbYnJyerpqbG5z41NTU+x588eVJ1dXVBy9oe7Zmf0+nU2rVrVVBQoMLCQg0dOlRTpkzRjh07QhE56CJp/dorUtfQNE3l5ubq0ksv1ahRo047LpLX0OocI20NKyoq1L17d8XFxWnBggXauHGjRowY4XNsJK6fP/OLtLWTpFdffVV79+5Vfn6+pfGRuIZAW2LsDhBpDMPwem+aZqttZxrva3u48Gd+Q4cO1dChQz3vs7OzdfjwYT366KOaOHFiUHOGSqStn78idQ1vueUWlZeXa+fOnWccG6lraHWOkbaGQ4cO1b59+3TkyBEVFBRo7ty5KikpOW25jLT182d+kbZ2hw8f1m233aYtW7YoPj7e8n6RtoZAWzjjbFFSUpKio6NbnX2tra1t9a/pFikpKT7Hx8TEqE+fPkHL2h7tmZ8vF198sT766KNAx7NFJK1fIIX7Gt56660qKirS9u3bNWDAgDbHRuoa+jNHX8J5Dbt06aLBgwcrKytL+fn5Gj16tB5//HGfYyNx/fyZny/hvHZ79+5VbW2txowZo5iYGMXExKikpESrVq1STEyMmpqaWu0TiWsItIXibFGXLl00ZswYbd261Wv71q1bdckll/jcJzs7u9X4LVu2KCsrS7GxsUHL2h7tmZ8vZWVlcjqdgY5ni0hav0AK1zU0TVO33HKLCgsLtW3bNqWnp59xn0hbw/bM0ZdwXUNfTNNUY2Ojz88ibf18aWt+voTz2k2ZMkUVFRXat2+f55WVlaWf/vSn2rdvn6Kjo1vt0xHWEPBiyy2JEerVV181Y2NjzXXr1pnvv/++efvtt5vdunUzq6qqTNM0zbvvvtu8/vrrPeM/+eQTs2vXruYdd9xhvv/+++a6devM2NhY809/+pNdU2iTv/NbsWKFuXHjRvPDDz809+/fb959992mJLOgoMCuKbTp66+/NsvKysyysjJTkvnYY4+ZZWVl5qeffmqaZuSvn2n6P8dIWsObbrrJdDgcZnFxselyuTyvo0ePesZE+hq2Z46RtIZLly41d+zYYVZWVprl5eVmXl6eGRUVZW7ZssU0zchfP3/nF0lrdzqnPlUj0tcQOBOKs5+eeuop89xzzzW7dOliXnTRRV6PiZo7d645adIkr/HFxcXmhRdeaHbp0sVMS0szV69eHeLE/vFnfr/5zW/M8847z4yPjzd79eplXnrppeYbb7xhQ2prWh79dOpr7ty5pml2jPXzd46RtIa+5iXJfP755z1jIn0N2zPHSFrDG264wfPfl759+5pTpkzxlErTjPz183d+kbR2p3NqcY70NQTOxDDNf1+lDwAAAOC0uMYZAAAAsIDiDAAAAFhAcQYAAAAsoDgDAAAAFlCcAQAAAAsozgAAAIAFFGcAAADAAoozAAAAYAHFGUDAGYah119/3e4YbSouLpZhGDpy5IjdUQAAEYLiDMCSefPmyTAMGYah2NhYJScn64orrtBzzz2n5uZmr7Eul0vTpk2zKak1l1xyiVwulxwOR1B/zo4dO3T11Verf//+EfEPCgDA6VGcAVh21VVXyeVyqaqqSm+++aYuu+wy3XbbbfrhD3+okydPesalpKQoLi7OxqRn1qVLF6WkpMgwjKD+nIaGBo0ePVpPPvlkUH8OACD4KM4ALIuLi1NKSorOOeccXXTRRcrLy9Of//xnvfnmm1q/fr1n3HfPrFZVVckwDL322muaMGGCEhISNHbsWH344Yfas2ePsrKy1L17d1111VX6xz/+4fXznn/+eQ0fPlzx8fEaNmyYnn76ac9nLcctLCzUZZddpq5du2r06NHatWuXZ8ynn36qq6++Wr169VK3bt00cuRIbdq0SZLvSzUKCgo0cuRIxcXFKS0tTb/73e+88qSlpWn58uW64YYb1KNHDw0cOFBr165t8+9s2rRpevjhhzVz5kx//qoBAGGI4gzgrFx++eUaPXq0CgsL2xx333336Ze//KXee+89xcTE6Nprr9WSJUv0+OOP6+2339ahQ4d07733esY/88wzWrZsmX71q1/pwIEDWr58ue655x698MILXsddtmyZFi9erH379mnIkCG69tprPWe/Fy5cqMbGRu3YsUMVFRX6zW9+o+7du/vMt3fvXv3nf/6nfvKTn6iiokL333+/7rnnHq9/EEjS7373O2VlZamsrEw333yzbrrpJn3wwQft+JsDAESaGLsDAIh8w4YNU3l5eZtjFi9erCuvvFKSdNttt+naa6/VW2+9pfHjx0uSbrzxRq+S+tBDD+l3v/ud50xtenq63n//ff3+97/X3LlzvY77gx/8QJL0wAMPaOTIkfr44481bNgwVVdXa9asWcrIyJAkDRo06LT5HnvsMU2ZMkX33HOPJGnIkCF6//339cgjj2jevHmecdOnT9fNN98sSbrrrru0YsUKFRcXa9iwYVb+qgAAEYwzzgDOmmmaZ7xWODMz0/Pn5ORkSfIU2pZttbW1kqR//OMfOnz4sG688UZ1797d83r44Yd16NCh0x7X6XRKkuc4ixYt0sMPP6zx48frvvvua7PcHzhwwFPiW4wfP14fffSRmpqafP48wzCUkpLi+XkAgI6N4gzgrB04cEDp6eltjomNjfX8uaVkn7qt5ekcLf/zmWee0b59+zyv/fv369133z3jcVv2nz9/vj755BNdf/31qqioUFZWlp544gmf+XyVf9M025zHqbkBAB0bxRnAWdm2bZsqKio0a9asgB0zOTlZ55xzjj755BMNHjzY63Wmgn6q1NRULViwQIWFhbrzzjv1zDPP+Bw3YsQI7dy502vbO++8oyFDhig6OrrdcwEAdBxc4wzAssbGRtXU1KipqUlffPGFNm/erPz8fP3whz/UnDlzAvqz7r//fi1atEiJiYmaNm2aGhsbVVpaqq+++kq5ubmWjnH77bdr2rRpGjJkiL766itt27ZNw4cP9zn2zjvv1NixY/XQQw9p9uzZ2rVrl5588kmvJ3m0x7/+9S99/PHHnveVlZXat2+fevfurYEDB57VsQEAoUVxBmDZ5s2b5XQ6FRMTo169emn06NFatWqV5s6dq6iowP4/sObPn6+uXbvqkUce0ZIlS9StWzdlZGTo9ttvt3yMpqYmLVy4UP/3f/+nxMREXXXVVVqxYoXPsRdddJFee+013XvvvXrooYfkdDr14IMPet0Y2B6lpaW67LLLPO9bSv/cuXNbPbEDABDeDNPXRXwAAAAAvHCNMwAAAGABxRkAAACwgOIMAAAAWEBxBgAAACygOAMAAAAWUJwBAAAACyjOAAAAgAUUZwAAAMACijMAAABgAcUZAAAAsIDiDAAAAFjw/wO7HazN18I5cgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from smt.sampling_methods import LHS\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "xlimits = np.array([[0.0, 4.0], [0.0, 3.0], [0.0, 3.0], [1.0, 5.0]])\n",
    "sampling = LHS(xlimits=xlimits, criterion='ese', random_state=1)\n",
    "\n",
    "num = 10\n",
    "x = sampling(num)\n",
    "\n",
    "\n",
    "### For visualization only\n",
    "\n",
    "intervals = []  \n",
    "subspace_bool = []\n",
    "for i in range(len(xlimits)):  \n",
    "    intervals.append(np.linspace(xlimits[i][0],xlimits[i][1],num+1))\n",
    "    subspace_bool.append([[intervals[i][j]<x[kk][i]<intervals[i][j+1] \n",
    "    for kk in range(len(x))] for j in range(len(intervals[i])-1)])\n",
    "    \n",
    "PLOT1 = 0 #First dimension to be depicted\n",
    "PLOT2 = 1 #Second dimension to be depicted\n",
    "\n",
    "fig, ax = plt.subplots(1)\n",
    "ax.plot(x[:, PLOT1], x[:, PLOT2], \".\")\n",
    "for i in range(len(intervals[0])):\n",
    "    ax.plot([intervals[PLOT1][i], intervals[PLOT1][i]], [intervals[PLOT2][0],intervals[PLOT2][-1]], linewidth=0.5,color='k')\n",
    "    ax.plot([intervals[PLOT1][0], intervals[PLOT1][-1]], [intervals[PLOT2][i],intervals[PLOT2][i]], linewidth=0.5,color='k')\n",
    "ax.set_xlabel(\"Dimension \" + str(PLOT1+1))\n",
    "ax.set_ylabel(\"Dimension \" + str(PLOT2+1))\n",
    "ax.legend(['Initial LHS'],bbox_to_anchor = (1.05, 0.6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Augmented Sampling - Basic Methodology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we expand our LHS sampling (10 points) by adding another 10 points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvkAAAGwCAYAAAAzA3pBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOwklEQVR4nO3de3wU9b3/8feQGwjshhhzIYQQEEJEglzLRRSlcmtBKuVotQIVPUARhBTUwPFuTc/xBnjBomig3tAGLZyiD1Au0QMWg6FBQUAghB8mhohkacAEwvz+SLO4JEAWdnazk9fz8dhH3NnvfPKenYz5ZPjOrGGapikAAAAAttEk0AEAAAAA+BZNPgAAAGAzNPkAAACAzdDkAwAAADZDkw8AAADYDE0+AAAAYDM0+QAAAIDNhAY6gL+dOnVK3377rVq2bCnDMAIdBwAA1INpmjp69Khat26tJk04RwmcT6Nr8r/99lslJiYGOgYAALgABw4cUJs2bQIdA2jwGl2T37JlS0nV/5NwOBwBTgMAAOrD5XIpMTHR/XscwLk1uia/ZoqOw+GgyQcAIMgw1RaoHya1AQAAADZDkw8AAADYDE0+AAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANgMTT4AAABgMzT5AAAAgM3Q5AMAAAA2Q5MPAAAA2ExAm/yFCxcqLS1NDodDDodD/fr10wcffHDOdTZs2KCePXuqadOmat++vV566SU/pQUAAACCQ0Cb/DZt2uhPf/qTcnNzlZubq+uvv1433nijvvrqqzrH79u3TyNGjNDAgQOVl5enOXPmaPr06crOzvZzcgAAAKDhMkzTNAMd4qeioqL05JNPauLEibVeu++++7RixQrt2LHDvWzy5Mn65z//qU2bNtWrvsvlktPpVFlZmRwOh89yAwAA6/D7G/BOg5mTX1VVpbffflvl5eXq169fnWM2bdqkIUOGeCwbOnSocnNzdeLEiTrXqaiokMvl8ngAAAAAdhYa6ADbtm1Tv3799OOPP6pFixZ67733dMUVV9Q5tri4WLGxsR7LYmNjdfLkSZWWlio+Pr7WOpmZmXrkkUcsyV6XxI4piktM8nnd/Tt3KCkl1ed1g7U2mf1Tm8z+qR2Mma2sTWb/1LYy87HvDp516i0A/wh4k5+SkqKtW7fqyJEjys7O1vjx47Vhw4azNvqGYXg8r5ltdObyGhkZGUpPT3c/d7lcSkxM9FH62uISk5SxcInP62ZOGW9J3WCtTWb/1Cazf2oHY2Yra5PZP7WtzJw1e5IldQHUX8Cb/PDwcF1++eWSpF69eunzzz/X/Pnz9ec//7nW2Li4OBUXF3ssKykpUWhoqC699NI660dERCgiIsL3wQGgHo6UV+qQq0KXOSIU2Tw80HEAAI1EwJv8M5mmqYqKijpf69evn1auXOmxbPXq1erVq5fCwsL8EQ8A6u2zXaV6d2OhTFMyDGls/7bq2yk60LEAAI1AQC+8nTNnjj755BMVFBRo27Ztmjt3rtavX6/bbrtNUvVUm3HjxrnHT548Wfv371d6erp27NihV199VYsXL9asWbMCtQkAUKcj5ZXuBl+STFN6d2OhjpRXBjYYAKBRCOiZ/O+++0633367ioqK5HQ6lZaWpg8//FA33HCDJKmoqEiFhYXu8cnJyVq1apVmzpypF154Qa1bt9aCBQs0ZsyYQG0CANTpkKtCZ96g2DSlUlcF03YAAJYLaJO/ePHic76elZVVa9m1116rL774wqJEAOAblzkiZBjyaPQNQ4p2cI0QAMB6DeY++QBgJ5HNwzW2f1vV3PirZk4+Z/EBAP7Q4C68BQC76NspWp0THCp1VSiau+sAAPyIJh8ALBTZPJzmHgDgd0zXAQAAAGyGJh8AAACwGZp8AAAAwGZo8gEAAACbockHAAAAbIYmHwAAALAZmnwAAADAZgzT/OmHrtufy+WS0+lUWVmZHA6Hz+vHJCQqKSXV53UPbM9Xnz59fF5XknI2blLHtO6W1LYqN5k9kfm0zZs3B92xEqzH9+78vKB7P6z6+QjGY8XKzIf27lJBQYFPa1r9+xuwGz4My8eSUlKVsXCJz+tmzZ6kFStW+LyuJPW+foglmSXrcpPZE5lPGzVqVNAdK8F6fGdOGR9074dVPx/BeKxYnRlAYDFdBwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABshiYfAAAAsBmafAAAAMBmaPIBAAAAm6HJB9BwlR2U9uVUfwUAAPXGh2EBaJi+WCqtvEcyT0lGE2nkfKnHuECnAgAgKHAmH0DDU3bwdIMvVX9dOYMz+gAA1BNNPoCG5/Ce0w1+DbNKOrw3MHkAAAgyNPkAGp6oDtVTdH7KCJGi2gcmDwAAQYYmH0DD40yonoNvhFQ/N0KkkfOqlwMAgPPiwlsADVOPcVKHwdVTdKLa0+ADAOAFmnwADZczgeYeAIALYJimaQY6hD+5XC45nU6VlZXJ4XD4vH5MQqKSUlJ9Xnd3fp46pnX3ed1grU1m/9Qms39qB2NmK2uT2T+1rcx8aO8uFRQU+LSm1b+/AbvhTL6PJaWkKmPhEp/XzZwy3pK6wVqbzP6pTWb/1A7GzFbWJrN/aluZOWv2JEvqAqg/LrwFAAAAbIYmHwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABshiYfAAAAsBmafAAAAMBmaPIBAAAAm6HJBwAAAGyGJh8AAACwGZp8NFpHyiu1u+iojpRXBjoKAACAT4UGOgAQCJ/tKtW7GwtlmpJhSGP7t1XfTtGBjgUAAOATnMlHo3OkvNLd4EuSaUrvbizkjD4AALANmnw0OodcFe4Gv4ZpSqWuisAEAgAA8DHDNM9sd+zN5XLJ6XSqrKxMDofD5/VjEhKVlJLq87oHtuerT58+Pq8rSTkbN6ljWndLaluV+2Iyn6gytetbV63lnVo7FBZiNMjM50Pm0zZv3hx0x0qwHt+78/OC7v2w6ucjGI8VKzMf2rtLBQUFPq1p9e9vwG6Yk+9jSSmpyli4xOd1s2ZP0ooVK3xeV5J6Xz/EksySdbkvNvO55uQ31MznQubTRo0aFXTHSrAe35lTxgfd+2HVz0cwHitWZwYQWAGdrpOZmanevXurZcuWiomJ0ejRo7Vz585zrrN+/XoZhlHr8fXXX/spNeygb6doPTD2Sv1+WEc9MPZKLroFAAC2EtAz+Rs2bNDUqVPVu3dvnTx5UnPnztWQIUO0fft2NW/e/Jzr7ty50+Of6y677DKr48JmIpuHK7J5eKBjAAAA+FxAm/wPP/zQ4/lrr72mmJgYbdmyRddcc805142JiVFkZKSF6QAAAIDg1KDurlNWViZJioqKOu/Y7t27Kz4+XoMHD9a6devOOq6iokIul8vjAQAAANhZg2nyTdNUenq6rr76al155ZVnHRcfH69FixYpOztby5cvV0pKigYPHqycnJw6x2dmZsrpdLofiYmJVm0CAAAA0CA0mLvr3H333crPz9enn356znEpKSlKSUlxP+/Xr58OHDigp556qs4pPhkZGUpPT3c/d7lcNPoAAACwtQZxJn/atGlasWKF1q1bpzZt2ni9ft++fbV79+46X4uIiJDD4fB4AI1K2UFpX071VwAA0CgE9Ey+aZqaNm2a3nvvPa1fv17JyckXVCcvL0/x8fE+TgfYwBdLpZX3SOYpyWgijZwv9RgX6FQAAMBiAW3yp06dqjfffFN/+9vf1LJlSxUXF0uSnE6nmjVrJql6us3Bgwe1dOlSSdK8efPUrl07denSRZWVlXr99deVnZ2t7OzsgG0H0CCVHTzd4EvVX1fOkDoMlpwJAY0GAACsFdAmf+HChZKkQYMGeSx/7bXXNGHCBElSUVGRCgsL3a9VVlZq1qxZOnjwoJo1a6YuXbro73//u0aMGOGv2EBwOLzndINfw6ySDu+lyQcAwOYCPl3nfLKysjye33vvvbr33nstSgTYSFSH6ik6P230jRApqn3gMgEAAL9oEBfeArCAM6F6Dr4RUv3cCJFGzuMsPgAAjUCDuYUmAAv0GFc9B//w3uoz+DT4AAA0CjT5gN05E2juAQBoZJiuAwAAANiMYdbn6lcbcblccjqdKisrs+SDsWISEpWUkurzurvz89QxrbvP6wZrbTL7pzaZ/VM7GDNbWZvM/qltZeZDe3epoKDApzWt/v0N2A3TdXwsKSVVGQuX+Lxu5pTxltQN1tpk9k9tMvundjBmPlftI+WVOuSq0GWOCEU2D/dZXV9gH1pfV5KyZk+ypC6A+qPJBwD4zGe7SvXuxkKZpmQY0tj+bdW3U3SgYwFAo8OcfACATxwpr3Q3+JJkmtK7Gwt1pLwysMEAoBGiyQcA+MQhV4XOvMrLNKVSV0VgAgFAI0aTDwDwicscETIMz2WGIUU7IgITCAAaMZp8AIBPRDYP19j+bd2Nfs2c/Au5+BYAcHG48BYA4DN9O0Wrc4JDpa4KRV/g3XUAABePJh8A4FORzcNp7gEgwJiuAwAAANgMTT4AAABgMzT5AAAAgM3Q5AMAAAA2Q5MPAAAA2AxNPgAAAGAz3EITAIBGpNmxYjU/WqCQqopARwFgIZp8AAAaiaQ9f1WPzQ/K0CnFFB9T0p6/an+HXwc6FgALGKZpmoEO4U8ul0tOp1NlZWVyOBw+rx+TkKiklFSf1z2wPV99+vTxeV1Jytm4SR3TultS26rcZPZE5tM2b94cdMdKsB7fu/Pzgu79sOrnIyiOlRM/SnvXuZ9uPlil3gkhOhTXX1UhERdf/ycO7d2lgoICn9a0+vc3YDecyfexpJRUZSxc4vO6WbMnacWKFT6vK0m9rx9iSWbJutxk9kTm00aNGhV0x0qwHt+ZU8YH3fth1c9HUBwr+3KkJf9wPx311jGt+M0lyrl+pkpjf3bx9X8ia/Ykn9YD4D0uvAUAoDGI6iAZnr/2TxlNVN4yKUCBAFiJJh8AgMbAmSCNnC8ZIZIkU1Je70d1/JK4wOYCYAmm6wAA0Fj0GCd1GCwd3qtD6/6Li24BG6PJBwCgMXEmSM4En19sC6BhYboOAAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANgMTT4AAABgMzT5AAAAgM3Q5AMAAAA2Q5MPAAAA2AxNPgAAAGAzNPkAAACAzdDkAwAAADZjmKZpBjqEP7lcLjmdTpWVlcnhcPi8fkxCopJSUn1ed3d+njqmdfd53WCtTWb/1Cazf2oHY2Yra5NZOlFlqvJklcJDQ1Tw1dagyPxTh/buUkFBgU9rWv37G7Cb0EAHsJuklFRlLFzi87qZU8ZbUjdYa5PZP7XJ7J/awZjZytqNPfNnu0r17sZCmaZ03JCif3yywWc+U9bsSZbUBVB/TNcBAKCBOFJe6W7wJck0pW8PH9eR8srABgMQdGjyAQBoIA65KlTXJNpSV4X/wwAIajT5AAA0EJc5ImQYtZdHOyL8HwZAUKPJBwCggYhsHq6x/du6G33DkFpHNVNk8/DABgMQdLjwFgCABqRvp2h1TnCo1FWhaEeEFv6DBh+A9wJ6Jj8zM1O9e/dWy5YtFRMTo9GjR2vnzp3nXW/Dhg3q2bOnmjZtqvbt2+ull17yQ1oAAPwjsnm4Lo9vyRl8ABcsoE3+hg0bNHXqVH322Wdas2aNTp48qSFDhqi8vPys6+zbt08jRozQwIEDlZeXpzlz5mj69OnKzs72Y3IAAACg4QrodJ0PP/zQ4/lrr72mmJgYbdmyRddcc02d67z00ktq27at5s2bJ0lKTU1Vbm6unnrqKY0ZM8bqyAAAAECD16AuvC0rK5MkRUVFnXXMpk2bNGTIEI9lQ4cOVW5urk6cOFFrfEVFhVwul8cDAAAAsLMG0+Sbpqn09HRdffXVuvLKK886rri4WLGxsR7LYmNjdfLkSZWWltYan5mZKafT6X4kJib6PDsAAADQkDSYJv/uu+9Wfn6+3nrrrfOONc64ibD5708OOXO5JGVkZKisrMz9OHDggG8CAwAAAA1Ug7iF5rRp07RixQrl5OSoTZs25xwbFxen4uJij2UlJSUKDQ3VpZdeWmt8RESEIiL4EBEAAAA0HgE9k2+apu6++24tX75ca9euVXJy8nnX6devn9asWeOxbPXq1erVq5fCwsKsigoAAAAEjYA2+VOnTtXrr7+uN998Uy1btlRxcbGKi4t1/Phx95iMjAyNGzfO/Xzy5Mnav3+/0tPTtWPHDr366qtavHixZs2aFYhNAAAAABqcgDb5CxcuVFlZmQYNGqT4+Hj3Y9myZe4xRUVFKiwsdD9PTk7WqlWrtH79el111VV67LHHtGDBAm6fCQAAAPybYdZctdpIuFwuOZ1OlZWVyeFw+Lx+TEKiklJSfV73wPZ89enTx+d1JSln4yZ1TOtuSW2rcpPZE5lP27x5c9AdK8F6fO/Ozwu698Oqn49gPFaszHxo7y4VFBT4tKbVv78Bu2kQF97aSVJKqjIWLvF53azZk7RixQqf15Wk3tcPsSSzZF1uMnsi82mjRo0KumMlWI/vzCnjg+79sOrnIxiPFaszAwisBnMLTQAIiLKD0r6c6q8AANgEZ/IBNF5fLJVW3iOZpySjiTRyvtRj3PnXAwCggeNMPoDGqezg6QZfqv66cgZn9AEAtkCTD6BxOrzndINfw6ySDu8NTB4AAHyIJh9A4xTVoXqKzk8ZIVJU+8DkAQDAh2jyATROzoTqOfhGSPVzI0QaOa96OQAAQY4LbwE0Xj3GSR0GV0/RiWpPgw8AsA2vzuQfP35cn376qbZv317rtR9//FFLly71WTAA8AtngpQ8kAYfAGAr9W7yd+3apdTUVF1zzTXq2rWrBg0apKKiIvfrZWVl+t3vfmdJSAAAAAD1V+8m/7777lPXrl1VUlKinTt3yuFwaMCAASosLLQyHwAAAAAv1bvJ37hxo5544glFR0fr8ssv14oVKzR8+HANHDhQe/dyyzkAAACgoaj3hbfHjx9XaKjn8BdeeEFNmjTRtddeqzfffNPn4QAAAAB4r95NfufOnZWbm6vU1FSP5c8995xM09SoUaN8Hg6BdaS8UodcFbrMEaHI5uGBjgMAAIB6qneT/6tf/UpvvfWWbr/99lqvPf/88zp16pReeukln4ZD4Hy2q1TvbiyUaUqGIY3t31Z9O0UHOhYAAADqwTBN0wx0CH9yuVxyOp0qKyuTw+Hwef2YhEQlpaSef6CXdufnqWNad5/Xrav2iSpTu7511RrXqbVDYSHGRdX2FX++Hw29rpW1yeyf2sGY2craZPZPbSszH9q7SwUFBT6tafXvb8Bu+DAsH0tKSVXGwiU+r5s5ZbwldeuqvbvoqBZ+uLvWuFuHddTl8S0vqrav+PP9aOh1raxNZv/UDsbMVtYms39qW5k5a/YkS+oCqD+vPgwLjcNljggZZ5ywNwwp2hERmEAAAADwCmfyUUtk83CN7d+21px8Lr4FAECqqqrSiRMnAh0DjUxYWJhCQkLqPZ4mH3Xq2ylanRMcKnVVKJq76wAAINM0VVxcrCNHjgQ6ChqpyMhIxcXFyThzykUdaPJxVpHNw2nuAQD4t5oGPyYmRpdcckm9Gi3AF0zT1LFjx1RSUiJJio+PP+86F9Tk79q1S+vXr1dJSYlOnTrl8dqDDz54ISUBAAAarKqqKneDf+mllwY6DhqhZs2aSZJKSkoUExNz3qk7Xjf5L7/8sqZMmaLo6Oha/1xgGAZNPgAAsJ2aOfiXXHJJgJOgMav5+Ttx4oTvm/zHH39cf/zjH3XfffddWDoAAIAgxRQdBJI3P39e30Lzhx9+0NixY71dDQAAAICfeN3kjx07VqtXr7YiCwAAABoAwzD0/vvvn3PMhAkTNHr06HrXLCgokGEY2rp160VlQ/14PV3n8ssv1wMPPKDPPvtMXbt2VVhYmMfr06dP91k4AAAAXJwJEyboyJEj523af6qoqEitWrWSVN2cJycnKy8vT1dddZV7zPz582Wapk+zDho0SFdddZXmzZtX5+uGYei9996r9cfFmdtYUlKiBx54QB988IG+++47tWrVSt26ddPDDz+sfv36+TRzQ+V1k79o0SK1aNFCGzZs0IYNGzxeMwyDJh8AACDIxcXFnXeM0+n0Q5ILM2bMGJ04cUJLlixR+/bt9d133+njjz/W4cOHAx3Nb7xu8vft22dFDgAAgEajqOy49pWWKzm6ueKdzfz6vQcNGqS0tDQ1bdpUr7zyisLDwzV58mQ9/PDD7jE/PWOenJwsSerevbsk6dprr9X69etrnT3/8MMP9fjjj+vLL79USEiI+vXrp/nz56tDhw5+3b4jR47o008/1fr163XttddKkpKSktSnTx+/5gg0r+fk/5Rpmj7/ZxoAAAA7W/Z5oQb8aa1uffkfGvCntVr2eaHfMyxZskTNmzfXP/7xD/3P//yPHn30Ua1Zs6bOsZs3b5YkffTRRyoqKtLy5cvrHFdeXq709HR9/vnn+vjjj9WkSRP96le/qvWZSlZr0aKFWrRooffff18VFRV+/d4NyQU1+UuXLlXXrl3VrFkzNWvWTGlpafrLX/7i62wAAAC2UlR2XBnLt+nUv8+RnjKlOcu/VFHZcb/mSEtL00MPPaSOHTtq3Lhx6tWrlz7++OM6x1522WWSpEsvvVRxcXGKioqqc9yYMWN00003qWPHjrrqqqu0ePFibdu2Tdu3b/dp9t/85jfuRr7m8cYbb7hfDw0NVVZWlpYsWaLIyEgNGDBAc+bMUX5+vk9zNHReN/nPPPOMpkyZohEjRuidd97RsmXLNGzYME2ePFnPPvusFRkBAABsYV9pubvBr1FlmiooPebXHGlpaR7P4+PjVVJSclE19+zZo1tvvVXt27eXw+FwT/MpLPTtv1Q8++yz2rp1q8dj1KhRHmPGjBmjb7/9VitWrNDQoUO1fv169ejRQ1lZWT7N0pB5PSf/ueee08KFCzVu3Dj3shtvvFFdunTRww8/rJkzZ/o0YLDZv3OHMqeM93ndA9vzlTV7ks/rStLu/DxLMkvW5SazJzJbX1ey7v0IxsxW1rby/eBYOc3KzIf27rKkrh0kRzdXE0MejX6IYahdtH8/SffMuyMahnHR02pGjhypxMREvfzyy2rdurVOnTqlK6+8UpWVlRdV90xxcXG6/PLLPZa1bNlSR44c8VjWtGlT3XDDDbrhhhv04IMP6s4779RDDz2kCRMm+DRPQ+V1k19UVKT+/fvXWt6/f38VFRX5JFQwS0pJVcbCJT6vmzV7klasWOHzupLU+/ohlmSWrMtNZk9kPm3UqFFBd6wE6/GdOWV80L0fVv18BOOxYnVm1C3e2UyZN3XVnOVfqso0FWIYeuKmK/1+8a03wsPDJUlVVVVnHfP9999rx44d+vOf/6yBAwdKkj799FO/5KuvK664wqvbiAa7C7pP/jvvvKM5c+Z4LF+2bJk6duzos2AAAAB2dHPvtrqm02UqKD2mdtGXNOgGX5JiYmLUrFkzffjhh2rTpo2aNm1a6/aZrVq10qWXXqpFixYpPj5ehYWFuv/++y/o+x06dKjWB2bFxcXV67aeUvUfHGPHjtUdd9yhtLQ0tWzZUrm5ufqf//kf3XjjjReUKRh53eQ/8sgjuvnmm5WTk6MBAwbIMAx9+umn+vjjj/XOO+9YkREAAMBW4p3NGnxzXyM0NFQLFizQo48+qgcffFADBw7U+vXrPcY0adJEb7/9tqZPn64rr7xSKSkpWrBggQYNGuT193vzzTf15ptveix76KGHPG7xeS4tWrTQz372Mz377LPas2ePTpw4ocTERN111121TlLbmddN/pgxY/SPf/xDzz77rN5//32ZpqkrrrhCmzdvdt8/FQAAAA3DmRebntmgS6o1jeXMW6TfeeeduvPOO89Z9+c//3mtO+n8tE67du3Oe+v1urKdK1ddWSIiIpSZmanMzMxz1rI7r5t8SerZs6def/11X2cBAAAA4AP1avJdLpccDof7v8+lZhzgoeygdHiPFNVBciYEOg0AAICt1avJb9WqlYqKihQTE6PIyEgZhlFrjGmaMgzjnFdeo5H6Yqm08h7JPCUZTaSR86Ue486/HgAAAC5IvZr8tWvXuj/dbN26dZYGgs2c+PF0gy9Vf105Q+owmDP6AAAAFqlXk3/ttdfW+d/AeZ0oP93g1zCrpMN7afIBAAAs0sTbFT788EOPDzd44YUXdNVVV+nWW2/VDz/84NNwsIGw5tVTdH7KCJGi2gcmDwAAQCPgdZM/e/Zs98W327ZtU3p6ukaMGKG9e/cqPT3d5wER5MKaVs/BN0Kqnxsh0sh5nMUHAACwkNe30Ny3b5+uuOIKSVJ2drZGjhypJ554Ql988YVGjBjh84CwgR7jqufgH95bfQafBh8AAMBSXp/JDw8P17FjxyRJH330kYYMGSJJioqKOu/tNc+Uk5OjkSNHqnXr1jIMo9YHMZxp/fr1Mgyj1uPrr7/2djPgb84EKXkgDT4AAIAfeN3kX3311UpPT9djjz2mzZs36xe/+IUkadeuXWrTpo1XtcrLy9WtWzc9//zzXq23c+dOFRUVuR8dO3b0an0AAADgbNq1a6d58+YFOsZF8brJf/755xUaGqq//vWvWrhwoRISqs/MfvDBBxo2bJhXtYYPH67HH39cN910k1frxcTEKC4uzv0ICQnxan0AAIDGZuPGjQoJCfG6XwsW/mzMH374YV111VVnfX3QoEGaMWNGreVZWVmKjIx0P6+qqlJmZqY6d+6sZs2aKSoqSn379tVrr7120Rm9npPftm1b/e///m+t5c8+++xFh6mv7t2768cff9QVV1yh//qv/9J111131rEVFRWqqKhwP/d2ShEAAIAdvPrqq5o2bZpeeeUVFRYWqm3btoGO1Og9/PDDWrRokZ5//nn16tVLLpdLubm5PrljpddNviSdOnVK33zzjUpKSnTqlOc90K+55pqLDnU28fHxWrRokXr27KmKigr95S9/0eDBg7V+/fqzft/MzEw98sgjlmU60/6dO5Q5ZbzP6+7Oz1Pv64f4vG5NbSsy19S2IjeZa9cms7V1a2pzfFtf2+r3g2PldF2rMh/au8uSurZSdlA6vEeK6uCX69XKy8v1zjvv6PPPP1dxcbGysrL04IMPul/PysrSjBkzdOTIEfey999/X7/61a9kmqZ72eOPP64FCxbo+PHjuvnmmxUdHa0PP/xQW7dulSRNmDBBR44cUZ8+fTR//nxVVFRo5syZmjt3rjIyMrR48WJdcsklevTRR3XHHXe46x48eFDp6elavXq1mjRpoquvvlrz589Xu3btPOpeffXVevrpp1VZWalbbrlF8+bNU1hYmAYNGqT9+/dr5syZmjlzpiS5c2/cuFH333+/Pv/8c0VHR+tXv/qVMjMz1bx5c0lSSUmJJk6cqI8++khxcXF6/PHHrdgFdVq5cqV+//vfa+zYse5l3bp180ltr5v8zz77TLfeeqv279/vsdMlyTAMVVVV+SRYXVJSUpSSkuJ+3q9fPx04cEBPPfXUWZv8jIwMj1t7ulwuJSYmWpYxKSVVGQuX+Lxu5pTxltQN1tpk9k9tMvundjBmtrI2mf1T28rMWbMnWVLXNr5YevrT4I0m1bea7jHO0m+5bNkydx/129/+VtOmTdMDDzwgwzDqXeONN97QH//4R7344osaMGCA3n77bT399NNKTk72GLd27Vq1adNGOTk5+r//+z9NnDhRmzZt0jXXXKN//OMfWrZsmSZPnqwbbrhBiYmJOnbsmK677joNHDhQOTk5Cg0N1eOPP65hw4YpPz9f4eHhkqR169YpPj5e69at0zfffKObb75ZV111le666y4tX75c3bp103/+53/qrrvucmfZtm2bhg4dqscee0yLFy/WoUOHdPfdd+vuu+92T4mZMGGCDhw4oLVr1yo8PFzTp09XSUmJD97184uLi9PatWv1+9//XpdddplPa3s9J3/y5Mnq1auXvvzySx0+fFg//PCD+3H48GGfhquPvn37avfu3Wd9PSIiQg6Hw+MBAAAQEGUHTzf4UvXXlTOql1to8eLF+u1vfytJGjZsmP71r3/p448/9qrGc889p4kTJ+p3v/udOnXqpAcffFBdu3atNS4qKkoLFixQSkqK7rjjDqWkpOjYsWOaM2eOOnbsqIyMDIWHh+v//u//JElvv/22mjRpoldeeUVdu3ZVamqqXnvtNRUWFmr9+vXuuq1atdLzzz+vzp0765e//KV+8YtfuLchKipKISEhatmypfuaTUl68skndeutt2rGjBnq2LGj+vfvrwULFmjp0qX68ccftWvXLn3wwQd65ZVX1K9fP/Xs2VOLFy/W8ePHL+Rt9vDiiy+qRYsWHo/Jkyd7jHnmmWd06NAhxcXFKS0tTZMnT9YHH3xw0d9buoAmf/fu3XriiSeUmpqqyMhIOZ1Oj4e/5eXlKT4+3u/fFwAAwGuH95xu8GuYVdWfJWORnTt3avPmzbrlllskSaGhobr55pv16quvel2nT58+HsvOfC5JXbp0UZMmp1vM2NhYjz8GQkJCdOmll7rPlm/ZskXffPONWrZs6W6Go6Ki9OOPP2rPnj0edX96s5X4+PjznnHfsmWLsrKyPBrtoUOH6tSpU9q3b5927Nih0NBQ9erVy71O586dPS6OvVC33Xabtm7d6vF49NFHPcZcccUV+vLLL/XZZ5/pd7/7nb777juNHDlSd95550V/f6+n6/zsZz/TN998o8svv/yiv/m//vUvffPNN+7n+/bt09atWxUVFaW2bdsqIyNDBw8e1NKlSyVJ8+bNU7t27dSlSxdVVlbq9ddfV3Z2trKzsy86CwAAgOWiOlRP0flpo2+EVH9YpEUWL16skydPuu+IKFXPVw8LC9MPP/ygVq1aqUmTJrWmYZ84caJWrTOn95y5jiSFhYXVWqeuZTXXdZ46dUo9e/bUG2+8UavWT6ewnKvG2Zw6dUqTJk3S9OnTa73Wtm1b7dy5013L15xOZ61+OSYmpta4Jk2aqHfv3urdu7dmzpyp119/Xbfffrvmzp1bayqUN7xu8qdNm6Y//OEPKi4uVteuXWu94WlpafWulZub63FnnJq58+PHj1dWVpaKiopUWFjofr2yslKzZs3SwYMH1axZM3Xp0kV///vf+aRdAAAQHJwJ1XPwV86oPoNvhEgj51l28e3Jkye1dOlSPf300+4PMK0xZswYvfHGG7r77rt12WWX6ejRoyovL3dfkFpzMW2NlJQUbd68Wbfffrt7WW5u7kVn7NGjh5YtW6aYmJiLmlYdHh5e69rQHj166KuvvjrryenU1FSdPHlSubm57n+V2Llzp8cFyP52xRVXSKq+WPpieN3kjxkzRpI8rog2DEOmaXp94e2gQYPq/AuwRlZWlsfze++9V/fee693gQEAABqSHuOkDoOrp+hEtbf07jr/+7//qx9++EETJ06sNa3617/+tRYvXqy7775bP/vZz3TJJZdozpw5mjZtmjZv3lyrD5s2bZruuusu9erVS/3799eyZcuUn5+v9u0v7l8hbrvtNj355JO68cYb9eijj6pNmzYqLCzU8uXLNXv27Hp/2Gq7du2Uk5OjW265RREREYqOjtZ9992nvn37aurUqbrrrrvUvHlz7dixQ2vWrNFzzz2nlJQUDRs2THfddZcWLVqk0NBQzZgxQ82aNTvv9zt+/HitP4RatGjh1WyXX//61xowYID69++vuLg47du3TxkZGerUqZM6d+5c7zp18XpO/r59+2o99u7d6/4KAACA83AmSMkDLb995uLFi/Xzn/+8zusmx4wZo61bt+qLL75QVFSUXn/9da1atUpdu3bVW2+9pYcffthj/G233aaMjAzNmjVLPXr00L59+zRhwgQ1bdr0ojJecsklysnJUdu2bXXTTTcpNTVVd9xxh44fP+7Vmf1HH31UBQUF6tChg3uaT1pamjZs2KDdu3dr4MCB6t69ux544AGP6zlfe+01JSYm6tprr9VNN92k//zP/6xzWs2Zdu3ape7du3s8vJ1LP3ToUK1cuVIjR45Up06dNH78eHXu3FmrV69WaOgF3enezeu1k5KSLuobAgAAwD9Wrlx51td69OjhMaNi9OjRGj16tMeYn96OUpIeeOABPfDAA+7nN9xwg8eZ6zPP/kvyuENOjYKCAo/ncXFxWrLk7Ld0ravumZ9u27dvX/3zn/+sNa53795avXr1WWvHxcXV+qDXn05JqsvDDz9c64+gn6prm6Xq23VOmDDB/fyuu+6q9R77itdn8iXpL3/5iwYMGKDWrVtr//79kqrf6L/97W8+DQcAAICG4dixY3rmmWf01Vdf6euvv9ZDDz2kjz76SOPHW/Oharg4Xjf5CxcuVHp6ukaMGKEjR4645+BHRkbW+osKAAAA9mAYhlatWqWBAweqZ8+eWrlypbKzs/Xzn/880NFQB6+n6zz33HN6+eWXNXr0aP3pT39yL+/Vq5dmzZrl03AAAABoGJo1a6aPPvoo0DFQTxd04W337t1rLY+IiLjoW/3gtCPlldpddFRHyisDHQUAAABBxusz+cnJydq6dWutC3A/+OAD9309cXE+21WqdzcWyjQlw5DG9m8b6EgAAEB1f/gT4C/e/Px53eTPnj1bU6dO1Y8//ijTNLV582a99dZbyszM1CuvvOJtOZzhSHmlu8GXJNOU3t1YKGcV/1MBACBQaj7889ixY/W6hzpghWPHjkmq/em/dfG6yf/d736nkydP6t5779WxY8d06623KiEhQfPnz9ctt9zifVp4OOSq0Jl/pJmmVHmy/h8yBgAAfCskJESRkZEqKSmRVH1vd8MwApwKjYVpmjp27JhKSkoUGRmpkJCQ865zQXfZr7mnZ2lpqU6dOlWvDwxA/VzmiJBhyKPRNwwpPPT8OxMAAFgnLi5OktyNPuBvkZGR7p/D8zHMRja5zOVyyel0qqyszKtPUauvmIREJaWkXlSNH8or9e3h4+7nraOa6V/7v1afPn0uNl6dcjZuUse02hdT+8KB7fmW5CazJzKftnnz5qA7Vqx6LyRr9+Hu/Lygez+s+vkIxmPFysyH9u6q9WFHF8vq39/eqKqq0okTJwKaAY1PWFhYvc7g1/D6TP7333+vBx98UOvWrVNJSYlOnTrl8frhw4e9LWkrSSmpylh49k9sq68j5ZUqdVUo2hGhyObhypo9SStWrPBBwtp6Xz/EJ5nrYlVuMnsi82mjRo0KumMlWI/vzCnjG+b7UXZQOrxHiuogORM8XrLq5yMYjxWrM9tZSEiIV80WEAheN/m//e1vtWfPHk2cOFGxsbHMR7NIZPNwRTYPD3QMAAguXyyVVt4jmacko4k0cr7UY1ygUwGA33nd5H/66af69NNP1a1bNyvyAABwYcoOnm7wpeqvK2dIHQbXOqMPAHbn9Ydhde7cWcePHz//QAAA/OnwntMNfg2zSjq8NzB5ACCAvG7yX3zxRc2dO1cbNmzQ999/L5fL5fEAACAgojpUT9H5KSNEimofmDwAEEBeT9eJjIxUWVmZrr/+eo/lpmnKMAxVVXE/dwBAADgTqufgr5xRfQbfCJFGzmOqDoBGyesm/7bbblN4eLjefPNNLrwFADQsPcZVz8E/vLf6DD4NPoBGyusm/8svv1ReXp5SUlKsyAMAwMVxJtDcA2j0vJ6T36tXLx04cMCKLAAAAAB8wOsz+dOmTdM999yj2bNnq2vXrgoLC/N4PS0tzWfhAAAAAHjP6yb/5ptvliTdcccd7mWGYXDhLQAAANBAeN3k79u3z4ocAAAAAHzE6yY/KSnJihwAAAAAfKReTf6KFSs0fPhwhYWFacWKFeccO2rUKJ8EAwAAAHBh6tXkjx49WsXFxYqJidHo0aPPOo45+QAAAEDg1avJP3XqVJ3/DQAAAKDhMUzTNAMdwp9cLpecTqfKysrkcDh8Xj8mIVFJKak+r7s7P08d07r7vG6w1iazf2qT2T+1gzGzlbXJ7J/aVmY+tHeXCgoKfFrT6t/fgN14deHtqVOnlJWVpeXLl6ugoECGYSg5OVm//vWvdfvtt8swDKtyBo2klFRlLFzi87qZU8ZbUjdYa5PZP7XJ7J/aDTnzkfJKHXJV6DJHhCKbh/u09tk05PfD33WtrG1l5qzZkyypC6D+6t3km6apUaNGadWqVerWrZu6du0q0zS1Y8cOTZgwQcuXL9f7779vYVQAgD99tqtU724slGlKhiGN7d9WfTtFBzoWAKAe6t3kZ2VlKScnRx9//LGuu+46j9fWrl2r0aNHa+nSpRo3bpzPQwIA/OtIeaW7wZck05Te3ViozgmOWmf0AQANT5P6Dnzrrbc0Z86cWg2+JF1//fW6//779cYbb/g0HAAgMA65KnTmFVumKZW6KgITCADglXo3+fn5+Ro2bNhZXx8+fLj++c9/+iQUACCwLnNE6MzLrAxDinZEBCYQAMAr9W7yDx8+rNjY2LO+Hhsbqx9++MEnoQAAgRXZPFxj+7d1N/o1c/KZqgMAwaHec/KrqqoUGnr24SEhITp58qRPQgEAAq9vp2h1TnCo1FWh6DrurgMAaLi8urvOhAkTFBFR9z/VVlQwTxMA7CayeTjNPQAEoXo3+ePHjz/vGO6sAwAAAARevZv81157zcocAAAAAHyk3hfeAgAAAAgONPkAAACAzdDkAwAAADZDkw8AAADYDE0+AAAAYDOGaZpmoEP4k8vlktPpVFlZmRwOh8/rxyQkKikl1ed1D2zPV58+fXxeV5JyNm5Sx7TultS2KjeZPZH5tM2bNwfdsRKsx/fu/Lygez989vNx4kfpRLkU1lwKaxqUx4qVmQ/t3aWCggKf1rT69zdgN/W+hSbqJyklVRkLl/i8btbsSVqxYoXP60pS7+uHWJJZsi43mT2R+bRRo0YF3bESrMd35pTxQfd++OTn44ul0sp7JPOUZDSRRs5X71k/Bt2xYvXxDSCwAjpdJycnRyNHjlTr1q1lGIbef//9866zYcMG9ezZU02bNlX79u310ksvWR8UAABJKjt4usGXqr+unKGQKj71HUDDEtAmv7y8XN26ddPzzz9fr/H79u3TiBEjNHDgQOXl5WnOnDmaPn26srOzLU4KAICkw3tON/g1zCqFnDwemDwAcBYBna4zfPhwDR8+vN7jX3rpJbVt21bz5s2TJKWmpio3N1dPPfWUxowZY1FKAAD+LapD9RSdnzb6RoiqQpsFLhMA1CGo7q6zadMmDRkyxGPZ0KFDlZubqxMnTtS5TkVFhVwul8cDAIAL4kyQRs6XjJDq50aINHKeqkIiApsLAM4QVBfeFhcXKzY21mNZbGysTp48qdLSUsXHx9daJzMzU4888oi/IgIA7K7HOKnDYOnwXimqfXXjr9cDnQoAPATVmXxJMgzD43nNHUDPXF4jIyNDZWVl7seBAwcszwgAsDlngpQ88N8NPgA0PEF1Jj8uLk7FxcUey0pKShQaGqpLL720znUiIiIUEcE/owIAAKDxCKoz+f369dOaNWs8lq1evVq9evVSWFhYgFIBAAAADUtAm/x//etf2rp1q7Zu3Sqp+haZW7duVWFhoaTqqTbjxo1zj588ebL279+v9PR07dixQ6+++qoWL16sWbNmBSI+AAAA0CAFdLpObm6urrvuOvfz9PR0SdL48eOVlZWloqIid8MvScnJyVq1apVmzpypF154Qa1bt9aCBQu4fSYAAADwEwFt8gcNGuS+cLYuWVlZtZZde+21+uKLLyxMBQAAAAS3oJqTDwAAAOD8aPIBAAAAm6HJBwAAAGyGJh8AAACwGcM815WvNuRyueR0OlVWViaHw+Hz+jEJiUpKSfV53d35eeqY1t3ndYO1Npn9U5vM/qkdjJmtrE1m/9S2MvOhvbtUUFDg05pW//4G7CaoPvE2GCSlpCpj4RKf182cMt6SusFam8z+qU1m/9QOxsxW1iazf2pbmTlr9iRL6gKoP6brAAAAADZDkw8AAADYDE0+AAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANgMTT4AAABgMzT5AAAAgM3Q5AMAAAA2Q5MPAAAA2AxNPgAAAGAzNPkALHekvFK7i47qSHlloKMAANAohAY6AAB7+2xXqd7dWCjTlAxDGtu/rfp2ig50LAAAbI0z+QAsc6S80t3gS5JpSu9uLOSMPgAAFqPJB2CZQ64Kd4NfwzSlUldFYAIBANBIGKZ55q9ge3O5XHI6nSorK5PD4fB5/ZiERCWlpPq87oHt+erTp4/P60pSzsZN6pjW3ZLaVuUms6eGmvlElald37pqLe/U2qHindssybx58+agO1aC9fjenZ8XdO+HVT8fjfH4PpdDe3epoKDApzWt/v0N2A1z8n0sKSVVGQuX+Lxu1uxJWrFihc/rSlLv64dYklmyLjeZPTXkzGebk29V5lGjRgXdsRKsx3fmlPFB935Y9fPRWI/vs8maPcmSugDqjyYfgKX6dopW5wSHSl0VinZEKLJ5eKAjAQBgezT5ACwX2Tyc5h4AAD/iwlvgfMoOSvtyqr8CAAAEAc7kA+fyxVJp5T2SeUoymkgj50s9xgU6FQAAwDlxJh84m7KDpxt8qfrryhmc0QcAAA0eTT5wNof3nG7wa5hV0uG9gckDAABQTzT5wNlEdaieovNTRogU1T4weQAAAOqJJh84G2dC9Rx8I6T6uREijZxXvRwAAKAB48Jb4Fx6jJM6DK6eohPVngYfAAAEBZp84HycCTT3AAAgqDBdBwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABshiYfAAAAsBnDNE0z0CH8yeVyyel0qqysTA6Hw+f1YxISlZSS6vO6B7bnq0+fPj6vK0k5GzepY1p3S2pblZvMnsh82ubNm4PuWAnW43t3fl7QvR9W/XwE47FiZeZDe3epoKDApzWt/v0N2A230PSxpJRUZSxc4vO6WbMnacWKFT6vK0m9rx9iSWbJutxk9kTm00aNGhV0x0qwHt+ZU8YH3fth1c9HMB4rVmcGEFhM1wEAAABshiYfAAAAsBmafAAAAMBmaPIBAAAAm6HJBwAAAGwm4E3+iy++qOTkZDVt2lQ9e/bUJ598ctax69evl2EYtR5ff/21HxMDAAAADVtAm/xly5ZpxowZmjt3rvLy8jRw4EANHz5chYWF51xv586dKioqcj86duzop8QAAABAwxfQJv+ZZ57RxIkTdeeddyo1NVXz5s1TYmKiFi5ceM71YmJiFBcX536EhIT4KTEAAADQ8AWsya+srNSWLVs0ZMgQj+VDhgzRxo0bz7lu9+7dFR8fr8GDB2vdunXnHFtRUSGXy+XxAAAAAOwsYE1+aWmpqqqqFBsb67E8NjZWxcXFda4THx+vRYsWKTs7W8uXL1dKSooGDx6snJycs36fzMxMOZ1O9yMxMdGn2+FXZQelfTnVXwEAAICzCA10AMMwPJ6bpllrWY2UlBSlpKS4n/fr108HDhzQU089pWuuuabOdTIyMpSenu5+7nK5grPRLzsgzbtSMk9JRhNp5Hypx7hApwIAAEADFLAz+dHR0QoJCal11r6kpKTW2f1z6du3r3bv3n3W1yMiIuRwODwewabZsWKp+MvqBl+q/rpyBmf0AQAAUKeANfnh4eHq2bOn1qxZ47F8zZo16t+/f73r5OXlKT4+3tfxGpTmRwtqLzSrpMN7/Z4FAAAADV9Ap+ukp6fr9ttvV69evdSvXz8tWrRIhYWFmjx5sqTqqTYHDx7U0qVLJUnz5s1Tu3bt1KVLF1VWVur1119Xdna2srOzA7kZlitv2a72QiNEimrv9ywAAABo+ALa5N988836/vvv9eijj6qoqEhXXnmlVq1apaSkJElSUVGRxz3zKysrNWvWLB08eFDNmjVTly5d9Pe//10jRowI1Cb4xfFL4qS4KyVjf/UZfCNEGjlPciYEOhoAAAAaIMM0TTPQIfzJ5XLJ6XSqrKzMkvn5MQmJSkpJ9Xnd3fl56tzlCoWcPK6q0GaqConwae2Oad19Vs8ftcnsn9pk9k/tYMxsZW0y+6e2lZkP7d2lgoICn9a0+vc3YDcBv7uO3SSlpCpj4RKf182cMl73WlC3prYVma2sTWb/1Cazf2oHY2Yra5PZP7WtzJw1e5IldQHUX0A/8RYAAACA79HkAwAAADZDkw8AAADYDE0+AAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANgMTT4AAABgMzT5AACgliPlldpddFRHyisDHQXABeATbwEAgIfPdpXq3Y2FMk3JMKSx/duqb6foQMcC4AXO5AMAALcj5ZXuBl+STFN6d2MhZ/SBIEOTDwAA3A65KtwNfg3TlEpdFYEJBOCC0OQDAAC3yxwRMgzPZYYhRTsiAhMIwAWhyQcAAG6RzcM1tn9bd6NfMyc/snl4YIMB8Iphmmf+o5y9uVwuOZ1OlZWVyeFw+Lx+TEKiklJSfV53d36eOqZ193ndYK1NZv/UJrN/agdjZitrk9k/tc9X90SVqcqTVQoPDVFYiHHWcXU5tHeXCgoKLjKhJ6t/fwN2w911fCwpJVUZC5f4vG7mlPGW1A3W2mT2T20y+6d2MGa2sjaZ/VPbysxZsydZUhdA/TFdBwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABshiYfAAAAsBmafKCBOFJeqd1FR/noeAAAcNG4hSbQAHy2q1TvbiyUaZ7+4Jm+naIDHQsAAAQpzuQDAXakvNLd4EuSaUrvbizkjD4AALhgNPlAgB1yVejMz502TanUVRGYQAAAIOjR5AMBdpkjQsYZnxhvGFK0IyIwgQAAQNCjyQcCLLJ5uMb2b+tu9Gvm5Ec2Dw9sMAAAELS48BZoAPp2ilbnBIdKXRWKdkTQ4AMAgItCkw80EJHNw2nuAQCATzBdBwAAALAZmnwAAADAZgzTPPPmffbmcrnkdDpVVlYmh8Ph8/oxCYlKSkn1ed0D2/PVp08fn9eVpJyNm9Qxrbslta3KTWZPZD5t8+bNQXesBOvxvTs/L+jeD6t+PoLxWLEy86G9u1RQUODTmlb//gbshjn5PpaUkqqMhUt8Xjdr9iStWLHC53Ulqff1QyzJLFmXm8yeyHzaqFGjgu5YCdbjO3PK+KB7P6z6+QjGY8XqzAACi+k6AAAAgM3Q5AMAAAA2Q5MPAAAA2AxNPgAAAGAzNPkAAACAzdDkAwAAADZDkw8AAADYDE0+AAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANhMwJv8F198UcnJyWratKl69uypTz755JzjN2zYoJ49e6pp06Zq3769XnrpJT8lBYAGpuygtC+n+isAAD8R0CZ/2bJlmjFjhubOnau8vDwNHDhQw4cPV2FhYZ3j9+3bpxEjRmjgwIHKy8vTnDlzNH36dGVnZ/s5OQAE2BdLpXlXSktGVn/9YmmgEwEAGpCANvnPPPOMJk6cqDvvvFOpqamaN2+eEhMTtXDhwjrHv/TSS2rbtq3mzZun1NRU3Xnnnbrjjjv01FNP+Tk5AARQ2UFp5T2Sear6uXlKWjmDM/oAALeANfmVlZXasmWLhgwZ4rF8yJAh2rhxY53rbNq0qdb4oUOHKjc3VydOnKhznYqKCrlcLo8HAAS1w3tON/g1zCrp8N7A5AEANDihgfrGpaWlqqqqUmxsrMfy2NhYFRcX17lOcXFxneNPnjyp0tJSxcfH11onMzNTjzzyiO+Cn8ex7w4qa/Ykn9fNz8/XqFGjfF5Xkg7t3WVJZsm63GT2RGbr60rWvR9eZz7xo7T32BkLDekfj0lhTT2WWrkPG8z70QBqB+OxYmXmPXv2WFIXQP0FrMmvYRiGx3PTNGstO9/4upbXyMjIUHp6uvu5y+VSYmLihcY9r6+++sqy2gDg9sXS6ik6ZpVkhEgj50k9xgU6FQCggQhYkx8dHa2QkJBaZ+1LSkpqna2vERcXV+f40NBQXXrppXWuExERoYiICN+EBoCGosc4qcPg6ik6Ue0lZ0KgEwEAGpCAzckPDw9Xz549tWbNGo/la9asUf/+/etcp1+/frXGr169Wr169VJYWJhlWQGgQXImSMkDafABALUE9O466enpeuWVV/Tqq69qx44dmjlzpgoLCzV58mRJ1VNtxo07/c/PkydP1v79+5Wenq4dO3bo1Vdf1eLFizVr1qxAbQIAAADQ4AR0Tv7NN9+s77//Xo8++qiKiop05ZVXatWqVUpKSpIkFRUVedwzPzk5WatWrdLMmTP1wgsvqHXr1lqwYIHGjBkTqE0AAAAAGhzDrLlytZFwuVxyOp0qKyuTw+EIdBwAAFAP/P4GvBPQ6ToAAAAAfI8mHwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABshiYfAAAAsBmafAAAAMBmaPIBAAAAm6HJBwAAAGwmNNAB/K3mA35dLleAkwAAgPqq+b1d83scwLk1uib/6NGjkqTExMQAJwEAAN46evSonE5noGMADZ5hNrI/iU+dOqVvv/1WLVu2lGEYPq3tcrmUmJioAwcOyOFw+LR2Q2D37ZPsv41sX/Cz+zayfcHPqm00TVNHjx5V69at1aQJs42B82l0Z/KbNGmiNm3aWPo9HA6Hbf/nLdl/+yT7byPbF/zsvo1sX/CzYhs5gw/UH38KAwAAADZDkw8AAADYDE2+D0VEROihhx5SREREoKNYwu7bJ9l/G9m+4Gf3bWT7gl9j2EYgGDS6C28BAAAAu+NMPgAAAGAzNPkAAACAzdDkAwAAADZDkw8AAADYDE2+l1588UUlJyeradOm6tmzpz755JNzjt+wYYN69uyppk2bqn379nrppZf8lPTCeLN969evl2EYtR5ff/21HxPXX05OjkaOHKnWrVvLMAy9//77510n2Paft9sYTPswMzNTvXv3VsuWLRUTE6PRo0dr586d510vmPbhhWxjMO3DhQsXKi0tzf0hSf369dMHH3xwznWCaf95u33BtO/qkpmZKcMwNGPGjHOOC6Z9CNgJTb4Xli1bphkzZmju3LnKy8vTwIEDNXz4cBUWFtY5ft++fRoxYoQGDhyovLw8zZkzR9OnT1d2drafk9ePt9tXY+fOnSoqKnI/Onbs6KfE3ikvL1e3bt30/PPP12t8sO0/yfttrBEM+3DDhg2aOnWqPvvsM61Zs0YnT57UkCFDVF5eftZ1gm0fXsg21giGfdimTRv96U9/Um5urnJzc3X99dfrxhtv1FdffVXn+GDbf95uX41g2Hdn+vzzz7Vo0SKlpaWdc1yw7UPAVkzUW58+fczJkyd7LOvcubN5//331zn+3nvvNTt37uyxbNKkSWbfvn0ty3gxvN2+devWmZLMH374wQ/pfEuS+d57751zTLDtvzPVZxuDeR+WlJSYkswNGzacdUyw78P6bGMw70PTNM1WrVqZr7zySp2vBfv+M81zb1+w7rujR4+aHTt2NNesWWNee+215j333HPWsXbYh0Cw4kx+PVVWVmrLli0aMmSIx/IhQ4Zo48aNda6zadOmWuOHDh2q3NxcnThxwrKsF+JCtq9G9+7dFR8fr8GDB2vdunVWxvSrYNp/FysY92FZWZkkKSoq6qxjgn0f1mcbawTbPqyqqtLbb7+t8vJy9evXr84xwbz/6rN9NYJt302dOlW/+MUv9POf//y8Y4N5HwLBjia/nkpLS1VVVaXY2FiP5bGxsSouLq5zneLi4jrHnzx5UqWlpZZlvRAXsn3x8fFatGiRsrOztXz5cqWkpGjw4MHKycnxR2TLBdP+u1DBug9N01R6erquvvpqXXnllWcdF8z7sL7bGGz7cNu2bWrRooUiIiI0efJkvffee7riiivqHBuM+8+b7Qu2fSdJb7/9trZs2aLMzMx6jQ/GfQjYRWigAwQbwzA8npumWWvZ+cbXtbyh8Gb7UlJSlJKS4n7er18/HThwQE899ZSuueYaS3P6S7DtP28F6z68++67lZ+fr08//fS8Y4N1H9Z3G4NtH6akpGjr1q06cuSIsrOzNX78eG3YsOGsjXCw7T9vti/Y9t2BAwd0zz33aPXq1WratGm91wu2fQjYBWfy6yk6OlohISG1zmqXlJTUOktRIy4urs7xoaGhuvTSSy3LeiEuZPvq0rdvX+3evdvX8QIimPafLzX0fTht2jStWLFC69atU5s2bc45Nlj3oTfbWJeGvA/Dw8N1+eWXq1evXsrMzFS3bt00f/78OscG4/7zZvvq0pD33ZYtW1RSUqKePXsqNDRUoaGh2rBhgxYsWKDQ0FBVVVXVWicY9yFgFzT59RQeHq6ePXtqzZo1HsvXrFmj/v3717lOv379ao1fvXq1evXqpbCwMMuyXogL2b665OXlKT4+3tfxAiKY9p8vNdR9aJqm7r77bi1fvlxr165VcnLyedcJtn14IdtYl4a6D+timqYqKirqfC3Y9l9dzrV9dWnI+27w4MHatm2btm7d6n706tVLt912m7Zu3aqQkJBa69hhHwJBKyCX+wapt99+2wwLCzMXL15sbt++3ZwxY4bZvHlzs6CgwDRN07z//vvN22+/3T1+79695iWXXGLOnDnT3L59u7l48WIzLCzM/Otf/xqoTTgnb7fv2WefNd977z1z165d5pdffmnef//9piQzOzs7UJtwTkePHjXz8vLMvLw8U5L5zDPPmHl5eeb+/ftN0wz+/Wea3m9jMO3DKVOmmE6n01y/fr1ZVFTkfhw7dsw9Jtj34YVsYzDtw4yMDDMnJ8fct2+fmZ+fb86ZM8ds0qSJuXr1atM0g3//ebt9wbTvzubMu+sE+z4E7IQm30svvPCCmZSUZIaHh5s9evTwuLXd+PHjzWuvvdZj/Pr1683u3bub4eHhZrt27cyFCxf6ObF3vNm+//7v/zY7dOhgNm3a1GzVqpV59dVXm3//+98DkLp+am5Xd+Zj/PjxpmnaY/95u43BtA/r2i5J5muvveYeE+z78EK2MZj24R133OH+/8tll11mDh482N0Am2bw7z9vty+Y9t3ZnNnkB/s+BOzEMM1/XwEDAAAAwBaYkw8AAADYDE0+AAAAYDM0+QAAAIDN0OQDAAAANkOTDwAAANgMTT4AAABgMzT5AAAAgM3Q5AMAAAA2Q5MPwOcMw9D7778f6BjntH79ehmGoSNHjgQ6CgAAPkeTD6BeJkyYIMMwZBiGwsLCFBsbqxtuuEGvvvqqTp065TG2qKhIw4cPD1DS+unfv7+KiorkdDot/T45OTkaOXKkWrduHRR//AAA7IEmH0C9DRs2TEVFRSooKNAHH3yg6667Tvfcc49++ctf6uTJk+5xcXFxioiICGDS8wsPD1dcXJwMw7D0+5SXl6tbt256/vnnLf0+AAD8FE0+gHqLiIhQXFycEhIS1KNHD82ZM0d/+9vf9MEHHygrK8s97qdnrAsKCmQYht555x0NHDhQzZo1U+/evbVr1y59/vnn6tWrl1q0aKFhw4bp0KFDHt/vtddeU2pqqpo2barOnTvrxRdfdL9WU3f58uW67rrrdMkll6hbt27atGmTe8z+/fs1cuRItWrVSs2bN1eXLl20atUqSXVP18nOzlaXLl0UERGhdu3a6emnn/bI065dOz3xxBO644471LJlS7Vt21aLFi0653s2fPhwPf7447rpppu8easBALgoNPkALsr111+vbt26afny5ecc99BDD+m//uu/9MUXXyg0NFS/+c1vdO+992r+/Pn65JNPtGfPHj344IPu8S+//LLmzp2rP/7xj9qxY4eeeOIJPfDAA1qyZIlH3blz52rWrFnaunWrOnXqpN/85jfuf1WYOnWqKioqlJOTo23btum///u/1aJFizrzbdmyRf/xH/+hW265Rdu2bdPDDz+sBx54wOOPF0l6+umn1atXL+Xl5en3v/+9pkyZoq+//voC3jkAAKwTGugAAIJf586dlZ+ff84xs2bN0tChQyVJ99xzj37zm9/o448/1oABAyRJEydO9GioH3vsMT399NPuM+DJycnavn27/vznP2v8+PEedX/xi19Ikh555BF16dJF33zzjTp37qzCwkKNGTNGXbt2lSS1b9/+rPmeeeYZDR48WA888IAkqVOnTtq+fbuefPJJTZgwwT1uxIgR+v3vfy9Juu+++/Tss89q/fr16ty5c33eKgAA/IIz+QAummma553bnpaW5v7v2NhYSXI33zXLSkpKJEmHDh3SgQMHNHHiRLVo0cL9ePzxx7Vnz56z1o2Pj5ckd53p06fr8ccf14ABA/TQQw+d8w+RHTt2uP/gqDFgwADt3r1bVVVVdX4/wzAUFxfn/n4AADQUNPkALtqOHTuUnJx8zjFhYWHu/675g+DMZTV36an5+vLLL2vr1q3ux5dffqnPPvvsvHVr1r/zzju1d+9e3X777dq2bZt69eql5557rs58df2hYprmObfjzNwAADQUNPkALsratWu1bds2jRkzxmc1Y2NjlZCQoL179+ryyy/3eJzvj4kzJSYmavLkyVq+fLn+8Ic/6OWXX65z3BVXXKFPP/3UY9nGjRvVqVMnhYSEXPC2AAAQCMzJB1BvFRUVKi4uVlVVlb777jt9+OGHyszM1C9/+UuNGzfOp9/r4Ycf1vTp0+VwODR8+HBVVFQoNzdXP/zwg9LT0+tVY8aMGRo+fLg6deqkH374QWvXrlVqamqdY//whz+od+/eeuyxx3TzzTdr06ZNev755z3u6HMh/vWvf+mbb75xP9+3b5+2bt2qqKgotW3b9qJqAwBwNjT5AOrtww8/VHx8vEJDQ9WqVSt169ZNCxYs0Pjx49WkiW//YfDOO+/UJZdcoieffFL33nuvmjdvrq5du2rGjBn1rlFVVaWpU6fq//2//yeHw6Fhw4bp2WefrXNsjx499M477+jBBx/UY489pvj4eD366KMeF91eiNzcXF133XXu5zV/oIwfP77WnXsAAPAVw6xr0ikAAACAoMWcfAAAAMBmaPIBAAAAm6HJBwAAAGyGJh8AAACwGZp8AAAAwGZo8gEAAACbockHAAAAbIYmHwAAALAZmnwAAADAZmjyAQAAAJuhyQcAAABs5v8DGdxzKEGE2AIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "new = num #Points to be added\n",
    "\n",
    "new_num = new+num\n",
    "x_new = sampling.expand_lhs(x,  new, method='basic')\n",
    "\n",
    "\n",
    "### For visualization only\n",
    "intervals = []  \n",
    "subspace_bool = []\n",
    "for i in range(len(xlimits)):  \n",
    "    intervals.append(np.linspace(xlimits[i][0],xlimits[i][1],new_num+1))\n",
    "    subspace_bool.append([[intervals[i][j]<x[kk][i]<intervals[i][j+1] \n",
    "    for kk in range(len(x))] for j in range(len(intervals[i])-1)])\n",
    "    \n",
    "    \n",
    "fig, ax = plt.subplots(1)\n",
    "ax.plot(x[:, PLOT1], x[:, PLOT2], \".\")\n",
    "ax.plot(x_new[len(x):, PLOT1], x_new[len(x):, PLOT2], \".\")\n",
    "for i in range(len(intervals[0])):\n",
    "    ax.plot([intervals[PLOT1][i], intervals[PLOT1][i]], [intervals[PLOT2][0],intervals[PLOT2][-1]], linewidth=0.5,color='k')\n",
    "    ax.plot([intervals[PLOT1][0], intervals[PLOT1][-1]], [intervals[PLOT2][i],intervals[PLOT2][i]], linewidth=0.5,color='k')\n",
    "ax.set_xlabel(\"Dimension \" + str(PLOT1+1))\n",
    "ax.set_ylabel(\"Dimension \" + str(PLOT2+1))\n",
    "\n",
    "i=PLOT1\n",
    "for j in range(len(subspace_bool[i])):\n",
    "    if True in subspace_bool[i][j]:\n",
    "        ax.add_patch(patches.Rectangle((intervals[PLOT1][j], intervals[PLOT2][0]), \n",
    "        (intervals[PLOT1][1]-intervals[PLOT1][0]), (intervals[PLOT2][-1]-intervals[PLOT2][0]),facecolor='lightblue'))\n",
    "\n",
    "i=PLOT2\n",
    "for j in range(len(subspace_bool[i])):\n",
    "    if True in subspace_bool[i][j]:\n",
    "        ax.add_patch(patches.Rectangle((intervals[PLOT1][0], intervals[PLOT2][j]), \n",
    "        (intervals[PLOT1][-1]-intervals[PLOT1][0]), (intervals[PLOT2][1]-intervals[PLOT2][0]),facecolor='lightblue'))\n",
    "        \n",
    "ax.legend(['Initial LHS','Augmented LHS'],bbox_to_anchor = (1.05, 0.6))\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Augmented Sampling - Optimization ESE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we expand our LHS sampling (10 points) by adding another 10 points. In this case, we try to optimize the augmented LHS using the Enhanced Stochastic Evolutionary algorithm (ESE)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvkAAAGwCAYAAAAzA3pBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABO3klEQVR4nO3de3wU9b3/8feQkHDdDTGEhBDCRQgRCHItF1GUyq0FqdSfVitQ0SOIIKSgBitea3qON8ALFkUD9YYWtNCiB5RLtGAxGBoUBARCOJgYKJKlARNI5vdHmsUlAbKws5sZXs/HYx9xZ7/zyXt2MuaT4TuzhmmapgAAAAA4Rr1QBwAAAAAQWDT5AAAAgMPQ5AMAAAAOQ5MPAAAAOAxNPgAAAOAwNPkAAACAw9DkAwAAAA4THuoAwVZRUaFvv/1WTZs2lWEYoY4DAABqwTRNHT16VC1btlS9epyjBM7lomvyv/32WyUmJoY6BgAAOA/79+9Xq1atQh0DqPMuuia/adOmkir/J+FyuUKcBgAA1IbH41FiYqL39ziAs7vomvyqKToul4smHwAAm2GqLVA7TGoDAAAAHIYmHwAAAHAYmnwAAADAYWjyAQAAAIehyQcAAAAchiYfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwmJA2+fPnz1dqaqpcLpdcLpf69eunDz744KzrrF+/Xj179lSDBg3Url07vfTSS0FKCwAAANhDSJv8Vq1a6Q9/+IOys7OVnZ2ta665Rtddd52++uqrGsfv3btXI0aM0MCBA5WTk6NZs2Zp6tSpWrp0aZCTAwAAAHWXYZqmGeoQPxYdHa0nn3xSEyZMqPbafffdp+XLl2v79u3eZRMnTtQ///lPbdy4sVb1PR6P3G63iouL5XK5ApYbAABYh9/fgH/qzJz88vJyvf322yopKVG/fv1qHLNx40YNGTLEZ9nQoUOVnZ2tEydO1LhOaWmpPB6PzwMAAABwsvBQB9i6dav69eunH374QU2aNNF7772nyy67rMaxhYWFatGihc+yFi1a6OTJkzp06JDi4+OrrZORkaFHHnnEkuw1SeyQrLjEpIDX3bdju5KSUwJe1661yRyc2mQOTm07ZrayNpmDU9vKzMe+O3DGqbcAgiPkTX5ycrK2bNmiI0eOaOnSpRo3bpzWr19/xkbfMAyf51WzjU5fXiU9PV1paWne5x6PR4mJiQFKX11cYpLS5y8KeN2MSeMsqWvX2mQOTm0yB6e2HTNbWZvMwaltZebMmXdaUhdA7YW8yY+IiNCll14qSerVq5c+//xzzZ07V3/84x+rjY2Li1NhYaHPsqKiIoWHh+uSSy6psX5kZKQiIyMDHxwAauFISZkOekrV3BWpqMYRoY4DALhIhLzJP51pmiotLa3xtX79+mnFihU+y1atWqVevXqpfv36wYgHALX22c5DendDvkxTMgzphv6t1bdjTKhjAQAuAiG98HbWrFn65JNPlJeXp61bt+qBBx7QunXrdMstt0iqnGozduxY7/iJEydq3759SktL0/bt2/Xqq69q4cKFmjFjRqg2AQBqdKSkzNvgS5JpSu9uyNeRkrLQBgMAXBRCeib/u+++06233qqCggK53W6lpqbqww8/1LXXXitJKigoUH5+vnd827ZttXLlSk2fPl0vvPCCWrZsqXnz5mnMmDGh2gQAqNFBT6lOv0GxaUqHPKVM2wEAWC6kTf7ChQvP+npmZma1ZVdddZW++OILixIBQGA0d0XKMOTT6BuGFOPiGiEAgPXqzH3yAcBJohpH6Ib+rVV146+qOfmcxQcABEOdu/AWAJyib8cYdUpw6ZCnVDHcXQcAEEQ0+QBgoajGETT3AICgY7oOAAAA4DA0+QAAAIDD0OQDAAAADkOTDwAAADgMTT4AAADgMDT5AAAAgMPQ5AMAAAAOY5jmjz903fk8Ho/cbreKi4vlcrkCXj82IVFJySkBr7t/W6769OkT8LqSlLVhozqkdrektlW5yeyLzKds2rTJdseKXY/vXbk5tns/rPr5sOOxYmXmg3t2Ki8vL6A1rf79DTgNH4YVYEnJKUqfvyjgdTNn3qnly5cHvK4k9b5miCWZJetyk9kXmU8ZNWqU7Y4Vux7fGZPG2e79sOrnw47HitWZAYQW03UAAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGFo8gEAAACHockHAAAAHIYmHwAAAHAYmnwER/EBaW9W5VcAAABYig/DgvW+WCytuEcyKySjnjRyrtRjbKhTAQAAOBZn8mGtEz+cavClyq8rpnFGHwAAwEI0+bDWiZJTDX4Vs1w6vCc0eQAAAC4CNPmwVv3GlVN0fswIk6LbhSYPAADARYAmH9aq36ByDr4RVvncCJNGzpHcCSGNBQAA4GRceAvr9RgrtR9cOUUnuh0NPgAAgMVo8hEc7gSaewAAgCAxTNM0Qx0imDwej9xut4qLi+VyuQJePzYhUUnJKQGvuys3Rx1Suwe8rl1rkzk4tckcnNp2zGxlbTIHp7aVmQ/u2am8vLyA1rT69zfgNJzJD7Ck5BSlz18U8LoZk8ZZUteutckcnNpkDk5tO2a2sjaZg1PbysyZM++0pC6A2uPCWwAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGFo8gEAAACHockHAAAAHIYmHwAAAHAYmnwAAADAYWjyAQAAAIehycdF60hJmXYVHNWRkrJQRwEAAAio8FAHAELhs52H9O6GfJmmZBjSDf1bq2/HmFDHAgAACAjO5OOic6SkzNvgS5JpSu9uyOeMPgAAcAyafFx0DnpKvQ1+FdOUDnlKQxMIAAAgwAzTPL3dcTaPxyO3263i4mK5XK6A149NSFRSckrA6+7flqs+ffoEvK4kZW3YqA6p3S2pbVXuC8l8otzUzm891ZZ3bOlS/TCjTmY+FzKfsmnTJtsdK3Y9vnfl5tju/bDq58OOx4qVmQ/u2am8vLyA1rT69zfgNMzJD7Ck5BSlz18U8LqZM+/U8uXLA15XknpfM8SSzJJ1uS8089nm5NfVzGdD5lNGjRplu2PFrsd3xqRxtns/rPr5sOOxYnVmAKEV0uk6GRkZ6t27t5o2barY2FiNHj1aO3bsOOs669atk2EY1R5ff/11kFLDCfp2jNGDN3TRXcM66MEbunDRLQAAcJSQnslfv369Jk+erN69e+vkyZN64IEHNGTIEG3btk2NGzc+67o7duzw+ee65s2bWx0XDhPVOEJRjSNCHQMAACDgQtrkf/jhhz7PX3vtNcXGxmrz5s268sorz7pubGysoqKiLEwHAAAA2FOdurtOcXGxJCk6OvqcY7t37674+HgNHjxYa9euPeO40tJSeTwenwcAABet4gPS3iyFlXNHMcDJ6kyTb5qm0tLSdMUVV6hLly5nHBcfH68FCxZo6dKlWrZsmZKTkzV48GBlZWXVOD4jI0Nut9v7SExMtGoTAACo275YLM3pIi0aqdjCDUra/edQJwJgkTpzd527775bubm5+vTTT886Ljk5WcnJyd7n/fr10/79+/XUU0/VOMUnPT1daWlp3ucej4dGHwBw8Sk+IK24RzIrvIu6fz5bRfFX6HijuBAGA2CFOnEmf8qUKVq+fLnWrl2rVq1a+b1+3759tWvXrhpfi4yMlMvl8nkAAHDRObzbp8GXpHpmhRof3ReiQACsFNIz+aZpasqUKXrvvfe0bt06tW3b9rzq5OTkKD4+PsDpAABwkOj2klHPp9GvMOqppGlSCEMBsEpIm/zJkyfrzTff1F/+8hc1bdpUhYWFkiS3262GDRtKqpxuc+DAAS1evFiSNGfOHLVp00adO3dWWVmZXn/9dS1dulRLly4N2XYAAFDnuROkkXOlFdMks1ympJzejzJVB3CokDb58+fPlyQNGjTIZ/lrr72m8ePHS5IKCgqUn5/vfa2srEwzZszQgQMH1LBhQ3Xu3Fl/+9vfNGLEiGDFBgDAnnqMldoPlg7v0cG1v9O+9r8MdSIAFgn5dJ1zyczM9Hl+77336t5777UoEQAADudOkNwJKg+LDHUSABaqExfeAgAAAAgcmnwAAADAYWjyAQAAAIehyQcAAAAcxjBrc/Wrg3g8HrndbhUXF1vywVixCYlKSk4JeN1duTnqkNo94HXtWpvMwalN5uDUtmNmK2uTOTi1rcx8cM9O5eXlBbSm1b+/AacJ6d11nCgpOUXp8xcFvG7GpHGW1LVrbTIHpzaZg1PbjpnPVvtISZkOekrV3BWpqMYRAasbCOxD6+tKUubMOy2pC6D2aPIBAAHz2c5DendDvkxTMgzphv6t1bdjTKhjAcBFhzn5AICAOFJS5m3wJck0pXc35OtISVlogwHARYgmHwAQEAc9pTr9Ki/TlA55SkMTCAAuYjT5AICAaO6KlGH4LjMMKcbFJ6sCQLDR5AMAAiKqcYRu6N/a2+hXzck/n4tvAQAXhgtvAQAB07djjDoluHTIU6qY87y7DgDgwtHkAwACKqpxBM09AIQY03UAAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGFo8gEAAACHockHAAAAHIYmHwAAAHAYmnwAAADAYQzTNM1Qhwgmj8cjt9ut4uJiuVyugNePTUhUUnJKwOvu35arPn36BLyuJGVt2KgOqd0tqW1VbjL7IvMpmzZtst2xYtfje1duju3eD6t+Pux4rFiZ+eCencrLywtoTat/fwNOwyfeBlhScorS5y8KeN3MmXdq+fLlAa8rSb2vGWJJZsm63GT2ReZTRo0aZbtjxa7Hd8akcbZ7P6z6+bDjsWJ1ZgChxXQdAAAAwGFo8gEAAACHockHAAAAHIYmHwAAAHAYmnzA6YoPSHuzKr8CAICLAnfXAZzsi8XSinsks0Iy6kkj50o9xoY6FQAAsBhn8gGnKj5wqsGXKr+umMYZfQAALgI0+YBTHd59qsGvYpZLh/eEJg8AAAgamnzAqaLbV07R+TEjTIpuF5o8AAAgaGjyAadyJ1TOwTfCKp8bYdLIOZXLAQCAo3HhLeBkPcZK7QdXTtGJbkeDDwDARYImH3A6dwLNPQAAFxmm6wAAAAAOQ5MPAAAAOAxNPgAAAOAwhmmaZqhDBJPH45Hb7VZxcbFcLlfA68cmJCopOSXgdXfl5qhDaveA17VrbTIHpzaZg1PbjpmtrE1m6US5qbKT5YoID1PeV1tskfnHDu7Zqby8vIDWtPr3N+A0XHgbYEnJKUqfvyjgdTMmjbOkrl1rkzk4tckcnNp2zGxl7Ys982c7D+ndDfkyTem4IcX88GSdz3y6zJl3WlIXQO0xXQcAgDriSEmZt8GXJNOUvj18XEdKykIbDIDt0OQDAFBHHPSUqqZJtIc8pcEPA8DWaPIBAKgjmrsiZRjVl8e4IoMfBoCt0eQDAFBHRDWO0A39W3sbfcOQWkY3VFTjiNAGA2A7XHgLAEAd0rdjjDoluHTIU6oYV6Tm/4MGH4D/QnomPyMjQ71791bTpk0VGxur0aNHa8eOHedcb/369erZs6caNGigdu3a6aWXXgpCWgAAgiOqcYQujW/KGXwA5y2kTf769es1efJkffbZZ1q9erVOnjypIUOGqKSk5Izr7N27VyNGjNDAgQOVk5OjWbNmaerUqVq6dGkQkwMAAAB1V0in63z44Yc+z1977TXFxsZq8+bNuvLKK2tc56WXXlLr1q01Z84cSVJKSoqys7P11FNPacyYMVZHBgAAAOq8OnXhbXFxsSQpOjr6jGM2btyoIUOG+CwbOnSosrOzdeLEiWrjS0tL5fF4fB4AAACAk9WZJt80TaWlpemKK65Qly5dzjiusLBQLVq08FnWokULnTx5UocOHao2PiMjQ2632/tITEwMeHYAAACgLqkzTf7dd9+t3NxcvfXWW+cca5x2E2HzP58ccvpySUpPT1dxcbH3sX///sAEBgAAAOqoOnELzSlTpmj58uXKyspSq1atzjo2Li5OhYWFPsuKiooUHh6uSy65pNr4yMhIRUbyISIAAAC4eIT0TL5pmrr77ru1bNkyrVmzRm3btj3nOv369dPq1at9lq1atUq9evVS/fr1rYoKwKEaHitUzHefqeGxwnMPBgDAJkLa5E+ePFmvv/663nzzTTVt2lSFhYUqLCzU8ePHvWPS09M1duxY7/OJEydq3759SktL0/bt2/Xqq69q4cKFmjFjRig2AYCNJe3+s4b95RpduWa8hv3lGiXt/nOoIwEAEBAhbfLnz5+v4uJiDRo0SPHx8d7HkiVLvGMKCgqUn5/vfd62bVutXLlS69at0+WXX67HHntM8+bN4/aZAPzS8FihemyaLUMVkiRDFer++WzO6AMAHMEwq65avUh4PB653W4VFxfL5XIFvH5sQqKSklMCXnf/tlz16dMn4HUlKWvDRnVI7W5Jbatyk9kXmU/ZtGlT7eoe+5e0f1P15Yl9pEbVr++RrHs/7Hp878rNsd37UeufDz/Z8VixMvPBPTuVl5cX0JpW//4GnKZOXHjrJEnJKUqfvyjgdTNn3qnly5cHvK4k9b5miCWZJetyk9kXmU8ZNWpU7eoWH5DmdJHMilPLjDBp2juSO6HGVax6P+x6fGdMGme796PWPx9+suOxYnVmAKFVZ26hCQBB5U6QRs6tbOylyq8j55yxwQcAwE44kw/g4tVjrNR+sHR4jxTdjgYfAOAYNPkALm7uBJp7AIDjMF0HAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAABzGryb/+PHj+vTTT7Vt27Zqr/3www9avHhxwIIBAAAAOD+1bvJ37typlJQUXXnlleratasGDRqkgoIC7+vFxcX6zW9+Y0lIAAAAALVX6yb/vvvuU9euXVVUVKQdO3bI5XJpwIABys/PtzIfAAAAAD/VusnfsGGDnnjiCcXExOjSSy/V8uXLNXz4cA0cOFB79uyxMiMAAAAAP9T6w7COHz+u8HDf4S+88ILq1aunq666Sm+++WbAwwEAAADwX62b/E6dOik7O1spKSk+y5977jmZpqlRo0YFPBxC60hJmQ56StXcFamoxhGhjgMAAIBaqnWT/4tf/EJvvfWWbr311mqvPf/886qoqNBLL70U0HAInc92HtK7G/JlmpJhSDf0b62+HWNCHQsAAAC1YJimaYY6RDB5PB653W4VFxfL5XIFvH5sQqKSklPOPdBPu3Jz1CG1e8Dr1lT7RLmpnd96qo3r2NKl+mHGBdUOlGC+H3W9rpW1yRyc2nbMbGVtMgentpWZD+7Zqby8vIDWtPr3N+A0tT6Tj9pJSk5R+vxFAa+bMWmcJXVrqr2r4Kjmf7ir2ribh3XQpfFNL6h2oATz/ajrda2sTebg1LZjZitrkzk4ta3MnDnzTkvqAqg9PvEW1TR3Rco47YS9YUgxrsjQBAIAAIBfOJOPaqIaR+iG/q2rzcnn4lsAAKTy8nKdOHEi1DFwkalfv77CwsJqPZ4mHzXq2zFGnRJcOuQpVQx31wEAQKZpqrCwUEeOHAl1FFykoqKiFBcXJ+P0KRc1oMnHGUU1jqC5BwDgP6oa/NjYWDVq1KhWjRYQCKZp6tixYyoqKpIkxcfHn3Od82ryd+7cqXXr1qmoqEgVFRU+r82ePft8SgIAANRZ5eXl3gb/kksuCXUcXIQaNmwoSSoqKlJsbOw5p+743eS//PLLmjRpkmJiYqr9c4FhGDT5AADAcarm4Ddq1CjESXAxq/r5O3HiROCb/Mcff1y///3vdd99951fOgAAAJtiig5CyZ+fP79vofn999/rhhtu8Hc1AAAAAEHid5N/ww03aNWqVVZkAQAAQB1gGIbef//9s44ZP368Ro8eXeuaeXl5MgxDW7ZsuaBsqB2/p+tceumlevDBB/XZZ5+pa9euql+/vs/rU6dODVg4AAAAXJjx48fryJEj52zaf6ygoEDNmjWTVNmct23bVjk5Obr88su9Y+bOnSvTNAOaddCgQbr88ss1Z86cGl83DEPvvfdetT8uTt/GoqIiPfjgg/rggw/03XffqVmzZurWrZsefvhh9evXL6CZ6yq/m/wFCxaoSZMmWr9+vdavX+/zmmEYNPkAAAA2FxcXd84xbrc7CEnOz5gxY3TixAktWrRI7dq103fffaePP/5Yhw8fDnW0oPG7yd+7d68VOQAAAC4aBcXHtfdQidrGNFa8u2FQv/egQYOUmpqqBg0a6JVXXlFERIQmTpyohx9+2Dvmx2fM27ZtK0nq3r27JOmqq67SunXrqp09//DDD/X444/ryy+/VFhYmPr166e5c+eqffv2Qd2+I0eO6NNPP9W6det01VVXSZKSkpLUp0+foOYINb/n5P+YaZoB/2caAAAAJ1vyeb4G/GGNbn75HxrwhzVa8nl+0DMsWrRIjRs31j/+8Q/9z//8jx599FGtXr26xrGbNm2SJH300UcqKCjQsmXLahxXUlKitLQ0ff755/r4449Vr149/eIXv6j2mUpWa9KkiZo0aaL3339fpaWlQf3edcl5NfmLFy9W165d1bBhQzVs2FCpqan605/+FOhsAAAAjlJQfFzpy7aq4j/nSCtMadayL1VQfDyoOVJTU/XQQw+pQ4cOGjt2rHr16qWPP/64xrHNmzeXJF1yySWKi4tTdHR0jePGjBmj66+/Xh06dNDll1+uhQsXauvWrdq2bVtAs//qV7/yNvJVjzfeeMP7enh4uDIzM7Vo0SJFRUVpwIABmjVrlnJzcwOao67zu8l/5plnNGnSJI0YMULvvPOOlixZomHDhmnixIl69tlnrcgIAADgCHsPlXgb/Crlpqm8Q8eCmiM1NdXneXx8vIqKii6o5u7du3XzzTerXbt2crlc3mk++fmB/ZeKZ599Vlu2bPF5jBo1ymfMmDFj9O2332r58uUaOnSo1q1bpx49eigzMzOgWeoyv+fkP/fcc5o/f77Gjh3rXXbdddepc+fOevjhhzV9+vSABrSbfTu2K2PSuIDX3b8tV5kz7wx4XUnalZtjSWbJutxk9kVm6+tK1r0fdsxsZW0r3w+OlVOszHxwz05L6jpB25jGqmfIp9EPMwy1iQnuJ+mefndEwzAueFrNyJEjlZiYqJdfflktW7ZURUWFunTporKysguqe7q4uDhdeumlPsuaNm2qI0eO+Cxr0KCBrr32Wl177bWaPXu2br/9dj300EMaP358QPPUVX43+QUFBerfv3+15f3791dBQUFAQtlZUnKK0ucvCnjdzJl3avny5QGvK0m9rxliSWbJutxk9kXmU0aNGmW7Y8Wux3fGpHG2ez+s+vmw47FidWbULN7dUBnXd9WsZV+q3DQVZhh64vouQb/41h8RERGSpPLy8jOO+de//qXt27frj3/8owYOHChJ+vTTT4OSr7Yuu+wyv24janfndZ/8d955R7NmzfJZvmTJEnXo0CFgwQAAAJzoxt6tdWXH5so7dExtYhrV6QZfkmJjY9WwYUN9+OGHatWqlRo0aFDt9pnNmjXTJZdcogULFig+Pl75+fm6//77z+v7HTx4sNoHZsXFxdXqtp5S5R8cN9xwg2677TalpqaqadOmys7O1v/8z//ouuuuO69MduR3k//II4/oxhtvVFZWlgYMGCDDMPTpp5/q448/1jvvvGNFRgAAAEeJdzes8819lfDwcM2bN0+PPvqoZs+erYEDB2rdunU+Y+rVq6e3335bU6dOVZcuXZScnKx58+Zp0KBBfn+/N998U2+++abPsoceesjnFp9n06RJE/3kJz/Rs88+q927d+vEiRNKTEzUHXfcUe0ktZP53eSPGTNG//jHP/Tss8/q/fffl2mauuyyy7Rp0ybv/VMBAABQN5x+senpDbqkatNYTr9F+u23367bb7/9rHV/+tOfVruTzo/rtGnT5py3Xq8p29ly1ZQlMjJSGRkZysjIOGstp/O7yZeknj176vXXXw90FgAAAAABUKsm3+PxyOVyef/7bKrGAQAAAAiNWjX5zZo1U0FBgWJjYxUVFSXDMKqNMU1ThmGc9cprAAAAANarVZO/Zs0a76ebrV271tJAAAAAAC5MrZr8q666qsb/BgAAAFD31PN3hQ8//NDnww1eeOEFXX755br55pv1/fffBzQcAAAAAP/53eTPnDnTe/Ht1q1blZaWphEjRmjPnj1KS0sLeEAAAAAA/vH7Fpp79+7VZZddJklaunSpRo4cqSeeeEJffPGFRowYEfCAAAAAAPzj95n8iIgIHTt2TJL00UcfaciQIZKk6Ojoc95e83RZWVkaOXKkWrZsKcMwqn0Qw+nWrVsnwzCqPb7++mt/NwOAHRQfkPZmVX4FAAC15neTf8UVVygtLU2PPfaYNm3apJ/97GeSpJ07d6pVq1Z+1SopKVG3bt30/PPP+7Xejh07VFBQ4H106NDBr/UB2MAXi6U5XaRFIyu/frE41IkAABeJNm3aaM6cOaGOcUH8bvKff/55hYeH689//rPmz5+vhIQESdIHH3ygYcOG+VVr+PDhevzxx3X99df7tV5sbKzi4uK8j7CwML/WB1DHFR+QVtwjmRWVz80KacU0zugDwAXYsGGDwsLC/O7X7CKYjfnDDz+syy+//IyvDxo0SNOmTau2PDMzU1FRUd7n5eXlysjIUKdOndSwYUNFR0erb9++eu211y44o99z8lu3bq2//vWv1ZY/++yzFxymtrp3764ffvhBl112mX73u9/p6quvPuPY0tJSlZaWep/7O6UIQAgc3n2qwa9ilkuH90juhNBkAgCbe/XVVzVlyhS98sorys/PV+vWrUMd6aL38MMPa8GCBXr++efVq1cveTweZWdnB+SOlX43+ZJUUVGhb775RkVFRaqo8P1FfOWVV15wqDOJj4/XggUL1LNnT5WWlupPf/qTBg8erHXr1p3x+2ZkZOiRRx6xLNPp9u3YroxJ4wJed1dujnpfMyTgdatqW5G5qrYVuclcvbaTMoeVlyq28JjPMlPSwbW/U3lY5HnXvVBWvR92zGxlbavfDycdKxda16rMB/fstKSuoxQfqDyhEd0+KCcvSkpK9M477+jzzz9XYWGhMjMzNXv2bO/rmZmZmjZtmo4cOeJd9v777+sXv/iFTNP0Lnv88cc1b948HT9+XDfeeKNiYmL04YcfasuWLZKk8ePH68iRI+rTp4/mzp2r0tJSTZ8+XQ888IDS09O1cOFCNWrUSI8++qhuu+02b90DBw4oLS1Nq1atUr169XTFFVdo7ty5atOmjU/dK664Qk8//bTKysp00003ac6cOapfv74GDRqkffv2afr06Zo+fbokeXNv2LBB999/vz7//HPFxMToF7/4hTIyMtS4cWNJUlFRkSZMmKCPPvpIcXFxevzxx63YBTVasWKF7rrrLt1www3eZd26dQtIbb+b/M8++0w333yz9u3b57PTJckwDJWXlwckWE2Sk5OVnJzsfd6vXz/t379fTz311Bmb/PT0dJ9be3o8HiUmJlqWMSk5RenzFwW8bsakcZbUtWttMgendigzJ+3+s7p/Plv1zApVGPWU0/tR7Wv/ywuueyGc+D7XxdpkDk5tKzNnzrzTkrqO8cXiU1MSjXrSyLlSj7GWfsslS5Z4+6hf//rXmjJlih588EEZhlHrGm+88YZ+//vf68UXX9SAAQP09ttv6+mnn1bbtm19xq1Zs0atWrVSVlaW/v73v2vChAnauHGjrrzySv3jH//QkiVLNHHiRF177bVKTEzUsWPHdPXVV2vgwIHKyspSeHi4Hn/8cQ0bNky5ubmKiIiQJK1du1bx8fFau3atvvnmG9144426/PLLdccdd2jZsmXq1q2b/uu//kt33HGHN8vWrVs1dOhQPfbYY1q4cKEOHjyou+++W3fffbd3Ssz48eO1f/9+rVmzRhEREZo6daqKiooC8K6fW1xcnNasWaO77rpLzZs3D2htv+fkT5w4Ub169dKXX36pw4cP6/vvv/c+Dh8+HNBwtdG3b1/t2rXrjK9HRkbK5XL5PADUffva/1L/O2qNsq5ZpP8dtaZWDT4A1HkhuuZo4cKF+vWvfy1JGjZsmP7973/r448/9qvGc889pwkTJug3v/mNOnbsqNmzZ6tr167VxkVHR2vevHlKTk7WbbfdpuTkZB07dkyzZs1Shw4dlJ6eroiICP3973+XJL399tuqV6+eXnnlFXXt2lUpKSl67bXXlJ+fr3Xr1nnrNmvWTM8//7w6deqkn//85/rZz37m3Ybo6GiFhYWpadOm3ms2JenJJ5/UzTffrGnTpqlDhw7q37+/5s2bp8WLF+uHH37Qzp079cEHH+iVV15Rv3791LNnTy1cuFDHjx8/n7fZx4svvqgmTZr4PCZOnOgz5plnntHBgwcVFxen1NRUTZw4UR988MEFf2/pPJr8Xbt26YknnlBKSoqioqLkdrt9HsGWk5Oj+Pj4oH9fANY73ihOh1r8RMcbxYU6CgAExtmuObLIjh07tGnTJt10002SpPDwcN1444169dVX/a7Tp08fn2WnP5ekzp07q169Uy1mixYtfP4YCAsL0yWXXOI9W75582Z98803atq0qbcZjo6O1g8//KDdu3f71P3xzVbi4+PPecZ98+bNyszM9Gm0hw4dqoqKCu3du1fbt29XeHi4evXq5V2nU6dOPhfHnq9bbrlFW7Zs8Xk8+uijPmMuu+wyffnll/rss8/0m9/8Rt99951Gjhyp22+//YK/v9/TdX7yk5/om2++0aWXXnrB3/zf//63vvnmG+/zvXv3asuWLYqOjlbr1q2Vnp6uAwcOaPHiylvnzZkzR23atFHnzp1VVlam119/XUuXLtXSpUsvOAsAAIDlottXTtH5caNvhEnR7Sz7lgsXLtTJkye9d0SUKuer169fX99//72aNWumevXqVZuGfeLEiWq1Tp/ec/o6klS/fv1q69S0rOq6zoqKCvXs2VNvvPFGtVo/nsJythpnUlFRoTvvvFNTp06t9lrr1q21Y8cOb61Ac7vd1frl2NjYauPq1aun3r17q3fv3po+fbpef/113XrrrXrggQeqTYXyh99N/pQpU/Tb3/5WhYWF6tq1a7U3PDU1tda1srOzfe6MUzV3fty4ccrMzFRBQYHy8/O9r5eVlWnGjBk6cOCAGjZsqM6dO+tvf/sbn7QLAADswZ1QOQd/xbTKM/hGmDRyjmUX3548eVKLFy/W008/7f0A0ypjxozRG2+8obvvvlvNmzfX0aNHVVJS4r0gtepi2irJycnatGmTbr31Vu+y7OzsC87Yo0cPLVmyRLGxsRc0rToiIqLataE9evTQV199dcaT0ykpKTp58qSys7O9/yqxY8cOnwuQg+2yyy6TVHmx9IXwu8kfM2aMJPlcEW0YhkzT9PvC20GDBtX4F2CVzMxMn+f33nuv7r33Xv8CAwAA1CU9xkrtB1dO0YluZ+nddf7617/q+++/14QJE6pNq/7lL3+phQsX6u6779ZPfvITNWrUSLNmzdKUKVO0adOman3YlClTdMcdd6hXr17q37+/lixZotzcXLVrd2H/CnHLLbfoySef1HXXXadHH31UrVq1Un5+vpYtW6aZM2fW+sNW27Rpo6ysLN10002KjIxUTEyM7rvvPvXt21eTJ0/WHXfcocaNG2v79u1avXq1nnvuOSUnJ2vYsGG64447tGDBAoWHh2vatGlq2LDhOb/f8ePHq/0h1KRJE79mu/zyl7/UgAED1L9/f8XFxWnv3r1KT09Xx44d1alTp1rXqYnfc/L37t1b7bFnzx7vVwAAAJyDO0FqO9Dy22cuXLhQP/3pT2u8bnLMmDHasmWLvvjiC0VHR+v111/XypUr1bVrV7311lt6+OGHfcbfcsstSk9P14wZM9SjRw/t3btX48ePV4MGDS4oY6NGjZSVlaXWrVvr+uuvV0pKim677TYdP37crzP7jz76qPLy8tS+fXvvNJ/U1FStX79eu3bt0sCBA9W9e3c9+OCDPtdzvvbaa0pMTNRVV12l66+/Xv/1X/9V47Sa0+3cuVPdu3f3efg7l37o0KFasWKFRo4cqY4dO2rcuHHq1KmTVq1apfDw87rTvZffayclJV3QNwQAAEBwrFix4oyv9ejRw2dGxejRozV69GifMT++HaUkPfjgg3rwwQe9z6+99lqfM9enn/2X5HOHnCp5eXk+z+Pi4rRo0Zlv6VpT3dM/3bZv37765z//WW1c7969tWrVqjPWjouLq/ZBrz+eklSThx9+uNofQT9W0zZLlbfrHD9+vPf5HXfcUe09DhS/z+RL0p/+9CcNGDBALVu21L59+yRVvtF/+ctfAhoOAAAAdcOxY8f0zDPP6KuvvtLXX3+thx56SB999JHGjbPmQ9VwYfxu8ufPn6+0tDSNGDFCR44c8c7Bj4qKqvYXFQAAAJzBMAytXLlSAwcOVM+ePbVixQotXbpUP/3pT0MdDTXwe7rOc889p5dfflmjR4/WH/7wB+/yXr16acaMGQENBwAAgLqhYcOG+uijj0IdA7V0Xhfedu/evdryyMjIC77VD045UlKmXQVHdaSkLNRRAAAAYDN+n8lv27attmzZUu0C3A8++MB7X09cmM92HtK7G/JlmpJhSDf0bx3qSAAAQDV/+BMQLP78/Pnd5M+cOVOTJ0/WDz/8INM0tWnTJr311lvKyMjQK6+84m85nOZISZm3wZck05Te3ZAvdzn/UwEAIFSqPvzz2LFjtbqHOmCFY8eOSar+6b818bvJ/81vfqOTJ0/q3nvv1bFjx3TzzTcrISFBc+fO1U033eR/Wvg46CnV6X+kmaZUdrL2HzIGAAACKywsTFFRUSoqKpJUeW93wzBCnAoXC9M0dezYMRUVFSkqKkphYWHnXOe87rJfdU/PQ4cOqaKiolYfGIDaae6KlGHIp9E3DCki/Nw7EwAQHA2PFarx0TyVNG2j443iQh0HQRIXV7mvqxp9INiioqK8P4fnYpgX2eQyj8cjt9ut4uJivz5FrbZiExKVlJxyQTW+LynTt4ePe5+3jG6of+/7Wn369LnQeDXK2rBRHVKrX0wdCPu35VqSm8y+yHzKpk2bbHesWPVeSNbuw125ObZ7PwLy81G8Xyr88tTzuC7K+ur/bHesWPmzcXDPzmofdnShrP797Y/y8nKdOHEipBlw8alfv36tzuBX8ftM/r/+9S/Nnj1ba9euVVFRkSoqKnxeP3z4sL8lHSUpOUXp88/8iW21daSkTIc8pYpxRSqqcYQyZ96p5cuXByBhdb2vGRKQzDWxKjeZfZH5lFGjRtnuWLHr8Z0xaZzt3o8L/vkoPiDN6SKZjU4tM/apr7rqXpsdK1Yf304WFhbmV7MFhILfTf6vf/1r7d69WxMmTFCLFi2Yj2aRqMYRimocEeoYAIAfO7xbMn1PbsksV9jJ4zWPB4AQ8bvJ//TTT/Xpp5+qW7duVuQBAKDuim4vGfV8G30jTOXh3G0FQN3i94dhderUScePc8YCAHARcidII+dKxn+mahhh0sg5Kg+LDG0uADiN32fyX3zxRd1///2aPXu2unTpUu0+naG+GAYAAEv1GCu1Hywd3iNFt6ts/PV6qFMBgA+/m/yoqCgVFxfrmmuu8VlumqYMw1B5OfdzBwA4nDvhP809ANRNfjf5t9xyiyIiIvTmm29y4S0AAABQB/nd5H/55ZfKyclRcnKyFXkAAAAAXCC/L7zt1auX9u/fb0UWAAAAAAHg95n8KVOm6J577tHMmTPVtWvXahfepqamBiwcAAAAAP/53eTfeOONkqTbbrvNu8wwDC68BQAAAOoIv5v8vXv3WpEDAAAAQID43eQnJSVZkQMAAABAgNSqyV++fLmGDx+u+vXra/ny5WcdO2rUqIAEAwAAAHB+atXkjx49WoWFhYqNjdXo0aPPOI45+QAAAEDo1arJr6ioqPG/AQAAANQ9hmmaZqhDBJPH45Hb7VZxcbFcLlfA68cmJCopOSXgdXfl5qhDaveA17VrbTIHpzaZg1PbjpmtrE3m4NS2MvPBPTuVl5cX0JpW//4GnMavC28rKiqUmZmpZcuWKS8vT4ZhqG3btvrlL3+pW2+9VYZhWJXTNpKSU5Q+f1HA62ZMGmdJXbvWJnNwapM5OLXrcuYjJWU66ClVc1ekohpHBLT2mdTl9yPYda2sbWXmzJl3WlIXQO3Vusk3TVOjRo3SypUr1a1bN3Xt2lWmaWr79u0aP368li1bpvfff9/CqACAYPps5yG9uyFfpikZhnRD/9bq2zEm1LEAALVQ6yY/MzNTWVlZ+vjjj3X11Vf7vLZmzRqNHj1aixcv1tixYwMeEgAQXEdKyrwNviSZpvTuhnx1SnBVO6MPAKh76tV24FtvvaVZs2ZVa/Al6ZprrtH999+vN954I6DhAAChcdBTqtOv2DJN6ZCnNDSBAAB+qXWTn5ubq2HDhp3x9eHDh+uf//xnQEIBAEKruStSp19mZRhSjCsyNIEAAH6pdZN/+PBhtWjR4oyvt2jRQt9//31AQgEAQiuqcYRu6N/a2+hXzclnqg4A2EOt5+SXl5crPPzMw8PCwnTy5MmAhAIAhF7fjjHqlODSIU+pYmq4uw4AoO7y6+4648ePV2Rkzf9UW1rKPE0AcJqoxhE09wBgQ7Vu8seNG3fOMdxZBwAAAAi9Wjf5r732mpU5AAAAAARIrS+8BQAAAGAPNPkAAACAw9DkAwAAAA5Dkw8AAAA4DE0+AAAA4DCGaZpmqEMEk8fjkdvtVnFxsVwuV8DrxyYkKik5JeB192/LVZ8+fQJeV5KyNmxUh9TultS2KjeZfZH5lE2bNtnuWLHr8b0rN8d274dVPx92PFaszHxwz07l5eUFtKbVv78Bp6n1LTRRO0nJKUqfvyjgdTNn3qnly5cHvK4k9b5miCWZJetyk9kXmU8ZNWqU7Y4Vux7fGZPG2e79sOrnw47HitWZAYRWSKfrZGVlaeTIkWrZsqUMw9D7779/znXWr1+vnj17qkGDBmrXrp1eeukl64MCAAAANhLSJr+kpETdunXT888/X6vxe/fu1YgRIzRw4EDl5ORo1qxZmjp1qpYuXWpxUgAAAMA+QjpdZ/jw4Ro+fHitx7/00ktq3bq15syZI0lKSUlRdna2nnrqKY0ZM8ailAAAAIC92OruOhs3btSQIUN8lg0dOlTZ2dk6ceJEjeuUlpbK4/H4PAAAAAAns1WTX1hYqBYtWvgsa9GihU6ePKlDhw7VuE5GRobcbrf3kZiYGIyoAAAAQMjYqsmXJMMwfJ5X3QH09OVV0tPTVVxc7H3s37/f8owAAABAKNnqFppxcXEqLCz0WVZUVKTw8HBdcsklNa4TGRmpyMjIYMQDAAAA6gRbncnv16+fVq9e7bNs1apV6tWrl+rXrx+iVAAAAEDdEtIm/9///re2bNmiLVu2SKq8ReaWLVuUn58vqXKqzdixY73jJ06cqH379iktLU3bt2/Xq6++qoULF2rGjBmhiA8AAADUSSGdrpOdna2rr77a+zwtLU2SNG7cOGVmZqqgoMDb8EtS27ZttXLlSk2fPl0vvPCCWrZsqXnz5nH7TAAAAOBHQtrkDxo0yHvhbE0yMzOrLbvqqqv0xRdfWJgKOE3xAenwbim6veROCHUaAACAc7LVhbdA0H2xWFpxj2RWSEY9aeRcqcfYc68HAAAQQra68BYIquIDpxp8qfLrimmVywEAAOowmnzgTA7vPtXgVzHLpcN7QpMHAACglmjygTOJbl85RefHjDApul1o8gAAANSSYZ7tylcH8ng8crvdKi4ulsvlCnj92IREJSWnBLzurtwcdUjtHvC6dq0drMwNSwrkPvK1DEmmpOKoTjreOP6C6waa3d9nO9S1srYdM1tZm8zBqW1l5oN7diovLy+gNa3+/Q04DRfeBlhScorS5y8KeN2MSeMsqWvX2sHM3PBYoRof3aeSpkk63iguYHUDyQnvc12va2VtO2a2sjaZg1PbysyZM++0pC6A2qPJB87heKO4C2ruAQAAgo05+QAAAIDD0OQDAAAADkOTDwAAADgMTT4AAADgMDT5AAAAgMPQ5AMAAAAOQ5MPAAAAOAxNPgAAAOAwNPkAAACAw9DkAwAAAA5Dkw/AckdKyrSr4KiOlJSFOgoAABeF8FAHAOBsn+08pHc35Ms0JcOQbujfWn07xoQ6FgAAjsaZfACWOVJS5m3wJck0pXc35HNGHwAAi9HkA7DMQU+pt8GvYprSIU9paAIBAHCRMEzz9F/BzubxeOR2u1VcXCyXyxXw+rEJiUpKTgl43f3bctWnT5+A15WkrA0b1SG1uyW1rcpNZl91NfOJclM7v/VUW96xpUuFO7ZaknnTpk22O1bsenzvys2x3fth1c/HxXh8n83BPTuVl5cX0JpW//4GnIY5+QGWlJyi9PmLAl43c+adWv6n+dLh3VJ0e8mdELDava8ZYklm6T+5ly8PeF0y+6rLmc80J9+qzKNGjbKkrmTdPrTqvZCs/bnLmDTOdu+HVT8fF+vxfSaZM++0pC6A2qPJt4vi/dKcLpJZIRn1pJFzpR5jQ50KOKe+HWPUKcGlQ55SxbgiFdU4ItSRAABwPObk20DDY4VS4ZeVDb5U+XXFNKn4QEhzAbUV1ThCl8Y3pcEHACBIaPJtoPHRvOoLzXLp8J6gZwEAAEDdR5NvAyVN21RfaIRJ0e2CngUAAAB1H02+DRxvFCfFdals7KXKryPnBPTiWwAAADgHF97ahTtRmvbXyik60e1o8AEAAHBGNPl24k6guQcAAMA5MV0HAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGFo8gEAAACHockHAAAAHMYwTdMMdYhg8ng8crvdKi4ulsvlCnj92IREJSWnBLzu/m256tOnT8DrSlLWho3qkNrdktpW5SazLzKfsmnTJtsdK3Y9vnfl5tju/bDq58OOx4qVmQ/u2am8vLyA1rT69zfgNHzibYAlJacoff6igNfNnHmnli9fHvC6ktT7miGWZJasy01mX2Q+ZdSoUbY7Vux6fGdMGme798Oqnw87HitWZwYQWkzXAQAAAByGJh8A7Kr4gLQ3q/IrAAA/wnQdALCjLxZLK+6RzArJqCeNnCv1GBvqVACAOoIz+QBgN8UHTjX4UuXXFdM4ow8A8KLJBwC7Obz7VINfxSyXDu8JTR4AQJ0T8ib/xRdfVNu2bdWgQQP17NlTn3zyyRnHrlu3ToZhVHt8/fXXQUwMACEW3b5yis6PGWFSdLvQ5AEA1DkhbfKXLFmiadOm6YEHHlBOTo4GDhyo4cOHKz8//6zr7dixQwUFBd5Hhw4dgpQYAOoAd0LlHHwjrPK5ESaNnFO5HAAAhfjC22eeeUYTJkzQ7bffLkmaM2eO/vd//1fz589XRkbGGdeLjY1VVFRUkFICQB3UY6zUfnDlFJ3odjT4AAAfITuTX1ZWps2bN2vIkCE+y4cMGaINGzacdd3u3bsrPj5egwcP1tq1a886trS0VB6Px+cBAI7gTpDaDqTBBwBUE7Im/9ChQyovL1eLFi18lrdo0UKFhYU1rhMfH68FCxZo6dKlWrZsmZKTkzV48GBlZWWd8ftkZGTI7XZ7H4mJiQHdDgAAAKCuCfl98g3D8Hlumma1ZVWSk5OVnJzsfd6vXz/t379fTz31lK688soa10lPT1daWpr3ucfjodEHAACAo4XsTH5MTIzCwsKqnbUvKiqqdnb/bPr27atdu3ad8fXIyEi5XC6fBwAAAOBkIWvyIyIi1LNnT61evdpn+erVq9W/f/9a18nJyVF8fHyg4wEAAAC2FdLpOmlpabr11lvVq1cv9evXTwsWLFB+fr4mTpwoqXKqzYEDB7R48WJJlXffadOmjTp37qyysjK9/vrrWrp0qZYuXRrKzQAAAADqlJA2+TfeeKP+9a9/6dFHH1VBQYG6dOmilStXKikpSZJUUFDgc8/8srIyzZgxQwcOHFDDhg3VuXNn/e1vf9OIESNCtQkAAABAnWOYpmmGOkQweTweud1uFRcXWzI/PzYhUUnJKQGvuys3Rx1Suwe8rl1rkzk4tckcnNp2zGxlbTIHp7aVmQ/u2am8vLyA1rT69zfgNCG/u47TJCWnKH3+ooDXzZg0zpK6dq1N5uDUJnNwatsxs5W1yRyc2lZmzpx5pyV1AdReyC68BQAAAGANmnwAAADAYWjyAQAAAIehyQcAAAAchiYfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAFRzpKRMuwqO6khJWaijADgPfOItAADw8dnOQ3p3Q75MUzIM6Yb+rdW3Y0yoYwHwA2fyAQCA15GSMm+DL0mmKb27IZ8z+oDN0OQDAACvg55Sb4NfxTSlQ57S0AQCcF5o8gEAgFdzV6QMw3eZYUgxrsjQBAJwXmjyAQCAV1TjCN3Qv7W30a+akx/VOCK0wQD4xTDN0/9Rztk8Ho/cbreKi4vlcrkCXj82IVFJySkBr7srN0cdUrsHvK5da5M5OLXJHJzadsxsZW0yB6f2ueqeKDdVdrJcEeFhqh9mnHFcTQ7u2am8vLwLTOjL6t/fgNNwd50AS0pOUfr8RQGvmzFpnCV17VqbzMGpTebg1LZjZitrkzk4ta3MnDnzTkvqAqg9pusAAAAADkOTDwAAADgMTT4AAADgMDT5AAAAgMPQ5AMAAAAOQ5MP1BFHSsq0q+AoHx0PAAAuGLfQBOqAz3Ye0rsb8mWapz54pm/HmFDHAgAANsWZfCDEjpSUeRt8STJN6d0N+ZzRBwAA540mHwixg55Snf6506YpHfKUhiYQAACwPZp8IMSauyJlnPaJ8YYhxbgiQxMIAADYHk0+EGJRjSN0Q//W3ka/ak5+VOOI0AYDAAC2xYW3QB3Qt2OMOiW4dMhTqhhXJA0+AAC4IDT5QB0R1TiC5h4AAAQE03UAAAAAh6HJBwAAABzGMM3Tb97nbB6PR263W8XFxXK5XAGvH5uQqKTklIDX3b8tV3369Al4XUnK2rBRHVK7W1Lbqtxk9kXmUzZt2mS7Y8Wux/eu3BzbvR9W/XzY8VixMvPBPTuVl5cX0JpW//4GnIY5+QGWlJyi9PmLAl43c+adWr58ecDrSlLva4ZYklmyLjeZfZH5lFGjRtnuWLHr8Z0xaZzt3g+rfj7seKxYnRlAaDFdBwAAAHAYmnwAAADAYWjyAQAAAIehyQcAAAAchiYfAAAAcBiafAAAAMBhaPIBAM5SfEDam1X5FQAuUtwnHwDgHF8sllbcI5kVklFPGjlX6jE21KkAIOg4kw8AcIbiA6cafKny64ppnNEHcFGiyQcAOMPh3aca/CpmuXR4T2jyAEAI0eQDAJwhun3lFJ0fM8Kk6HahyQMAIRTyJv/FF19U27Zt1aBBA/Xs2VOffPLJWcevX79ePXv2VIMGDdSuXTu99NJLQUoKAKjT3AmVc/CNsMrnRpg0ck7lcgC4yIT0wtslS5Zo2rRpevHFFzVgwAD98Y9/1PDhw7Vt2za1bt262vi9e/dqxIgRuuOOO/T666/r73//u+666y41b95cY8aMCcEWAADqlB5jpfaDK6foRLejwQdw0QrpmfxnnnlGEyZM0O23366UlBTNmTNHiYmJmj9/fo3jX3rpJbVu3Vpz5sxRSkqKbr/9dt1222166qmngpwcAFBnuROktgNp8AFc1ELW5JeVlWnz5s0aMmSIz/IhQ4Zow4YNNa6zcePGauOHDh2q7OxsnThxosZ1SktL5fF4fB4AAACAk4Vsus6hQ4dUXl6uFi1a+Cxv0aKFCgsLa1ynsLCwxvEnT57UoUOHFB8fX22djIwMPfLII4ELfg7HvjugzJl3Brxubm6uRo0aFfC6knRwz05LMkvW5SazLzJbX1ey7v2wY2Yra1v5fnCsnGJl5t27d1tSF0DthfzDsAzD8Hlumma1ZecaX9PyKunp6UpLS/M+93g8SkxMPN+45/TVV19ZVhsAAACojZA1+TExMQoLC6t21r6oqKja2foqcXFxNY4PDw/XJZdcUuM6kZGRioyMDExoAAAAwAZCNic/IiJCPXv21OrVq32Wr169Wv37969xnX79+lUbv2rVKvXq1Uv169e3LCsAAABgJyG9u05aWppeeeUVvfrqq9q+fbumT5+u/Px8TZw4UVLlVJuxY8d6x0+cOFH79u1TWlqatm/frldffVULFy7UjBkzQrUJAAAAQJ0T0jn5N954o/71r3/p0UcfVUFBgbp06aKVK1cqKSlJklRQUKD8/Hzv+LZt22rlypWaPn26XnjhBbVs2VLz5s3jHvkAAADAjxhm1ZWrFwmPxyO3263i4mK5XK5QxwEAALXA72/APyGdrgMAAAAg8GjyAQAAAIehyQcAAAAchiYfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGHCQx0g2Ko+4Nfj8YQ4CQAAqK2q39tVv8cBnN1F1+QfPXpUkpSYmBjiJAAAwF9Hjx6V2+0OdQygzjPMi+xP4oqKCn377bdq2rSpDMMIaG2Px6PExETt379fLpcroLXrAqdvn+T8bWT77M/p28j22Z9V22iapo4ePaqWLVuqXj1mGwPnctGdya9Xr55atWpl6fdwuVyO/Z+35Pztk5y/jWyf/Tl9G9k++7NiGzmDD9QefwoDAAAADkOTDwAAADgMTX4ARUZG6qGHHlJkZGSoo1jC6dsnOX8b2T77c/o2sn32dzFsI2AHF92FtwAAAIDTcSYfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJt9PL774otq2basGDRqoZ8+e+uSTT846fv369erZs6caNGigdu3a6aWXXgpS0vPjz/atW7dOhmFUe3z99ddBTFx7WVlZGjlypFq2bCnDMPT++++fcx277T9/t9FO+zAjI0O9e/dW06ZNFRsbq9GjR2vHjh3nXM9O+/B8ttFO+3D+/PlKTU31fkhSv3799MEHH5x1HTvtP3+3z077riYZGRkyDEPTpk076zg77UPASWjy/bBkyRJNmzZNDzzwgHJycjRw4EANHz5c+fn5NY7fu3evRowYoYEDByonJ0ezZs3S1KlTtXTp0iAnrx1/t6/Kjh07VFBQ4H106NAhSIn9U1JSom7duun555+v1Xi77T/J/22sYod9uH79ek2ePFmfffaZVq9erZMnT2rIkCEqKSk54zp224fns41V7LAPW7VqpT/84Q/Kzs5Wdna2rrnmGl133XX66quvahxvt/3n7/ZVscO+O93nn3+uBQsWKDU19azj7LYPAUcxUWt9+vQxJ06c6LOsU6dO5v3331/j+Hvvvdfs1KmTz7I777zT7Nu3r2UZL4S/27d27VpTkvn9998HIV1gSTLfe++9s46x2/47XW220c77sKioyJRkrl+//oxj7L4Pa7ONdt6HpmmazZo1M1955ZUaX7P7/jPNs2+fXffd0aNHzQ4dOpirV682r7rqKvOee+4541gn7EPArjiTX0tlZWXavHmzhgwZ4rN8yJAh2rBhQ43rbNy4sdr4oUOHKjs7WydOnLAs6/k4n+2r0r17d8XHx2vw4MFau3atlTGDyk7770LZcR8WFxdLkqKjo884xu77sDbbWMVu+7C8vFxvv/22SkpK1K9fvxrH2Hn/1Wb7qtht302ePFk/+9nP9NOf/vScY+28DwG7o8mvpUOHDqm8vFwtWrTwWd6iRQsVFhbWuE5hYWGN40+ePKlDhw5ZlvV8nM/2xcfHa8GCBVq6dKmWLVum5ORkDR48WFlZWcGIbDk77b/zZdd9aJqm0tLSdMUVV6hLly5nHGfnfVjbbbTbPty6dauaNGmiyMhITZw4Ue+9954uu+yyGsfacf/5s31223eS9Pbbb2vz5s3KyMio1Xg77kPAKcJDHcBuDMPweW6aZrVl5xpf0/K6wp/tS05OVnJysvd5v379tH//fj311FO68sorLc0ZLHbbf/6y6z68++67lZubq08//fScY+26D2u7jXbbh8nJydqyZYuOHDmipUuXaty4cVq/fv0ZG2G77T9/ts9u+27//v265557tGrVKjVo0KDW69ltHwJOwZn8WoqJiVFYWFi1s9pFRUXVzlJUiYuLq3F8eHi4LrnkEsuyno/z2b6a9O3bV7t27Qp0vJCw0/4LpLq+D6dMmaLly5dr7dq1atWq1VnH2nUf+rONNanL+zAiIkKXXnqpevXqpYyMDHXr1k1z586tcawd958/21eTurzvNm/erKKiIvXs2VPh4eEKDw/X+vXrNW/ePIWHh6u8vLzaOnbch4BT0OTXUkREhHr27KnVq1f7LF+9erX69+9f4zr9+vWrNn7VqlXq1auX6tevb1nW83E+21eTnJwcxcfHBzpeSNhp/wVSXd2Hpmnq7rvv1rJly7RmzRq1bdv2nOvYbR+ezzbWpK7uw5qYpqnS0tIaX7Pb/qvJ2bavJnV53w0ePFhbt27Vli1bvI9evXrplltu0ZYtWxQWFlZtHSfsQ8C2QnK5r029/fbbZv369c2FCxea27ZtM6dNm2Y2btzYzMvLM03TNO+//37z1ltv9Y7fs2eP2ahRI3P69Onmtm3bzIULF5r169c3//znP4dqE87K3+179tlnzffee8/cuXOn+eWXX5r333+/KclcunRpqDbhrI4ePWrm5OSYOTk5piTzmWeeMXNycsx9+/aZpmn//Wea/m+jnfbhpEmTTLfbba5bt84sKCjwPo4dO+YdY/d9eD7baKd9mJ6ebmZlZZl79+41c3NzzVmzZpn16tUzV61aZZqm/fefv9tnp313JqffXcfu+xBwEpp8P73wwgtmUlKSGRERYfbo0cPn1nbjxo0zr7rqKp/x69atM7t3725GRESYbdq0MefPnx/kxP7xZ/v++7//22zfvr3ZoEEDs1mzZuYVV1xh/u1vfwtB6tqpul3d6Y9x48aZpumM/efvNtppH9a0XZLM1157zTvG7vvwfLbRTvvwtttu8/7/pXnz5ubgwYO9DbBp2n//+bt9dtp3Z3J6k2/3fQg4iWGa/7kCBgAAAIAjMCcfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGFo8gEAAACHockHAAAAHIYmH0DAGYah999/P9QxzmrdunUyDENHjhwJdRQAAAKOJh9ArYwfP16GYcgwDNWvX18tWrTQtddeq1dffVUVFRU+YwsKCjR8+PAQJa2d/v37q6CgQG6329Lvk5WVpZEjR6ply5a2+OMHAOAMNPkAam3YsGEqKChQXl6ePvjgA1199dW655579POf/1wnT570jouLi1NkZGQIk55bRESE4uLiZBiGpd+npKRE3bp10/PPP2/p9wEA4Mdo8gHUWmRkpOLi4pSQkKAePXpo1qxZ+stf/qIPPvhAmZmZ3nE/PmOdl5cnwzD0zjvvaODAgWrYsKF69+6tnTt36vPPP1evXr3UpEkTDRs2TAcPHvT5fq+99ppSUlLUoEEDderUSS+++KL3taq6y5Yt09VXX61GjRqpW7du2rhxo3fMvn37NHLkSDVr1kyNGzdW586dtXLlSkk1T9dZunSpOnfurMjISLVp00ZPP/20T542bdroiSee0G233aamTZuqdevWWrBgwVnfs+HDh+vxxx/X9ddf789bDQDABaHJB3BBrrnmGnXr1k3Lli0767iHHnpIv/vd7/TFF18oPDxcv/rVr3Tvvfdq7ty5+uSTT7R7927Nnj3bO/7ll1/WAw88oN///vfavn27nnjiCT344INatGiRT90HHnhAM2bM0JYtW9SxY0f96le/8v6rwuTJk1VaWqqsrCxt3bpV//3f/60mTZrUmG/z5s36f//v/+mmm27S1q1b9fDDD+vBBx/0+eNFkp5++mn16tVLOTk5uuuuuzRp0iR9/fXX5/HOAQBgnfBQBwBgf506dVJubu5Zx8yYMUNDhw6VJN1zzz361a9+pY8//lgDBgyQJE2YMMGnoX7sscf09NNPe8+At23bVtu2bdMf//hHjRs3zqfuz372M0nSI488os6dO+ubb75Rp06dlJ+frzFjxqhr166SpHbt2p0x3zPPPKPBgwfrwQcflCR17NhR27Zt05NPPqnx48d7x40YMUJ33XWXJOm+++7Ts88+q3Xr1qlTp061easAAAgKzuQDuGCmaZ5zbntqaqr3v1u0aCFJ3ua7allRUZEk6eDBg9q/f78mTJigJk2aeB+PP/64du/efca68fHxkuStM3XqVD3++OMaMGCAHnroobP+IbJ9+3bvHxxVBgwYoF27dqm8vLzG72cYhuLi4rzfDwCAuoImH8AF2759u9q2bXvWMfXr1/f+d9UfBKcvq7pLT9XXl19+WVu2bPE+vvzyS3322WfnrFu1/u233649e/bo1ltv1datW9WrVy8999xzNear6Q8V0zTPuh2n5wYAoK6gyQdwQdasWaOtW7dqzJgxAavZokULJSQkaM+ePbr00kt9Huf6Y+J0iYmJmjhxopYtW6bf/va3evnll2scd9lll+nTTz/1WbZhwwZ17NhRYWFh570tAACEAnPyAdRaaWmpCgsLVV5eru+++04ffvihMjIy9POf/1xjx44N6Pd6+OGHNXXqVLlcLg0fPlylpaXKzs7W999/r7S0tFrVmDZtmoYPH66OHTvq+++/15o1a5SSklLj2N/+9rfq3bu3HnvsMd14443auHGjnn/+eZ87+pyPf//73/rmm2+8z/fu3astW7YoOjparVu3vqDaAACcCU0+gFr78MMPFR8fr/DwcDVr1kzdunXTvHnzNG7cONWrF9h/GLz99tvVqFEjPfnkk7r33nvVuHFjde3aVdOmTat1jfLyck2ePFn/93//J5fLpWHDhunZZ5+tcWyPHj30zjvvaPbs2XrssccUHx+vRx991Oei2/ORnZ2tq6++2vu86g+UcePGVbtzDwAAgWKYNU06BQAAAGBbzMkHAAAAHIYmHwAAAHAYmnwAAADAYWjyAQAAAIehyQcAAAAchiYfAAAAcBiafAAAAMBhaPIBAAAAh6HJBwAAAByGJh8AAABwGJp8AAAAwGH+P0tTeNqWbIbbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "new = num #Points to be added\n",
    "\n",
    "new_num = new+num\n",
    "x_new = sampling.expand_lhs(x,  new, method='ese')\n",
    "\n",
    "\n",
    "### For visualization only\n",
    "intervals = []  \n",
    "subspace_bool = []\n",
    "for i in range(len(xlimits)):  \n",
    "    intervals.append(np.linspace(xlimits[i][0],xlimits[i][1],new_num+1))\n",
    "    subspace_bool.append([[intervals[i][j]<x[kk][i]<intervals[i][j+1] \n",
    "    for kk in range(len(x))] for j in range(len(intervals[i])-1)])\n",
    "    \n",
    "    \n",
    "fig, ax = plt.subplots(1)\n",
    "ax.plot(x[:, PLOT1], x[:, PLOT2], \".\")\n",
    "ax.plot(x_new[len(x):, PLOT1], x_new[len(x):, PLOT2], \".\")\n",
    "for i in range(len(intervals[0])):\n",
    "    ax.plot([intervals[PLOT1][i], intervals[PLOT1][i]], [intervals[PLOT2][0],intervals[PLOT2][-1]], linewidth=0.5,color='k')\n",
    "    ax.plot([intervals[PLOT1][0], intervals[PLOT1][-1]], [intervals[PLOT2][i],intervals[PLOT2][i]], linewidth=0.5,color='k')\n",
    "ax.set_xlabel(\"Dimension \" + str(PLOT1+1))\n",
    "ax.set_ylabel(\"Dimension \" + str(PLOT2+1))\n",
    "\n",
    "i=PLOT1\n",
    "for j in range(len(subspace_bool[i])):\n",
    "    if True in subspace_bool[i][j]:\n",
    "        ax.add_patch(patches.Rectangle((intervals[PLOT1][j], intervals[PLOT2][0]), \n",
    "        (intervals[PLOT1][1]-intervals[PLOT1][0]), (intervals[PLOT2][-1]-intervals[PLOT2][0]),facecolor='lightblue'))\n",
    "\n",
    "i=PLOT2\n",
    "for j in range(len(subspace_bool[i])):\n",
    "    if True in subspace_bool[i][j]:\n",
    "        ax.add_patch(patches.Rectangle((intervals[PLOT1][0], intervals[PLOT2][j]), \n",
    "        (intervals[PLOT1][-1]-intervals[PLOT1][0]), (intervals[PLOT2][1]-intervals[PLOT2][0]),facecolor='lightblue'))\n",
    "       \n",
    "ax.legend(['Initial LHS','Augmented LHS'],bbox_to_anchor = (1.05, 0.6))\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Augmented Sampling for mixed integer variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from smt.utils.kriging import XSpecs\n",
    "from smt.surrogate_models import XType, XRole\n",
    "from smt.applications.mixed_integer import MixedIntegerSamplingMethod\n",
    "from smt.utils.mixed_integer import (\n",
    "    cast_to_mixed_integer, unfold_with_enum_mask, fold_with_enum_index, encode_with_enum_index \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlimits = [\n",
    "    [\"6,7\", \"3,7\", \"4,6\", \"3,4\"], \n",
    "    [0, 1],  \n",
    "    [0, 100],  \n",
    "]\n",
    "xroles = [\n",
    "    XRole.META,\n",
    "    XRole.NEUTRAL,\n",
    "    XRole.NEUTRAL,\n",
    "]\n",
    "xtypes = [\n",
    "    (XType.ENUM, 4),\n",
    "    XType.ORD,\n",
    "    XType.FLOAT,\n",
    "]\n",
    "xspecs = XSpecs(xtypes=xtypes, xlimits=xlimits, xroles=xroles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling = MixedIntegerSamplingMethod(LHS, xspecs, criterion=\"ese\", random_state=42)\n",
    "n_doe = 5\n",
    "Xt = sampling(n_doe)\n",
    "n_add = 10\n",
    "x_new = sampling.expand_lhs(Xt, n_add, method=\"ese\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial DOE \n",
      " [[ 2.          0.         67.32723687]\n",
      " [ 1.          1.         50.49512863]\n",
      " [ 2.          1.         39.39819704]\n",
      " [ 0.          0.         80.92900825]\n",
      " [ 3.          0.          3.11989041]]\n"
     ]
    }
   ],
   "source": [
    "print('Initial DOE \\n', Xt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expand DOE \n",
      " [[ 0.          0.         51.20707681]\n",
      " [ 0.          1.         29.85512515]\n",
      " [ 2.          1.         69.37809697]\n",
      " [ 2.          0.         15.3675912 ]\n",
      " [ 1.          1.         95.12575685]\n",
      " [ 1.          0.         47.66738419]\n",
      " [ 1.          0.         70.928392  ]\n",
      " [ 3.          0.         86.7619275 ]\n",
      " [ 3.          1.          1.61037191]\n",
      " [ 2.          1.         38.59810412]]\n"
     ]
    }
   ],
   "source": [
    "print('Expand DOE \\n', x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial DOE in the initial space (labels) from integer encoding: \n",
      "Doe point  i=0 =['4,6', 0, 67.32723686587384] \n",
      "\n",
      "Doe point  i=1 =['3,7', 1, 50.49512863264476] \n",
      "\n",
      "Doe point  i=2 =['4,6', 1, 39.39819704323989] \n",
      "\n",
      "Doe point  i=3 =['6,7', 0, 80.92900825439996] \n",
      "\n",
      "Doe point  i=4 =['3,4', 0, 3.1198904067240534] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "#to give the initial doe in the initial space from integer encoding\n",
    "print('Initial DOE in the initial space (labels) from integer encoding: ')\n",
    "\n",
    "doe_in_enum_space = []\n",
    "for i in range(n_doe):\n",
    "    doe_in_enum_space.append( cast_to_mixed_integer(xspecs, Xt[i]))\n",
    "    print(\"Doe point  i={} ={}\".format((i),doe_in_enum_space[i]  ),'\\n')   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}