tutorial/SMT_ExpandedLHS.ipynb

{
 "cells": [
  {
   "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:** July 2021"
   ]
  },
  {
   "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": [
    "<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": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x18c183aee20>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAEGCAYAAACjAHa5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhRklEQVR4nO3df3Ad9Xnv8fdHkk1NoLFjK9jBWLYAAzYEwVEcUdJeJdN2gFJQe2lDSn40vR4GkkxLaG9D01xE/uhMcnNv0utAcBiHBhrnR+cm1RhjQklBUNKrkCOQMWCTKC5OhO3YGGNDcGyL89w/zpoqQj+Oj7Tas9LnNfMdn7P73X2eXevo0e7Z3a8iAjMzM8uHuqwTMDMzs8q5cJuZmeWIC7eZmVmOuHCbmZnliAu3mZlZjjRkncDxWrBgQSxdujTrNMzMcqW3t/eFiGjMOg+buNwV7qVLl1IsFrNOw8wsVyTtyDoHmxw+VW5mZpYjLtxmZmY54sJtZmaWIy7cZmZmOeLCbWZmliOpFW5JvybpMUmbJT0t6dMj9JGkNZL6JT0p6cK08jEzM5sO0rwd7DDwnoh4RdIs4FFJ90VEz5A+lwJnJu2dwO3Jv2ZmZjaC1I64o+yV5O2spA0fQ/RK4O6kbw8wV9KitHIyMzPLu1QfwCKpHugFzgBui4gfDOtyKvCzIe8Hkmm7hq3nWuBagCVLllSdz+rVq1m8eHHVy09EX18fLS0tMyZulrG9zdM/bpaxs9zmgYEB1q1bl0lsqyERkXoD5gIPAecOm34v8K4h7/8VKIy1rkKhENXq7OysetmJyiq2t3lmxJ5pcbOMnddtBooxBb/v3dJvU3JVeUS8BHQDlwybNQCcNuT9YmDnVORklje9O/Zz20P99O7Yn3UqZpah1E6VS2oEjkbES5LmAL8NfHZYtw3AxyR9k/JFaQciYhdm9it6d+znmnU9HBksMbuhjvWr2yg0zcs6LTPLQJrfcS8C7kq+564D/ikiNkq6DiAi1gKbgMuAfuBV4MMp5mOWWz3b93FksEQp4OhgiZ7t+1y4zWao1Ap3RDwJXDDC9LVDXgfw0bRyMJsu2prnM7uhjqODJWY11NHWPD/rlMwsI7kb1tNsJio0zWP96jZ6tu+jrXm+j7bNZjAXbrOcKDTNc8E2Mz+r3MzMLE9cuM3MzHLEhdvMzCxHXLjNzMxyxIXbzMwsR1y4zczMcsSF28zMLEdcuM3MzHLEhdvMzCxHXLjNzMxyxIXbzMwsR1y4zczMckTlkTXzo7W1NYrFYlXLdnR00NLSMrkJVai7u5v29vYZEzfL2N7m6R83y9hZbnNfXx9dXV1VLSupNyJaJzcjy0RE5KoVCoWoVmdnZ9XLTlRWsb3NMyP2TIubZey8bjNQjBr4He428eZT5WZmZjniwm1mZpYjLtxmZmY54sJtZmaWIy7cZmZmOeLCba/r3bGf2x7qp3fH/qxTMTOzUTRknYDVht4d+7lmXQ9HBkvMbqhj/eo2Ck3zsk7LzMyG8RG3AdCzfR9HBkuUAo4OlujZvi/rlMzMbASpFW5Jp0l6SNJWSU9L+osR+rRLOiCpL2k3p5WPja2teT6zG+qoF8xqqKOteX7WKZmZ2QjSPFU+CPxlRDwu6WSgV9IDEfHMsH7/FhGXp5iHVaDQNI/1q9vo2b6Ptub5Pk1uZlajUivcEbEL2JW8flnSVuBUYHjhthpRaJrngm1mVuOm5DtuSUuBC4AfjDD7IkmbJd0naeUoy18rqSipuHfv3jRTNTMzq2mpF25JJwHfBm6IiIPDZj8ONEXE+cAXga6R1hERd0REa0S0NjY2ppqvmZlZLUu1cEuaRblor4+I7wyfHxEHI+KV5PUmYJakBWnmZGZmlmdpXlUu4CvA1oj4/Ch9Fib9kLQqycf3IZmZmY0izavKLwY+AGyR1JdM+ySwBCAi1gJXAddLGgQOAVdHRKSYk5mZWa6leVX5o4DG6XMrcGtaOZiZmU03fnKamZlZjihvZ6ZbW1ujWCxWtWxHRwctLS2Tm1CFuru7aW9vnzFxs4ztbc427s6XDjGw/xCL583hbXPnTGnsqZDlz1dfXx9dXV1VLSupNyJaJzcjy0RE5KoVCoWoVmdnZ9XLTlRWsb3NMyN2rcQtPvdinPWpTbHspo1x1qc2RfG5F6cs9lTJ688XUIwa+B3uNvHmU+VmNmk8WI1Z+ly4zWzSeLAas/R5PG4zmzQerMYsfS7cZjapPFiNWbp8qtzMzCxHXLjNzMxyxIXbzMwsR1y4zczMcsSF28zMLEdcuM3MzHLEhdvMzCxHXLjNzMxyxIXbzMwsR1y4zczMcsSF28zMLEdcuM3MzHLEhdvMzCxHXLjNzMxyRBGRdQ7HpbW1NYrFYlXLdnR00NLSMrkJVai7u5v29vYZEzfL2N7m6R83i9g7XzrEwP5D/PzZXn7/kt+ZsrhD9fX10dXVVdWyknojonVyM7JMRESuWqFQiGp1dnZWvexEZRXb2zwzYs+0uFMdu/jci3HWpzbFsps2xvzf/JMoPvfilMUeaiLbDBSjBn6Hu028+VS5mdk4erbv48hgiVLAa6WgZ/u+rFOyGSy1wi3pNEkPSdoq6WlJfzFCH0laI6lf0pOSLkwrHzOzarU1z2d2Qx31gvo60dY8P+uUbAZrSHHdg8BfRsTjkk4GeiU9EBHPDOlzKXBm0t4J3J78a2ZWMwpN81i/uo2e7fvYfGgxhaZ5WadkM1hqR9wRsSsiHk9evwxsBU4d1u1K4O7kK5geYK6kRWnlZGZWrULTPD767jN429w5WadiM9yUfMctaSlwAfCDYbNOBX425P0AbyzuSLpWUlFSce/evanlaWZmVutSL9ySTgK+DdwQEQeHzx5hkTfcnxYRd0REa0S0NjY2ppGmmZlZLqRauCXNoly010fEd0boMgCcNuT9YmBnmjmZmZnlWZpXlQv4CrA1Ij4/SrcNwAeTq8vbgAMRsSutnMzMzPJu1KvKJf068DeUj4Lvi4ivD5n3pYj4yDjrvhj4ALBFUl8y7ZPAEoCIWAtsAi4D+oFXgQ9XtxlmZmYzw1i3g/0D8GPKp7r/TNJ/Bf4kIg4DbeOtOCIeZeTvsIf2CeCjladrZmY2s411qvz0iLgpIroi4grgceBBSX7ygJmZWUbGOuI+QVJdRJQAIuLvJA0AjwAnTUl2ZmZm9ivGOuK+B3jP0AkRcRfwl8CRNJMyMzOzkY16xB0Rfz3K9O9SfkSpmZmZTTGPDmZmZpYjLtxmZmY54sKdgd4d+7ntoX56d+zPOhUzM8sZlW+lHqeT9BvAUoZ8Jx4Rd6eX1uhaW1ujWCxWtWxHRwctLS2Tm1CFuru7aW9vZ+dLh/jO4wO8Vgrq68QfXrg41dGGjsXNQlaxvc3TP26WsbPc5r6+Prq6uqpaVlJvRLRObkaWiYgYswH/CPw78CXgi0lbM95yabVCoRDV6uzsrHrZiToW+9YHfxzLbtoYTZ/YGM03bYxbH/zxlMTNQlaxvc3TP26WsfO6zUAxMvq97Ta5baz7uI9pBVYk//E2QW3N85ndUMfRwRKzGupoa/bzbMzMrHKVFO6ngIWAB/+YBIWmeaxf3UbP9n20Nc+n0DQv65TMzKrW29v71oaGhnXAufi6qclQAp4aHBxcXSgU9ozUoZLCvQB4RtJjwOFjE6P8GFSrQqFpngu2mU0LDQ0N6xYuXHhOY2Pj/rq6Op+ZnaBSqaS9e/eu2L179zpgxDpbSeG+ZVKzMjOz6eRcF+3JU1dXF42NjQd279597mh9xi3cEfGwpFOAdySTHouIEQ/fzcxsxqlz0Z5cyf4c9WuHcb+PkPTHwGPAHwF/DPxA0lWTlqGZmZlVrJILCf4WeEdEfCgiPgisAv5HummZmZlV5sQTT7xgvD7vfe97m3p7e38N4Kabblo4dN4FF1xwdrUxRpp+4403vu3mm28+Zei0U0899bxdu3Y1AHziE59YeMYZZ6xcvnz5irPPPnvFgw8++Kbx4g9VyXfcdcNOje/DVw6amVmOfOtb39px7PWaNWsWfeYzn9l97P0TTzyxbary+N73vvem+++/f+6WLVuemTNnTuzatavh8OHDOp51VFKAvyvpfkl/KulPgXuBTdUkbGZm9mj/C2/67He3LXy0/4XjOtIcz8aNG09etWrVWZdccknzsmXLVl5xxRXLSqUSAKtWrTrrkUceOfEjH/nIqYcPH647++yzV1xxxRXL4D+Pmg8cOFB30UUXLV+xYsU5y5cvX/G1r31t7mTmB/D888/Pestb3jI4Z86cAFi0aNHg0qVLjx7POsYt3BHx34E7gLcD5wN3RMQnqknYzMxmtkf7X3jTn331h8u//PBPTv2zr/5w+WQX761bt8657bbbftbf3//0T3/60xMeeOCBk4bO/9KXvvT8CSecUNq2bdszGzZs+I+h80488cTSvffe2//MM89sffjhh3/0yU9+cvGxwn+81q5de8rZZ5+94ljbs2fPLICOjo6DO3funL106dJz3//+9y+59957TxpvXcNVdMo7Ir4dETdGxMcj4p+PN4iZmRnA9/tfOHnwtVJdKWDwtVLd9/tfOHky13/eeef94vTTTz9aX1/PypUrX/3JT34yu9JlS6WSbrjhhsXLly9f8e53v3v5nj17Zg8MDFTylfIbXHfddT/ftm3bM8faW9/61qMAb37zm0tPPfXUM7feeuuOxsbGwQ996EOnr1mz5rgeoTlq4Zb0aPLvy5IODmkvSzpYzYaYmdnMdvEZC15uqK8r1Qka6utKF5+x4OXJXP8JJ5zw+q1p9fX1DA4OVvz98Ze//OW37Nu3r2HLli1bt23b9sz8+fOPHjp0aNKv6WpoaODyyy9/+Qtf+MLOz33ucz/t6uo6ridyjfqXRES8K/l3Uv8aMjOzmetdZyz4xZ1/+o4ffb//hZMvPmPBy+86Y8EvpjqHhoaGOHz4sIYWeYADBw7UL1iw4OgJJ5wQ99xzz8k7d+6s+Gi9Ups3bz6hrq6O88477zDAE088MWfx4sVHjmcd454CkHQ6MBARhyW1U/6u++6IeOn4UzYzs5nuXWcs+EUWBfuYa665Zu8555yz4txzz3116Pfcq1evfvHSSy8949xzzz1n5cqVry5btuyX463rl7/8Zd0pp5zy9mPvr7/++p+P1f/gwYP1f/7nf77k4MGD9fX19bF06dLDd911146xlhmuknP33wZaJZ0BfAXYAHwduOx4ApmZmaXh1VdffQLg8ssvf/nyyy9//dT73Xff/dNjrx977LFnj72+/fbbnweeH778okWLBvv6+ka8NexYn+FKpVJvJTk+//zzW47FmOjtZ5Wcuy9FxCDwB8DfR8THgUXjLSTpTkl7JD01yvx2SQck9SXt5uNL3czMbOap5Ij7qKT3AR8Cfj+ZNquC5b4K3ArcPUaff4uIyytYl5mZmVHZEfeHgYuAv4uI/5C0DPjaeAtFxCPAixPMz8zMalupVCod15O/bGzJ/hz1BnJFpDeoi6SlwMaIeMPwZMmFbt8GBoCdwF9FxNOjrOda4FqAJUuWFHbsOK7v8V/X0dFBS0tLVctOVHd3N+3t7TMmbpaxvc3TP26WsbPc5r6+Prq6uqpaVlJvRLRObkawefPmDQsXLlzR2Nh4wKOETVwyHvebd+/e/cz5559f3Xjcki6mPCZ3U9JfQERE8wTzexxoiohXJF0GdAFnjtQxIu6g/PQ2Wltbq/7BaGlp4ZZbbql28Qm55ZZbMomdVdwsY3ubp3/cLGNnvc21ZnBwcPXu3bvXJeNHexyLiSsBTw0ODq4erUMl33F/Bfg40Au8NkmJEREHh7zeJOlLkhZExAuTFcPMzNJVKBT2ACMeGVo6KincByLivskOLGkh8POICEmrKP+ltm+y45iZmU0nlRTuhyR9DvgOcPjYxIh4fKyFJH0DaAcWSBoAOkmuRo+ItcBVwPWSBoFDwNWR5hfuZmZm00Alhfudyb9DL2oI4D1jLRQR7xtn/q2UbxczMzOzCo1buCPi3VORiJmZmY1v3CsAJZ0i6SuS7kver5D039JPLXu9O/Zz20P99O7Yn3UqZmZmQGWnyr8K/APwt8n7HwHfony1+bTVu2M/16zr4chgidkNdaxf3Uah6bhGXjMzM5t0ldxztyAi/onkKS7Jc8sn7bawWtWzfR9HBkuUAo4OlujZ7gvezcwse5UU7l9Imk/5gjQktQEHUs2qBrQ1z2d2Qx31glkNdbQ1z886JTMzs4pOld9IeSjP0yV9H2ikfCvXtFZomsf61W30bN9HW/N8nyY3M7OaUMlV5Y9L+i/AWZQfd/psRBxNPbMaUGia54JtZmY1pZJnldcDlwFLk/6/K4mI+HzKuZmZmdkwlZwqvwf4JbCFMYYZMzMzs/RVUrgXR8TbU8/EzMzMxlXJVeX3Sfrd1DMxMzOzcVVyxN0D/LOkOuAo/zke96+nmpmZmZm9QSWF+38DFwFbPHqXmZlZtio5Vf5j4CkXbTMzs+xpvHos6atAM3Afvzoedya3g7W2tkaxWKxq2Y6ODlpaWiY3oQp1d3fT3t4+Y+JmGdvbnJ6dLx1iYP8hFs+bw9vmzvG+nmJ9fX10dXVVtayk3ohoHb+n1byIGLMBnSO18ZZLqxUKhahWZ2dn1ctOVFaxvc0zI/ZUxC0+92Kc9alNseymjXHWpzZF8bkXva9zFBsoRka/t90mt1Xy5LRPp/mHg5nlgwfeMasNoxZuSX8fETdIuodkgJGhIuKKVDMzs5pybOCdo4Ol1wfeuefhrLMym3nGOuL+x+Tf/zUViZhZbRtp4J17sk7KbAYatXBHRG/y78OSGpPXe6cqMTOrPR54xyx7o94OprJbJL0AbAN+JGmvpJunLj0zMzMbaqz7uG8ALgbeERHzI2Ie8E7gYkkfn4rkzMzM7FeNVbg/CLwvIv7j2ISI2A68P5lnZmZmU2yswj0rIl4YPjH5nntWeimZmZnZaMYq3EeqnAeApDsl7ZH01CjzJWmNpH5JT0q6cLx1mpmZzXRjFe7zJR0cob0MnFfBur8KXDLG/EuBM5N2LXB7pUmbmZnNVGPdDlY/kRVHxCOSlo7R5Urg7ogIoEfSXEmLImLXROKamZlNZ5WMDpaWU4GfDXk/kEx7A0nXSipKKu7d61vJzcxs5sqycGuEaSMOVRYRd0REa0S0NjY2ppyWmZlZ7cqycA8Apw15vxjYmVEuZmZmuZBl4d4AfDC5urwNOODvt83MzMY27rCe1ZL0DaAdWCBpgPI43rMAImItsAm4DOgHXgU+nFYuZmZm00VqhTsi3jfO/AA+mlZ8MzOz6SjLU+VmZmZ2nFQ+8M2P1tbWKBaLVS3b0dFBS0vL5CZUoe7ubtrb22dM3Cxje5unf9wsY2e5zX19fXR1dVW1rKTeiGid3IwsExGRq1YoFKJanZ2dVS87UVnF9jbPjNgzLW6WsfO6zUAxauB3uNvEm0+Vm5mZ5YgLt5mZWY64cJuZmeWIC7eZmVmOuHCbmZnliAu3zXi9O/Zz20P99O7Yn3UqZmbjSu3JaWZ50LtjP9es6+HIYInZDXWsX91GoWle1mmZmY3KR9w2o/Vs38eRwRKlgKODJXq278s6JTOzMblw24zW1jyf2Q111AtmNdTR1jw/65TMzMbkU+U2oxWa5rF+dRs92/fR1jzfp8nNrOa5cNuMV2ia54JtZrnhU+VmZmY54sJtZmaWIy7cZmZmOeLCbWZmliMu3GZmZjniwm1mZpYjioisczgura2tUSwWq1q2o6ODlpaWyU2oQt3d3bS3t8+YuFnG9jZP/7hZxs5ym/v6+ujq6qpqWUm9EdE6uRlZJiIiV61QKES1Ojs7q152orKK7W2eGbFnWtwsY+d1m4Fi1MDvcLeJN58qNzMzyxEXbjMzsxxJtXBLukTSs5L6Jd00wvx2SQck9SXt5jTzMTMzy7vUnlUuqR64DfgdYAD4oaQNEfHMsK7/FhGXp5WHmZnZdJLmEfcqoD8itkfEEeCbwJUpxjMzM5v20izcpwI/G/J+IJk23EWSNku6T9LKkVYk6VpJRUnFvXv3ppGrmZlZLqRZuDXCtOE3jT8ONEXE+cAXga6RVhQRd0REa0S0NjY2Tm6WZmZmOZJm4R4AThvyfjGwc2iHiDgYEa8krzcBsyQtSDEnMzOzXEuzcP8QOFPSMkmzgauBDUM7SFooScnrVUk++1LMyczMLNdSu6o8IgYlfQy4H6gH7oyIpyVdl8xfC1wFXC9pEDgEXB0R+XoGq5mZ2RRKrXDD66e/Nw2btnbI61uBW9PMwczMbDrxk9PMzMxyxIXbzMwsR1y4zczMcsSF28wsQ7079nPbQ/307tifdSqWE6lenGZmZqPr3bGfa9b1cGSwxOyGOtavbqPQNC/rtKzG+YjbzCwjPdv3cWSwRCng6GCJnu1+jIWNz4XbzCwjbc3zmd1QR71gVkMdbc3zs07JckB5e95Ja2trFIvFqpbt6OigpaVlchOqUHd3N+3t7TMmbpaxvc3TP26WsSc77s6XDjGw/xCL583hbXPnjNm3r6+Prq6uquJI6o2I1qoWttoSEblqhUIhqtXZ2Vn1shOVVWxv88yIPdPiZhk7r9sMFKMGfoe7Tbz5VLmZmVmOuHCbmZnliAu3TSnfs2pmNjG+j9umjO9ZNTObOB9x25TxPatmZhPnwm1TxvesmplNnE+V25QpNM1j/eo2erbvo615vk+Tm5lVwYXbplShaZ4LtpnZBPhUuZmZWY64cJuZmeWIC7eZmVmOuHCbmZnliAu3mZlZjrhwm5mZ5YgLt5mZWY6kWrglXSLpWUn9km4aYb4krUnmPynpwjTzMTMzy7vUCrekeuA24FJgBfA+SSuGdbsUODNp1wK3p5WPmZnZdJDmEfcqoD8itkfEEeCbwJXD+lwJ3B1lPcBcSYtSzMnMzCzXFBHprFi6CrgkIlYn7z8AvDMiPjakz0bgMxHxaPL+X4FPRERx2LqupXxEzpIlSwo7duyoKqfVq1ezePHiqpadqL6+PlpaWmZM3Cxje5unf9wsY2e5zQMDA6xbt66qZSX1RkTrJKdkGUjzWeUaYdrwvxIq6UNE3AHcAdDa2lr1XxrV/sCbmZnVijRPlQ8Apw15vxjYWUUfMzMzS6RZuH8InClpmaTZwNXAhmF9NgAfTK4ubwMORMSuFHMyMzPLtdROlUfEoKSPAfcD9cCdEfG0pOuS+WuBTcBlQD/wKvDhtPIxMzObDlIdjzsiNlEuzkOnrR3yOoCPppmDmZnZdOInp5mZmeWIC7eZmVmOuHCbmZnliAu3mZlZjqT25LS0SNoLVPfoNFgAvDCJ6UyWWs0Lajc353V8nNfxmY55NUVE42QmY9nIXeGeCEnFWnzkX63mBbWbm/M6Ps7r+Dgvq2U+VW5mZpYjLtxmZmY5MtMK9x1ZJzCKWs0Lajc353V8nNfxcV5Ws2bUd9xmZmZ5N9OOuM3MzHLNhdvMzCxHpmXhlnSJpGcl9Uu6aYT5krQmmf+kpAtrJK92SQck9SXt5inK605JeyQ9Ncr8rPbXeHlN+f6SdJqkhyRtlfS0pL8Yoc+U768K88pif/2apMckbU7y+vQIfbLYX5XklcnnMYldL+kJSRtHmJfJ59FqSERMq0Z5CNGfAM3AbGAzsGJYn8uA+wABbcAPaiSvdmBjBvvst4ALgadGmT/l+6vCvKZ8fwGLgAuT1ycDP6qRn69K8spifwk4KXk9C/gB0FYD+6uSvDL5PCaxbwS+PlL8rD6PbrXTpuMR9yqgPyK2R8QR4JvAlcP6XAncHWU9wFxJi2ogr0xExCPAi2N0yWJ/VZLXlIuIXRHxePL6ZWArcOqwblO+vyrMa8ol++CV5O2spA2/IjaL/VVJXpmQtBj4PWDdKF0y+Txa7ZiOhftU4GdD3g/wxl9glfTJIi+Ai5LTd/dJWplyTpXKYn9VKrP9JWkpcAHlo7WhMt1fY+QFGeyv5LRvH7AHeCAiamJ/VZAXZPPz9ffAXwOlUebX8ufRpsB0LNwaYdrwv6Qr6TPZKon5OOXnCZ8PfBHoSjmnSmWxvyqR2f6SdBLwbeCGiDg4fPYIi0zJ/honr0z2V0S8FhEtwGJglaRzh3XJZH9VkNeU7y9JlwN7IqJ3rG4jTKuFz6NNkelYuAeA04a8XwzsrKLPlOcVEQePnb6LiE3ALEkLUs6rElnsr3Fltb8kzaJcHNdHxHdG6JLJ/hovr6x/viLiJaAbuGTYrEx/vkbLK6P9dTFwhaTnKH+d9h5JXxvWpyY/jzZ1pmPh/iFwpqRlkmYDVwMbhvXZAHwwuTqzDTgQEbuyzkvSQklKXq+i/P+zL+W8KpHF/hpXFvsrifcVYGtEfH6UblO+vyrJK6P91ShpbvJ6DvDbwLZh3bLYX+PmlcX+ioi/iYjFEbGU8u+IByPi/cO61eTn0aZOQ9YJTLaIGJT0MeB+yldy3xkRT0u6Lpm/FthE+crMfuBV4MM1ktdVwPWSBoFDwNURkfopMEnfoHwF7QJJA0An5Yt1MttfFeaVxf66GPgAsCX5fhTgk8CSIXllsb8qySuL/bUIuEtSPeXC908RsTHrz2OFeWXyeRxJDewvqyF+5KmZmVmOTMdT5WZmZtOWC7eZmVmOuHCbmZnliAu3mZlZjrhwm5mZ5YgLt007kl5TeTSnp5PHVd4oqS6Z1yppTUZ5/fskreePkm0rSWqdjHWaWX74djCbdiS9EhEnJa/fSnmUpe9HRGe2mU0OSedQfo71l4G/iohiximZ2RTyEbdNaxGxB7gW+FjypKl2JWMcS7pF0l2S/kXSc5L+UNL/lLRF0neTR4giqSDpYUm9ku5XMhKTpG5Jn1V5XOcfSfrNZPrKZFqfyuMln5lMfyX5V5I+J+mpJNZ7k+ntyTr/r6RtktYfe3LXsG3aGhHPTsX+M7Pa48Jt015EbKf8s/7WEWafTnkIxSuBrwEPRcR5lJ+U9XtJ8f4icFVEFIA7gb8bsnxDRKwCbqD8ZDeA64D/kwxg0Ur52dJD/SHQApxP+VGbn9N/Dst4QbKuFZTHbr+4mm02s+lr2j3y1GwUI42oBHBfRByVtIXyo2i/m0zfAiwFzgLOBR5IDn7rgaHPhT42mEdv0h/g/wF/q/K4yt+JiB8Pi/ku4BsR8Rrwc0kPA+8ADgKPRcQAQPLo0qXAo8e5rWY2jfmI26Y9Sc3Aa5THXR7uMEBElICjQ55FXaL8h62ApyOiJWnnRcTvDl8+WX9Dsq6vA1dQPmq/X9J7hqc0RrqHh7x+fZ1mZse4cNu0JqkRWAvcWuUAEc8CjZIuStY3S9LKcWI2A9sjYg3lkZzePqzLI8B7JdUn+f0W8FgVuZnZDOTCbdPRnGO3gwHfA/4F+HQ1K4qII5RHifqspM1AH/Ab4yz2XuCp5FT32cDdw+b/M/AksBl4EPjriNhdaU6S/iAZLe0i4F5J91e6rJnln28HMzMzyxEfcZuZmeWIC7eZmVmOuHCbmZnliAu3mZlZjrhwm5mZ5YgLt5mZWY64cJuZmeXI/wcO5mzd95fQmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAEGCAYAAAAJ9c30AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt2klEQVR4nO3de3iU9Zk38O89CeeoGJMQhOAgAQMBEYjxAOwCtS0q6tazbZW663LBq69a9bWnXZ16XX1rX9RadOWwngAtQm2lSBUrq9Qg1ZgIhEA4BAhLgBACkgOHhMnc7x/zhI1xkpnMM8/88mS+n+uaa07Pb+57fp3Um+dw/0RVQURERNQRj+kEiIiIqOtjwUBERERhsWAgIiKisFgwEBERUVgsGIiIiCisZNMJdFZaWpp6vV7TaRARuUpxcXGNqqabzoPcy3UFg9frRVFRkek0iIhcRUT2mc6B3I2HJIiIiCgsFgxEREQUFgsGIiIiCosFAxEREYXFgoGIiIjCcqxgEJHeIlIoIptFZKuI/DLENiIi80SkXERKRGS8U/kQERFR9Jy8rLIRwDRVbRCRHgDWi8j7qvpZq22uBTDcul0BYL51T0RERF2IY3sYNKjBetrDurVdS/smAEusbT8D0F9EBjqVExEREUXH0cZNIpIEoBhANoD/UNXP22wyCMD+Vs8rrdcOtfmcWQBmAcCQIUOizudbt34faZnR1yN7y7Zi6MjchBlrMrYbx5qM7caxJmO79Tufc7oWL7/8ctSxiexwtGBQ1WYAl4lIfwDviMhoVS1ttYmEGhbicxYBWAQAeXl533g/UmmZA3HH/34s2uFY/sIzUY9341iTsd041mRsN441Gdut37lk2cKo4xLZFZerJFT1OIB1AKa3easSQFar54MBHIxHTkRuU1HdgLUlVaiobgi/MRFRjDl5lUS6tWcBItIHwDUAtrfZbBWAe6yrJa4EUKuqh0BEX1NR3YD5a3bh/S8PYv6aXSwaiCjunDwkMRDAYus8Bg+AFaq6WkRmA4CqLgDwHoDrAJQDOAngXgfzIXKt8qoG+AMKVcAfUJRXNcCbkWI6LSJKII4VDKpaAmBciNcXtHqsAO53Kgei7iI7MwXJHoE/oEj2CLIzWSwQUXy5bnlrokTkzUjBnOnDUV7VgOzMFO5dIKK4Y8FA5BLeDBYKRGQO15IgIiKisFgwEBERUVgJdUhib9lWLH/hmajHHy4tjrpxSmnhBiDK2Kbimoztxrh2x7txvkz+vtw413ZjN1bujjoukV0JVTAMHZlrq0NbybKF8Pl8UY0tq6m31d3NRFyTsd0YFwB8Pp/r8rb1G7HZ9dBO3m6c61jEJjKFhySIiIgoLBYMREREFBYLBqJEsb8QKHg2eE9E1EkJdQ4DUcLaXwgsvhFobgKSegIzVwFZ+aazIiIX4R4GokRQURAsFrQ5eF9RYDojInIZFgxEicA7ObhnQZKC997JpjMiIpfhIQmiRJCVHzwMUVEQLBZ4OIKIOokFA1GiyMpnoUBEUUuogsFup8fSwg0oq6mPeqydbnom4pqM7ca4dse7db5M/r7cNtd2Y7PTI5mUUAWD3U6PtrrauXGsydhuHGsythvHmozt0u/MTo9kEk96JCIiorBYMBAREVFYLBiIiIgoLBYMREREFBYLBiIiIgqLBQOdVVHdgLUlVaiobjCdChERdTEJdVklta+iugHz1+yCP6BI9gjmTB8Ob0aK6bSIiKiL4B6GCKXWbMQFR4qRWrPRdCqOKK9qgD+gUAX8AUV5FfcyEBHR/3BsD4OIZAFYAiATQADAIlX9XZttpgD4M4C91kt/UtWnnMop2k6PfU5WYcjeldi2txHr/r0IGHsXcN6gTn2Gne5uh0uLo27YEmnc4yeaUFt+DAFVeESwuyoVR//WMy6xQ3FjXLvj3Thfdjs92snbjXNtNzY7PZJJTh6S8AN4VFW/FJFzABSLyIequq3NdgWqOsPBPM6KttPjiK0LMaqkJ55SP3z/2BOYlgNMfrRTn1FWU2+ru5vP54tqbGfifre6AeVVDcjOTDl7OCJesdtyY1wA8Pl8rsvbTly7XQ/t5O3GuY5FbCJTHCsYVPUQgEPW43oRKQMwCEDbgqHLqxmQj0BSTwCnu/XSwN6MFJ63QEREIcXlHAYR8QIYB+DzEG9fJSKbReR9EcltZ/wsESkSkaIjR444mWpIx9LGYf2014Chk4NLBHPFPyIiSjCOFwwikgLgjwAeVtW6Nm9/CeAiVR0L4AUAK0N9hqouUtU8Vc1LT093NN/2HEsbB1x0NYsFIiJKSI4WDCLSA8Fi4U1V/VPb91W1TlUbrMfvAeghImlO5kRERESd51jBICIC4BUAZar6XDvbZFrbQUTyrXyOOpUTERERRcfJqyQmArgbwBYR2WS99nMAQwBAVRcAuBXAHBHxAzgF4E5VVQdzIiIioig4eZXEegASZpsXAbzoVA5EREQUG+z0SERERGEl1FoS0XZ6bFFauAFlNfVRj7XTTc9EXJOx3RjX7ni3zlfbscdPNOFYQxNSU3qif7+eYcfbydttc203Njs9kkkJVTBE2+nxLDtd7dw41mRsN441GbuLjG29iNnxSBYx6yJ5x22szfHs9Egm8ZAEEcUMFzEj6r5YMBBRzGRnpiDZIxABkj2C7Ey2GifqLhLqkAQROcubkYI504d/YxEzInI/FgxEFFNcxIyoe+IhCSIiIgqLBQMRERGFxYKBiIiIwmLBQERERGEl1EmPdjs9Hi4tjrpxip3ubqbimoztxrh2x7txvkz+vtw413Zjs9MjmZRQBYPdTo8lyxbC5/NFNbaspt5WdzcTcU3GdmNcAPD5fK7L29ZvxGbXQzt5u3GuYxGbyBQekiB32F8IFDwbvCciorhLqD0M5FL7C4HFNwLNTUBST2DmKiAr33RWREQJhXsYqOurKAgWC9ocvK8oMJ0REVHCYcFAXZ93cnDPgiQF772TTWdERJRweEiCur6s/OBhiIqCYLHAwxFERHHHgoHcISufhQIRkUE8JEFERERhsWAgIiKisBLqkITdTo+lhRtQVlMf9Vg73fRMxDUZ241x7Y5363yZ/H3F8zsfP9GEYw1NOLZro+3fCDs9khslVMFgt9Ojra52bhxrMrYbx5qM7caxJmN3cmxFdQPmr9kFf0BRW9WAK+6YDW9GSlxit8ZOj2QSD0kQEYVRXtUAf0ChCgRUUV7VYDolorhzrGAQkSwR+VhEykRkq4g8FGIbEZF5IlIuIiUiMt6pfIiIopWdmYJkj0AE8IggOzPKvQtELubkIQk/gEdV9UsROQdAsYh8qKrbWm1zLYDh1u0KAPOteyKiLsObkYI504ejvKoBu6tSoz8cQeRiju1hUNVDqvql9bgeQBmAQW02uwnAEg36DEB/ERnoVE5ERNHyZqTgmksz0b9fT9OpEBkRl3MYRMQLYByAz9u8NQjA/lbPK/HNogIiMktEikSk6MiRI47lSURERKE5XjCISAqAPwJ4WFXr2r4dYoh+4wXVRaqap6p56enpTqRJREREHXC0YBCRHggWC2+q6p9CbFIJIKvV88EADjqZExEREXWeYyc9iogAeAVAmao+185mqwA8ICJvIXiyY62qHnIqJ7uNmw6XFkd9HbSdZi2m4pqM7ca4dse7cb5M/r7cONd2Y7NxE5nUbsEgIucC+BmC/+p/X1V/3+q9l1T1f4X57IkA7gawRUQ2Wa/9HMAQAFDVBQDeA3AdgHIAJwHcG93XiIzdxk0lyxbC5/NFNbaspt5WsxYTcU3GdmNcAPD5fK7L29ZvxGbjJjt5u3GuYxGbyJSO9jC8BmAXgocU/llEbgHwfVVtBHBluA9W1fUIfY5C620UwP2Rp0tkyP5CLq9NRAmto4JhmKreYj1eKSK/APCRiNwYh7yIuo79hcDiG4HmJiCpJzBzFYsGIko4HZ302EtEzr6vqr8CsAjAJwAucDoxoi6joiBYLGhz8L6iwHRGRERx11HB8C6Aaa1fUNXFAB4F0ORkUkRdindycM+CJAXvvZNNZ0REFHftHpJQ1cfbeX0Ngq2ciRJDVn7wMATPYSCiBJZQy1sTRS0rn4UCESU0Lm9NREREYbFgMKCiugFrS6pQUd1gOhUiIqKIRHRIQkSuBuBtvb2qLnEoJ8fY7fRYWrgBZTX1UY/FC8/g+IkmfFF+DAFVeERweXZq2NXvYhE3WqZiuzGu3fFunS+Tvy+3zbXd2Oz0SCaFLRhEZCmAYQA2AWi2XlYArisY7HZ6tNXVzhq7tqQKO788CFVABBg2/kJcc2mm43GjZiq2G8eajO3GsSZju/Q7s9MjmRTJHoY8AKOsroxkU3ZmCpI9An9AkewRZGemmE6JiIgorEgKhlIAmQAcWxQqkXgzUjBn+nCUVzUgOzMF3gwWDETkXsXFxRnJyckvAxgNnhfnZgEApX6//74JEyZUh9ogkoIhDcA2ESkE0NjyoqqyRXSUvBksFIioe0hOTn45MzNzZHp6+lcej4d7ol0qEAjIkSNHRlVVVb0MIOR/3yMpGHwxzYqIiLqT0SwW3M/j8Wh6enptVVXV6Pa2CVswqOrfRGQAgMutlwpVNeTuCiIiSjgeFgvdg/W/Y7uHlcIebxKR2wEUArgNwO0APheRW2OWIREREXV5kZyg8gsAl6vqTFW9B0A+gH93Ni0iIqLI9O3bd1y4be64446LiouLewPAT3/6069dyz5u3LicaGOEev2RRx658IknnhjQ+rVBgwaNOXToUDIA/OQnP8nMzs7OHTFixKicnJxRH330Ub9w8buCSM5h8LQ5BHEUPBOWiIhcZPny5ftaHs+bN2/g008/XdXyfOPGjdvjlcfatWv7ffDBB/23bNmyrU+fPnro0KHkxsZGiVd8OyIpGNaIyAcAllnP7wDwnnMpOcdup8fDpcVRN06x093NVFyTsd0Y1+54N86Xyd+XG+fabuzu0ulxfXlNv0/La86ZmJ1WPyk77USsPnf16tXnPPXUUxempqae2bFjR58xY8acXLly5V6Px4P8/PxLnnnmmf1vvfXW+Y2NjZ6cnJxRI0aMOLVq1aq9ffv2HXfy5MmNtbW1nunTp2fX1tYm+f1+eeKJJw7+8Ic/PB6r/ADgwIEDPVJTU/19+vRRABg4cKA/lp/vpEhOevw/InILgIkABMAiVX3H8cwcYLfTY8myhfD5fFGNLaupt9XdzURck7HdGBcAfD6f6/K29Rux2fXQTt5unOtYxHa79eU1/f759S9G+JsDnlfW7w28+qPLd8ayaCgrK+uzadOmPV6v98yECRNyPvzww5Tvfve7Zxfueemllw68/vrrGdu3b9/Wdmzfvn0Df/nLX8pTU1MDhw4dSr7iiityvv/97x/3eDq/U33BggUDVqxYcUHL8+rq6h4A8E//9E91v/71ry/0er2jJ02aVHfXXXcdu/76612xsFBEs6Cqf1TVR1T1x24tFoiIyLxPy2vO8TcHPAEF/M0Bz6flNefE8vPHjBlzYtiwYWeSkpKQm5t7cvfu3R0v1tNKIBCQhx9+ePCIESNGTZ06dUR1dXXPysrKiNZcamv27NmHt2/fvq3llpGRcQYAzjvvvEBpaem2F198cV96erp/5syZw+bNm3dBuM/rCtotGERkvXVfLyJ1rW71IlIXvxSJiKi7mJidVp+c5Al4BEhO8gQmZqdFv5JXCL169Tp7iWdSUhL8fn/E5wcsXLgw9ejRo8lbtmwp2759+7YLLrjgzKlTp2J+zl5ycjJmzJhR/9vf/vbg3Llz/3vlypXnxzqGE9qtnFR1knUf0+qPiIgS16TstBOv/ujynU6cwxCp5ORkbWxslNbFBQDU1tYmpaWlnenVq5e+++675xw8eDDivROR2rx5cy+Px4MxY8Y0AsDGjRv7DB48uCnWcZwQyWqVwwBUqmqjiEwBcCmAJap63NnUiIhiaH8hUFEAeCebziThTcpOO2GiUGjxgx/84MjIkSNHjR49+uSqVav2trx+3333Hbv22muzR48ePTI3N/fk0KFDT4f7rNOnT3sGDBhwacvzOXPmHO5o+7q6uqQHH3xwSF1dXVJSUpJ6vd7GxYsX7+toTFcRybGZPwLIE5FsAK8AWAXg9wCuczIxIqKY2V8ILL4RaG4Cknqiz8lvmc6IYujkyZMbAWDGjBn1M2bMOHuIY8mSJf/d8riwsHBHy+P58+cfAHCg7fiBAwf6N23aFPISy5Zt2goEAsWR5HjgwIEtLTHieRlnLEVybCagqn4A3wPwvKr+GMDAcINE5FURqRaR0nbenyIitSKyybo90bnUiYgiVFEQLBa0GWhuQt8TB8KPIaKviWQPwxkRuQvATAA3WK/1iGDc6wBeBLCkg20KVHVGBJ9FRBQ972QgqefZPQwn+w0ynRGR60Syh+FeAFcB+JWq7hWRoQDeCDdIVT8BcMxmfkRE9mXlAzNXAdN+AcxchVN9M8OPIaKviaRx0zYAD7Z6vhfA0zGKf5WIbAZwEMBjqro11EYiMgvALAAYMmRI1MHsdnosLdyAsprorgCy203PRFyTsd0Y1+54t86Xyd9XdGOXG/+NJHqnR3KnSK6SmAjAB+Aia3sBoKp6sc3YXwK4SFUbROQ6ACsBDA+1oaouArAIAPLy8qJeRtVup0dbXe3cONZkbDeONRnbjWNNxnbpd+4OnR7JvSI5h+EVAD8GUAygOVaBVbWu1eP3ROQlEUlT1ZpYxSAiIqLYiOQchlpVfV9Vq1X1aMvNbmARyRQRsR7nW7nY/lwiIko8S5Ys6S8iEzZu3NjbdC6hbNiwoc/y5cvP6+y4/Pz8Sz755JO+kby+evXqc6ZOnZrd+rVbbrnF+9prr50PAMuWLTtv5MiRoy655JJRw4YNy507d25aZ3KJZA/DxyIyF8CfADS2vKiqX3Y0SESWAZgCIE1EKgE8CevqClVdAOBWAHNExA/gFIA7VTXqww1ERJS43nrrrdTx48c3LF26NHXcuHEHTefTVlFRUd+ioqJ+d9xxR62J+I2NjfLQQw9d9Pe//71s2LBhZ06dOiU7d+7sVCfLSPYwXAEgD8D/BfCsdQt7xo6q3qWqA1W1h6oOVtVXVHWBVSxAVV9U1VxVHauqV6rqhs4kTkRELrVnXT+s9WViz7p+sfi42tpaT1FRUcprr71W8c4775xdl6Htv7jvueeeIS0LPS1fvvy8oUOH5k6YMOGSH/3oR1kt2z3yyCMX3nzzzd6JEycOHzRo0JjFixf3nz179uARI0aMmjx58vDGxkYBgIKCgr6XX375Jbm5uSMnTZo0fN++fT2A4L/858yZM2jMmDEjvV7v6DVr1qScPn1afv3rX1/47rvvnp+TkzPqP//zP8+vq6vz3Hbbbd7Ro0ePHDly5Kg33nijPwA0NDTIjBkzLh4xYsSo66+//uLTp09HvBZGR44fP+7x+/0yYMAAPwD06dNHx44d2xhuXGthCwZVnRriNi3apImIKIHtWdcPb94+Ap/+bhDevH1ELIqGN998s/+UKVNqL7300sb+/fs3r1+//hu78Fs7efKkPPTQQxe9//77u4qLi3ccPXr0a3vb9+3b1+ujjz4qf/vtt8tnz549dNq0aXU7d+7c1rt378CKFSvOa2xslAcffHDIn//8591bt24tmzlzZs1jjz12trmH3++XLVu2lP3mN7/Z/9RTT13Yu3dv/dnPfnbwhhtu+Gr79u3b/vVf//Wrn//85wOnTp1aV1paWlZQULDj3/7t3wbX1dV5nnnmmYw+ffoEdu7cue2JJ544tG3btk7NT1FRUUpOTs6oltvatWv7A8CAAQOav/3tbx8fMmTIpTfccMPQ+fPnpzY3d+60xLAFg4gMEJFXROR96/koEfmXTkVxqYrqBqwtqUJFtSuWKici6vr2rDsHgTMeaAAInPFgzzrbCxyuWLEi9a677voKAG655ZZjS5cuTe1o+02bNvXOyspqzMnJaQKAO++882s9g6655praXr16aX5+/qnm5ma59dZb6wAgNzf31N69e3uWlJT02rVrV59p06aNyMnJGTV37tyBBw8ePNvQ8LbbbvsKAK6++uoTlZWVIXf7r1u37tzf/va3A3NyckZNmjTpksbGRikvL++5fv36lLvvvvsoAFxxxRWnRowYcbIzc5GXl9fQelnta6655njLe8uXL9+3Zs2anXl5eSfmzZuXefvtt3s789mRnMPwOoDXAPzCer4TwHIEr57otiqqGzB/zS74A4pkj2DO9JBXfBJRHKTWbETa4ULUDMjHsbRxptMhOy6eUo+/vxRA4IwHnh4BXDzF1vLWVVVVSZ999tm5O3fu7PPAAw+gublZRETnz59f2aNHDw0EAme3bTmcEO50uZZVLJOSkpCcnKweT/Df1h6PB36/X1RVsrOzT7W37kTv3r0VCC5j3dzcHPKQgqri7bffLg91WMC6HsAR+fn5p/Lz80/NmjXrWHZ29hgAFZGOjeQchjRVXQEgAADWuhIxu7yyqyqvaoA/oFAF/AFFeRX3MhCZkFqzEZM+uhejSn6HSR/di9SakGsAkVtcPOUEfrBiJyY+dAA/WLETF0+xtWrl0qVLz7/55puPHjx4cMuBAwe2VFVVlQwePLjpr3/9a8qwYcMay8vL+5w6dUqOHj2atH79+nMBYOzYsaf379/fa8eOHT0BYPny5R3ukWjr0ksvPX3s2LHktWvX9gOChUhRUVGHV2ece+65zQ0NDWf/mzt16tS6Z599dkBLQfPpp5/2AYBJkyY1vPHGG6kA8MUXX/TeuXNnh4dXIlVbW+tZvXr12b05n3/+eZ8LL7ywU8tqR7KH4YSIXABAAUBErgRg5CxPuzrT6fH4iSbUlh9DQBUeEeyuSkXj3i1RN06x093tcGmxkbgmY7sxrt3xbpyvuPy+9m3Aur11CP6b5TSw53ngoqtdOdd2Y3ebTo8XTzlht1Bo8Yc//OGCxx9//FDr12666aavli5dmvrmm2/+9w033PDVyJEjc4cOHXo6Nzf3JACkpKToc889t2/69OnDU1NT/ePGjetULr1799a33npr94MPPjikvr4+qbm5WebMmXM4Ly+v3eWwr7322vpnnnlmYE5OzqhHH3300NNPP31w1qxZQ3JyckapqgwePLjx448/Ln/ssceq77zzzqEjRowYlZube3LMmDHt5va9731veHJysgLA+PHjG+6///4j7W0bCAQwd+7cAQ888MBFvXv3DvTt2zfwyiuv7G1v+1Ak3K4ZERkP4AUAowGUAkgHcKuqlnQmUKzk5eVpUVFRVGPveODRTnVYq6huQHlVA7IzU+DNSEHJsoXw+Xxxid2aqbgmY7sxLgD4fD7X5W0n7nKbXQ8jyrvN0tSYuQrIynflXJuMLSLFqpoX1eAObN68uWLs2LGua7hXW1vrOe+88wKBQAD33HPPkOHDh59+8sknq03nZdrmzZvTxo4d6w31XiRrSXwpIv8I4BIE20LvUNUzsU2xa/JmBAsFIjKoZeGoioLgqpNZ+aYzom7g+eefT1u2bFnamTNnJDc39+QjjzziuqIn3iJZSyIJwHUAvNb23xERqOpzDudGRBSUlc9CgWLqySefrOYehc6J5ByGdwGcBrAF1omPRERElkAgEBCPx8NOvS4XCAQEHfx3PpKCYbCqXhq7lIiIqBspPXLkyKj09PRaFg3uFQgE5MiRI+cheK5iSJEUDO+LyHdU9a+xS42IiLoDv99/X1VV1ctVVVWjEdml+tQ1BQCU+v3++9rbIJKC4TMA74iIB8AZBE98VFU9NzY5EhGRW02YMKEawI2m8yDnRVIwPAvgKgBbuJokERFRYopk99EuAKUsFoiIiBJXJHsYDgFYZy0+dbbntRsvq+xMp8dQSgs3oKwmurbndrvpmYhrMrYb49od74b5On6iCccampCa0hP9+/U0/vty21zbjd1tOj2SK0VSMOy1bj2tm2sNHZlrqysd7HS1c+NYk7HdONZk7DiMbb0g2/FWC7JxvuI33k5LaiK7Iun0+Mt4JEJEXRsXZCNKbO0WDCLyvKo+LCLvwlp4qjVV5VmxRAkkOzMFyR45u+R7dmYKjppOiojipqM9DEut++gPUBJRt+HNSMGc6cO/tiDb56aTIqK4abdgUNVi6/5vIpJuPW536Uwi6v64IBtR4mr3skoJ8olIDYDtAHaKyBEReSJ+6REREVFX0FEfhocBTARwuapeoKrnA7gCwEQR+XE8kiMiIqKuoaOC4R4Ad6nq3pYXVHUPgB9a7xEREVGC6Khg6KGqNW1ftM5j6OFcSkRERNTVdHSVRFOU7wEARORVADMAVKvq6BDvC4DfAbgOwEkAP1LVL8N9rh12Oz0eLi2OunGKne5upuJ+LXbtAeD4PqD/RcB5gxyP7ca5tjvejfMVs99XnMea/I2w0yO5VUcFw1gRqQvxugDoHcFnvw7gRQBL2nn/WgDDrdsVAOZb946x2+mxZNlC+Hy+qMaW1dTb6u5mIu7Z2P9yHbD4RqB/E5BUAsxcBWTlOxrbjXMNAD6fz3V52/qN2Ox6aCdvN851LGITmdLuIQlVTVLVc0PczlHVsIckVPUTAMc62OQmAEs06DMA/UVkYOe/AjmuogBobgK0OXhfUWA6IyIiirNIVqt0yiAA+1s9r7Re+wYRmSUiRSJSdOQIW0HEnXcykNQTkKTgvXey6YyIiCjOIll8yikS4rWQS2ir6iIAiwAgLy+Py2zHW1Z+8DBERUGwWIjgcAQREXUvJguGSgBZrZ4PBnDQUC4UTlY+CwUiogRm8pDEKgD3WB0lrwRQq6qHDOZDRERE7XBsD4OILAMwBUCaiFQCeBJW/wZVXQDgPQQvqSxH8LLKe53KhYiIiOxxrGBQ1bvCvK8A7ncqPhEREcWOyUMSRERE5BImT3qMO7udHksLN6Cspj7qsXa66ZmIazK2G+PaHe/W+TL5+3LbXNuNzU6PZFJCFQx2Oz3a6mrnxrEmY7txrMnYbhxrMrZLvzM7PZJJPCRBREREYbFgICIiorBYMBAREVFYLBiIiIgoLBYMREREFFZCXSVBFEpFdQPKqxqQnZliOhUioi6LBQMltIrqBsxfswv+gCLZI7j4RJPplIiIuqSEKhjsNm46XFoc9XXQdpq1mIprMna84u453IAjh1o14anfbeta9+4+X7EcC9jL241zbTc2GzeRSQlVMNht3FSybCF8Pl9UY8tq6m01azER12TseMVtu4dh6om/RR0XAHw+X7eer2+w2cTITt5unOtYxCYyJaEKBqK2vBkpmDN9+NlzGOr+q9B0SkREXRILBkp43owUeDOCJzyWGM6FiKir4mWVREREFBYLBiIiIgqLhySIiOJhfyFQUYA+J6tMZ0IUFRYMRERO218ILL4RaG7CkL2NSK35AY6ljTOdFVGn8JAEEZHTKgqA5iZAmyEaQNphXo1D7sOCgYjIad7JQFJPQJKg4kHNgHzTGRF1mqiq6Rw6JS8vT4uKiqIam/+t6Rg6Mjfq2IdLizFlypSoxv7hvQ8wOv9qV8U1GduNcQFg3bp1rsvbTtzSwg3Gfl+um+vaA8DxfVi56RAumfidqGI3Vu7GypUroxorIsWqmhfVYCIk2DkM7PTonthujAu4s/sgOz12ju1Ojw88yk6P5Eo8JEFERERhsWAgIiKisBwtGERkuojsEJFyEflpiPeniEitiGyybk84mQ8RERFFx7FzGEQkCcB/APg2gEoAX4jIKlXd1mbTAlWd4VQeREREZJ+TexjyAZSr6h5VbQLwFoCbHIxHREREDnGyYBgEYH+r55XWa21dJSKbReR9EQl5zaOIzBKRIhEpOnLkiBO5EhERUQecLBgkxGttmz58CeAiVR0L4AUAK0N9kKouUtU8Vc1LT0+PbZZusb8QKHg2eE9ERBRnTvZhqASQ1er5YAAHW2+gqnWtHr8nIi+JSJqq1jiYl/vUHjjbhx5JPYGZq4AsdoojIqL4cbJg+ALAcBEZCuAAgDsBfL/1BiKSCeCwqqqI5CO4x+OoUwntLduK5S88E/X40sINKKupj3osooy95+MPoenHIVAoTuHIlsdxNH2C43Fbxpv4zm6Ma3e8W+fL5O/LbXNtN3Zj5e6o4xLZ5VjBoKp+EXkAwAcAkgC8qqpbRWS29f4CALcCmCMifgCnANypDvaqttvp0VZXOxtjV52swr8PXQtP4AwCnh5YP+3xyFe6s9mJz9R3duVYk7HdONZkbJd+Z3Z6JJMcbQ2tqu8BeK/NawtaPX4RwItO5tAdnOqbifXTXkPa4ULUDMjnsrhERBR3CbWWhJsdSxvHQoGIiIxha2giIiIKiwUDERERhcWCgYjIoIrqBqwtqUJFdYPpVIg6xHMYiIgMqahuwPw1u+APKJI9gjnTh8ObkWI6LaKQuIeBiMiQ8qoG+AMKVcAfUJRXcS8DdV0sGIiIDMnOTEGyRyACJHsE2Zncu0BdlzjYJ8kReXl5WlRUFNXY/G9Nx9CRIde3ikhp4QaMzr86YcaajO3GsSZju3Gsydhd6TsfP9GEYw1NSE3pif79enY4trFyN1auXBlVXBEpVtW8qAYTIcHOYXBrp0e3dqVLuLEmY7txrMnYLv3O7PRIJvGQBBEREYXFgoGIiIjCYsFAccVrzomI3CmhzmEgs3jNORGRe3EPA8UNrzknInIvFgwUN7zmnIjIvXhIguLGm5GCOdOHo7yqAdmZKTwcQUTkIiwYKK68GSwUiIjcKKEKhr1lW7H8hWeiHn+4tDjqximlhRuAKGObimsythvj2h3vxvky+fty41zbjd1YuTvquER2JVTBYLfTY8myhfD5fFGNLaupt9XdzURck7HdGBcAfD6f6/K29Rux2fXQTt5unOtYxCYyhSc9EhERUVgsGIiIiCgsFgxEREQUFgsGIiIiCosFA1G87C8ECp4N3hMRuYyjBYOITBeRHSJSLiI/DfG+iMg86/0SERnvZD5ExuwvBBbfCHz0q+A9iwYichnHCgYRSQLwHwCuBTAKwF0iMqrNZtcCGG7dZgGY71Q+REZVFADNTYA2B+8rCkxnRETUKU7uYcgHUK6qe1S1CcBbAG5qs81NAJZo0GcA+ovIQAdzIjLDOxlI6glIUvDeO9l0RkREnSKq6swHi9wKYLqq3mc9vxvAFar6QKttVgN4WlXXW8//C8BPVLWozWfNQnAPBIYMGTJh3759UeV03333YfDgwVGNBYBNmzbhsssuS5ixJmO7cWzY8bUHgOP7gP4XAecNimlsN441Gdut37myshIvv/xyVGNFpFhV86IaTARnOz1KiNfaVieRbANVXQRgEQDk5eVFXeFE+4dGRESU6Jw8JFEJIKvV88EADkaxDRERERnmZMHwBYDhIjJURHoCuBPAqjbbrAJwj3W1xJUAalX1kIM5ERERURQcOyShqn4ReQDABwCSALyqqltFZLb1/gIA7wG4DkA5gJMA7nUqHyIiIoqeo6tVqup7CBYFrV9b0OqxArjfyRyIiIjIPnZ6JCIiorBYMBAREVFYLBiIiIgoLBYMREREFJZjnR6dIiJHAETX6hFIA1ATw3RipavmBXTd3JhX5zCvzumOeV2kqumxTIYSi+sKBjtEpKgrtkbtqnkBXTc35tU5zKtzmBfRN/GQBBEREYXFgoGIiIjCSrSCYZHpBNrRVfMCum5uzKtzmFfnMC+iNhLqHAYiIiKKTqLtYSAiIqIosGAgIiKisLplwSAi00Vkh4iUi8hPQ7wvIjLPer9ERMZ3kbymiEitiGyybk/EKa9XRaRaRErbed/UfIXLK+7zJSJZIvKxiJSJyFYReSjENnGfrwjzMjFfvUWkUEQ2W3n9MsQ2JuYrkryM/D1asZNEZKOIrA7xnpG/RyKoare6IbiU9m4AFwPoCWAzgFFttrkOwPsABMCVAD7vInlNAbDawJz9A4DxAErbeT/u8xVhXnGfLwADAYy3Hp8DYGcX+X1FkpeJ+RIAKdbjHgA+B3BlF5ivSPIy8vdoxX4EwO9DxTf198gbb91xD0M+gHJV3aOqTQDeAnBTm21uArBEgz4D0F9EBnaBvIxQ1U8AHOtgExPzFUlecaeqh1T1S+txPYAyAIPabBb3+Yowr7iz5qDBetrDurU909rEfEWSlxEiMhjA9QBebmcTI3+PRN2xYBgEYH+r55X45v9xRrKNibwA4CprN+n7IpLrcE6RMjFfkTI2XyLiBTAOwX+dtmZ0vjrICzAwX9bu9U0AqgF8qKpdYr4iyAsw8/t6HsDjAALtvN+V/x6pG+uOBYOEeK3tvxwi2SbWIon5JYL93scCeAHASodzipSJ+YqEsfkSkRQAfwTwsKrWtX07xJC4zFeYvIzMl6o2q+plAAYDyBeR0W02MTJfEeQV9/kSkRkAqlW1uKPNQrzWFf4eqZvrjgVDJYCsVs8HAzgYxTZxz0tV61p2k6rqewB6iEiaw3lFwsR8hWVqvkSkB4L/UX5TVf8UYhMj8xUuL9O/L1U9DmAdgOlt3jL6+2ovL0PzNRHAjSJSgeBhy2ki8kabbbrk3yN1f92xYPgCwHARGSoiPQHcCWBVm21WAbjHOtv4SgC1qnrIdF4ikikiYj3OR/B/n6MO5xUJE/MVlon5suK9AqBMVZ9rZ7O4z1ckeRmar3QR6W897gPgGgDb22xmYr7C5mVivlT1Z6o6WFW9CP5/xEeq+sM2m3XJv0fq/pJNJxBrquoXkQcAfIDglQmvqupWEZltvb8AwHsInmlcDuAkgHu7SF63ApgjIn4ApwDcqaqO72oUkWUInhGeJiKVAJ5E8CQwY/MVYV4m5msigLsBbLGOfwPAzwEMaZWXifmKJC8T8zUQwGIRSULwP7grVHW16b/HCPMy8vcYSheYLyK2hiYiIqLwuuMhCSIiIooxFgxEREQUFgsGIiIiCosFAxEREYXFgoGIiIjCYsFA3Y6INEtwdcGtVlvfR0TEY72XJyLzDOW1IUafc5v13QIikheLzyQiCoeXVVK3IyINqppiPc5AcNW/T1X1SbOZxYaIjERwnYGFAB5T1SLDKRFRAuAeBurWVLUawCwAD1id8aaIyGoAEBGfiCwWkb+KSIWI3Cwi/09EtojIGqvVMkRkgoj8TUSKReQDsVYGFJF1IvIbESkUkZ0iMtl6Pdd6bZOIlIjIcOv1ButeRGSuiJRase6wXp9ifebbIrJdRN5s6TTY5juVqeqOeMwfEVELFgzU7anqHgR/6xkh3h6G4FLCNwF4A8DHqjoGwc5+11tFwwsAblXVCQBeBfCrVuOTVTUfwMMIdqIEgNkAfmctbJSHYO//1m4GcBmAsQi2JJ4r/7M88Tjrs0YBuBjBDo5ERMZ1u9bQRO0ItcIfALyvqmdEZAuCLbvXWK9vAeAFcAmA0QA+tP6xnwSgdd/+lkWeiq3tAeDvAH4hIoMB/ElVd7WJOQnAMlVtBnBYRP4G4HIAdQAKVbUSAKwWz14A6zv5XYmIYo57GKjbE5GLATQDqA7xdiMAqGoAwJlWawUEECyoBcBWVb3Muo1R1e+0HW99frL1Wb8HcCOCeyk+EJFpbVPqIN3GVo/PfiYRkWksGKhbE5F0AAsAvBjlwkE7AKSLyFXW5/UQkdwwMS8GsEdV5yG4suClbTb5BMAdIpJk5fcPAAqjyI2IKG5YMFB31KflskoAawH8FcAvo/kgVW1CcNXC34jIZgCbAFwdZtgdAEqtQwo5AJa0ef8dACUANgP4CMDjqloVaU4i8j1r9c6rAPxFRD6IdCwRUbR4WSURERGFxT0MREREFBYLBiIiIgqLBQMRERGFxYKBiIiIwmLBQERERGGxYCAiIqKwWDAQERFRWP8fMjUM12aPyZsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAEGCAYAAAAJ9c30AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtt0lEQVR4nO3de3iU9bUv8O+ahEsgKsYkBCE4SMBAQARivAD7gLW7qKi73q1V6j5uHjh61KrH3vbRqc/TU7tB26Jb0O0N0CLUVooUsFJFg1RjwiUEwiVIKAFCCEguBBIms84fM2Gn4yQzmTfv/ObNfD/PM8/c3t+slZ8TWXkv6yeqCiIiIqLOuEwnQERERPGPBQMRERGFxYKBiIiIwmLBQERERGGxYCAiIqKwkk0n0FXp6enqdrtNp0FE5CglJSW1qpphOg9yLscVDG63G8XFxabTICJyFBHZbzoHcjYekiAiIqKwWDAQERFRWCwYiIiIKCwWDERERBQWCwYiIiIKy7aCQUT6ikiRiGwVke0i8vMQ24iIzBeRChEpFZEJduVDRERE0bPzsspmANeoaqOI9AKwQUTWqOrn7ba5DsCIwO0KAAsC90RERBRHbNvDoH6Ngae9ArfgtbRvBrA4sO3nAAaIyCC7ciIiIqLo2Nq4SUSSAJQAyAHwn6r6RdAmgwEcaPe8KvDa4aDPmQVgFgAMHTo06ny+ddv3kJ4VfT2yr3w7ho3KS5ixJmM7cazJ2E4cazK2U3/mc07X4dVXX406NpEVthYMqtoK4DIRGQDgPREZo6pl7TaRUMNCfM4rAF4BgPz8/G+8H6n0rEG4838/Ee1wLHthXtTjnTjWZGwnjjUZ24ljTcZ26s9cuvTlqOMSWRWTqyRU9QSA9QCmB71VBSC73fMhAA7FIicip6msacS60mpU1jSG35iIqJvZeZVERmDPAkQkBcC1AHYGbbYSwH2BqyWuBFCnqodBRP+gsqYRC9buwZpNh7Bg7R4WDUQUc3YekhgEYFHgPAYXgOWqukpEZgOAqi4EsBrA9QAqADQBuN/GfIgcq6K6EV6fQhXw+hQV1Y1wZ6aaTouIEohtBYOqlgIYH+L1he0eK4AH7cqBqKfIyUpFskvg9SmSXYKcLBYLRBRbjlvemigRuTNTMWf6CFRUNyInK5V7F4go5lgwEDmEO5OFAhGZw7UkiIiIKCwWDERERBRWQh2S2Fe+HctemBf1+CNlJVE3Tikr2ghEGdtUXJOxnRjX6ngnzpfJ75cT59pq7OaqvVHHJbIqoQqGYaPyLHVoK136MjweT1Rjy2sbLHV3MxHXZGwnxgUAj8fjuLwtfUcsdj20krcT57o7YhOZwkMSREREFBYLBiIiIgqLBQMRERGFxYKBiIiIwmLBQJE5UAQUPue/JyKihJNQV0lQlA4UAYtuAlpbgKTewMyVQHaB6ayIiCiGuIeBwqss9BcL2uq/ryw0nREREcUYCwYKzz3Fv2dBkvz37immMyIiohhLqEMSVjs9lhVtRHltQ9RjrXTTMxG3feyUpm+h38mDaOo/GKd+tQzAMltjO3GurY536nx1x/fLxFiT3xF2eiQnSqiCwWqnR0td7Zw41mRsJ441GduJY03GdujPzE6PZBIPSRAREVFYLBiIiIgoLBYMREREFBYLBiIiIgqLBQMRERGFxYKBzqqsacS60mpU1jSaToWIiOJMQl1WSR2rrGnEgrV74PUpkl2COdNHwJ2ZajotIiKKE9zDQACAiupGeH0KVcDrU1RUcy8DERH9N9v2MIhINoDFALIA+AC8oqq/DdpmKoA/AdgXeOmPqvqMXTlZ7fR4pKwk6sYpVrq7xSLuiZMtqKs4Dp8qXCLYW52GY5/07tE/c3fHtTreifNltdOjlbydONdWY7PTI5lk5yEJL4DHVXWTiJwDoEREPlTVHUHbFarqDBvzOMtqp8fSpS/D4/FENba8tsFSd7dYxP1OTSMqqhuRk5V69nBET/+ZuzMuAHg8HsflbSWu1a6HVvJ24lx3R2wiU2wrGFT1MIDDgccNIlIOYDCA4IKB4oQ7M5XnLRARUUgxOYdBRNwAxgP4IsTbV4nIVhFZIyJ5HYyfJSLFIlJ89OhRO1MlIiKiEGwvGEQkFcAfADyqqvVBb28CcJGqjgPwAoAVoT5DVV9R1XxVzc/IyLA1XyIiIvomWwsGEekFf7Hwtqr+Mfh9Va1X1cbA49UAeolIup05ERERUdfZeZWEAHgNQLmqPt/BNlkAjqiqikgB/AXMMbtyIiIy5kARUFmIlKZq05kQRcXOqyQmAbgXwDYR2RJ47acAhgKAqi4EcBuAOSLiBXAKwF2qqjbmREQUeweKgEU3Aa0tGLqvGWm19+B4+njTWRF1iZ1XSWwAIGG2eRHAi3blQEQUFyoLgdYWQFsh6kP6kSIWDOQ47PRIRGQ39xQgqTcgSVBxoXZggemMiLosodaSsNrpsaxoI8prG6Iea6Wbnom4JmM7Ma7V8U6dr+CxJ0624HhjC9JSe2NA/95hx1vJ20lzndL0LfQ7eRAbTxzDrqV/BfDXLn8GOz2SSQlVMFjt9Gipq50Tx5qM7cSxJmPHydj2i5idiGQRszjJO2ZjARyzMJ6dHskkHpIgom7DRcyIei4WDETUbXKyUpHsEogAyS5BThZbjRP1FAl1SIKI7OXOTMWc6SO+sYgZETkfCwYi6lZcxIyoZ+IhCSIiIgqLBQMRERGFxYKBiIiIwmLBQERERGEl1EmPVjs9HikribpxipXOcqbimoztxLhWxztxvkx+v5w411Zjs9MjmZRQBYPVTo+lS1+Gx+OJamx5bYOl7m4m4pqM7cS4AODxeByXt6XviMWuh1byduJcd0dsIlN4SMIpDhQBhc/574mIiGIsofYwOFbdQWDRTf7lcZN6AzNXAtlc7Y6IiGKHexic4MR+f7Ggrf77ykLTGRERUYJhweAEAy7y71mQJP+9e4rpjIiIKMHwkIQTnDfYfxiistBfLPBwBBERxRgLBqfILmChQERExvCQBBEREYXFgoGIiIjCSqhDElY7PZYVbUR5bUPUY6100zMR12RsJ8a1Ot6p82Xy+xXLn/nEyRYcb2zB8T2bLX9H2OmRnCihCgarnR4tdbVz4liTsZ041mRsJ441GbuLYytrGrFg7R54fYq66kZccedsuDNTYxK7PXZ6JJN4SIKIKIyK6kZ4fQpVwKeKiupG0ykRxZxtBYOIZIvIxyJSLiLbReSRENuIiMwXkQoRKRWRCXblQ0QUrZysVCS7BCKASwQ5WVHuXSByMDsPSXgBPK6qm0TkHAAlIvKhqu5ot811AEYEblcAWBC4JyKKG+7MVMyZPgIV1Y3YW50W/eEIIgezbQ+Dqh5W1U2Bxw0AygEMDtrsZgCL1e9zAANEZJBdORERRcudmYprL83CgP69TadCZERMzmEQETeA8QC+CHprMIAD7Z5X4ZtFBURklogUi0jx0aNHbcuTiIiIQrO9YBCRVAB/APCoqtYHvx1iiH7jBdVXVDVfVfMzMjLsSJOIiIg6YWvBICK94C8W3lbVP4bYpApAdrvnQwAcsjMnIiIi6jrbTnoUEQHwGoByVX2+g81WAnhIRN6B/2THOlU9bFdOVhs3HSkrifo6aCvNWkzFNRnbiXHDjq876F+qfMBF/gXFgjhxvkx+v5z6HWHjJnKqDgsGETkXwE/g/6t/jar+rt17L6nq/wrz2ZMA3Atgm4hsCbz2UwBDAUBVFwJYDeB6ABUAmgDcH92PERmrjZtKl74Mj8cT1djy2gZLzVpMxDUZ24lxAcDj8YQef6AIWHQTMKAFSCr1rz4atJiYE+fLauMmK3l3ONcRMPkdsRqbyJTO9jC8AWAP/IcU/lVEbgXwPVVtBnBluA9W1Q0IfY5C+20UwIORp0vkUJWFQGsLoK3++8pCrj5KRI7S2TkMw1X1x6q6QlVvArAJwEcickGMciPqOdxTgKTegCT5791TTGdERNQlne1h6CMiLlX1AYCq/kJEqgB8CoBdS4i6IrvAfxiistBfLHDvAhE5TGcFw/sArgGwru0FVV0kIkcAvGB3YkQ9TnYBCwUicqwOCwZVfbKD19fC38qZiIiIEgRXqyQiIqKwWDAQERFRWCwYDKisacS60mpU1jSaToWIiCgiEXV6FJGrAbjbb6+qi23KyTZWOz2WFW1EeW1D1GPxwjycONmCLyuOw6cKlwguz0kLu/pdd8SNlqnYToxrdbxT58vk98tpc201Njs9kklhCwYRWQJgOIAtAFoDLysAxxUMVjs9WupqFxi7rrQauzcdgiogAgyfcCGuvTTL9rhRMxXbiWNNxnbiWJOxHfozs9MjmRTJHoZ8AKMDXRnJopysVCS7BF6fItklyMliSwsiIop/kRQMZQCyANi2KFQicWemYs70EaiobkROVircmSwYiMi5SkpKMpOTk18FMAY8L87JfADKvF7vAxMnTqwJtUEkBUM6gB0iUgSgue3FQLtoioI7k4UCEfUMycnJr2ZlZY3KyMj42uVycU+0Q/l8Pjl69Ojo6urqVwGE/Pc9koLB061ZERFRTzKGxYLzuVwuzcjIqKuurh7T0TZhCwZV/UREBgK4PPBSkaqG3F1BREQJx8VioWcI/Hfs8LBS2ONNInIHgCIAtwO4A8AXInJbt2VIREREcS+SE1R+BuByVZ2pqvcBKADwf+1Ni4iIKDL9+vUbH26bO++886KSkpK+APDjH//4H65lHz9+fG60MUK9/thjj1341FNPDWz/2uDBg8cePnw4GQB+9KMfZeXk5OSNHDlydG5u7uiPPvqof7j48SCScxhcQYcgjoFnwhKRw6TVbkb6kSLUDuSKoYlo2bJl+9sez58/f9Czzz5b3fZ88+bNO2OVx7p16/p/8MEHA7Zt27YjJSVFDx8+nNzc3Cyxim9FJAXDWhH5AMDSwPM7Aay2LyX7WO30eKSsJOrGKVa6u5mKazK2E+NaHe/E+TL5/erS2LqDwNalgK8VcCVh19fDsCyqqNa/I+z0CGyoqO3/WUXtOZNy0hsm56Sf7K7PXbVq1TnPPPPMhWlpaWd27dqVMnbs2KYVK1bsc7lcKCgouGTevHkH3nnnnfObm5tdubm5o0eOHHlq5cqV+/r16ze+qalpc11dnWv69Ok5dXV1SV6vV5566qlD3//+9090V34AcPDgwV5paWnelJQUBYBBgwZ5u/Pz7RTJSY//R0RuBTAJgAB4RVXfsz0zG1jt9Fi69GV4PJ6oxpbXNljq7mYirsnYTowLAB6Px3F5W/qOWOx6aCXvLs114XPAgN6AtgKShJodF+BaQ98Rq/+dnW5DRW3/f33zy5HeVp/rtQ37fK//4PLd3Vk0lJeXp2zZsuUrt9t9ZuLEibkffvhh6ne+852zC/e89NJLB998883MnTt37gge269fP9+f//znirS0NN/hw4eTr7jiitzvfe97J1yuru9UX7hw4cDly5df0Pa8pqamFwD8y7/8S/0vf/nLC91u95jJkyfX33333cdvuOEGRywsFNEsqOofVPUxVf2hU4sFIkpg7ilAUm9AkoCk3mjqP9h0Rgnrs4rac7ytPpdPAW+rz/VZRe053fn5Y8eOPTl8+PAzSUlJyMvLa9q7d2/ni/W04/P55NFHHx0ycuTI0dOmTRtZU1PTu6qqKqI1l4LNnj37yM6dO3e03TIzM88AwHnnnecrKyvb8eKLL+7PyMjwzpw5c/j8+fMvCPd58aDDgkFENgTuG0Skvt2tQUTqY5ciEZFF2QXAzJXANT8DZq7EqX5h1m8h20zKSW9ITnL5XAIkJ7l8k3LSo1/JK4Q+ffqcvcQzKSkJXq834vMDXn755bRjx44lb9u2rXznzp07LrjggjOnTp3q9nP2kpOTMWPGjIZf//rXh+bOnfv3FStWnN/dMezQYeWkqpMD991a/RERGZFd4L8BQNRnMJBVk3PST77+g8t323EOQ6SSk5O1ublZ2hcXAFBXV5eUnp5+pk+fPvr++++fc+jQoYj3TkRq69atfVwuF8aOHdsMAJs3b04ZMmRIS3fHsUMkq1UOB1Clqs0iMhXApQAWq+oJe1MjIqKeaHJO+kkThUKbe+655+ioUaNGjxkzpmnlypX72l5/4IEHjl933XU5Y8aMGZWXl9c0bNiw0+E+6/Tp066BAwde2vZ8zpw5Rzrbvr6+Punhhx8eWl9fn5SUlKRut7t50aJF+zsbEy8iOTbzBwD5IpID4DUAKwH8DsD1diZGREQUiaamps0AMGPGjIYZM2acPcSxePHiv7c9Lioq2tX2eMGCBQcBHAweP2jQIO+WLVtCXmLZtk0wn89XEkmOBw8e3NYWI5aXcXanSI7N+FTVC+C7AH6jqj8EMCjcIBF5XURqRKSsg/enikidiGwJ3J7qWupEREQUK5HsYTgjIncDmAngxsBrvSIY9yaAFwEs7mSbQlWdEcFnERERkUGR7GG4H8BVAH6hqvtEZBiAt8INUtVPARy3mB8RERHFgUgaN+0A8HC75/sAPNtN8a8Ska0ADgF4QlW3h9pIRGYBmAUAQ4cOjTqY1U6PZUUbUV4b3RVAVrvpmYhrMrYT41od79T5Mvn9ctpcW43dUzo9kjNFcpXEJAAeABcFthcAqqoXW4y9CcBFqtooItcDWAFgRKgNVfUVAK8AQH5+ftTLqFrt9Gipq50Tx5qM7cSxJmM7cazJ2A79mXtCp0dyrkjOYXgNwA8BlABo7a7Aqlrf7vFqEXlJRNJVtba7YhAREVH3iOQchjpVXaOqNap6rO1mNbCIZImIBB4XBHKx/LlERJR4Fi9ePEBEJm7evLmv6VxC2bhxY8qyZcvO6+q4goKCSz799NN+kby+atWqc6ZNm5bT/rVbb73V/cYbb5wPAEuXLj1v1KhRoy+55JLRw4cPz5s7d256V3KJZA/DxyIyF8AfATS3vaiqmzobJCJLAUwFkC4iVQCeRuDqClVdCOA2AHNExAvgFIC7VDXqww1ERJS43nnnnbQJEyY0LlmyJG38+PGHTOcTrLi4uF9xcXH/O++8s85E/ObmZnnkkUcu+tvf/lY+fPjwM6dOnZLdu3d3qZNlJHsYrgCQD+D/AXgucAt7xo6q3q2qg1S1l6oOUdXXVHVhoFiAqr6oqnmqOk5Vr1TVjV1JnIiIHOqr9f2xzpOFr9b3746Pq6urcxUXF6e+8cYble+9997ZdRmC/+K+7777hrYt9LRs2bLzhg0bljdx4sRLfvCDH2S3bffYY49deMstt7gnTZo0YvDgwWMXLVo0YPbs2UNGjhw5esqUKSOam5sFAAoLC/tdfvnll+Tl5Y2aPHnyiP379/cC/H/5z5kzZ/DYsWNHud3uMWvXrk09ffq0/PKXv7zw/fffPz83N3f0f/3Xf51fX1/vuv32291jxowZNWrUqNFvvfXWAABobGyUGTNmXDxy5MjRN9xww8WnT5+OeC2Mzpw4ccLl9Xpl4MCBXgBISUnRcePGNYcb117YgkFVp4W4XRNt0kRElMC+Wt8fb98xEp/9djDevmNkdxQNb7/99oCpU6fWXXrppc0DBgxo3bBhwzd24bfX1NQkjzzyyEVr1qzZU1JSsuvYsWP/sLd9//79fT766KOKd999t2L27NnDrrnmmvrdu3fv6Nu3r2/58uXnNTc3y8MPPzz0T3/6097t27eXz5w5s/aJJ544uwSq1+uVbdu2lf/qV7868Mwzz1zYt29f/clPfnLoxhtv/Hrnzp07/u3f/u3rn/70p4OmTZtWX1ZWVl5YWLjr3//934fU19e75s2bl5mSkuLbvXv3jqeeeurwjh07ujQ/xcXFqbm5uaPbbuvWrRsAAAMHDmz99re/fWLo0KGX3njjjcMWLFiQ1tratdMSwxYMIjJQRF4TkTWB56NF5H92KYpDVdY0Yl1pNSprHLFUORFR/Ptq/TnwnXFBfYDvjAtfrbe8wOHy5cvT7r777q8B4NZbbz2+ZMmStM6237JlS9/s7Ozm3NzcFgC46667/qFn0LXXXlvXp08fLSgoONXa2iq33XZbPQDk5eWd2rdvX+/S0tI+e/bsSbnmmmtG5ubmjp47d+6gQ4cOnW1oePvtt38NAFdfffXJqqqqkLv9169ff+6vf/3rQbm5uaMnT558SXNzs1RUVPTesGFD6r333nsMAK644opTI0eObOrKXOTn5ze2X1b72muvPdH23rJly/avXbt2d35+/sn58+dn3XHHHe6ufHYk5zC8CeANAD8LPN8N/1Jvr3UlkNNU1jRiwdo98PoUyS7BnOkhr/gkIqKuuHhqA/72kg++My64evlw8VRLy1tXV1cnff755+fu3r075aGHHkJra6uIiC5YsKCqV69e6vP5zm7bdjgh3OlybatYJiUlITk5WV0u/9/WLpcLXq9XVFVycnJOdbTuRN++fRXwL2Pd2toa8pCCquLdd9+tCHVYIHA9gC0KCgpOFRQUnJo1a9bxnJycsQAqIx0byTkM6aq6HIAPAALrSnTb5ZXxqqK6EV6fQhXw+hQV1dzLQERk2cVTT+Ke5bsx6ZGDuGf5blw81dKqlUuWLDn/lltuOXbo0KFtBw8e3FZdXV06ZMiQlr/85S+pw4cPb66oqEg5deqUHDt2LGnDhg3nAsC4ceNOHzhwoM+uXbt6A8CyZcs63SMR7NJLLz19/Pjx5HXr1vUH/IVIcXFxp1dnnHvuua2NjY1n/82dNm1a/XPPPTewraD57LPPUgBg8uTJjW+99VYaAHz55Zd9d+/e3enhlUjV1dW5Vq1adXZvzhdffJFy4YUXdmlZ7Uj2MJwUkQsAKACIyJUAjJzlaVVXOj2eONmCuorj8KnCJYK91Wlo3rct6sYpVrq7HSkrMRLXZOxvxK07CJzYDwy4CDhvcMcDuztuDMc78Tti8vvlxLm2GrvHdHq8eOpJq4VCm9///vcXPPnkk4fbv3bzzTd/vWTJkrS333777zfeeOPXo0aNyhs2bNjpvLy8JgBITU3V559/fv/06dNHpKWlecePH9+lXPr27avvvPPO3ocffnhoQ0NDUmtrq8yZM+dIfn5+h8thX3fddQ3z5s0blJubO/rxxx8//Oyzzx6aNWvW0Nzc3NGqKkOGDGn++OOPK5544omau+66a9jIkSNH5+XlNY0dO7bD3L773e+OSE5OVgCYMGFC44MPPni0o219Ph/mzp078KGHHrqob9++vn79+vlee+21fR1tH4qE2zUjIhMAvABgDIAyABkAblPV0q4E6i75+flaXFwc1dg7H3q8Sx3WKmsaUVHdiJysVLgzU1G69GV4PJ6YxG7PVFyTsf8h7oEiYNFNQGsLkNQbmLkSyC6wP24UPB6P+fmKYdxlFrseWsnbiXNtMraIlKhqflSDO7F169bKcePGOa7hXl1dneu8887z+Xw+3HfffUNHjBhx+umnn64xnZdpW7duTR83bpw71HuRrCWxSUT+B4BL4G8LvUtVz3RvivHJnekvFMiwykJ/saCt/vvKwk4LBiKicH7zm9+kL126NP3MmTOSl5fX9Nhjjzmu6Im1SNaSSAJwPQB3YPt/FhGo6vM250bk557i37PQtofBPcV0RkTkcE8//XQN9yh0TSTnMLwP4DSAbQic+EgUU9kF/sMQlYX+YoF7F4jiic/n84nL5WKnXofz+XyCTv6dj6RgGKKql3ZfSkRRyC5goUAUn8qOHj06OiMjo45Fg3P5fD45evToefCfqxhSJAXDGhH5Z1X9S/elRkREPYHX632gurr61erq6jGI7FJ9ik8+AGVer/eBjjaIpGD4HMB7IuICcAb+Ex9VVc/tnhyJiMipJk6cWAPgJtN5kP0iKRieA3AVgG1cTZKIiCgxRbL7aA+AMhYLREREiSuSPQyHAawPLD51tue1Ey+r7Eqnx1DKijaivDa6tudWu+mZiGsythPjWh3vhPk6cbIFxxtbkJbaGwP69zb+/XLaXFuN3WM6PZIjRVIw7AvcegdujjVsVJ6lrnSw0tXOiWNNxnbiWJOxYzC2/YJsJ9otyMb5it14Ky2piayKpNPjz2ORCBHFNy7IRpTYOiwYROQ3qvqoiLyPwMJT7akqz4olSiA5WalIdsnZJd9zslJxzHRSRBQzne1hWBK4j/4AJRH1GO7MVMyZPuIfFmT7wnRSRBQzHRYMqloSuP9ERDICjztcOpOIej4uyEaUuDq8rFL8PCJSC2AngN0iclREnopdekRERBQPOuvD8CiASQAuV9ULVPV8AFcAmCQiP4xFckRERBQfOisY7gNwt6rua3tBVb8C8P3Ae0RERJQgOisYeqlqbfCLgfMYetmXEhEREcWbzq6SaInyPQCAiLwOYAaAGlUdE+J9AfBbANcDaALwA1XdFO5zrbDa6fFIWUnUjVOsdHczFddkbCfGtTreifNl8vvlxLm2GpudHsmkzgqGcSJSH+J1AdA3gs9+E8CLABZ38P51AEYEblcAWBC4t43VTo+lS1+Gx+OJamx5bYOl7m4m4pqM7cS4AODxeByXt6XviMWuh1byduJcd0dsIlM6u6wyycoHq+qnIuLuZJObASwOLGr1uYgMEJFBqnrYSlwiIiLqfpGsVmmXwQAOtHteFXjtG0RklogUi0jx0aNsBUFERBRrJgsGCfFayCW0VfUVVc1X1fyMjAyb0yIiIqJgJguGKgDZ7Z4PAXDIUC5ERETUCZMFw0oA9wU6Sl4JoI7nLxAREcWnsMtbR0tElgKYCiBdRKoAPI1A/wZVXQhgNfyXVFbAf1nl/Xbl0h3SajcD+zcCB4qA7ALT6RAREcWUbQWDqt4d5n0F8KBd8btTWu1mTP7ofqzfVw8sugmYuZJFAxERJRSThyQcI/1IEVytLQB8QGsLUFloOiUiIqKYsm0PQzyKttNjSlM1hu5rwSeVrfCtb8bf92/CqWWPd+kzrHbTK69tiHqslU58pmI7Ma7V8U6dL5PfL6fNtdXY7PRIJiVUwWCl02Na7T2onf8fuObhJ3E8fXzXP8BKRzxTY03GduJYk7GdONZkbIf+zOz0SCbxkESEjqePx7GMidEVC0RERA7HgoGIiIjCYsFAREREYbFgICIiorBYMBAREVFYCXWVBFEolTWNqKhuRE5WqulUiIjiFgsGSmiVNY1YsHYPvD5Fsktw8ckW0ykREcWlhCoYom3c1OZIWUnU10FbadZiKq7J2LGK+9WRRhw93K4JT8NeS9e69/T56s6xgLW8nTjXVmOzcROZlFAFg5XGTYC/aYrH44lqbHltg6VmLSbimowdq7jBeximnfwk6rgA4PF4evR8fYPFJkZW8nbiXHdHbCJTEqpgIArmzkzFnOkjzp7DUP/XItMpERHFJRYMlPDcmalwZ/pPeCw1nAsRUbziZZVEieJAEVD4nP+eiKiLuIeBKBEcKAIW3eRfnj2pNzBzJZBdYDorInIQ7mEgSgSVhf5iQVv995WFpjMiIodhwUCUCNxT/HsWJMl/755iOiMichgekiBKBNkF/sMQlYX+YoGHI4ioi1gwECWK7AIWCkQUNVFV0zl0SX5+vhYXF0c1tuBb0zFsVF7UsY+UlWDq1KlRjf396g8wpuBqR8U1GduJcQFg/fr1jsvbStyyoo3Gvl9OnGursZur9mLFihVRjRWRElXNj2owERJsDwM7PTonthPjAs7sPshOj13DTo+UqHjSIxEREYXFgoGIiIjCsrVgEJHpIrJLRCpE5Mch3p8qInUisiVwe8rOfIiIiCg6tp3DICJJAP4TwLcBVAH4UkRWquqOoE0LVXWGXXkQUQ9woIiXhBIZZudJjwUAKlT1KwAQkXcA3AwguGAgIuoY21oTxQU7D0kMBnCg3fOqwGvBrhKRrSKyRkRCXvMoIrNEpFhEio8ePWpHrkQUr9jWmigu2FkwSIjXgps+bAJwkaqOA/ACgBWhPkhVX1HVfFXNz8jI6N4siSi+sa01UVyw85BEFYDsds+HADjUfgNVrW/3eLWIvCQi6apaa2NeROQkbGtNFBfsLBi+BDBCRIYBOAjgLgDfa7+BiGQBOKKqKiIF8O/xOGZXQvvKt2PZC/OiHl9WtBHltQ1Rj0WUsU3FNRnbiXGtjnfqfMX2+7UscHPmXFuN3Vy1N+q4RFbZVjCoqldEHgLwAYAkAK+r6nYRmR14fyGA2wDMEREvgFMA7lIbe1Vb7fRoqaudE8eajO3EsSZjO3GsydgO/ZnZ6ZFMsrU1tKquBrA66LWF7R6/COBFO3MgIiIi69jpkYiIiMJiwUBERERhsWAgIiKisFgwEBEZVFnTiHWl1aisaTSdClGnbD3pkYiIOlZZ04gFa/fA61MkuwRzpo+AOzPVdFpEIXEPAxGRIRXVjfD6FKqA16eoqOZeBopfLBiIiAzJyUpFsksgAiS7BDlZ3LtA8Uts7JNki/z8fC0uLo5qbMG3pmPYqJDrW0WkrGgjxhRcnTBjTcZ24liTsZ041mTsePqZT5xswfHGFqSl9saA/r07HdtctRcrVqyIKq6IlKhqflSDiZBg5zCw06ODYjtxrMnYThxrMrZDf2Z2eiSTeEiCiIiIwmLBQERERGGxYKCY4jXnRETOlFDnMJBZvOaciMi5uIeBYobXnBMRORcLBooZXnNORORcPCRBMePOTMWc6SNQUd2InKxUHo4gInIQFgwUU+5MFgpERE6UUAXDvvLtWPbCvKjHHykribpxSlnRRiDK2KbimoztxLhWxztxvkx+v5w411ZjN1ftjToukVUJVTBY7fRYuvRleDyeqMaW1zZY6u5mIq7J2HEX90ARUFkIuKcA2QUdjvd4PPGVt81xrXY9tJK3E+e6O2ITmZJQBQNRVA4UAYtuAlpbgKTewMyVnRYNREQ9Ea+SIAqnstBfLGir/76y0HRGREQxx4KBKBz3FP+eBUny37unmM6IiCjmeEiCKJzsAv9hiAjOYSAi6qlYMBBFIruAhQIRJTRbD0mIyHQR2SUiFSLy4xDvi4jMD7xfKiIT7MyHiIiIomNbwSAiSQD+E8B1AEYDuFtERgdtdh2AEYHbLAAL7MqHiIiIomfnHoYCABWq+pWqtgB4B8DNQdvcDGCx+n0OYICIDLIxJyIiIoqCqKo9HyxyG4DpqvpA4Pm9AK5Q1YfabbMKwLOquiHw/K8AfqSqxUGfNQv+PRAYOnToxP3790eV0wMPPIAhQ4ZENRYAtmzZgssuuyxhxpqM7cSxJmM7cazJ2E79mauqqvDqq69GNVZESlQ1P6rBRLD3pEcJ8VpwdRLJNlDVVwC8AgD5+flRVzjR/qIRERElOjsPSVQByG73fAiAQ1FsQ0RERIbZWTB8CWCEiAwTkd4A7gKwMmiblQDuC1wtcSWAOlU9bGNOREREFAXbDkmoqldEHgLwAYAkAK+r6nYRmR14fyGA1QCuB1ABoAnA/XblQ0RERNGztXGTqq6Gvyho/9rCdo8VwIN25kBERETWcS0JIiIiCosFAxEREYXFgoGIiIjCYsFAREREYdnW6dEuInIUQHStHoF0ALXdmE53ide8gPjNjXl1DfPqmp6Y10WqmtGdyVBicVzBYIWIFMdja9R4zQuI39yYV9cwr65hXkTfxEMSREREFBYLBiIiIgor0QqGV0wn0IF4zQuI39yYV9cwr65hXkRBEuocBiIiIopOou1hICIioiiwYCAiIqKwemTBICLTRWSXiFSIyI9DvC8iMj/wfqmITIiTvKaKSJ2IbAncnopRXq+LSI2IlHXwvqn5CpdXzOdLRLJF5GMRKReR7SLySIhtYj5fEeZlYr76ikiRiGwN5PXzENuYmK9I8jLy+xiInSQim0VkVYj3jPw+EkFVe9QN/qW09wK4GEBvAFsBjA7a5noAawAIgCsBfBEneU0FsMrAnP0TgAkAyjp4P+bzFWFeMZ8vAIMATAg8PgfA7jj5fkWSl4n5EgCpgce9AHwB4Mo4mK9I8jLy+xiI/RiA34WKb+r3kTfeeuIehgIAFar6laq2AHgHwM1B29wMYLH6fQ5ggIgMioO8jFDVTwEc72QTE/MVSV4xp6qHVXVT4HEDgHIAg4M2i/l8RZhXzAXmoDHwtFfgFnymtYn5iiQvI0RkCIAbALzawSZGfh+JemLBMBjAgXbPq/DN/3FGso2JvADgqsBu0jUikmdzTpEyMV+RMjZfIuIGMB7+v07bMzpfneQFGJivwO71LQBqAHyoqnExXxHkBZj5fv0GwJMAfB28H8+/j9SD9cSCQUK8FvyXQyTbdLdIYm6Cv9/7OAAvAFhhc06RMjFfkTA2XyKSCuAPAB5V1frgt0MMicl8hcnLyHypaquqXgZgCIACERkTtImR+Yogr5jPl4jMAFCjqiWdbRbitXj4faQericWDFUAsts9HwLgUBTbxDwvVa1v202qqqsB9BKRdJvzioSJ+QrL1HyJSC/4/1F+W1X/GGITI/MVLi/T3y9VPQFgPYDpQW8Z/X51lJeh+ZoE4CYRqYT/sOU1IvJW0DZx+ftIPV9PLBi+BDBCRIaJSG8AdwFYGbTNSgD3Bc42vhJAnaoeNp2XiGSJiAQeF8D/3+eYzXlFwsR8hWVivgLxXgNQrqrPd7BZzOcrkrwMzVeGiAwIPE4BcC2AnUGbmZivsHmZmC9V/YmqDlFVN/z/j/hIVb8ftFlc/j5Sz5dsOoHupqpeEXkIwAfwX5nwuqpuF5HZgfcXAlgN/5nGFQCaANwfJ3ndBmCOiHgBnAJwl6ravqtRRJbCf0Z4uohUAXga/pPAjM1XhHmZmK9JAO4FsC1w/BsAfgpgaLu8TMxXJHmZmK9BABaJSBL8/+AuV9VVpn8fI8zLyO9jKHEwX0RsDU1ERETh9cRDEkRERNTNWDAQERFRWCwYiIiIKCwWDERERBQWCwYiIiIKiwUD9Tgi0ir+1QW3B9r6PiYirsB7+SIy31BeG7vpc24P/Gw+Ecnvjs8kIgqHl1VSjyMijaqaGnicCf+qf5+p6tNmM+seIjIK/nUGXgbwhKoWG06JiBIA9zBQj6aqNQBmAXgo0BlvqoisAgAR8YjIIhH5i4hUisgtIvIfIrJNRNYGWi1DRCaKyCciUiIiH0hgZUARWS8ivxKRIhHZLSJTAq/nBV7bIiKlIjIi8Hpj4F5EZK6IlAVi3Rl4fWrgM98VkZ0i8nZbp8Ggn6lcVXfFYv6IiNqwYKAeT1W/gv+7nhni7eHwLyV8M4C3AHysqmPh7+x3Q6BoeAHAbao6EcDrAH7RbnyyqhYAeBT+TpQAMBvAbwMLG+XD3/u/vVsAXAZgHPwtiefKfy9PPD7wWaMBXAx/B0ciIuN6XGtoog6EWuEPANao6hkR2QZ/y+61gde3AXADuATAGAAfBv7YTwLQvm9/2yJPJYHtAeBvAH4mIkMA/FFV9wTFnAxgqaq2AjgiIp8AuBxAPYAiVa0CgECLZzeADV38WYmIuh33MFCPJyIXA2gFUBPi7WYAUFUfgDPt1grwwV9QC4DtqnpZ4DZWVf85eHzg85MDn/U7ADfBv5fiAxG5JjilTtJtbvf47GcSEZnGgoF6NBHJALAQwItRLhy0C0CGiFwV+LxeIpIXJubFAL5S1fnwryx4adAmnwK4U0SSAvn9E4CiKHIjIooZFgzUE6W0XVYJYB2AvwD4eTQfpKot8K9a+CsR2QpgC4Crwwy7E0BZ4JBCLoDFQe+/B6AUwFYAHwF4UlWrI81JRL4bWL3zKgB/FpEPIh1LRBQtXlZJREREYXEPAxEREYXFgoGIiIjCYsFAREREYbFgICIiorBYMBAREVFYLBiIiIgoLBYMREREFNb/B0hmCtrpXtMMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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() "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}