{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate D-efficient designs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains example code from the article [Two-level designs to estimate all main effects and two-factor interactions](https://doi.org/10.1080/00401706.2016.1142903) by Eendebak, P. T. and Schoen, E. D.\n",
    "This example shows how to generate D-efficient designs with a user-specified optimization function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import oapackage\n",
    "import oapackage.Doptim\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the class of designs to generate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We generate D-efficient designs with 40 rows and 7 columns\n",
      "\n"
     ]
    }
   ],
   "source": [
    "run_size = 40\n",
    "number_of_factors = 7\n",
    "factor_levels = 2\n",
    "strength = 0\n",
    "arrayclass = oapackage.arraydata_t(factor_levels, run_size, strength, number_of_factors)\n",
    "print(\"We generate D-efficient designs with %d rows and %d columns\\n\" % (run_size, number_of_factors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a single D-efficient design using $\\alpha=(1,2,0)$ as the parameters for the optimization function. For details on this parameter and its corresponding optimization function, see [Two-Level Designs to Estimate All Main Effects and Two-Factor Interactions](https://doi.org/10.1080/00401706.2016.1142903)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doptim: optimization class 40.2-2-2-2-2-2-2\n",
      "Doptimize: iteration 0/30\n",
      "Doptimize: iteration 29/30\n",
      "Doptim: done (11 arrays, 0.4 [s])\n"
     ]
    }
   ],
   "source": [
    "alpha = [1, 2, 0]\n",
    "scores, design_efficiencies, designs, ngenerated = oapackage.Doptim.Doptimize(\n",
    "    arrayclass, nrestarts=30, optimfunc=alpha, selectpareto=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Generated 11 designs, the efficiencies for these designs are:\n",
      "array 0: D-efficiency 0.8815, Ds-efficiency 0.9807\n",
      "array 1: D-efficiency 0.9076, Ds-efficiency 0.9628\n",
      "array 2: D-efficiency 0.8670, Ds-efficiency 0.9827\n",
      "array 3: D-efficiency 0.8945, Ds-efficiency 0.9669\n",
      "array 4: D-efficiency 0.9027, Ds-efficiency 0.9540\n",
      "array 5: D-efficiency 0.8851, Ds-efficiency 0.9549\n",
      "array 6: D-efficiency 0.8737, Ds-efficiency 0.9581\n",
      "array 7: D-efficiency 0.9036, Ds-efficiency 0.9400\n",
      "array 8: D-efficiency 0.8614, Ds-efficiency 0.9595\n",
      "array 9: D-efficiency 0.8897, Ds-efficiency 0.9418\n",
      "array 10: D-efficiency 0.9046, Ds-efficiency 0.9203\n",
      "\n",
      "The design with the highest D-efficiency (0.9076) is:\n",
      "\n",
      "0010011001011001010101011101101100011101\n",
      "0011010111100101111100010111001100100010\n",
      "1011011111001010100010111001000010011010\n",
      "0101100001100111110011010001001100111110\n",
      "0011100010001101001011101101110100010001\n",
      "0011110000011011101001110011101010101000\n",
      "1000110011101000101110001001101111001101\n"
     ]
    }
   ],
   "source": [
    "print(\"\\nGenerated %d designs, the efficiencies for these designs are:\" % len(designs))\n",
    "for ii, d in enumerate(designs):\n",
    "    dd = d.Defficiencies()\n",
    "    print(\"array %d: D-efficiency %.4f, Ds-efficiency %.4f\" % (ii, dd[0], dd[1]))\n",
    "\n",
    "D = [d.Defficiency() for d in designs]\n",
    "best = np.argmax(D)\n",
    "print(f\"\\nThe design with the highest D-efficiency ({D[best]:.4f}) is:\\n\")\n",
    "\n",
    "designs[best].transposed().showarraycompact()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimizing with a different optimization target leads to different D-efficient designs. Below we compare the sets of designs generated with\n",
    "optimization target [1,0,0] and [1,2,0]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doptim: optimization class 40.2-2-2-2-2-2-2\n",
      "Doptimize: iteration 0/30\n",
      "Doptimize: iteration 29/30\n",
      "Doptim: done (13 arrays, 0.4 [s])\n"
     ]
    }
   ],
   "source": [
    "scores0, design_efficiencies0, designs0, _ = oapackage.Doptim.Doptimize(\n",
    "    arrayclass, nrestarts=30, optimfunc=[1, 0, 0], selectpareto=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pareto: 23 optimal values, 24 objects\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEICAYAAABxiqLiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlcTfn/B/DXvSVJIrJUV5SSFu3Zl5JEDAnZlWUylpmMZb4zYyxjGGZkMsRgBlnCMNl+hoZB2SVkqyxtKjQtlva69fn9YdxxdZfu7Xbvrd7Px6PHwz2fc895Fe678/mc8/lwGGMMhBBCiIy4qg5ACCGkbqICQgghRC5UQAghhMiFCgghhBC5UAEhhBAiFyoghBBC5EIFhBBCiFyUXkAiIyNhaWkJc3NzrFmzpkp7WloaPDw8YGdnBzc3N2RkZAjanj59ikGDBsHKygrW1tZITU1VYnJCCCHv4yjzQcKKigp07twZZ86cAY/Hg6urK/bv3w9ra2vBPmPGjMGwYcPg7++Pc+fOYefOndizZw8AwM3NDYsXL4anpycKCgrA5XKho6OjrPiEEELeo6nMk8XExMDc3BxmZmYAgHHjxuHYsWNCBSQ+Ph4hISEAAHd3d/j4+Ai28/l8eHp6AgB0dXWlns/AwAAdO3ZU8HdBCCH1W2pqKnJycqTup9QCkpmZifbt2wte83g8XL9+XWgfe3t7REREICgoCEeOHEF+fj5yc3Px6NEjtGjRAr6+vkhJScHAgQOxZs0aaGhoCL1/27Zt2LZtGwCgadOmiI2Nrf1vjBBC6hEXF5dq7afUMRBRvWUcDkfodXBwMKKjo+Ho6Ijo6GgYGxtDU1MTfD4fFy9eRHBwMG7cuIHk5GSEhYVVOV5gYCBiY2MRGxuL1q1b19a3QgghDZ5SCwiPx0N6errgdUZGBoyMjIT2MTIywuHDh3H79m2sWrUKANC8eXPweDw4OjrCzMwMmpqa8PHxwa1bt5QZnxBCyHuUWkBcXV3x+PFjpKSkoKysDAcOHMDw4cOF9snJyUFlZSUAYPXq1Zg2bZrgvS9fvkR2djYA4Ny5c0JjJ4QQQpRLqQVEU1MToaGh8PLygpWVFfz8/GBjY4OlS5fi+PHjAICoqChYWlqic+fOyMrKwuLFiwEAGhoaCA4OhoeHB7p27QrGGD7++GNlxieEEPIepd7Gq2wuLi40iE4IITKq7mcnPYlOCCFELlRACCGEyIUKCCGkRiorK7FiRTQ6d94IE5MQDBu2D48f56o6FlECpT5ISAipX/75pwCWlpvw6lWJYFt6+hv8+edjBAd7YsGCXipMR2obXYEQQuTWv3+YUPF438KFZ5CUlKfkRESZqIAQQuSSlJSHxETJXVXz5kUqKQ1RBSoghBC5nDz5WOo+d+5kKSEJURUqIIQQuejra0vdp0mTRkpIQlSFCgghRC4TJnSFpqbkj5ApU+yUlIaoAhUQQohcuFwu5s51Fduur6+Nr77qo8RERNmogBBC5BYSMhiff94DjRoJf5R06WKAxMQ54HLpI6Y+o+dACCE18tNPXggO9sT+/Q/w8mUxhg61gKmpvqpjESWgAkIIqTEul4uJE7uqOgZRMrq+JIQQIhcqIIQQQuRCBYQQQohcaAyE1CsPHvyD4mI+rKwM0LSplqrjEFKvUQEh9UJERDyWLDmPhIQcAEDz5o3x8cdOWLXKA1paGipOR0j9RAWE1HkHDtzHhAkReH9x5tevSxEcfBVPnrzEkSNjVReOkHqMxkBInVZZyfDll38LFY/3HT2aiCtX0pUbipAGggoIqdOuXctAWtprifscOHBfSWkIaViogJA6LT+/tBr7lCkhCSENDxUQUqd17dpW6oywTk7tlJSGkIaFCgip04yMmmHUKCux7c2bN8aUKfZKTERIw0EFhNR5v/wyFC4uRlW2N2umhcOHx6J5c+kLHxFCZEe38ZI6T1+/Ca5enY4jRxIQEZGA4mI+evQwxvTpTmjTpqmq4xFSb1EBIfWCpiYXY8bYYMwYG1VHIaTBoC4sQgghcqECQgghRC5UQAghhMhF6QUkMjISlpaWMDc3x5o1a6q0p6WlwcPDA3Z2dnBzc0NGRoZQ+5s3b2BsbIy5c+cqK3K9VFrKR3j4XXz99VkEB19BRsYbVUcihNQxSi0gFRUVmDNnDk6dOoX4+Hjs378f8fHxQvssXLgQU6ZMwd27d7F06VJ89dVXQu1LlixB//79lRm73rl4MQ0mJusxadIRrF59CYsWnUHHjuvxzTfnVB2NEFKHKLWAxMTEwNzcHGZmZtDS0sK4ceNw7NgxoX3i4+Ph4eEBAHB3dxdqv3nzJrKysjBo0CBlxq5XMjLeYOjQffjnn0Kh7RUVDKtWXcTWrbEqSkYIqWuUWkAyMzPRvn17wWsej4fMzEyhfezt7REREQEAOHLkCPLz85Gbm4vKykosWLAAa9eulXiObdu2wcXFBS4uLsjOzlb8N1HHbdkSK3FuqODgq2DiprYlhJD3KLWAiPpg4nA4Qq+Dg4MRHR0NR0dHREdHw9jYGJqamti8eTO8vb2FCpAogYGBiI2NRWxsLFq3bq3Q/PVBVFSqxPYnT/KQnk7jIYQQ6ZT6ICGPx0N6+n9rM2RkZMDISHgKCiMjIxw+fBgAUFBQgIiICDRv3hxXr17FxYsXsXnzZhQUFKCsrAy6uroiB+KJeNImHqzuPoQQotQC4urqisePHyMlJQXGxsY4cOAA9u3bJ7RPTk4OWrZsCS6Xi9WrV2PatGkAgPDwcME+YWFhiI2NpeIhh6FDLRAdnSa23cnJEEZGzZSYiBBSVyn1V01NTU2EhobCy8sLVlZW8PPzg42NDZYuXYrjx48DAKKiomBpaYnOnTsjKysLixcvVmbEem/6dCcYG4svEN9801eJaQghdRmH1eMRUxcXF8TG0l1FH3r4MAcTJx7GzZvPBdtatWqCtWs9MXWqowqTEULUQXU/O2kyxQbI0tIAsbGBuH49A/Hx2dDXb4IhQ8zRuDH9cyCEVB99YjRg3bvz0L07T9UxCCF1FN1uQwghRC5UQAghhMiFurBIrTp7NhmhoTdw924WWrTQxoQJtvj4Y2fo6TVWyPEvXXqKjRtjcPPmMzRtqgU/P2vMmuWKli2bKOT4hBDxqICQWrNiRTSWLYsS2nbr1nPs3BmHqKgAGBjo1Oj4GzdeR1BQJN6/j/Du3Sxs334bFy5MBY+nV6PjE0Ikoy4sUitiY59VKR7vPHiQjUWLztTo+E+e5GHevL8g6ib0lJRXmDv3ZI2OTwiRjgoIqRXbtt2U2H7gwH28fl0i9/F/++0WKivFP8J04sQjPHuWL/fxCSHSUReWBCUlfBw+nICkpDy0a6cLPz8bNG+urepYdcKTJ3kS20tK+MjMzJf75ynt+BUVDCkpL2laFkJqERUQMU6deozJk48gN7dYsO3zz//CTz95ITDQWYXJ6oa2bXUltnO5HLRuLf8YSNu2TaXu066d5AyEkJqhLiwR4uOz4et7UKh4AEBhYTk++eQETp16rKJkdUdAgL3E9qFDLdC6tfQiIP74DhLb+/QxQadOLeU+PiFEOiogIqxffw0lJXyRbYwBP/54RcmJ6p5Bgzph3DhbkW2tWjXBjz961uj4rq7G+OQT0VeCzZppYf16rxodnxAiHRUQEc6eTZHYHhWVCj6/Uklp6iYOh4O9e0ciJMQL5uZvrwS0tTUREOCA69dnoEsXgxqfY/PmodiyZSisrd8uHKalpYFx42xx9ep0ODsbSXk3IaSmaDZeEczNNyAp6aXYdi6Xg7Kyb6ChQfW3uoqLy9G4sSa4XI70neU8vpaWBv2dEKIA1f3slOl/28SJE3Hx4kW5Q9UV3t4WEtsHDepEH1QyatKkkcKLR2rqK8yadQItW/4APb016NlzO/buvavQcxBCxJPpU/Dq1atwc3ODtbU1NmzYgFevXtVWLpUKCuqOZs20RLZpaHDw5Ze9lZyIfCgxMQcuLtuwZctNvHxZAj6/EjduPMPkyUcwb16kquMR0iDIVECSk5Nx8uRJdOnSBQsXLoSxsTGmTp2Ka9eu1VY+lejUqSVOnZqIjh1bCG1v3VoH+/aNQv/+HVUTjAh89tmpKnfJvfPzz9dx40amkhMR0vDIPQby4sUL/Prrr9i+fTvS09NhZ2eHmTNnYtKkSdDVVY/772u6ImFlJcPp00lISsqDoWEzDB1qQYsuKVhRUTm2bIlFWFgcXrwogJmZPgIDnREQ4CC2yysj4w3atw+ReNzRo61w6JBfbUQmpN6r7mdnjQfRnz17hgkTJuDChQsAAF1dXcycORPLly9H06by3+evCLSkrXorLCzDwIF7cO1aRpW2MWOsceDAaJFF5Nq1dPTsuUPisdu318PTp58rLCshDUmtDKK/79y5c/Dz84OpqSnu3buHzz//HFeuXMGnn36KLVu2YMqUKfIemkiRn1+KZcvOo337EHC536Jjx/VYufICiorKVR1NJj/+eFlk8QCAQ4fisW/fPZFt2tqNpB5b3HM8hBDFkak/Jjc3Fzt37sS2bduQlJQEZ2dnbNmyBePHj4e29ts5jXr06IGuXbti+vTptRK4ocvPL4W7+y7cvPlcsC0t7TWWLDmPU6ee4O+/J6NJE+kfsOrg119vSW2fNMmuynYzM32px/5w/IoQongyFRBjY2NwuVyMHTsW4eHhcHV1Fblfly5d0KZNG4UEJMJCQq4JFY/3XbmSjs2bb2DBgl5KTiW7srIKPH9eIHGflBTRz+Lo6TXGgAGmOHdO/AOfH3/sVKN81fXPP4XIzi4Ej6dHE22SBkemLqxVq1YhMzMTO3fuFFs8AMDBwQEpKZKf5iby2bkzrkbt6kJLS0PqglKSFoQKDvaEjo7oKy0Hh3Yir1wUKSEhG8OG7YOh4TrY2v6Cdu3Wwd//KLKzC2v1vISoE5kKyIIFC6CvL737gNQeaWtcZGbWnTUwpk6VPCHitGmOYtscHQ0RFeWP/v07CLZpa2ti6lQHnDs3pVa78Z48yUPfvjvx55+PBWuSlJTwsXv3HfTvH4Y3b0pr7dyEqBOZurA+//xz5OTkYM+ePVXaJk+ejLZt2yI4OFhh4UhVHTu2wKNHuWLbTU3rTt//V1/1wcmTj/HgQXaVtoEDzeDvL3lGX1dXY0RFBeDp09fIzS2Cqak+WrSofjdSXl4xIiLikZdXDBubNhgyxLxaMwx8990Fsc+gJCTkYNu2m1i4UP27EQmpKZmuQI4fP45BgwaJbPPy8sLRo0cVEoqIJ61vX1l9/4qgr98EFy9Oxddf94GRUTNwuRyYm7fEjz8OxIkT49GokUa1jmNi0hyOjoYyFY+1ay+Dx/sJgYEn8OWXZ/HRR/thbr4RsbHPJL6voqISBw8+kLiPuLvHCKlvZLoCyczMRPv27UW28Xg8ZGbS07+17dNPu+HUqSciB5C9vS3w8cd1a7Erff0mWLXKA6tWeSjtnLt338EXX/xdZXtq6it4ee3FgwezxS5GVVpaIfUW4devqQuLNAwyXYHo6+vjyZMnItuePHmCZs1o+dDa1rixJk6dmohNm7zh5GSINm2awtXVCFu3DsOxY+OgqUmTPEqzZs0lsW15ecUS13PX0WkkmJ5eHDu7tnJnI6QukenTZuDAgVi1ahWysrKEtmdlZeH777+Hp2fNFgki1aOlpYHZs11x82YgsrIWIibmYwQGOlPxqIaMjDdISMiRuM/p00kS2+fMEX8HIgDMnSu5nZD6QqYurO+++w6urq6wsLDAsGHDBN1WJ06cQOPGjbFy5craykmIQnCqMaM8R8pOn33WHTExmdi//36VthUr3ODhYSZnOkLqFpkKSMeOHXHjxg0sXboUZ86cQW5uLgwMDDBy5Eh8++236NChg/SDEKJCxsZ6sLVtg/v3/xG7z+DBnSQeg8vlIDzcF1OnOiAs7A5evCiAufnbSSBpJUTSkCh9RcLIyEgEBQWhoqICM2bMwJdffinUnpaWhmnTpiE7OxstW7bE3r17wePxEBcXh1mzZuHNmzfQ0NDA4sWLMXbsWInnoskUiSj799/DhAmHRba1bq2D+/dno00b1U4ESogq1fpkivKoqKjAnDlzcOrUKcTHx2P//v2Ij48X2mfhwoWYMmUK7t69i6VLl+Krr74CAOjo6GD37t148OABIiMjMW/evHq7oBWpXePHd8XPPw+Grq7womGdO7fCmTOTqXgQUk0yL24RHR2N/fv34+nTpygpKRFq43A4OHv2rNj3xsTEwNzcHGZmb/uIx40bh2PHjsHa2lqwT3x8PEJC3q714O7uDh8fHwBA586dBfsYGRmhTZs2yM7ORosWdefBOaI+PvusOwICHHDsWCLy8ophbd0aAweaSR3/IIT8R6YrkK1bt8Ld3R0RERF49eoVGGNCX5WVlRLf/+FzJKKeHbG3t0dERAQA4MiRI8jPz0durvCT1zExMSgrK0OnTlX7qrdt2wYXFxe4uLggO7vqE86EvKOn1xiTJ9sjKKgHPD07UfEgREYyXYGsW7cOEyZMwI4dO6ClJXrNcElEDbd8+J82ODgYc+fORVhYGPr16wdjY2Noav4X8/nz55g8eTJ27doFLrdq/QsMDERgYCCAt/14hBBCaofMT6JPnTpVruIBvL3iSE9PF7zOyMiAkZHwXStGRkY4fPjtAGdBQQEiIiLQvHlzAMCbN28wdOhQrFy5Ej169JArAyGEEMWQqQvL2dkZycnJcp/M1dUVjx8/RkpKCsrKynDgwAEMHz5caJ+cnBxBV9jq1asxbdo0AEBZWRlGjhyJKVOmYMyYMXJnIIQQohgyFZANGzZg/fr1gvXPZaWpqYnQ0FB4eXnBysoKfn5+sLGxwdKlS3H8+HEAQFRUFCwtLdG5c2dkZWVh8eLFAICDBw/iwoULCAsLg4ODAxwcHBAXVzfWviCEkPpIpudA2rdvjzdv3qCgoAA6OjpV1gbhcDhIS0tTeEh50XMghBAiu+p+dso0BuLh4UF3qhBCCAEgYwEJCwurpRiEEELqGpq+lRBCiFxkLiC3b9+Gr68vDAwMoKmpiVu3bgEAvv76a0RGRio8ICGEEPUkUwG5dOkSevbsicTEREyYMEHoyXMul4stW7YoPCAhhBD1JFMB+fLLL+Hl5YUHDx7gp59+EmpzcnISXI0QQgip/2QaRL916xYOHz4MDodT5W4sAwMDmnuKEEIaEJmuQLS1tVFUVCSy7fnz54IpRwghhNR/MhWQPn36YP369aioqBBse3clsn37dgwYMECx6QghhKgtmddE7927N+zt7TF69GhwOBzs2rUL8+fPx82bN3Hjxo3aykkIIUTNyHQFYm9vjwsXLqBt27ZYtWoVGGMIDQ0F8HahKUtLy1oJSQghRP3IvCKhk5MTzp49i5KSEuTl5aFFixbQ0dGpjWyEEELUmMwF5B1tbe0qa3kQQghpOKQWkBUrVmDGjBkwMjLCihUrJO7L4XCwZMkShYUjhBB1Usnng6sp9+/d9Y7U6dy5XC6uXbuGbt26iVxCVuhgHI7QHVqqRtO5E0JqqqK8HNd//hmxv/yCl8nJ0G7RAnaTJ6Pv119Dt107VcerFQqbzv396Ure/zMhhNR3lRUVODhqFB793/8JtpW8eoWYjRvx8PhxTL9yBc0acFc+zcZLCCFiPDh4UKh4vO91WhrOL12q5ETqRaYCcuLECcFtux/atGkTTp48qZBQhBCiDu5IWQPp/v794JeWKieMGpKpgHz33XcoLCwU2VZcXIzvvvtOIaEIIUQd5D9/LrG9vKgIpa9fKymN+pGpgCQmJsLJyUlkm4ODAxISEhQSihBSc6yyEiWvX6NSjW5sqWv0zcwktmvr60NbX19JadSPTAWksrISBQUFItvy8/NRXl6ukFCEEPmVvH6NM198gbVt2uCHFi2w1sAAp4KCUJSbq+podY5zYKDEdoeAAGg0aqSkNOpH5qlMwsPDRbaFh4fDzs5OIaEIIfIpzc/HLnd3XFm7FsX/FoySV68Qs2EDdvbpg+K8PBUnfKv45Us8v3ULr9LSVB1FIgtvbzjPnCmyzdDJCf2XLVNyIvUiUwFZsGABDh8+jDFjxuD06dOIj4/HmTNnMGbMGBw5cgSLFi2qrZyEkGq4sWkTXty+LbItJzERV9atU3IiYcUvX+LY1Kn4ycgI25yd8XPHjtjZty8yrl9XaS5Jhm3ZgjGHDsHUwwPNjI3RzsEBg376CQEXLkC7gS9hIfVBwg9t3LgRixcvFgymM8agq6uL1atXY86cObUSUl70ICFpaDZZWSEnMVFsezNjY8zPyFBiov+UFxdjR+/eIgtcIx0dBFy4ACNnZ8Wft6gISWfOgF9cDOPu3aFvaqrwc9Q3CnuQ8EOffvopAgICcOXKFeTm5sLAwAC9evWCrq6uXEEJIYpT+M8/ktuzspSUpKq7e/eKvToqLypC1NKlmPDnnwo95+W1a3Fx1SrBnVIcLheWw4dj+I4daNKAB78VRa5JXZo1awYvLy9FZyGE1FBLCwtkSugOatW5sxLTCHtw4IDE9ieRkSh5/Vph3ULX1q/H3198IbSNVVYi8ehRFLx4gWmXL4MjZXomIpnUAnLhwgU4OTlBV1cXFy5ckHrAfv36KSQYIUR2zjNnSiwg4gaElaH0zRuJ7ayyEmUFBQopIBVlZbj4/fdi2zOuXcOTyEhYeHvX+FwNmdQC4ubmJphM0c3NTbCE7YcYY2o3mSIhDY2Dvz+Sz5zB/f37q7R1/ugjuM6erYJUb7VzdMQzCf3quoaG0G3bViHnSr96FUXZ2RL3STx2jApIDUktIOfPn4e1tbXgz4QQ9cXhcuEbHo4uPj649dtveJWaCj0eD47TpsF2/HhwNTRUls11zhzc3r4dTMykrC6ffKKwqdIrysqk7lNZjX2IZFL/to4dOwYzMzPo6uqCw+EIurMIIeqJw+HAxs8PNn5+qo4ipJ29PYb+8gv+nDWrShHp4uODPl99pbBzGbu6opGODsqLisTu06F/f4Wdr6GSOoL0888/4/m/88G4u7sjPj6+RieMjIyEpaUlzM3NsWbNmirtaWlp8PDwgJ2dHdzc3JDx3i2Hu3btgoWFBSwsLLBr164a5SCEKJ9zYCDmJCai54IFMB8yBHaTJmHSX3/B7/BhhT7Rrd2iBRymTRPb3szYGLbjxinsfA0Wk6JNmzZs7969jDHGOBwOi4mJkfYWsfh8PjMzM2NJSUmstLSU2dnZsQcPHgjtM3r0aBYWFsYYY+zs2bNs0qRJjDHGcnNzmampKcvNzWV5eXnM1NSU5eXlSTyfs7Oz3FkJIXVbeUkJOzh6NFsOCH2FdOjAsu7fV3U8tVbdz06pXViDBg3C1KlTsXjxYgCAj48PGjduLHJfDoeDpKQksceKiYmBubk5zP6doGzcuHE4duyYYIwFAOLj4xESEgLg7RWPj48PAOCvv/6Cp6cnWrZsCQDw9PREZGQkxo8fL+1bIIQ0QJqNG2PMoUN4fvs2EiIiUF5cDF6PHuji49Og569SJKkFZOvWrejatSsSEhKwe/du2NnZwcDAQK6TZWZmon379oLXPB4P1z+45dDe3h4REREICgrCkSNHkJ+fj9zcXJHvzczMrHKObdu2Ydu2bQCAbCl3YRBC6j9DR0cYOjqqOka9JLWA6Ojo4It/H8bZtWsXvv32W3Tr1k2ukzERs6Z8eFtwcHAw5s6di7CwMPTr1w/GxsbQ1NSs1nsBIDAwEIH/zqDp4uIiV05CCCHSSR1Eb9myJW7dugUA8Pf3R+vWreU+GY/HQ3p6uuB1RkYGjD5YT9jIyAiHDx/G7du3sWrVKgBA8+bNq/VeQgghyiO1gBQWFqL03yUbd+/eXaNuIVdXVzx+/BgpKSkoKyvDgQMHMHz4cKF9cnJyUPnvLX6rV6/GtH/vpPDy8sLp06fx8uVLvHz5EqdPn6bpVAghRIWkdmF16NABv/76K0pLS8EYw+3bt1FSUiJ2f0lTmWhqaiI0NBReXl6oqKjAtGnTYGNjg6VLl8LFxQXDhw9HVFQUvvrqK3A4HPTr1w+bNm0C8PZKaMmSJXB1dQUALF26VDCgTgghRPmkTue+Y8cOzJw5U3BVIA5Tw6lMaDp3QgiRncKmc582bRqGDBmCR48ewd3dHRs2bICVlZVCQhJCSH1Smp8PzcaNoaGlpeooSlGtiWcMDQ1haGgIf39/DB06FKa0IAshhAB42/sS+8svuL5hA3IfPgRXUxOWI0ag35IlaGdvr+p4tUrmFQnfKSgoQG5uLoyMjNBITR/KoS4sQkhtOzFrFm5u2VJle6OmTTHl7FnwundXQaqaqe5np8yrqZw4cQJOTk5o3rw5OnXqhHv37gEAZsyYgX379smelBBC6qjMGzdEFg8AKC8sRGRQkJITKZdMBeTo0aMYMWIEDAwM8MMPPwgNrJuamtIEh4SQBuWOlM+8zOvXkfPwoZLSKJ9MBeTbb7/F1KlTcfr0acybN0+ozdbWFvfv31doOEIIUWfS1qAHVLsOfW2TqYAkJCRg7NixAKpOI6Kvr4/c3FzFJSOEEDXX0sJCYjuHy4V+p05KSqN8MhUQPT095OTkiGxLTU2t0TQnhBBS1zjNmAEOV/zHaOdhw6BnbKzERMolUwHx9PTE6tWr8erVK8E2DoeD0tJShIaGYsiQIQoPSAgh6krf1BTemzYBIiZ2bfGurR6TaQHiVatWoVu3brC0tIS3tzc4HA7WrFmDu3fv4vXr1zh69Ght5SSEELXk8sknaOfoiJiNG/Hsxg1o6erCeswYOM+ciSb6+qqOV6tkKiAdO3bErVu3sGzZMvz111/Q0NDAhQsXMHjwYKxYsYJmxyWEyOzZzZtIPHIE/NJSmPTpg87DhoGroaHqWDLhde9eJ5/3qCm5HySsC+hBQkLUV3lxMSLGjcPD48eFtre0sMDEkyfR0txcRclIrT1I+L7du3fj5cuXNTkEIaSBOjl3bpXiAQB5jx9j7+DBqOTzVZCKyELuAlJRUYGpU6ciJSVFkXkIIQ1AQVYW7u3dK7b9ZVISEmlMVe3V6AqkHvd+EUJqUcbVq6goK5O4T2p0tJLSEHnVqIAQQog8uNWYgFVDTSdpJf+Ru4DMQl5DAAAgAElEQVRwOBx06NABjRs3VmQeQkgD0NHNDY319CTuY/nBctdE/ch0G29JSQliY2Px/PlzcDgc7N27F53q8WP6hJDaodW0KXouXIiopUtFtnfo3x8d3dyUG4rIrFoFpLS0FF988YXQ2ujA26sQbW1tzJo1C99//z20GsgqXISQmuv3zTdglZW4GhyMsoICAG/njrLy9cVHv/2m4nSkOqQWEMYYhg0bhnPnzmHEiBHw9vaGiYkJGGNIT0/HiRMnEBISgvj4eJw8eVIZmQkh9QCHw4HbsmXoOX8+ks+cAb+0FO179kSLjh1r5XwZ167h8alTYJWV6OTpiQ79+tXKeRoSqQ8SHjp0COPHj8ehQ4cwcuRIkfscPnwYfn5+OHjwIHx9fWslqDzoQUJCSMmrVzg4ahRSzp0T2s7r0QPjjh9HU5oEtorqfnZKLSC+vr7Q1taWutrg+PHjUVZWhoiICNmS1iJRP4Ty8nJkZGSgpKRERamIKmlra4PH46ntMsxE8fZ/9BEenTghss2kb19MvXBByYnUX3ULiNQurNu3b2PlypVSDzRs2DB888031UunQhkZGWjWrBk6duxYZU0TUr8xxpCbm4uMjAyYmpqqOg5RguyEBLHFAwCeXryIjOvXG+Q8Voog9Tbe7OxsmJiYSD2QiYkJ/qnG6lyqVlJSglatWlHxaIA4HA5atWpFV58NSGpUlPR9zp+v/SD1lNQCUlRUVK1nPbS0tOrMf0wqHg0X/d03LFxN6TeaVmcfIlq1fnKZmZlITk6WuE9GRoZCAhFCiKJYDBkCjoYGWEWF2H06f/SREhPVL9V6En306NGwsLCQ+DVmzJjazlpvcDgcLFiwQPA6ODgYy5cvV12gavr+++9r9P6wsDDMnTtX5vc9e/YMo0ePrtG5ScOkx+PBYepUse02fn4wsLRUYiLFKCsowKMTJ5B49CgKsrJUlkPqFcjOnTuVkUPtlZaWKmzalsaNG+Pw4cP46quvYGBgoJBjKkJFRQU0JCzk8/333+Prr79WYqK3jIyM8Mcffyj9vKR+GLp5M7iamri9fTsqy8sBABwNDdhNnIihW7aoOJ1sGGOIXrECV9etQ1l+PoC384rZTZ4M79BQNGrSRKl5pBYQf39/ZeRQWy9evMDhw4eRnZ2N1q1bw9fXF+3atavRMTU1NREYGIiQkBCsWrVKqC0tLQ3Tpk0TnG/nzp0wMTFBQEAA9PT0EBsbixcvXuDHH38U+Vt5QEAAtLW18eDBA2RlZeGnn37CsGHDkJqaismTJ6OwsBAAEBoail69eiEqKgrffvstDA0NERcXh/j4eOzduxcbNmxAWVkZunfvjs2bN2Px4sUoLi6Gg4MDbGxsEB4ejp9++gk7duwAAMyYMQPz5s2rkmfnzp1YvXo1DA0N0blzZ0ERzs7OxieffIKnT58CANavX4/evXsjOjoaQUFBAN5eqV24cAG5ubkYNmwY7t+/j6KiIgQEBCAxMRFWVlZITU3Fpk2b4OLiAl1dXQQFBeHEiRNo0qQJjh07hrZt29bo74rUfRqNGmHYL7/AbflyJP/9N1hlJUzd3aHH46k6msyili3Dhe++E9pWWV6OuB07UJKXh7FHjig3EKvHnJ2dq2yLj4+X6RibNm1iy5cvF3xt2rSpxrmaNm3KXr9+zTp06MBevXrF1q5dy5YtW8YYY2zYsGEsLCyMMcbY9u3b2YgRIxhjjPn7+7PRo0eziooK9uDBA9apUyeRx/b392deXl6soqKCPXr0iBkbG7Pi4mJWWFjIiouLGWOMPXr0SPCzOX/+PNPR0WHJycmMsbc/n2HDhrGysjLGGGOzZs1iu3btEuR+JzY2ltna2rKCggKWn5/PrK2t2a1bt4SyPHv2jLVv3579888/rLS0lPXq1YvNmTOHMcbY+PHj2cWLFxljjKWlpbEuXboIvv9Lly4xxhjLz89n5eXlLCUlhdnY2DDGGFu7di0LDAxkjDF27949pqGhwW7cuMEYYwwAO378OGOMsUWLFrHvvvtO5M9I1n8DhKiD4lev2KqmTdlyQOzX89u3FXIuUZ+dotDtBxKUlpYiOztbaFt2drZCurP09PQwZcoUbNiwAU3eu+y8evUqDh8+DACYPHkyvvjiC0Gbj48PuFwurK2tkSWh39PPzw9cLhcWFhYwMzNDYmIiTE1NMXfuXMTFxUFDQwOPHj0S7N+tWzfBcxFnz57FzZs34erqCgAoLi5GmzZtqpzj0qVLGDlyJJo2bQrg7QOnFy9ehKOjo2Cf69evw83NDa3/fdJ37NixgvP+/fffiI+PF+z75s0b5Ofno3fv3pg/fz4mTpwIX19f8D74LfHSpUuCKxRbW1vY2dkJ2rS0tDBs2DAAgLOzM86cOSP2Z0RIXZN85gzK/+1BECfhyBG0c3BQUiIVrAcSGRkJS0tLmJubY82aNVXanz59Cnd3dzg6OsLOzk4wv1Z5eTn8/f3RtWtXWFlZYfXq1bWetXHjxoIPv3dat26tsLGQefPmYfv27YJuJVHev+30/fOyfycQWLx4MRwcHODw3j+aD29V5XA4CAkJQdu2bXHnzh3Exsai7L3FfN4VgXfH9ff3R1xcHOLi4vDw4UORA/ysmouJibtttrKyElevXhWcJzMzE82aNcOXX36J3377DcXFxejRowcSExOrfd5GjRoJzqehoQE+LYlK6hF+aanUfSqqsY8iKbWAVFRUYM6cOTh16hTi4+Oxf/9+od9CAWDlypXw8/PD7du3ceDAAcyePRvA2zm5SktLce/ePdy8eRNbt25FampqrWf29fUVFJF3YyCK0rJlS/j5+WH79u2Cbb169cKBAwcAAOHh4ejTp4/EY6xatUrwIfzOoUOHUFlZiaSkJCQnJ8PS0hKvX7+GoaEhuFwu9uzZgwoxtzV6eHjgjz/+EDwUmpeXh7S0NABvP6DL/x2E7NevH44ePYqioiIUFhbiyJEj6Nu3r9CxunfvjqioKOTm5qK8vByHDh0StA0aNAihoaGC1+/yJyUloWvXrvjf//4HFxeXKgWkT58+OHjwIAAgPj4e9+7dk/jzIaS+aN+rFzhcyR/ZJlI+LxRNqV1YMTExMDc3h5mZGQBg3LhxOHbsGKytrQX7cDgcvHnzBgDw+vVrGBkZCbYXFhaCz+ejuLgYWlpa0JOyII0itGvXDrNnz1boXVjvW7BggdAH6YYNGzBt2jSsXbtWMIguK0tLS/Tv3x9ZWVnYsmULtLW1MXv2bIwaNQqHDh2Cu7u70FXH+6ytrbFy5UoMGjQIlZWVaNSoETZt2oQOHTogMDAQdnZ2cHJyQnh4OAICAtCtWzcAbwfR3+++AgBDQ0MsX74cPXv2hKGhIZycnASFa8OGDZgzZw7s7OzA5/PRr18/bNmyBevXr8f58+ehoaEBa2trDBkyBM+fPxccc/bs2fD394ednZ3gKrV58+Yy/4wIqWv0TU3RxccHCf92cX/IoEsXWHh7KzeUQkZcqunQoUNs+vTpgte7d+8WDKq+8+zZM2Zra8uMjY1ZixYtWGxsLGOMsbKyMjZ27FhmYGDAdHR02NatW0WeY+vWrczZ2Zk5OzszExOTKu31fQDV39+fHTp0SNUxag2fzxfcDPDkyRPWoUMHVlpaKtMx6vu/AVJ/Fb96xXb2719l8HyDhQXLffJEYedRy0F0JqL/+sM+8v379yMgIAALFizA1atXMXnyZNy/fx8xMTHQ0NDAs2fP8PLlS/Tt2xcDBw4UXM28ExgYiMDAQABvZ5Qk9UtRURHc3d1RXl4Oxhh++eUXWsiMNBjazZsjICoKqVFRePh//4fK8nJ0dHOD5fDhKpmSRaln5PF4SE9PF7zOyMgQdFG9s337dkRGRgIAevbsiZKSEuTk5GDfvn0YPHgwGjVqhDZt2qB3796IjY2tUkAaurCwMFVHqFXNmjWjNV5Ig9fRzU0tlvxV6iC6q6srHj9+jJSUFJSVleHAgQMYPny40D4mJiY4e/YsACAhIQElJSVo3bo1TExMcO7cOTDGUFhYiGvXrqFLly7KjE8IIeQ9Si0gmpqaCA0NhZeXF6ysrODn5wcbGxssXboUx48fBwCsW7cOv/76K+zt7TF+/HiEhYWBw+Fgzpw5KCgogK2tLVxdXTF16lShZwAIIYQol9QVCesyUatqJSQkwMrKSkWJiDqgfwOESFbdFQmV/iAhIYSQ+oEKCCGEAKh4b3YGUj1UQKrh7t0s/PXXEzx6lKvqKIQQBSrKzcXphQvxo4EBVjZujJ94PER9+y3Ki4tVHa1OoAIiQUxMJlxctsHefgsGDw6HpWUo+vcPQ2JiTo2PnZGRgREjRsDCwgKdOnVCUFCQ0PxUorx69QqbN28W2tarVy+Zzivr/uJs3boVhoaGcHBwgL29PcaMGYOUlBSFHJsQZSjKzcXOvn1xdd06FOe+/eUwPzMT0cuXY++gQeDXkSW6VYkKiBjx8dkYOHA3bt58LrT9woU0uLmF4dmzfLmPzRiDr68vfHx88PjxYzx69AgFBQVYvHixxPeJKiBXrlyR6dyy7i/O3bt3sWLFCsTFxeHOnTvw8PCAr69vtSdZJETVLq1Zg5yEBJFtTy9dwq3fflNyorqHCogYq1dfQn6+6CuCrKxCbNhwXe5jnzt3Dtra2pj671KbGhoaCAkJwY4dOxAfH48uXboI5nsaPXo0ioqKAABffvklkpKS4ODggEWLFgEAdHV1kZqaii5dumDGjBmwtbXFxIkT8ffff6N3796wsLBATEyM4Ny6uroAgC1btghm8TU1NYW7uzsAYO/evejWrRscHBwwc+ZMsZMu3rt3D7a2toLXn3zyCV68eCH0oCgh6uyOlIdu4+r5Q7mKQAVEjIiIeIntf/whuV2SBw8ewNnZWWibnp4eTExMwOfz8fDhQwQGBuLu3bvQ09MTXHWsWbMGnTp1QlxcHNauXSv0/idPniAoKAh3795FYmIi9u3bh0uXLiE4OFjkWuaffPIJ4uLicOPGDfB4PMyfPx8JCQn4/fffcfnyZcG6IeHh4SK/h/v378PGxkZoW5MmTfDy5Uu5fy6EKEtlRQWKciR3RRc8fy6xnVABEYkxhpISyWtJFBfLv9YEY0zkOhnvtrdv3x69e/cGAEyaNAmXLl2SekxTU1N07doVXC4XNjY28PDwAIfDQdeuXSVOex8UFIQBAwbgo48+ElpMysHBAWfPnkVycnKV96Snp6NZs2ZCsyGXl5fj+fPnNLUMqRO4GhpobmIicR/9Tp2UlKbuogIiAofDQbduxhL36dFD/vWUbWxsqjyk8+bNG6Snp0NDQ0PkglDSvD/VPJfLFbzmcrliF1YKCwtDWloali1bBqD6i0ndvXu3ytXHzp07MWDAADRr1kxqVkLUgdO/k66K4yylnVABEWvevB5i2zgcICiou9zH9vDwQFFREXbv3g3g7UJbCxYsQEBAAHR0dPD06VNcvXoVwNvZid8tKtWsWTPk58s/eP++mzdvIjg4GHv37gX330VqJC0m9b4Pxz9Onz6N1atXIzg4WCHZCFGGXgsWiJ2Q0HbcOHSdMEG5geogKiBijBtniyVL+uHDX/41NDjYuHEI+vXrIPexORwOjhw5gkOHDsHCwgKdO3eGtra2YKzCysoKu3btgp2dHfLy8jBr1iwAQKtWrdC7d2/Y2toKBtHlFRoairy8PLi7u8PBwQEzZswQWkzKzs4Onp6eQos5vXPv3j2Eh4fD2dkZTk5O2LVrFyIjI2l6EFKnaGprY2JkJIZu2QLj7t2h1749OvTvD9/wcPiGh0td/Y/QXFhSPXyYg7CwODx7VgBT0xaYNs0RJia1twJeamoqhg0bhvv379faORo6mguLEMmqOxeW8lcgqWMsLQ2wevVAVccghBC1Q9doaqZjx4509UEIqROogBBCCJELdWERQkgdVl5cjPv79yP577/B4XBgPmQIrMeMgeZ7t/bXFioghBBSR+U+fow9Awfi9dOngm339u3Dhe++w+QzZ6Q+LFlT1IVFCCF1EGMMB0aMECoe7+Q+eoQ/xo6t9QxUQAghpA5KPnNG7GzCAJBx7Roy35tItTZQASGEkDoo88YNhexTE1RACCGkDtJq2lT6Pv8u31BbqIBUQ9bdu3jy11/IffRI1VEIIQQAYOXrC46Ghth2TW1tWH70Ua1moAIiQWZMDLa5uGCLvT3CBw9GqKUlwvr3R05iYo2PTUvavpWeng53d3dYWVnBxsYGP//8c62ch5D6prmJCVznzBHb3mvRIjRp2bJWM1ABESM7Ph67Bw7E85s3hbanXbiAMDc35D97JvexG9qStlFRUQgICBB5HE1NTaxbtw4JCQm4du0aNm3ahPj4/xbroqVzCRFvcEgI3L79VqhQNG3bFp5r18J9xYpaPz8VEDEurV6NMjFTpxdmZeH6hg1yH5uWtP2PoaEhnJycALydrt7KygqZmZkKPw8h9RGHy0X/pUsxPzMT0y5fxvSrV/H506fotXChUs5PBUSM+IgIye1//CH3sWlJW9FSU1Nx+/ZtdO/+31ortHQuIdJpamujfa9e4PXoAQ0tLeWdV2lnqkMYY+CXlEjch19cXKPjy7Kk7YYNG7BQym8U75a0BSD3krahoaGCJW0BoLi4GG3atKnynuouadu9e3eUlpaioKAAeXl5cHBwAAD88MMP8PLyEjpmQUEBRo0ahfXr1wuOS0vnEmVLOX8eDw4eRFl+PgydnOAQEFDr4wh1GRUQETgcDoy7dUPm9eti9+H1EL9ioTQ2NjaI+OAKR5VL2oaGhgL4b0nb1atXSzxXdZe0vf7vzy8qKgphYWEICwsTebzy8nKMGjUKEydOhK+vb7XPk56ejhUrVqB58+YYPHgwBg6kafeJfPilpTg4ahQe//mnYNu98HBELV+OsYcPw4z+bYlEXVhi9Jg3T3wjh4PuQUFyH5uWtP0PYwzTp0+HlZUV5s+fL9N5EhMToaWlhc8++4yKB6mRv//3P6Hi8U5Zfj5+HzkShf/+nyDCqICIYTtuHPotWYIP17TlaGhgyMaN6NCvn9zHpiVt/3P58mXs2bMH586dEwzqnzx5slrn8fT0xKeffoq5c+cKDbwTIouyggLc3rFDYvut7duVmKjuUPqStpGRkQgKCkJFRQVmzJiBL7/8Uqj96dOn8Pf3x6tXr1BRUYE1a9bA29sbwNsujZkzZ+LNmzfgcrm4ceMGtLW1xZ5LEUva5jx8iLiwMBQ8e4YWpqZwnDatVme4pCVtq+9///sfKioqUFJSgpCQEDRq1Kha76Mlbcn7MmNi8Nt7N26I0sXHB2OPHFFSItVTyyVtKyoqMGfOHJw5cwY8Hg+urq4YPnw4rK2tBfusXLkSfn5+mDVrFuLj4+Ht7Y3U1FTw+XxMmjQJe/bsgb29PXJzc6v9gVETBpaWGChlTICoxg8//KDqCKQeaKSjo5B9GiKldmHFxMTA3NwcZmZm0NLSwrhx43Ds2DGhfTgcDt68eQMAeP36NYyMjAC87f+2s7ODvb09gLfdORoSHuOvq2hJW0KUq42tLQy6dJG4j9Xo0UpKU7cotYBkZmaiffv2gtc8Hq9K3/Xy5cuxd+9e8Hg8eHt7Y+PGjQCAR48egcPhwMvLC05OTvjxxx+VGZ0QUo8NWLWqynjnO7wePWA5fLiSE9UNSi0gooZbPrxFdf/+/QgICEBGRgZOnjyJyZMno7KyEnw+H5cuXUJ4eDguXbqEI0eO4OzZs1WOt23bNri4uMDFxQXZ2dm19r0QQuoPK19fjP79d+h36iTYpqGlBbvJkzExMhLcetjboQhKHQPh8XhCU1BkZGQIuqje2b59OyIjIwEAPXv2RElJCXJycsDj8dC/f38YGBgAALy9vXHr1i14eHgIvT8wMBCBgYEA3g4EiSLuQT5S/9EcWkQcmzFjYD16NJ7FxqIsPx9tbG3RVMSDtOQ/Sr0CcXV1xePHj5GSkoKysjIcOHAAwz+4NDQxMRFcWSQkJKCkpAStW7eGl5cX7t69i6KiIvD5fERHRwsNvleXtrY2cnNz6YOkAWKMITc3V+Kde6Rh43A4MHZ1hemAAVQ8qkGpVyCampoIDQ2Fl5cXKioqMG3aNNjY2GDp0qVwcXHB8OHDsW7dOnz88ccICQkBh8NBWFgYOBwO9PX1MX/+fLi6uoLD4cDb2xtDhw6VOQOPx0NGRgZ1bzVQ2tra4PF4qo5BSL2g9OdAlKm69zITQgj5T3U/O+lJdEIIIXKhAkIIIUQuVEAIIYTIpV6PgRgYGKBjx46qjgEAyM7ORuvWrVUdQyR1zgaodz51zgZQvppQ52xA7eZLTU1FTk6O1P3qdQFRJ+o8oK/O2QD1zqfO2QDKVxPqnA1Qj3zUhUUIIUQuVEAIIYTIRWP58uXLVR2ioXB2dlZ1BLHUORug3vnUORtA+WpCnbMBqs9HYyCEEELkQl1YhBBC5EIFhBBCiFyogNRQZGQkLC0tYW5ujjVr1lRpf/r0Kdzd3eHo6Ag7OzucPHmySruuri6Cg4PVLt/du3fRs2dP2NjYoGvXrigpKVGLbOXl5fD390fXrl1hZWWF1bW05LC0fGlpafDw8ICdnR3c3NyQkZEhaNu1axcsLCxgYWGBXbt2qVW+uLg4wd+rnZ0dfv/9d7XJ9s6bN29gbGyMuXPnKjxbTfM9ffoUgwYNgpWVFaytrZGamqo22b744gvY2NjAysoKn332We3POs6I3Ph8PjMzM2NJSUmstLSU2dnZsQcPHgjt8/HHH7PNmzczxhh78OAB69Chg1C7r68vGz16NFu7dq1a5SsvL2ddu3ZlcXFxjDHGcnJyGJ/PV4ts4eHhbOzYsYwxxgoLC1mHDh1YSkqKwrJVN9/o0aNZWFgYY4yxs2fPskmTJjHGGMvNzWWmpqYsNzeX5eXlMVNTU5aXl6c2+R4+fMgePXrEGGMsMzOTtWvXjr18+VItsr3z2WefsfHjx7M5c+YoLJei8vXv35+dPn2aMcZYfn4+KywsVItsly9fZr169WJ8Pp/x+XzWo0cPdv78eYVlE4WuQGqgJmu8A8DRo0dhZmYGGxsbtctX22vQ1yQbh8NBYWEh+Hw+iouLoaWlBT09PYVlq26++Ph4wYJm7u7ugva//voLnp6eaNmyJfT19eHp6SlYJE0d8nXu3BkWFhYAACMjI7Rp00ahyxvUJBsA3Lx5E1lZWRg0aJDCMikqX3x8PPh8Pjw9PQEAurq60NHRUYtsHA4HJSUlKCsrQ2lpKcrLy9G2bVuFZROFCkgN1GSN98LCQvzwww9YtmyZWuar7TXoa5Jt9OjRaNq0KQwNDWFiYoKFCxeiZcuWSs9nb2+PiIgIAMCRI0eQn5+P3Nzcar1XlfneFxMTg7KyMnR6bylXVWarrKzEggULsHbtWoXlUWS+R48eoUWLFvD19YWjoyMWLVqEiooKtcjWs2dPuLu7w9DQEIaGhvDy8oKVlZXCsolCBaQGWA3WeF+2bBk+//xz6OrqqmW+6q5Br4psMTEx0NDQwLNnz5CSkoJ169YhOTlZYdmqmy84OBjR0dFwdHREdHQ0jI2NoampWa33qjLfO8+fP8fkyZOxc+dOcLmK+yioSbbNmzfD29tb6ENU0WqSj8/n4+LFiwgODsaNGzeQnJyMsLAwtcj25MkTJCQkICMjA5mZmTh37hwuXLigsGyiKHVFwvqmJmu8X79+HX/88Qe++OILvHr1ClwuF9ra2godNFTGGvSqyLZv3z4MHjwYjRo1Qps2bdC7d2/ExsbCzMxMIdmqm8/IyAiHDx8GABQUFCAiIgLNmzcHj8dDVFSU0Hvd3NwUlq2m+YC3g9RDhw7FypUr0aNHD7XJdvXqVVy8eBGbN29GQUEBysrKoKurK3IwWRX5eDweHB0dBf/WfHx8cO3aNUyfPl3l2bZt24YePXoIfikdMmQIrl27hn79+ikkm0i1OsJSz5WXlzNTU1OWnJwsGPC6f/++0D6DBw9mO3fuZIwxFh8fzwwNDVllZaXQPsuWLauVQfSa5MvLy2OOjo6ssLCQlZeXMw8PD3bixAm1yLZmzRoWEBDAKisrWUFBAbOysmJ37txRWLbq5svOzmYVFRWMMca+/vprtmTJEsbY20H0jh07sry8PJaXl8c6duzIcnNz1SZfaWkpGzBgAAsJCVFoJkVke9/OnTtrZRC9Jvn4fD6zs7Nj//zzD2OMsYCAABYaGqoW2Q4cOMA8PDxYeXk5KysrYwMGDGDHjx9XWDZRqIDU0J9//sksLCyYmZkZW7lyJWOMsSVLlrBjx44xxt7ePdSrVy9mZ2fH7O3t2V9//VXlGLVVQGqab8+ePcza2prZ2NiwRYsWqU22/Px8Nnr0aGZtbc2srKzYjz/+qPBs1cl36NAhZm5uziwsLNj06dNZSUmJ4L3bt29nnTp1Yp06dWI7duxQq3x79uxhmpqazN7eXvB1+/Zttcj2vtoqIDXNd/r0ada1a1dma2vL/P39WWlpqVpk4/P5LDAwkHXp0oVZWVmxzz//XKG5RKGpTAghhMiFBtEJIYTIhQoIIYQQuVABIYQQIhcqIIQQQuRCBYQQQohcqIAQQgiRCxUQQgghcqECQuqkX3/9FRwOR/DVuHFjtG/fHn5+frU+/8/7duzYAQsLC2hpaaFFixZity1fvlzm+bDkeQ8hykQPEpI6ac6cOdixYwfOnz8PACgtLUVSUhJ27tyJy5cvIzg4GPPnz6/VDM+ePYOJiQkmTpyIjz/+GNra2jAyMqqyzcXFBRkZGcjIyJBp3il53kOIMlEBIXVS7969UVZWhhs3bghtZ4zBw8MDFy9exMOHDxU6weKHoqOj4ebmhrNnz2LAgAFitxFSX1EXFqlzGGO4e/euYLGr93E4HCxatAh8Ph8nTpyQ+xx37tzB8OHDoa+vjyZNmqB37964ePGioD0gIEAww66Hhwc4HI7YbYDo7q/ibkgAAAUFSURBVKg7d+5g5MiRaNWqFZo0aQJLS0uh5XnFdWFJy/b+ex8/foyhQ4dCV1cXHTp0wIoVK1BZWVmtHH/88Qc4HA7u3LlTJYObmxt69uxZ7Z8nqZ+ogJA65/HjxygoKBBZQADA3NwcAORexOnWrVvo1asX8vLy8OuvvyIiIgKtWrXCwIEDcfPmTQDAkiVLsGHDBgDApk2bcPXqVbHbRImJiUHPnj2RlJSEkJAQ/Pnnn5g/f36VtcHlyfa+kSNHYsCAATh69Ch8fHywbNkyoTXaJeXw8fGBkZERtm7dKnTMhw8fIjo6GjNnzqz+D5XUT7U+XSMhCvb7778zACw6Olpk+/379xkAtmLFCrmOP2DAANalSxehWVb5fD7r0qULGzFihGDbmTNnGAChdadFbWPs7YzL7/9369u3L+PxeBLX0/7wPbJke/feD2cCtrW1ZZ6entXOsWzZMqanp8cKCgoE2z7//HPWokULVlRUJDY7aRjoCoTUObdv3wYA2NnZiWx/8uQJAMDS0lLsMRhj4PP5Ql8AUFxcjOjoaIwZMwZcLlfQxhjDwIEDFXKHV1FRES5fvoyJEyfKtJ62PNmGDh0q9NrW1hZPnz6tdo7AwEAUFRVh//79AICSkhLs2rULU6ZMQZMmTaqdndRPVEBInRMXFwcTExPBLbIfOnXqFBo1aiRx9cTo6Gg0atRI6AsA8vLyUFFRge+++65Ke2hoKF6+fFllDEFW747B4/Fkep882T5cK75x48YoKSmpdg4jIyOMGDECW7ZsAQAcOnQIeXl51H1FANCStqQOiouLg6urq8i25ORk7NmzBxMnTkSrVq0AAOfPn8fixYtRUlKCgoICLFy4EOPHj69yBxcAtGjRAlwuF3PmzMGUKVNEnqOm64fr6+uDy+XKPEaj6GzVzTF79mx4eHjg5s2b2Lp1K/r27Qtra2uZspP6iQoIqVOysrLw4sULkQPo6enpGD58OPT09ITW0J4wYQJiY2NhbGyMyspKvHnzBs2aNYOLi0uVYzRt2hR9+/bFnTt34OTkVONiIYqOjg769OmDvXv3YunSpdXuClJ0turmGDBgAKysrDB//nxcvnwZ4eHhNTovqT+ogJA65d34R6NGjXDt2jXw+XxkZWXh/PnzCAsLg7GxMc6dO4e2bdsK3mNoaIjZs2dj3Lhx+Oijj8R2fb3z008/oV+/fvDy8sL06dNhaGiInJwc3Lp1CxUVFULFSV7BwcHo378/evbsiQULFoDH4yE5ORlxcXHYuHGj0rJVN8cnn3yCoKAgGBgYYNSoUXJ/36R+oTEQUqfExcUBAJYtW4aePXvCy8sL8+fPR2ZmJkJDQ3H//n1YWVkJvScmJgZz585FVFQULC0tUVRUJPEcTk5OuHHjBlq1aoXPPvsMgwYNQlBQEO7du4d+/fop5PtwdXXF5cuX0b59e3z66afw9vbG2rVrpY6LKDpbdXOMGTMGwNvnXxo3bizzeUj9RE+ik3otMTERnTt3BpfLRVpaGpycnPDs2TP6EJTRr7/+ipkzZ+LRo0eC52wIoS4sUq+FhITg/PnzaNq0KZo0aYKDBw9S8ZBBfHw8kpKSsGzZMvj4+FDxIELoCoQQIpabmxuuXLmCXr16Yd++fTAyMlJ1JKJGqIAQQgiRCw2iE0IIkQsVEEIIIXKhAkIIIUQuVEAIIYTIhQoIIYQQuVABIYQQIhcqIIQQQuTy/5KudC9In4lVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def combineEfficiencyData(lst):\n",
    "    data = np.zeros((0, 4))\n",
    "\n",
    "    for jj, dds in enumerate(lst):\n",
    "        dds_index = np.hstack((dds, jj * np.ones((len(dds), 1))))\n",
    "        data = np.vstack((data, dds_index))\n",
    "    return data\n",
    "\n",
    "\n",
    "design_efficiencies_combined = combineEfficiencyData([design_efficiencies0, design_efficiencies])\n",
    "plot_handles = oapackage.generateDscatter(\n",
    "    design_efficiencies_combined, ndata=3, lbls=[\"Optimize $D$\", \"Optimize $D+2D_s$\"], verbose=0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}