mogp example in Python

mogp-example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.integrate import solve_ivp\n",
    "import mogp_emulator\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t, y, c):\n",
    "    \"Compute RHS of system of differential equations, returning vector derivative\"\n",
    "\n",
    "    # check inputs and extract\n",
    "\n",
    "    assert len(y) == 4\n",
    "    assert c >= 0.\n",
    "\n",
    "    vx = y[0]\n",
    "    vy = y[1]\n",
    "\n",
    "    # calculate derivatives\n",
    "\n",
    "    dydt = np.zeros(4)\n",
    "\n",
    "    dydt[0] = -c*vx*np.sqrt(vx**2 + vy**2)\n",
    "    dydt[1] = -9.8 - c*vy*np.sqrt(vx**2 + vy**2)\n",
    "    dydt[2] = vx\n",
    "    dydt[3] = vy\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event(t, y, c):\n",
    "    \"event to trigger end of integration\"\n",
    "\n",
    "    assert len(y) == 4\n",
    "    assert c >= 0.\n",
    "\n",
    "    return y[3]\n",
    "\n",
    "event.terminal = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sol(x):\n",
    "    \"simulator to solve ODE system for projectile motion with drag. returns distance projectile travels\"\n",
    "\n",
    "    # unpack values\n",
    "\n",
    "    assert len(x) == 2\n",
    "    assert x[1] > 0.\n",
    "\n",
    "    c = 10.**x[0]\n",
    "    v0 = x[1]\n",
    "\n",
    "    # set initial conditions\n",
    "\n",
    "    y0 = np.zeros(4)\n",
    "\n",
    "    y0[0] = v0/np.sqrt(2.)\n",
    "    y0[1] = v0/np.sqrt(2.)\n",
    "    y0[3] = 2.\n",
    "\n",
    "    # run simulation\n",
    "\n",
    "    results = solve_ivp(f, (0., 1.e8), y0, events=event, args = (c,))\n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulator(x):\n",
    "    return sol(x).y_events[0][0][2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = sol([np.log10(0.01),100.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/X0lEQVR4nO3dd3xV9fnA8c9DAgkrhBEgJEDYyIYgbmW4B6BQ1Fr36K+KdbZqnW1t66y17oWrDhQHTlARHLUCCXtvSJiBMMJIQpLn98f5pl7TjBvIzbk393m/XueVe/ZzT+59zrnf8z3fr6gqxhhjokc9vwMwxhhTuyzxG2NMlLHEb4wxUcYSvzHGRBlL/MYYE2Us8RtjTJSxxB/BRGSviHQO5XZF5BURuf8Qt6Mi0rWCeReJyBeHE6cx5tBY4q9lIrJORA645LrVJdYmh7ItVW2iqmsOM54ZInJVTW+3Kqr6hqqeGhBHhSeJwyEi94nIv6pYZryIZIhIgYi8Us78RiLytIhsF5HdIvJtwLxhIjLdTV9XxX4aiMgk9xlQERlaZv7n7nNROhSKyMJytnOSW//+gGlxIvKYiGwSkZ0u3voB82eISH7AtpcHzDtLRL4XkV0iskVEXhSRpgHzW4jIRBHZ4Y7BGyKSUCamG0RkrYjsE5GlItI94PgsdNveISIfiEhKmbgniMget++by2z3KhFZ5WKeIiLtyqz7rPse5YrIx2W2vbfMUCwiT7h5ae4YBs6/u7L/X11iid8f56hqE2AQMBi4q+wCIhJb61FFr03A/cCECuY/D7QAjnB/bwqYt8+t97sg9/U98CtgS9kZqnqGO+k2cZ+PH4B3A5dxyfxxYGaZ1W/H+yz1AbrjfbbKfq7GB2y/R8D0Znjvv517jynAwwHz7weaA52ALkAb4L6AmK4CrgTOApoAZwPb3ewlwGmqmui2vxJ4JmDb9wHdgI7AMOD3InK62+5Q4K/AKLzjvhZ4K2DdG4BjgH5u2zuBJ0pnljmWbYEDlDmeQGLAcn8mWqiqDbU4AOuAkwPGHwY+ca8VuA7vy7HWTbsaWAXkAh8B7QLWVaCrex0HPAJsALYCzwINA5YdBcwD9gCrgdOBvwDFQD6wF3iynO2+AtwfsJ2z3XZ24SWmfpW8VwX+z72fXcBTgLh5lwHfu9ffumX3uTjOB1oBn7j1coHvgHoV7OdxIMu9t0zgBDf9dKAQOOi2O7+K/839wCtlpvV0202oYt2TgXXV+BxkA0MrmZ/m/jdpZabfDjxUzv8lA/hFwPgvgayA8RnAVUHGdh6wMGD8c+DagPHrgKnudT137EcEsd044G/AkoBpm4BTA8b/DLztXj8CPBUwr537nHRx488ADwXMPwtYXsG+LwXWBHz+0ty2Ymv6Ox4Jg13x+0hE2gNnAnMDJo8GjgJ6ichwvC/KOCAZWA+8XcHmHsC70hsAdMW7arvH7WcI8BreVWkicCJekroTL6GWXgmOryLegXhXt78GWgLPAR+JSFwlq50NHIl3VTYOOK3sAqp6onvZ38UxEbgFLzkm4V1h/gHvi1qe2e59twDeBN4VkXhVnYJ3xTjRbbd/Ze+vAkPwjvsfXTHHQhEZcwjbqa5LgO9UdV3pBBHpCFwB/KmCdaTM61QRaRYw7W/uPfy7bDFTGScCiwPGnwLOFpHmItIcGIN3MgBIdUMfEclyxT1/FJH/5hYR6SAiu/CuuG/FO3HhtpUMzA/Y13ygdyXvCbxfNQAvAceJSDsRaQRcFBBXWZcCr6nL+gHWi0i2iLwsIq0qWLfOscTvjw/dF+F74Bu85FTqb6qaq6oH8D7IE1R1jqoWAHcAx4hIWuDGRESAa4Cb3Lp5bpsXuEWudNv5UlVLVHWjqi47hLivAZ5T1ZmqWqyqrwIFwNGVrPOAqu5S1Q3AdLwEHYyDeEmho6oeVNXvyvnSAqCq/1LVHapapKqP4l1Z9ihv2UOQipdoduNdcY4HXhWRI2po+xW5BO+qPtA/gbtVdW85y08BbhCRJBFpC/zWTW/k/t4GdMa7IHge+FhEupTdiIicgpck7wmYPAdoAOxwQzHwtJuX6v6eCvTFK665EO8zB4CqblCvqKcVXvFT6Wev9N7W7oB97QZK7y9MAcaJSD8Raehi0oD3tBLv18ZGvF9lR1DOSdGdME8CXg2YvB3vgqQjkO72+UbZdesqS/z+GK2qiaraUVWvdUm+VFbA63Z4V5sAuC/8Drwvb6AkvC9DpruJtgvvS5Pk5rfHK945XB2BW0r34fbT3sVZkcCy7P389GWvysN4RVxfiMgaEbm9ogVF5FZ3Q3G3i6kZXpKpCQfwTkL3q2qhqn6DdwI7tfLVDp2IHI9XJj0pYNo5QFP3a6g8f8H75TgPrwjuQxf3VgB3ss5T1QJ3wv433q/NwP0ejfeLaayqrgiY9Q6wAi85JuB9lkpvmJd+dh9yJ/h1eL8Ef7ZtF0MuXvKd7O5hlZ7AAm8UJwB5bvmvgHuB9/CKSNe5edlu2afwTvItgcbA+5R/xX8xXrHi2oBY9qpqhrtY2Ip3Qj818KZ2XWaJP/wEXtVuwku2AIhIY7wP+cYy62zH+wL2dieURFVtpt5NLfBOJv9zdVfO/qqSBfwlYB+JqtpIVd+qcs1qcknqFlXtDIwEbhaREWWXE5ETgN/jFSM1d1eWu/mpWOBwm59dUF54h7nNqlwKvF/myn4EMNjVfNmCdx/kRhGZDKCqB1R1vKqmuGO2A8hU1ZIK9qEEFKO4YryPgCtUdVqZZQfg/dLb52J6lp8S+3K8+yiBx6Sy4xMLtMa7Z7IT2AwEFsH1J6CYSVWfUtVuqtoG7wQQCywKiOsV9yu3AO/G7pByimwu4edX++UpjTkqcmJUvMkI9hZwuYgMcOXofwVmBpb7Argv9wvAYyLSGkBEUkSktDz9JbedESJSz83r6eZtxSsCCMYLwP+JyFHiaSxeVcCauEr6WRwicraIdHXFWLvxihfKS2JNgSIgB4gVkXv4+RXkViAtsMy5LBGJFZF4IAaIEZH4gFpV3+LdML/DLXccXnHGVLduPbdufW9U4kWkQSX7inPLAzRwywcm4IZ4J7FXyqx6Nz/dwxmAl6RfAC5366W4sm5xV+53410tIyKJInJa6fsSkYvwyvGnuPl93OvrVfXjcsKeDVwlIg1dfNfgToiquh+YiFcbp6mIpLr5n7htnyciPdxxSgL+Dsx1V//g3Xu6y90/6IlXmeEVt268iPRx76kDXhHV4+6EURrXJSLSTLzaTtcCm1S1tEYRInIs3i/ksrWjjgqIqyVeMdoMVQ0sdqq7avNOsg3/W6unzLz/1qYJmPZ/eD+tc/G+TKnlLQ/E450Y1uCVdy4Ffhuw7Ll4X9Y8vCKU09z0Y/B+xu8E/lnOdl/h57VHTsf7wu3Cu1p7F68Iosr3E7gtAmr1BLzPzW674/CqTK7Dq+mTjVe2Xd4+YvBuOO9x6/8+8Bjj/UL63r2/ORVs4z4Xa+BwX8D83sB/XCxLgHMD5g0tZ90ZAfMXAxeV+f+XXT4tYP6FeMV7UsXnqOz/5US37f14V+GB+0xy/7M8d3x/BE4JmP8y3kl1b8CwOGB+J+BjvF8RuXgniW4B8xPwKh3k4f0qvKc0fuB6vGqY+/CK/d7Gu29Tum5cwP9vK3BzwLxEvM9s6bp/A2IC5rfEK5ff5t7X98CQMsfpOeD1co7fhQFxbcY7AbX1Oz/U1lD6zzERxl3BFuN9iTb4HY8xJnJYUU/k6oNX//5/HgQyxpjKWOKPQOLVI58O3KaqhX7HY4yJLFbUY4wxUcau+I0xJspERENgrVq10rS0NL/DMMaYiJKZmbldVZPKTo+IxJ+WlkZGRobfYRhjTEQRkfXlTbeiHmOMiTKW+I0xJspY4jfGmChjid8YY6KMJX5jjIkylviNMSbKWOI3xpgoY4nfmAiTuX4n72Vmc6Cw2O9QTISKiAe4jDGQk1fA3z5fyvtzvA7Y/vrZUq46oTO/OroDTePr+xydiSSW+I0Jc0XFJfzrx/U8+sUKCopKGD+sK8d0aclz367hwSnLeGbGKi4/rhOXH5dGYqMKO/8y5r8s8RsTxmavy+XuDxexbEseJ3ZP4o8je9OpVWMAjuvaigXZu3jy61U8Pm0lL363hl8d05Grju9MUtM4nyM34SxkzTKLSA+8vjhLdcbrku01Nz0Nr6u4cfpTH5rlGjx4sFpbPSaa5OQV8MDny3hvTjYpiQ25++xenNa7DQHd8/7Msi17eGr6aj5dsIn6MfW4cEgHfn1SZ5KbNazlyE04EZFMVR38P9Nroz1+EYkBNgJHAdcBuar6gIjcDjRX1dsqW98Sv4kW/y3W+XIF+QeLuebEzlw3rCuNGgT343xNzl6embGaD+ZuRATGpqfym5O60qFloxBHbsKR34n/VOBeVT1ORJYDQ1V1s4gk43VM3aOy9S3xm2iQsS6XuycvZunmPZzQrRV/HNmbzklNDmlb2Tv389w3a5iYkUVxiTKqfzuuHdaFrq2b1nDUJpz5nfgnAHNU9UkR2aWqiW66ADtLx8uscw1wDUCHDh3S168vt3VRYyLe9r1esc6kzGzaNYvn7rN7cXqfthUW61TH1j35vPDtGt6YuYH8omLO6NOW64Z1pXe7ZjUQuQl3viV+EWkAbAJ6q+rWwMTv5u9U1eaVbcOu+E1dVFRcwhszN/DIF8vJP1jM1Sd0Zvzw4It1qmPH3gIm/Hstr/2wnryCIkb0bM11w7syqEOlXz0T4SpK/LVRq+cMvKv9rW58q4gkBxT1bKuFGIwJK5nrc7n7w8UsccU6943sTZdDLNYJRssmcfzutJ5cc2IXXvthHRP+vZbznv6B47q2ZPywbhzduUWN/MIwkaE2rvjfBqaq6stu/GFgR8DN3Raq+vvKtmFX/Kau2L63gAc/X8a7mdkku2KdM2qoWKc69hUU8ebMDTz/3Rpy8gpI79ic8cO7MrR7kp0A6hBfinpEpDGwAeisqrvdtJbAO0AHYD1edc7cyrZjid9EuuIS5c2Z63l46nIOHCzmqhM6M35YVxrH+fsoTf7BYt7NyOLZb9awcdcB+qQkMH5YN07t1YZ69ewEEOl8vbl7uCzxm0iWuX4n90xexOJNezi+ayv+OCq0xTqHorCohA/nbuTpGatYt2M/3ds04bphXTmrbzKxMdakV6SyxG9MLduxt4AHpyzjnQx/i3Wqo6i4hE8Xbuap6atYsXUvaS0b8ZuhXTh3YCoNYu0EEGks8RtTSwKLdfYXesU61w/3v1inOkpKlC+WbOWp6atYuHE37ZrF839DuzBucHvi68f4HZ4JkiV+Y2rBnA1esc6ijXs4rmtL/jiyd0Q/NKWqfLMihye/XkXG+p0kNY3j6hM6cdFRHSPqRBatLPEbE0I79hbw0JTlTMzIom1CPHedfQRn9U0O62Kd6lBVZq7N5cmvV/H9qu0kNqrPlcd14pJj02jW0JqEDleW+I0JgeIS5c1ZG3hk6nL2FRRx5Qmd+O3wbnX6anjOhp089fUqpi3bRtO4WC45tiNXHNeJlk2sRdBwY4nfmBo2d8NO7nbFOsd2acmfRkV2sU51Ld60m6enr+azRZuJj43hoqM6cPWJnWmTEO93aMaxxG9MDcndV8hDU5bx9uws2iTEcddZvTi7X90p1qmuVdvyeHrGaibP20SMCOOOTOXXJ3ahfQtrEdRvlviNOUzFJcrbszfw0BRXrHN8J64f0Y0mdbhYpzo27NjPM9+sZlJmFqowemAK1w7tcsgtjJrDZ4nfmMMwL2sX90xexILs3RzT2SvW6dYmeop1qmPz7gM8980a3pq1gYPFJZzVrx3XDetCz7YJfocWdSzxG3MIcvcV8vBUr1inddM47jyrF+dEcbFOdeTkFfDS92t5/T/r2FdYzCm92jB+WFf6t0/0O7SoYYnfmGooLlEmzs7ioanL2JtfxBXHd+K3VqxzSHbtL+SVH9bx8r/XsfvAQU7snsT4YV0Z0qmF36HVeZb4jQnSfFesMz97N0d3bsGfRvWhuxXrHLa9BUX868f1vPjdGrbvLWRIpxZcP7wrx3dtZb+gQsQSvzFV2LmvkIemLuft2RtIahLHnWcdwcj+7Swp1bADhcW8PXsDz32zhi178unfPpHrh3VlxBGt7VjXMEv8xlSgpER52xXr5OUXccVxafx2RDeaxtsTqaFUUFTMe5kbeeabVWTlHqBn26aMH96VM/okE2NNQtcIS/zGlCOwWOeoTl6xTo+2VqxTm4qKS/ho/iaemr6K1Tn76JzUmGuHdmXUgHbUtyahD4slfmMC7NxXyMNfLOetWRto1SSOu6xYx3fFJcrUxVt44utVLN28h9TmDfnN0C6MTU8lLtZaBD0UlviNwSvWeScjiwenLGNPfhGXH5vGDSdbsU44UVW+XraNJ75exbysXbRJiOOaE7vwyyEdaNjATgDVYYnfRL2F2bu5a/Ii5mftYkinFvzZinXCmqryw+odPPH1Sn5ck0vLxg248oROXHx0RztRB8mvPncTgReBPoACVwDLgYlAGrAOr8/dnZVtxxK/ORy79hfy8NTlvOmKde488whGDbBinUiSsS6XJ6evYsbyHBLiY7nsuE5cfUInOwFUwa/E/yrwnaq+KCINgEbAH4BcVX1ARG4HmqvqbZVtxxK/ORQlJcq7mVk88LlXrHPZsWncaMU6EW1h9m6enL6SqYu3MmpAOx6/YKDfIYW1Wk/8ItIMmAd01oCdiMhyYKiqbhaRZGCGqvaobFuW+E11Ldq4m7s+XMS8rF0MSWvBn0b3trZi6pB7Jy/irdlZzP7DyTRrZCfyilSU+ENZV6oTkAO8LCJzReRFEWkMtFHVzW6ZLUCbEMZgotCL363hnCe/J3vnAR47vz8Tf320Jf065heD21NYVMLHCzb5HUpECmXijwUGAc+o6kBgH3B74ALul0C5PzlE5BoRyRCRjJycnBCGaeqSdzKyuP/TpZzeuy1f33oS5w5MtbL8Oqh3uwR6tGnKpMxsv0OJSKFM/NlAtqrOdOOT8E4EW10RD+7vtvJWVtXnVXWwqg5OSkoKYZimrvhyyVZuf28BJ3RrxeMXDCTByvLrLBFhbHoq87J2sWrbXr/DiTghS/yqugXIEpHS8vsRwBLgI+BSN+1SYHKoYjDRY+aaHVz35hz6piby7K/SaRBrT3zWdaMGtiOmnvDeHLvqr65QfzuuB94QkQXAAOCvwAPAKSKyEjjZjRtzyJZs2sNVr2bQvnlDXr7syDrd0bn5Seum8ZzUPYn352RTXBL+zyOFk5B+Q1R1HvA/d5Txrv6NOWzrd+zjkgmzaBIfy+tXHkWLxg38DsnUorHpqXy9bBvfr9rOSd2tSDhY9nvYRKxteflc/NIsiktKeP3KIbRLbOh3SKaWjTiiNc0a1rebvNVkid9EpN0HDnLphNls31vAy5cPoWtra3ohGsXFxjBqQDu+WLyF3QcO+h1OxLDEbyJO/sFirn4tg1Xb8nj2V+kMsD5co9qYQakUFJXw6YLNVS9sAEv8JsIUFZcw/s25zF6Xy9/HDeBEK9eNev1Sm9GtdRMmZWb5HUrEsMRvIoaqcsf7C/lq6Vb+OLI35/Rv53dIJgyU1umfs2EXq3OsTn8wLPGbiPHAlGW8m5nNDSO6cckxaX6HY8LIuQNTqCfwvtXpD4olfhMRnv92Nc99s4aLj+7IjSd38zscE2ZaJ8RzYvck3p+z0er0B8ESvwl772Zk8dfPlnF2v2TuG9nb2t4x5Rqbnsrm3fn8sHq736GEPUv8Jqx9tWQrt7+/kBO6teLv4wYQU8+SvinfyUe0ISE+1ur0B8ESvwlbs9bmct2bc+jTLsHa3zFViq8fw8gB7Zi6eAt78q1Of2Xsm2TC0pJNe7jy1dmkNG/Iy5cPsfZ3TFDGDEol/2AJn1md/kpVmfhF5BgReUpEFohIjohsEJHPROQ618uWMTVqw479XPryLJrEWfs7pnoGtE+kS1JjK+6pQqWJX0Q+B64CpgKnA8lAL+AuIB6YLCIjQx2kiR7b8vK5eMJMDhaX8NoVQ0ix9ndMNXh1+tuTsX4na7fv8zucsFXVFf/Fqnqlqn6kqptUtUhV96rqHFV9VFWHAj/UQpwmCuzJP8hlE2azbU8BEy47km5trP0dU31Wp79qlSZ+Vf1ZvSgRSRCRFqVDecsYcyjyDxZz9asZrNiax7MXpzOoQ3O/QzIRqm2zeI7vlsR7mdmUWJ3+cgV1c1dEfi0iW4AFQKYbMkIZmIkeRcUl/Patucxal8uj4/pbu+rmsI1NT2XT7nz+s2aH36GEpWCrStwK9LGre1PTVJU7P1jEF0u2ct85vRg1IMXvkEwdcGqvNjR1dfqP69rK73DCTrDVOVcD+0MZiIlOD05ZzsSMLH47vCuXHdfJ73BMHRFfP4Zz+rfj80WbybM6/f8j2MR/B/CDiDwnIv8sHUIZmKn7Xvh2Dc9+s5qLjurATad09zscU8eU1un/fOEWv0MJO8EW9TwHfA0sBEqC3biIrAPygGKgSFUHu5vCE4E0YB0wTlV3Bh+yqQvey8zmL58t5cy+bfnTqD7W/o6pcYM6JNK5lVenf9yR7f0OJ6wEm/jrq+rNh7iPYWXuDdwOTFPVB0Tkdjd+2yFu20SgaUu38vv3FnB811Y8dr61v2NCQ0QYk57Kw1OXs37HPjq2bOx3SGEj2KKez0XkGhFJLlud8xCMAl51r18FRh/idkwEmr0ul2vfmEPvdgk8e3E6cbExfodk6rDzBqUgAu/N2eh3KGEl2MR/Ia6cn+pV51TgCxHJFJFr3LQ2qlrakMYWoE15K7oTTYaIZOTk5AQZpglnSzfv4YpXXPs7lx1JE2t/x4RYcrOGHN+1ldXpLyOoxK+qncoZOgex6vGqOgg4A7hORE4ss13FOzmUt8/nVXWwqg5OSrJ63ZEuK3c/l06YReMGsbx2xRBaNonzOyQTJcamp7Jx1wF+XGt1+ktV1VbP8VXMTxCRPhXNV9WN7u824ANgCLBVRJLd+snAtuoGbSJLTl4BF780k4KiEl67cgipzRv5HZKJIqf2akvTOGunP1BVV/xjROQHEblHRM4SkSEicqKIXCEirwOfAOW2oiUijUWkaelr4FRgEfARcKlb7FJgco28ExOW9uQf5LKXZ7HVtb/T3drfMbWsYYMYzu6fzJRFW9hXUOR3OGGh0kJWVb3J3cQdA/wCr3XOA8BS4DlV/b6S1dsAH7hqerHAm6o6RURmA++IyJXAemDc4b8NE47yDxZzzWsZLN+Sx4uXDia9o7W/Y/wxZlAqb83K4rOFm/nFYKvaWeXdNVXNBV5wQ9BUdQ3Qv5zpO4AR1dmWiTxFxSXc8PZcflyTy+MXDGBoj9Z+h2SiWHrH5qS1bMSkzGxL/FgPXCYEStvfmbp4K/da+zsmDHjt9Kcyc20uWbnW+owlflPjHp7qtb9z/fCuXG7t75gwce6gVFen327yWuI3NerF79bw9IzV/PKoDtxs7e+YMJKS2JBju7TkvTlWpz/Y9vgbicjdIvKCG+8mImeHNjQTad6fk839ny7ljD5t+bO1v2PC0Nj0VLJyDzBrXa7fofgq2Cv+l4EC4Bg3vhG4PyQRmYj09bKt/G7SAo7t0pJ/XGDt75jwdFrvtjSJi+W9KK/TH2zi76KqDwEHAVR1P2DfbANAhmt/p1dyAs9fMtja3zFhq1GDWM7qm8ynCzdHdZ3+YBN/oYg0xDWvICJd8H4BmCi3bIvX/k67Zg155XJrf8eEvzHpqewvLGbKouhtpz/YxH8vMAVoLyJvANOA34csKhMRsnL3c8lLs2jYIIbXrrT2d0xkODKtOR1aNIrqJhyCbaTtS+A84DLgLWCwqs4IXVgm3G3f+1P7O69feZS1v2MiRmmd/v+s2UH2zuis0x9srZ5z8XrQ+lRVPwGKRGR0SCMzYSvPtb+zZU8+Ey4bbO3vmIhz3iDvocL3o7Sd/qCLelR1d+mIqu7CK/4xUcZrfyeTZZvzeOZX6aR3PNT+eIzxT2rzRhzTuSWTMrPxWoePLsEm/vKWs7t4Uaa4RLnx7Xn8Z80OHvlFf4ZZ+zsmgo1NT2VD7n5mr4u+Lr+DTfwZIvJ3Eenihr/j9cJlooSqcteHC5myeAv3nN2L0QOt/R0T2c7o25bGDWKisk5/sIn/eqAQmOiGAuC6UAVlws8jXyznrVlZXDesC1ccb+3vmMjXqEEsZ7o6/fsLo6tOf7C1evap6u2lXSGq6h2qui/UwZnw8NL3a3lq+mouHNKeW0/t4Xc4xtSYMemp7C0oYuri6KrTH1Q5vYh0B24F0gLXUdXhoQnLhIsP5mbz50+WcHrvttw/uq+1v2PqlCFpLWjfoiGTMrM5d2Cq3+HUmmBv0L4LPAu8CBSHLhwTTqYv28bv3l3AMZ2t/R1TN9WrJ4wZlMrj01aycdcBUhLL7Um2zgm2jL9IVZ9R1Vmqmlk6hDQy46vM9bn85o1MeiY35flL0omvb+3vmLppzKBUVOGDKGqnP9jE/7GIXCsiySLSonQIZkURiRGRuSLyiRvvJCIzRWSViEwUkQaHHL0JieVb8rj85dkkN2vIK5cPoWl8fb9DMiZk2rdoxFGdWkRVnf5gE/+lwO+AH/CqcWYCGUGuewNe5+ylHgQeU9WuwE7gyiC3Y2pBVu5+Lpkwk/j6Mbx2xRBaWfs7JgqMTU9l3Y79ZK6Pjjr9wdbq6VTO0Lmq9UQkFTgL794A4t0ZHA5Mcou8Cow+pMhNjcs/WMzVr2VwoLCY164cQvsW1v6OiQ5n9k2mUYOYqOmWMeiuF0Wkj4iME5FLSocgVvsHXiueJW68JbBLVUsrzWYD5T4JJCLXiEiGiGTk5OQEG6Y5DPdOXsyyLXn888KB9Gyb4Hc4xtSaxnGxnNEnmU/mb+ZAYd2vvxJsI233Ak+4YRjwEDCyinXOBrYd6k1gVX2+9LmBpKSkQ9mEqYZJmdlMzMhi/LCuDLWmGEwUGpOeQl5BEV8sqft1+oO94h8LjAC2qOrlQH+gWRXrHAeMFJF1wNt4RTyPA4kiUlqNNBWvG0fjo+Vb8rjrw4Uc3bkFN57cze9wjPHF0Z1akpLYMCra6Q828R9Q1RK85pgTgG1A+8pWcE/3pqpqGnAB8LWqXgRMxzuRgHfTePIhRW5qxL6CIq59I5MmcfX55wUDiY0JuvTPmDqlXj1hTHoq36/azubdB/wOJ6Sq00hbIvACXo2eOcB/DnGftwE3i8gqvDL/lw5xO+YwqSp/+GAha7fv458XDqB1QrzfIRnjqzGDUlCt++30B/Xkrqpe614+KyJTgARVXRDsTlxvXTPc6zXAkOqFaULhrVlZTJ63iVtO6c6xXVr5HY4xvuvYsjFD0lrwXmY21w7tUmebKAn25u600tequk5VFwROM5Fn0cbd3PfxYk7snsR1w7r6HY4xYWNseiprtu9jzoZdfocSMpUmfhGJd0/othKR5gFP7aZRQTVME/725B/kujfn0KJRAx4b15961gaPMf91Zr9kGtav23X6q7ri/zVemX5PfnpiNxPvhuyToQ3NhIKqctukBWTvPMCTvxxIS3sy15ifaRIXy+l92vLx/E3kH6ybdforTfyq+riqdgJuVdXOAU/t9ldVS/wR6JUf1vH5oi3cdnoPBqdZf7nGlGdseip5+UV8sWSr36GERLC1eraISFMAEblLRN4XkUEhjMuEwNwNO/nrZ0s5+YjWXH1ClS1uGBO1juncknbN4utsnf5gE//dqponIscDJ+NVwXwmdGGZmrZrfyHj35xLm4R4Hv3FgDpbW8GYmvDfOv0rc9iyO9/vcGpcsIm/tKDrLOB5Vf0UsOaUI0RJiXLLO/PZlpfPU78cRLNG1syyMVU5b1AqJQofzK17dfqDTfwbReQ54HzgMxGJq8a6xmfPf7eGacu2cddZvejfPtHvcIyJCJ1aNWZwx+ZMysyqc+30B5u8xwFTgdNUdRfQAq99fhPmZq3N5eGpyzmrbzKXHNPR73CMiShj01NZnbOPeVm7/A6lRlVVj7+0bd54vCdvd7h6/QUE3xGL8cn2vQVc/9Yc2jdvyANjrKN0Y6rrzH7JxNevV+fq9Fd1xf+m+1va41ZgXX5L/GGsuES58e157Nx/kKcvSrfuE405BAnx9Tmtd1s+mle36vRXVY//bPe3U5l6/EH1wGX88+TXq/h+1Xb+NLI3vdpZpyrGHKqx6ansyS/iq6V1p05/pY20VVVXX1Xn1Gw4pib8e9V2/jFtBecNTOH8IyttPdsYU4Vju7Qi2dXpP7tfO7/DqRFVtc75qPsbDwwG5gMC9MMr6jkmdKGZQ7F1Tz43vD2XrklNuP/cPlaub8xhiqknnDcohWdmrGbbnvw60Xx5VUU9w1R1GLAZGOS6QkwHBmI9Z4WdouISrn9rLvsKinn6okE0ahBUq9vGmCrUtTr9wVbn7KGqC0tHVHURcERoQjKH6u9frmDW2lz+el4furVp6nc4xtQZXZKaMKhDIpMys+tEnf5gE/8CEXlRRIa64QUg6I5YTOhNX7aNp2es5sIh7Tl3YKrf4RhT54xNb8/KbXtZkL3b71AOW7CJ/3JgMXCDG5a4aSYMbNx1gJvemccRyQnce05vv8Mxpk46q18ycbF1o05/UIlfVfNV9TFVPdcNj6lq3Wu5KAIVFpUw/s05FBUrT180iPj6MX6HZEyd1KxhfU7t3ZbJ8zZRUBTZdfpD1t6O671rlojMF5HFIvJHN72TiMwUkVUiMlFErLG3w/DA58uYu2EXD47pR6dWjf0Ox5g6bWx6KrsPHGTa0m1+h3JYQtnQWgEwXFX7AwOA00XkaOBB4DFV7QrsBK4MYQx12pRFm5nw77VcdmwaZ/VL9jscY+q847u2ok1CHO9FeDv9IUv86tnrRuu7QYHhwCQ3/VVgdKhiqMvW79jH795dQP/UZtxxZk+/wzEmKnh1+lOZsSKHbXmRW9pd1ZO7H+Ml63Kp6sgq1o/Ba9enK/AUsBrYpapFbpFsKui0XUSuAa4B6NChQ2W7iTr5B4u59o05iMCTvxxEXKyV6xtTW8YMSuWZGauZPHcTV58YmS3XVPWEzyOHs3FVLQYGiEgi8AFep+3Brvs88DzA4MGDI7/ibA368ydLWLxpDy9eMpj2LRr5HY4xUaVr6yYMaO/V6b/qhE4R+XR8pYlfVb+piZ2o6i4RmY7XxEOiiMS6q/5U7Angapk8byNvzNzAr0/qzMm92vgdjjFRaWx6Knd9uIjFm/bQJ6WZ3+FUW1Bl/CLSTUQmicgSEVlTOlSxTpK70kdEGgKnAEuB6cBYt9ilwORDjj7KrNq2lzveX8iRac259dQefodjTNQ6p187GsTWi9jO2IO9ufsyXufqRcAw4DXgX1WskwxMF5EFwGzgS1X9BLgNuFlEVgEt8TpuN1U4UFjMtW9kEl8/hicuHET9GOv50hi/NGtUn1N6teHDeRsjsk5/sK14NVTVaSIiqroeuE9EMoF7KlpBVRfgNeZWdvoaYMghRRvF7p68iJXb9vLq5UNo2yzyWwc0JtKNTU/l0wWbmb5sG6f3iazq1MFeNhaISD1gpYiMF5FzgSYhjMsEeCcji0mZ2Vw/vBsndk/yOxxjDHBC11a0bhrHpMzIu00ZbOK/AWgE/BZIB36FVz5vQmzp5j3c/eEiju3SkhtGdPM7HGOMExtTj3MHpTB9+TZy8gr8Dqdagm2rZ7aq7lXVbFW9XFXHqOqPoQ4u2u0tKOK6N+aQ0LA+j18wkJh6kVdtzJi6bOygVIpLlMnzIuuqP9haPV+W1tBx481FZGrIojKoKne8v5B1O/bxxIUDSWoa53dIxpgyurVpSv/UZhFXuyfYop5WqrqrdERVdwKtQxKRAeBfMzfw8fxN3HJqD47u3NLvcIwxFRibnsqyLXks3hQ57fQHm/hLROS/7SaISEcqacrBHJ6F2bv588dLGNojid+c1MXvcIwxlTinfzsaxERWnf5gE/+dwPci8rqI/Av4FrgjdGFFr90HDnLtm5m0atKAx8YNoJ6V6xsT1hIbNeDkXq2ZPG8ThUUlfocTlGBv7k4BBgETgbeBdFW1Mv4apqr87t35bN6VzxO/HETzxtZVgTGRYGx6Krn7Cpm+PDLa6a808YtIT/d3ENAB2OSGDm6aqUEvfb+WL5Zs5fYzepLesbnf4RhjgnRityRaNYmcdvqrenL3ZrymkR8tZ15p2/qmBmSu38kDny/j1F5tuPL4Tn6HY4yphtiYepw3KIUJ369lx94CWjYJ71p4lV7xq+o17uUZqjoscADODH140SF3XyHj35xDcmI8D/+if0Q282pMtBszKJWiEmXyvE1+h1KlYG/u/hDkNFNNJSXKze/MY8feQp7+ZTrNGtb3OyRjzCHo0bYpfVMio05/VWX8bUUkHWgoIgNFZJAbhuI14WAO0zPfrGbG8hzuPqcXfVMjr11vY8xPxqansmTzHpZs2uN3KJWqqoz/NOAyvA5THgVKyyD2AH8IXVjR4cc1O3j0i+Wc078dvzrKupc0JtKN7N+O+z9dwntzsunVrpff4VSoqh64XhWR14ELVfWNWoopKuTkFXD9W3NJa9mYv53X18r1jakDmjduwIiebfhw7kZuP6Nn2PabUWVUqloC3FQLsUSN4hLlhrfnsufAQZ7+1SCaxAXbLYIxJtyNTU9lx75CZizP8TuUCgV7OvpKRG4VkfYi0qJ0CGlkddjj01byw+od/Hl0H3q2TfA7HGNMDTqpRxKtmjQI6zr9wV5qnu/+XhcwTYHONRtO3fftihye+HolY9NTGTe4vd/hGGNqWP2YeowekMKr/1lH7r5CWoThE/jBNtnQqZyh0qTvfh1Mdx20LxaRG9z0Fq6Z55Xub9Q8orpldz43TpxH99ZN+fOoPn6HY4wJkTHpqRwsVj4K03b6g22Pv5GI3CUiz7vxbiJydhWrFQG3qGov4GjgOhHpBdwOTFPVbsA0N17nHSwu4fq35pB/sJinLhpEwwYxfodkjAmRI5IT6N0ugUlzwrO4J9gy/peBQuBYN74RuL+yFVR1s6rOca/zgKVACjAKeNUt9iowunohR6ZHvljO7HU7+dt5fena2rorNqauG5ueyqKNe1i2Jfzq9Aeb+Luo6kPAQQBV3c9PdfqrJCJpwEBgJtBGVTe7WVuANhWsc42IZIhIRk5O+N4dD8ZXS7by3DdruOioDowakOJ3OMaYWjBqQAr1YyQsb/IGm/gLRaQhrvMVEekCBNW7sIg0Ad4DblTVn536VFWpoEMXVX1eVQer6uCkpKQgwww/Wbn7ueXd+fRul8DdZ4fvAx3GmJrVonEDhvVozQdzN3GwOLza6Q828d8HTAHai8gbeGXzv69qJRGpj5f031DV993krSKS7OYnA5HRgPUhKCwqYfybcygpUZ6+aBDx9a1c35hoMjY9le17C/h2RXiVWlTVVs9TInKcqn4BnIfXfMNbwGBVnVHFugK8BCxV1b8HzPoIuNS9vhSYfGihh78HpyxjfvZuHv5FPzq2bOx3OMaYWjasZ2taNm7Ae2F2k7eqK/4VwCMisg64Ddikqp+o6vYgtn0ccDEwXETmueFM4AHgFBFZCZzsxuuc71bm8NL3a7nkmI6c3ifZ73CMMT6oH1OPUQNS+GrJNnbuK/Q7nP+qqj3+x1X1GOAkYAcwQUSWici9ItK9inW/V1VR1X6qOsANn6nqDlUdoardVPVkVc2twfcTFnbuK+TWd+fTtXUT/nDmEX6HY4zx0Zj0FAqLS/h4Qfi00x/sA1zrVfVBVR0IXIhXBXNpKAOLVKrKHe8vJHdfIf84f4CV6xsT5Xq3a8YRyQlh1U5/sA9wxYrIOe7G7ufAcrwyf1PGu5nZTFm8hVtO7UGfFGtf3xjj3eRdkL2bFVvz/A4FqPrm7ikiMgHIBq4GPsWr03+BqtbZm7KHav2Offzxo8Uc3bkFV59gzRgZYzyjBrQjtl741Omv6or/DrwuFo9Q1ZGq+qaq7quFuCJOUXEJN02cR716wqPjBhBTz9rXN8Z4WjWJY2iP1rw/dyNFYVCnv6qbu8NV9UVV3VlbAUWqp6avZs6GXfzl3L6kJDb0OxxjTJgZm55KTl4B360MplJkaIVn9zARZs6Gnfzz65WMHtCOkf3b+R2OMSYMDe/ZmuaN6odFw22W+A/TvoIibpo4j7YJ8fxptDW1bIwpX4NYr07/l4u3snv/QV9jscR/mP708RI25O7n7+P6kxBf3+9wjDFhbGx6KoXFJXzkc51+S/yHYcqiLUzMyOI3J3XhqM4t/Q7HGBPmerdLoGfbpr7X6bfEf4i27snnjvcX0CclgRtPrvQhZmOMAUBEGJueyvysXaza5l+dfkv8h6CkRLn13fkcOFjMP84fSINYO4zGmOCMGpBCTD1hUqZ/3TJaxjoEr/5nHd+t3M6dZ/Wy3rSMMdWS1DSOod2T+GBuNsUl5XZHEnKW+KtpxdY8/vb5Mob3bM2vjurgdzjGmAg0Nj2VrXsK+H6VP3X6LfFXQ0FRMb99ay5N42J5cEw/vC4HjDGmeoYf0ZrERvV9u8lrib8aHv1iBcu25PHQ2H4kNY3zOxxjTISKi41hZP92TF28hd0Har9OvyX+IP2wajsvfOd1mD7iiHL7hzfGmKCNTU+lsKiET3yo02+JPwi79x/k5nfm06lVY+46yzpMN8Ycvr4pzejepokvLXZa4q+CqvKHDxeyfW8Bj58/kIYNrGMVY8zhK63TP2fDLlbn7K3VfYcs8YvIBBHZJiKLAqa1EJEvRWSl+9s8VPuvKR/M3cinCzZz0ynd6ZtqHasYY2rO6AEp1BNq/ao/lFf8rwCnl5l2OzBNVbsB09x42MrK3c89kxczJK0F/3dSF7/DMcbUMa0T4jmpexLvz9lYq3X6Q5b4VfVboGxH6qOAV93rV/H67g1LxSXKTRPnIcCj4/pbxyrGmJAYm96eLXvy+WF17dXpr+0y/jaqutm93gKEbfWYZ2asImP9Tv40ujftWzTyOxxjTB014ojWNI2P5bOFm6teuIb4dnNXVRWo8LeNiFwjIhkikpGTk1OLkcH8rF3846uVnNO/HaMHpNTqvo0x0SW+fgwD2icyP2t3re2zthP/VhFJBnB/t1W0oKo+r6qDVXVwUlJSrQW4v9DrWKV10zjuH9XHns41xoRcv9RmrNiaR/7B4lrZX20n/o+AS93rS4HJtbz/Kv35k6Ws3bGPR8b1p1kj61jFGBN6fVMSKSpRlm7eUyv7C2V1zreA/wA9RCRbRK4EHgBOEZGVwMluPGx8uWQrb83awDUndObYLq38DscYEyX6uariCzfWTnFPbKg2rKoXVjBrRKj2eTi25eVz23sL6JWcwM2nWscqxpjak9wsnlZNGrAgu3YSvz25i/d07u8nLWBfQRGPXzCAuFh7OtcYU3tEhL4pzVhoib/2vP7jemYsz+EPZx5BtzZN/Q7HGBOF+qYmsnJbHvsLi0K+r6hP/Ku25fGXT5dyUvckLjmmo9/hGGOiVL+UZpQoLNkU+hu8UZ34C4tKuOHteTSOi+XhX1jHKsYY/5S2BVYb5fwhu7kbCf7+5QoWb9rD8xen07ppvN/hGGOiWJuEeNokxNVKzZ6oveL/z+odPPftai4c0p5Te7f1OxxjjKFvSiILsneFfD9Rmfh3HzjILe/MI61lY+4+2zpWMcaEh36pzVizfR95+aHtjjEqE//dHy5ia14Bj50/gEYNorq0yxgTRvqmNkMVFof4Bm/UJf7J8zby0fxN3DiiGwPaJ/odjjHG/FffFPcEb4hv8EZV4s/euZ+7PlhEesfm/GaodaxijAkvrZrEkZLYkAUhvsEbNYm/uES5+Z35KPCP8wcQGxM1b90YE0H6pCSwMMQ3eKMm+7343Rpmrc3l3nN6Wccqxpiw1S81kXU79rN7f+hu8EZF4l+6eQ+PfrGC03q3YWx6qt/hGGNMhUrL+RdtCl1xT51P/AVFxdw0cR4JDevz13P72tO5xpiwVpr4Q/kEb52vy/jYlytZtiWPly4dTMsmcX6HY4wxlWreuAHtWzRk4cZdIdtHnb7in70ul+e+Xc0FR7ZnxBFh26+7Mcb8TL+UxJA23VCnE/+Dny8jtXlD7rKnc40xEaRvajOycg+wc19hSLZfp4t6nr9kMNvy8mkSV6ffpjGmjumX8lNXjCd2T6rx7dfpK/4WjRvQs22C32EYY0y19E4JbR+8viR+ETldRJaLyCoRud2PGIwxJlw1a1ifTq0ah6ylzlpP/CISAzwFnAH0Ai4UESuEN8aYAKHsg9ePK/4hwCpVXaOqhcDbwCgf4jDGmLDVL7UZm3bnk5NXUOPb9iPxpwBZAePZbtrPiMg1IpIhIhk5OTm1FpwxxoSDI9NacFbfZA4UFtf4tsP25q6qPq+qg1V1cFJSzd/VNsaYcNa/fSJPXTSIDi1rvm0xPxL/RqB9wHiqm2aMMaYW+JH4ZwPdRKSTiDQALgA+8iEOY4yJSrX+ZJOqFonIeGAqEANMUNXFtR2HMcZEK18eaVXVz4DP/Ni3McZEu7C9uWuMMSY0LPEbY0yUscRvjDFRxhK/McZEGVFVv2OokojkAOsPcfVWwPYaDKemWFzBC8eYwOKqrnCMKxxjgpqLq6Oq/s8TsBGR+A+HiGSo6mC/4yjL4gpeOMYEFld1hWNc4RgThD4uK+oxxpgoY4nfGGOiTDQk/uf9DqACFlfwwjEmsLiqKxzjCseYIMRx1fkyfmOMMT8XDVf8xhhjAljiN8aYKFOnE384dOouIu1FZLqILBGRxSJyg5veQkS+FJGV7m9zn+KLEZG5IvKJG+8kIjPdMZvoms6u7ZgSRWSSiCwTkaUicozfx0tEbnL/v0Ui8paIxPt1rERkgohsE5FFAdPKPT7i+aeLcYGIDKrFmB52/8MFIvKBiCQGzLvDxbRcRE4LRUwVxRUw7xYRURFp5cZr5VhVFpeIXO+O2WIReShges0eL1WtkwNek8+rgc5AA2A+0MuHOJKBQe51U2AFXifzDwG3u+m3Aw/6dJxuBt4EPnHj7wAXuNfPAr/xIaZXgavc6wZAop/HC69r0LVAw4BjdJlfxwo4ERgELAqYVu7xAc4EPgcEOBqYWYsxnQrEutcPBsTUy30f44BO7nsaU1txuent8ZqGXw+0qs1jVcnxGgZ8BcS58dahOl4h/5D6NQDHAFMDxu8A7giDuCYDpwDLgWQ3LRlY7kMsqcA0YDjwifvAbw/4sv7sGNZSTM1ckpUy0307XvzUT3QLvKbMPwFO8/NYAWllkka5xwd4DriwvOVCHVOZeecCb7jXP/suugR8TG0dKzdtEtAfWBeQ+GvtWFXwP3wHOLmc5Wr8eNXlop6gOnWvTSKSBgwEZgJtVHWzm7UFaONDSP8Afg+UuPGWwC5VLXLjfhyzTkAO8LIrgnpRRBrj4/FS1Y3AI8AGYDOwG8jE/2MVqKLjEy7fgyvwrqbB55hEZBSwUVXnl5nl97HqDpzgig+/EZEjQxVXXU78YUVEmgDvATeq6p7Aeeqdxmu1Xq2InA1sU9XM2txvEGLxfgI/o6oDgX14RRf/VdvHy5WXj8I7KbUDGgOn19b+q8uPz1NlROROoAh4IwxiaQT8AbjH71jKEYv3q/Jo4HfAOyIiodhRXU78YdOpu4jUx0v6b6jq+27yVhFJdvOTgW21HNZxwEgRWQe8jVfc8ziQKCKlPbP5ccyygWxVnenGJ+GdCPw8XicDa1U1R1UPAu/jHT+/j1Wgio6Pr98DEbkMOBu4yJ2Q/I6pC94JfL777KcCc0Skrc9xgffZf189s/B+ibcKRVx1OfGHRafu7oz9ErBUVf8eMOsj4FL3+lK8sv9ao6p3qGqqqqbhHZuvVfUiYDow1se4tgBZItLDTRoBLMHf47UBOFpEGrn/Z2lMvh6rMio6Ph8Bl7gaK0cDuwOKhEJKRE7HK0ocqar7y8R6gYjEiUgnoBswqzZiUtWFqtpaVdPcZz8br/LFFnw8Vs6HeDd4EZHueBUbthOK4xWqGxfhMODdpV+Bdxf8Tp9iOB7vZ/cCYJ4bzsQrT58GrMS7k9/Cx+M0lJ9q9XR2H6pVwLu4Gga1HM8AIMMdsw+B5n4fL+CPwDJgEfA6Xg0LX44V8BbevYaDeInryoqOD94N+6fcd2AhMLgWY1qFVzZd+rl/NmD5O11My4EzavNYlZm/jp9u7tbKsarkeDUA/uU+Y3OA4aE6XtZkgzHGRJm6XNRjjDGmHJb4jTEmyljiN8aYKGOJ3xhjoowlfmOMiTKW+E3IiMjeMuOXiciT1dzGSKnBllXFa/nz2iCX3RvEMveJyK3u9Z9E5ORKlh0tIr2Cj7Zmuf1X64lVEflKfGo51oSOJX4TtkQkVlU/UtUHanCziUBQib+6VPUeVf2qkkVG47W06JffA09Xc53XCdHxMv6xxG98ISJpIvK1a/d8moh0cNNfEZFnRWQm8FDgrwQRmRcwHBCRk8Rrh/5Dt50fRaSfW/Y+1+b5DBFZIyK/dbt+AOjitvGwiDRx+58jIgtdA15VxX6niKwQke+BHgHTXxGRse71A+L1wbBARB4RkWOBkcDDbt9dRORqEZktIvNF5D3Xjkzpdv4pIj+42McG7OM2F+d8EXnATesiIlNEJFNEvhORnuXE3B0oUNXtAft4xh2zNSIy1B2vpSLySsCqHwEXBvdfNRGjNp40tCE6B6CYn57anIfX9MGTbt7HwKXu9RXAh+71K3jNHse48ctK1wnY7jnAd0B94AngXjd9ODDPvb4P+AHvCdtWwA63fBo/bwo3Fkhwr1vhPW1a+mDj3nLeUzreU52NgAS3/K0BsY/Fe4p2ecB2EgPnB2yrZcDr+4HrA5Z7F+/CrBewyk0/w72nRm689OncaUA39/oovOY3ysZ9OfBowPgreG00CV4DdHuAvm6fmcCAgGVXBsZqQ+QPpQ1MGRMKB1R1QOmIa7BrsBs9BjjPvX4dryORUu+qanF5GxSRbsDDwDBVPSgixwNjAFT1axFpKSIJbvFPVbUAKBCRbZTflLMAfxWRE/EaxUpxy22p4D2dAHygru0ZESmv/afdQD7wkng9m31Swbb6iMj9eMVPTfDaWS/1oaqWAEtEpDTuk4GXS/etqrnitfp6LPCu/NSQY1w5+0rGa+460MeqqiKyENiqqgvde1qMd4Kc55bbhtcq6Y4K3oeJMJb4TTjaV95El+TeAa7W4BrPKgh4XUz5n/eLgCQg3Z1I1gHx1Qv351S1SESG4DXmNhYYj/drpKxXgNGqOt+dFIdWEHtlTfPWw+sXYEAVYR3A6+QmUOk+Ssrsr4SfH6t4t76pI6yM3/jlB7xWQcFLvt8Fsc4EvCvewGW/c+sjIkOB7Vqmv4My8vC6wCzVDK9fgoMiMgzoWEUM3wKjRaShiDTFK3b6GXeCaqaqnwE34fX0VN6+mwKbxWu2+6Iq9gvwJXB5wL2AFu69rhWRX7hpIiL9y1l3KdA1iH2UfS8CtMVrzMzUEXbFb/xyPV4vW7/DK4K4vLKFRaQj3tVzdxG5wk2+Cq8sf4KILAD281PTxOVS1R0i8m/xOrn+HK8v2I9dcUcGXgucla0/R0Qm4vWBug2v+e+ymgKTRSQe72r9Zjf9beAFd6N5LHA3Xm9sOe5v03K2FbjvKSIyAMgQkULgM7xORS4CnhGRu/DuY7zt4gv0LfCoiIiqVqdlxnTgR/2ppzFTB1jrnMZECRF5HK9cv7Iqp+Wt85GqTgtdZKa2WVGPMdHjr3i1kapjkSX9useu+I0xJsrYFb8xxkQZS/zGGBNlLPEbY0yUscRvjDFRxhK/McZEmf8HQX7WhI0nC8kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s.y[2,:],s.y[3,:])\n",
    "plt.ylabel(\"Vertical distance (m)\")\n",
    "plt.xlabel(\"Horizontal distance (m)\")\n",
    "plt.title(\"Projectile hits at \" + str(s.y_events[0][0][2]))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_results(inputs, predictions):\n",
    "    \"convenience function for printing out results and computing mean square error\"\n",
    "\n",
    "    print(\"Target Point                   Predicted mean            Actual Value\")\n",
    "    print(\"------------------------------------------------------------------------------\")\n",
    "\n",
    "    error = 0.\n",
    "\n",
    "    for pp, m in zip(inputs, predictions):\n",
    "        trueval = simulator(pp)\n",
    "        print(\"{}      {}       {}\".format(pp, m, simulator(pp)))\n",
    "        error += (trueval - m)**2\n",
    "\n",
    "    print(\"Mean squared error: {}\".format(np.sqrt(error)/len(predictions)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "lhd = mogp_emulator.LatinHypercubeDesign([(-5., 1.), (0., 1000.)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 308 ms, sys: 0 ns, total: 308 ms\n",
      "Wall time: 307 ms\n"
     ]
    }
   ],
   "source": [
    "n_simulations = 50\n",
    "simulation_points = lhd.sample(n_simulations)\n",
    "%time simulation_output = np.array([simulator(p) for p in simulation_points])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAARxElEQVR4nO3de6xlZX3G8e/jcI1gATmdTFA6eEmt2jrYI/UWo1grYhu0sVbSKG1tx1ZptLVG1KTVxKZotdpbNGOlYuMNUSPeHZVqjQoOOsAgKoiYMhmZ8UKVJrWCv/6xX3Rz5lz2HGbtfc6830+yc9Z+11p7/d5Z+zxnzbvWXjtVhSSpH3ebdQGSpOky+CWpMwa/JHXG4Jekzhj8ktSZw2ZdwCROPPHE2rx586zLkKR15YorrvhOVc0tbF8Xwb9582Z27Ngx6zIkaV1J8q3F2h3qkaTOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4JekzqyLT+7eFZvP+9DMtn3j+U+e2bYlaSke8UtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpM4MFf5Kjklye5Mok1yR5RWt/S5JvJtnZHluGqkGStL8h7875I+D0qro1yeHAZ5N8pM17UVVdPOC2JUlLGCz4q6qAW9vTw9ujhtqeJGkyg47xJ9mQZCewF9heVZe1WX+T5Kokr0ty5BLrbk2yI8mOffv2DVmmJHVl0OCvqturagtwL+C0JA8GXgI8AHgYcALw4iXW3VZV81U1Pzc3N2SZktSVqVzVU1W3AJcCZ1TVnhr5EfBvwGnTqEGSNDLkVT1zSY5r00cDTwC+mmRTawvwFGDXUDVIkvY35FU9m4ALk2xg9Afmoqr6YJJPJZkDAuwE/mTAGiRJCwx5Vc9VwKmLtJ8+1DYlSSvzk7uS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM4Y/JLUGYNfkjoz5JetH5Xk8iRXJrkmySta+ylJLktyfZJ3JTliqBokSfsb8oj/R8DpVfUQYAtwRpKHA68CXldV9wO+Dzx7wBokSQsMFvw1cmt7enh7FHA6cHFrvxB4ylA1SJL2N+gYf5INSXYCe4HtwDeAW6rqtrbITcBJS6y7NcmOJDv27ds3ZJmS1JVBg7+qbq+qLcC9gNOABxzAutuqar6q5ufm5oYqUZK6M5WreqrqFuBS4BHAcUkOa7PuBeyeRg2SpJEhr+qZS3Jcmz4aeAJwLaM/AE9ri50DvH+oGiRJ+zts5UVWbRNwYZINjP7AXFRVH0zyFeCdSV4JfBl484A1SJIWGCz4q+oq4NRF2m9gNN4vSZoBP7krSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzQ37Z+r2TXJrkK0muSfL81v7yJLuT7GyPM4eqQZK0vyG/bP024IVV9aUkxwJXJNne5r2uql4z4LYlSUsY8svW9wB72vQPk1wLnDTU9iRJk5nKGH+SzcCpwGWt6dwkVyW5IMnx06hBkjQyePAnOQZ4D/CCqvoB8AbgvsAWRv8jeO0S621NsiPJjn379g1dpiR1Y9DgT3I4o9B/W1W9F6Cqbq6q26vqJ8CbgNMWW7eqtlXVfFXNz83NDVmmJHVlyKt6ArwZuLaq/n6sfdPYYk8Fdg1VgyRpf0Ne1fMo4JnA1Ul2traXAmcn2QIUcCPwnAFrkCQtMORVPZ8FssisDw+1TUnSyvzkriR1xuCXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOTBT8SR41SZskae2b9Ij/nyZskyStccvepC3JI4BHAnNJ/mJs1j2ADUMWJkkaxkp35zwCOKYtd+xY+w+Apw1VlCRpOMsGf1V9Gvh0krdU1bemVJMkaUCT3o//yCTbgM3j61TV6UMUJUkazqTB/27gjcC/ArcPV44kaWiTBv9tVfWGQSuRJE3FpJdzfiDJc5NsSnLCHY9BK5MkDWLSI/5z2s8XjbUVcJ+lVkhyb+CtwMa27Laq+of2B+NdjM4X3Ag8vaq+f2BlS5JWa6Lgr6pTVvHatwEvrKovJTkWuCLJduD3gU9W1flJzgPOA168iteXJK3CRMGf5FmLtVfVW5dap6r2AHva9A+TXAucBJwFPLYtdiHwHxj8kjQ1kw71PGxs+ijg8cCXGA3lrCjJZuBU4DJgY/ujAPBtRkNBi62zFdgKcPLJJ09YpiRpJZMO9fzZ+PMkxwHvnGTdJMcA7wFeUFU/SDL+upWkltjmNmAbwPz8/KLLSJIO3Gpvy/w/wIrj/kkOZxT6b6uq97bmm5NsavM3AXtXWYMkaRUmHeP/AKMrc2B0c7ZfAi5aYZ0Abwauraq/H5t1CaOrhM5vP99/gDVLku6CScf4XzM2fRvwraq6aYV1HgU8E7g6yc7W9lJGgX9RkmcD3wKePnm5kqS7atIx/k8n2cjPTvJeN8E6nwWyxOzHT1aeJOlgm/QbuJ4OXA78DqMj9MuSeFtmSVqHJh3qeRnwsKraC5BkDvgEcPFQhUmShjHpVT13uyP0m+8ewLqSpDVk0iP+jyb5GPCO9vx3gQ8PU5IkaUgrfefu/Rh90vZFSX4beHSb9XngbUMXJ0k6+FY64n898BKA9gGs9wIk+eU277cGrE2SNICVxuk3VtXVCxtb2+ZBKpIkDWql4D9umXlHH8Q6JElTslLw70jyxwsbk/wRcMUwJUmShrTSGP8LgPcl+T1+FvTzwBHAUwesS5I0kGWDv6puBh6Z5HHAg1vzh6rqU4NXJkkaxKT36rkUuHTgWiRJU+CnbyWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnBgv+JBck2Ztk11jby5PsTrKzPc4cavuSpMUNecT/FuCMRdpfV1Vb2sN7+kvSlA0W/FX1GeB7Q72+JGl1ZjHGf26Sq9pQ0PFLLZRka5IdSXbs27dvmvVJ0iFt2sH/BuC+wBZgD/DapRasqm1VNV9V83Nzc1MqT5IOfVMN/qq6uapur6qfAG8CTpvm9iVJUw7+JJvGnj4V2LXUspKkYUx0d87VSPIO4LHAiUluAv4aeGySLUABNwLPGWr7kqTFDRb8VXX2Is1vHmp7kqTJ+MldSeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdGSz4k1yQZG+SXWNtJyTZnuS69vP4obYvSVrckEf8bwHOWNB2HvDJqro/8Mn2XJI0RYMFf1V9BvjeguazgAvb9IXAU4baviRpcdMe499YVXva9LeBjUstmGRrkh1Jduzbt2861UlSB2Z2creqCqhl5m+rqvmqmp+bm5tiZZJ0aJt28N+cZBNA+7l3ytuXpO5NO/gvAc5p0+cA75/y9iWpe0NezvkO4PPALya5KcmzgfOBJyS5Dvj19lySNEWHDfXCVXX2ErMeP9Q2JUkr85O7ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6M9hXLy4nyY3AD4Hbgduqan4WdUhSj2YS/M3jquo7M9y+JHXJoR5J6sysgr+Ajye5IsnWGdUgSV2a1VDPo6tqd5KfB7Yn+WpVfWZ8gfYHYSvAySefPIsaJemQNJMj/qra3X7uBd4HnLbIMtuqar6q5ufm5qZdoiQdsqYe/EnunuTYO6aB3wB2TbsOSerVLIZ6NgLvS3LH9t9eVR+dQR2S1KWpB39V3QA8ZNrblSSNeDmnJHXG4Jekzhj8ktQZg1+SOmPwS1JnZnmTNg1k83kfmsl2bzz/yTPZrqQD4xG/JHXG4Jekzhj8ktQZg1+SOmPwS1JnvKpnQLO6umZWeusveCWT1ieP+CWpMwa/JHXG4Jekzhj8ktQZT+5Kd8EsT2jP6sSyfZ6uIfrsEb8kdcbgl6TOzCT4k5yR5GtJrk9y3ixqkKReTT34k2wA/gV4EvBA4OwkD5x2HZLUq1kc8Z8GXF9VN1TV/wHvBM6aQR2S1KVZXNVzEvBfY89vAn5t4UJJtgJb29Nbk3xtlds7EfjOKtddS+zH2jLzfuRVd/klZt6HA7VEn9ddP5awaD/u4n7+hcUa1+zlnFW1Ddh2V18nyY6qmj8IJc2U/VhbDoV+HAp9APuxGrMY6tkN3Hvs+b1amyRpCmYR/F8E7p/klCRHAM8ALplBHZLUpakP9VTVbUnOBT4GbAAuqKprBtzkXR4uWiPsx9pyKPTjUOgD2I8Dlqqa1rYkSWuAn9yVpM4Y/JLUmUM6+Nf6rSGS3Jjk6iQ7k+xobSck2Z7kuvbz+NaeJP/Y+nJVkoeOvc45bfnrkpwzhbovSLI3ya6xtoNWd5Jfbf8u17d1M8V+vDzJ7rZPdiY5c2zeS1pNX0vyxLH2Rd9n7QKGy1r7u9rFDAe7D/dOcmmSryS5JsnzW/u62h/L9GO97Y+jklye5MrWj1cst+0kR7bn17f5m1fbvwNSVYfkg9GJ428A9wGOAK4EHjjruhbUeCNw4oK2VwPntenzgFe16TOBjwABHg5c1tpPAG5oP49v08cPXPdjgIcCu4aoG7i8LZu27pOm2I+XA3+5yLIPbO+hI4FT2ntrw3LvM+Ai4Blt+o3Anw7Qh03AQ9v0scDXW63ran8s04/1tj8CHNOmDwcua/92i24beC7wxjb9DOBdq+3fgTwO5SP+9XpriLOAC9v0hcBTxtrfWiNfAI5Lsgl4IrC9qr5XVd8HtgNnDFlgVX0G+N4Qdbd596iqL9ToN+CtY681jX4s5SzgnVX1o6r6JnA9o/fYou+zdlR8OnBxW3/83+Sgqao9VfWlNv1D4FpGn45fV/tjmX4sZa3uj6qqW9vTw9ujltn2+H66GHh8q/WA+negdR7Kwb/YrSGWeyPNQgEfT3JFRreoANhYVXva9LeBjW16qf6slX4erLpPatML26fp3DYMcsEdQyQceD/uCdxSVbctaB9MGyY4ldFR5rrdHwv6AetsfyTZkGQnsJfRH9BvLLPtn9bb5v93q3XQ3/dDOfjXg0dX1UMZ3an0eUkeMz6zHWGtu+tt12vdzRuA+wJbgD3Aa2dazYSSHAO8B3hBVf1gfN562h+L9GPd7Y+qur2qtjC6K8FpwANmW9H+DuXgX/O3hqiq3e3nXuB9jN4kN7f/XtN+7m2LL9WftdLPg1X37ja9sH0qqurm9ov7E+BNjPYJHHg/vstoGOWwBe0HXZLDGYXl26rqva153e2PxfqxHvfHHarqFuBS4BHLbPun9bb5P9dqHfb3/WCf3FgrD0afSr6B0YmRO06CPGjWdY3Vd3fg2LHpzzEam/877nxS7tVt+snc+aTc5a39BOCbjE7IHd+mT5hC/Zu580nRg1Y3+59MPHOK/dg0Nv3njMZZAR7EnU+23cDoRNuS7zPg3dz5hN5zB6g/jMbdX7+gfV3tj2X6sd72xxxwXJs+GvhP4DeX2jbwPO58cvei1fbvgOoc6hdqLTwYXcHwdUZjbC+bdT0LartP22lXAtfcUR+j8b1PAtcBnxj75QujL7D5BnA1MD/2Wn/I6OTP9cAfTKH2dzD6b/ePGY0xPvtg1g3MA7vaOv9M+4T5lPrx763OqxjdQ2o8eF7WavoaY1e2LPU+a/v48ta/dwNHDtCHRzMaxrkK2NkeZ663/bFMP9bb/vgV4Mut3l3AXy23beCo9vz6Nv8+q+3fgTy8ZYMkdeZQHuOXJC3C4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfilVUjysHbjsKOS3L3de/3Bs65LmoQf4JJWKckrGX3y8mjgpqr62xmXJE3E4JdWqX2L0heB/wUeWVW3z7gkaSIO9Uird0/gGEbfGHXUjGuRJuYRv7RKSS5h9A1IpzC6edi5My5JmshhKy8iaaEkzwJ+XFVvT7IB+FyS06vqU7OuTVqJR/yS1BnH+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6sz/A5z6Ou+6fvl/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(simulation_output)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp = mogp_emulator.GaussianProcess(simulation_points, simulation_output)\n",
    "gp = mogp_emulator.fit_GP_MAP(gp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 12.6 ms, total: 12.6 ms\n",
      "Wall time: 1.81 ms\n"
     ]
    }
   ],
   "source": [
    "%time gp_output = gp.predict(simulation_points);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fbdacc14f70>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXqklEQVR4nO3df4zc9X3n8eebxSZ7+dE1YWvB2j6TxHJlxJ1NRkBFVOWoig09yQ6KUtCpWDlUVxeQklPPqt1EB01SQWol6aFLqMiBYno5DCWOsVLSrY8g5dQ7G9YxwRjOZcOP2ouD3RhDqlhg+973x3yWjNe769nd2fnheT6k0X7n/f1+Zz4fz3pf8/1+Pt+ZyEwkSd3tvFY3QJLUeoaBJMkwkCQZBpIkDANJEnB+qxswXRdddFEuXry41c2QpI6ye/fuf8rM/rH1jg2DxYsXMzQ01OpmSFJHiYhXx6t7mkiSZBhIkgwDSRKGgSQJw0CSRAfPJpKkbrJtzwibBvfz2rHjXNLXy/qVS1mzYqBhj28YSFKb27ZnhI1b93L8xCkARo4dZ+PWvQANCwRPE0lSm9s0uP/dIBh1/MQpNg3ub9hzGAaS1OZeO3Z8SvXpOGsYRMR7IuKpiPhJROyLiD8t9UsjYldEDEfEwxExt9QvKPeHy/rFNY+1sdT3R8TKmvqqUhuOiA0N650knQMu6eudUn066jkyeBu4NjP/NbAcWBURVwNfAb6emR8B3gBuLdvfCrxR6l8v2xERy4CbgMuAVcA3I6InInqAbwDXA8uAm8u2kiRg/cql9M7pOa3WO6eH9SuXNuw5zhoGWfXP5e6cckvgWuDRUt8MrCnLq8t9yvrfjogo9S2Z+XZmvgwMA1eW23BmvpSZ7wBbyraSJKqDxHfdeDkDfb0EMNDXy103Xt782UTl3ftu4CNU38X/FDiWmSfLJgeB0VYNAAcAMvNkRLwJfLDUd9Y8bO0+B8bUr5qgHeuAdQCLFi2qp+mSdE5Ys2KgoX/8x6prADkzT2XmcmAB1XfyvzFrLZq8HfdlZiUzK/39Z3wCqyRpmqY0mygzjwFPAr8J9EXE6JHFAmCkLI8ACwHK+l8Dfl5bH7PPRHVJUpPUM5uoPyL6ynIv8DvAC1RD4ZNls7XAY2V5e7lPWf/DzMxSv6nMNroUWAI8BTwNLCmzk+ZSHWTe3oC+SZLqVM+YwcXA5jJucB7wSGZ+PyKeB7ZExJeBPcD9Zfv7gb+KiGHgKNU/7mTmvoh4BHgeOAnclpmnACLidmAQ6AEeyMx9DeuhJOmsovqmvfNUKpX0m84kaWoiYndmVsbWvQJZkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRRRxhExMKIeDIino+IfRHx2VK/MyJGIuKZcruhZp+NETEcEfsjYmVNfVWpDUfEhpr6pRGxq9Qfjoi5je6oJGli9RwZnAT+KDOXAVcDt0XEsrLu65m5vNweByjrbgIuA1YB34yInojoAb4BXA8sA26ueZyvlMf6CPAGcGuD+idJqsNZwyAzD2Xmj8vyL4AXgIFJdlkNbMnMtzPzZWAYuLLchjPzpcx8B9gCrI6IAK4FHi37bwbWTLM/kqRpmNKYQUQsBlYAu0rp9oh4NiIeiIh5pTYAHKjZ7WCpTVT/IHAsM0+OqY/3/OsiYigiho4cOTKVpkuSJlF3GETE+4DvAp/LzLeAe4EPA8uBQ8BXZ6OBtTLzvsysZGalv79/tp9OkrrG+fVsFBFzqAbBdzJzK0Bmvl6z/lvA98vdEWBhze4LSo0J6j8H+iLi/HJ0ULu9JKkJ6plNFMD9wAuZ+bWa+sU1m30CeK4sbwduiogLIuJSYAnwFPA0sKTMHJpLdZB5e2Ym8CTwybL/WuCxmXVLkjQV9RwZXAP8PrA3Ip4ptT+hOhtoOZDAK8AfAmTmvoh4BHie6kyk2zLzFEBE3A4MAj3AA5m5rzzeHwNbIuLLwB6q4SNJapKovjHvPJVKJYeGhlrdDEnqKBGxOzMrY+tegSxJMgwkSYaBJAnDQJKEYSBJos6LziSp0bbtGWHT4H5eO3acS/p6Wb9yKWtWTPaxZ5pNhoGkptu2Z4SNW/dy/MQpAEaOHWfj1r0ABkKLeJpIUtNtGtz/bhCMOn7iFJsG97eoRTIMJDXda8eOT6mu2WcYSGq6S/p6p1TX7DMMJDXd+pVL6Z3Tc1qtd04P61cubVGL5ACypKYbHSR2NlH7MAwktcSaFQP+8W8jniaSJBkGkiTDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSSJOsIgIhZGxJMR8XxE7IuIz5b6hRGxIyJeLD/nlXpExD0RMRwRz0bEFTWPtbZs/2JErK2pfzQi9pZ97omImI3OSpLGV8+RwUngjzJzGXA1cFtELAM2AE9k5hLgiXIf4HpgSbmtA+6FangAdwBXAVcCd4wGSNnmD2r2WzXzrkmS6nXWMMjMQ5n547L8C+AFYABYDWwum20G1pTl1cCDWbUT6IuIi4GVwI7MPJqZbwA7gFVl3Qcyc2dmJvBgzWNJkppgSmMGEbEYWAHsAuZn5qGy6mfA/LI8AByo2e1gqU1WPzhOXZLUJHWHQUS8D/gu8LnMfKt2XXlHnw1u23htWBcRQxExdOTIkdl+OknqGnWFQUTMoRoE38nMraX8ejnFQ/l5uNRHgIU1uy8otcnqC8apnyEz78vMSmZW+vv762m6JKkO9cwmCuB+4IXM/FrNqu3A6IygtcBjNfVbyqyiq4E3y+mkQeC6iJhXBo6vAwbLurci4uryXLfUPJYkqQnq+aaza4DfB/ZGxDOl9ifA3cAjEXEr8CrwqbLuceAGYBj4JfBpgMw8GhFfAp4u230xM4+W5c8A3wZ6gR+UmySpSaJ6ur/zVCqVHBoaanUzJKmjRMTuzKyMrXsFsiTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJIEnN/qBkjnkm17Rtg0uJ/Xjh3nkr5e1q9cypoVA61ulnRWhoHUINv2jLBx616OnzgFwMix42zcuhfAQFDbO+tpooh4ICIOR8RzNbU7I2IkIp4ptxtq1m2MiOGI2B8RK2vqq0ptOCI21NQvjYhdpf5wRMxtZAelZtk0uP/dIBh1/MQpNg3ub1GLpPrVM2bwbWDVOPWvZ+bycnscICKWATcBl5V9vhkRPRHRA3wDuB5YBtxctgX4SnmsjwBvALfOpENSq7x27PiU6lI7OWsYZOaPgKN1Pt5qYEtmvp2ZLwPDwJXlNpyZL2XmO8AWYHVEBHAt8GjZfzOwZmpdkNrDJX29U6pL7WQms4luj4hny2mkeaU2AByo2eZgqU1U/yBwLDNPjqmPKyLWRcRQRAwdOXJkBk2XGm/9yqX0zuk5rdY7p4f1K5e2qEVS/aYbBvcCHwaWA4eArzaqQZPJzPsys5KZlf7+/mY8pVS3NSsGuOvGyxno6yWAgb5e7rrxcgeP1RGmNZsoM18fXY6IbwHfL3dHgIU1my4oNSao/xzoi4jzy9FB7fZSx1mzYsA//upI0zoyiIiLa+5+AhidabQduCkiLoiIS4ElwFPA08CSMnNoLtVB5u2ZmcCTwCfL/muBx6bTJknS9J31yCAiHgI+DlwUEQeBO4CPR8RyIIFXgD8EyMx9EfEI8DxwErgtM0+Vx7kdGAR6gAcyc195ij8GtkTEl4E9wP2N6pwkqT5RfXPeeSqVSg4NDbW6GZLUUSJid2ZWxtb9bCJJkh9Hoe7mZwlJVYaBupafJST9iqeJ1LX8LCHpVzwyUFepPS000dQJP0tI3cgwUNcYe1poIn6WkLqRp4nUNcY7LTSWnyWkbuWRgbrGZKd/ApxNpK5mGKhrXNLXy8g4gTDQ18vfb7i2BS2S2oenidQ1/IhpaWIeGahrjJ7+8SIz6UyGgc4J9V5J7EdMS+MzDNTxvJJYmjnHDNTxvJJYmjnDQB1voimjXkks1c8wUMeb6IphrySW6mcYqOM5ZVSaOQeQ1fGcMirNnGGgc4JTRqWZMQzUdvz2Man5DAO1Fa8ZkFrDAWS1Fa8ZkFrDMFBb8ZoBqTXOGgYR8UBEHI6I52pqF0bEjoh4sfycV+oREfdExHBEPBsRV9Tss7Zs/2JErK2pfzQi9pZ97omIaHQn1Tm8ZkBqjXqODL4NrBpT2wA8kZlLgCfKfYDrgSXltg64F6rhAdwBXAVcCdwxGiBlmz+o2W/sc6mLeM2A1BpnDYPM/BFwdEx5NbC5LG8G1tTUH8yqnUBfRFwMrAR2ZObRzHwD2AGsKus+kJk7MzOBB2seS11ozYoB7rrxcgb6egmqXzxz142XO3gszbLpziaan5mHyvLPgPlleQA4ULPdwVKbrH5wnPq4ImId1SMOFi1aNM2mq915zYDUfDMeQC7v6LMBbannue7LzEpmVvr7+5vxlJLUFaYbBq+XUzyUn4dLfQRYWLPdglKbrL5gnLokqYmmGwbbgdEZQWuBx2rqt5RZRVcDb5bTSYPAdRExrwwcXwcMlnVvRcTVZRbRLTWPJUlqkrOOGUTEQ8DHgYsi4iDVWUF3A49ExK3Aq8CnyuaPAzcAw8AvgU8DZObRiPgS8HTZ7ouZOToo/RmqM5Z6gR+UmySpiaJ6yr/zVCqVHBoaanUzJKmjRMTuzKyMrXsFsiTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScKvvdQE/B5iqbsYBjqD30MsdR9PE+kMfg+x1H0MA53B7yGWuo9hoDP4PcRS9zEMuti2PSNcc/cPuXTD33DN3T9k257qV0n4PcRS93EAuUvVM0jsbCKpexgGXaR2uuh5EZwa8/Hlo4PEo99B7B9/qXsYBl1i7JHA2CAY5SCx1J0cM+gS400XHY+DxFJ3Mgy6RD3v+B0klrqXYdAlJnrH3xNBAAN9vdx14+WOE0hdyjGDLrF+5dLTxgygeiRgAEgCw6BrOF1U0mQMgy7idFFJE3HMQJI0szCIiFciYm9EPBMRQ6V2YUTsiIgXy895pR4RcU9EDEfEsxFxRc3jrC3bvxgRa2fWJUnSVDXiyODfZObyzKyU+xuAJzJzCfBEuQ9wPbCk3NYB90I1PIA7gKuAK4E7RgNEktQcs3GaaDWwuSxvBtbU1B/Mqp1AX0RcDKwEdmTm0cx8A9gBrJqFdkmSJjDTMEjg7yJid0SsK7X5mXmoLP8MmF+WB4ADNfseLLWJ6pKkJpnpbKKPZeZIRPw6sCMi/m/tyszMiBj/Q3CmoQTOOoBFixY16mElqevN6MggM0fKz8PA96ie83+9nP6h/DxcNh8BFtbsvqDUJqqP93z3ZWYlMyv9/f0zabokqca0wyAi3hsR7x9dBq4DngO2A6MzgtYCj5Xl7cAtZVbR1cCb5XTSIHBdRMwrA8fXlZokqUlmcppoPvC9iBh9nP+RmX8bEU8Dj0TErcCrwKfK9o8DNwDDwC+BTwNk5tGI+BLwdNnui5l5dAbtkiRNUeQEn2vf7iqVSg4NDbW6GZLUUSJid82lAO/yCmRJkp9N1Cpf2LaXh3Yd4FQmPRHcfNVCvrzm8lY3S1KXMgxa4Avb9vLfd/7ju/dPZb5730CQ1AqeJmqBh3YdmFJdkmabYdACE30Z/UR1SZpthkEL9FSn49Zdl6TZZhi0wM1XLZxSXZJmmwPILTA6SOxsIkntwovOGmzbnhG/Z1hS25roojOPDBpk254R7ty+j2PHT7xbGzl2nI1b9wIYCJLammMGDbBtzwgbt+49LQhGHT9xik2D+1vQKkmqn2HQAJsG93P8xKkJ17927HgTWyNJU2cYNMDZ/thf0tfbpJZI0vQYBg0w2R/73jk9rF+5tImtkaSpMwwaYP3KpfTO6TmjPu9fzOGuGy938FhS23M2UQOM/rF3SqmkTmUYTGIq1wysWTHgH39JHcswmMC/+9b/4e9/+qtv3/SaAUnnMsOgxuiRwMgEs4NGrxkwDCSdawyDYuyRwES8ZkDSucjZRNQfBOA1A5LOTV19ZPCFbXv5zs5/ZCof1ec1A5LORV0bBlf92Q5e/8U7U9rnmg9f6HiBpHNS14XBVE4JjfL7BiSd69omDCJiFfBfgB7gv2Xm3Y1+jsUb/mZK258X8LVPLfdoQNI5ry0GkCOiB/gGcD2wDLg5IpY18jmmGgQXnH+eQSCpa7TLkcGVwHBmvgQQEVuA1cDzzW7Ie+f28Gef8POEJHWXdgmDAeBAzf2DwFVjN4qIdcA6gEWLFjW8EfPfP5ddn/+dhj+uJLW7tjhNVK/MvC8zK5lZ6e/vb+hjL/n19xoEkrpWuxwZjAALa+4vKLWm+Ivfc2xAUndrlyODp4ElEXFpRMwFbgK2N/IJXrn7dyesGwSSul1bHBlk5smIuB0YpDq19IHM3Nfo55koECSp27VFGABk5uPA461uhyR1o3Y5TSRJaiHDQJJkGEiSDANJEhCZU/k0//YREUeAV6e5+0XAPzWwOa1iP9rHudAHsB/tZjb68S8z84yrdjs2DGYiIoYys9LqdsyU/Wgf50IfwH60m2b2w9NEkiTDQJLUvWFwX6sb0CD2o32cC30A+9FumtaPrhwzkCSdrluPDCRJNQwDSVJ3hUFErIqI/RExHBEbWt2e8UTEKxGxNyKeiYihUrswInZExIvl57xSj4i4p/Tn2Yi4ouZx1pbtX4yItU1o9wMRcTginqupNazdEfHR8u8yXPaNJvbjzogYKa/JMxFxQ826jaVN+yNiZU193N+18jHtu0r94fKR7Y3uw8KIeDIino+IfRHx2VLvqNdjkn502uvxnoh4KiJ+Uvrxp5M9d0RcUO4Pl/WLp9u/KcnMrrhR/WjsnwIfAuYCPwGWtbpd47TzFeCiMbU/BzaU5Q3AV8ryDcAPgACuBnaV+oXAS+XnvLI8b5bb/VvAFcBzs9Fu4KmybZR9r29iP+4E/tM42y4rv0cXAJeW36+eyX7XgEeAm8ryXwL/YRb6cDFwRVl+P/APpa0d9XpM0o9Oez0CeF9ZngPsKv924z438BngL8vyTcDD0+3fVG7ddGRwJTCcmS9l5jvAFmB1i9tUr9XA5rK8GVhTU38wq3YCfRFxMbAS2JGZRzPzDWAHsGo2G5iZPwKOzka7y7oPZObOrP6veLDmsZrRj4msBrZk5tuZ+TIwTPX3bNzftfLu+Vrg0bJ/7b9Jw2Tmocz8cVn+BfAC1e8Z76jXY5J+TKRdX4/MzH8ud+eUW07y3LWv06PAb5e2Tql/U21nN4XBAHCg5v5BJv/FapUE/i4idkfEulKbn5mHyvLPgPlleaI+tUtfG9XugbI8tt5Mt5dTKA+Mnl5h6v34IHAsM0+Oqc+acophBdV3ox37eozpB3TY6xERPRHxDHCYaqj+dJLnfre9Zf2bpa2z+v+9m8KgU3wsM68Argdui4jfql1Z3ol13HzgTm13cS/wYWA5cAj4aktbU6eIeB/wXeBzmflW7bpOej3G6UfHvR6ZeSozl1P9fvcrgd9obYvO1E1hMAIsrLm/oNTaSmaOlJ+Hge9R/cV5vRyaU34eLptP1Kd26Wuj2j1SlsfWmyIzXy//mf8f8C2qrwlMvR8/p3oK5vwx9YaLiDlU/4B+JzO3lnLHvR7j9aMTX49RmXkMeBL4zUme+932lvW/Vto6u//fGz1Y0q43ql/x+RLVgZfRQZbLWt2uMW18L/D+muX/TfVc/yZOH/j787L8u5w+8PdUqV8IvEx10G9eWb6wCe1fzOkDrw1rN2cOWN7QxH5cXLP8H6metwW4jNMH9F6iOpg34e8a8NecPmj4mVlof1A9j/8XY+od9XpM0o9Oez36gb6y3Av8L+DfTvTcwG2cPoD8yHT7N6V2ztZ/qHa8UZ018Q9Uz9d9vtXtGad9Hyov5E+AfaNtpHq+8AngReB/1vyHDOAbpT97gUrNY/17qgNMw8Cnm9D2h6gesp+ges7y1ka2G6gAz5V9/ivl6vkm9eOvSjufBbaP+WP0+dKm/dTMqJnod628xk+V/v01cMEs9OFjVE8BPQs8U243dNrrMUk/Ou31+FfAntLe54D/PNlzA+8p94fL+g9Nt39TuflxFJKkrhozkCRNwDCQJBkGkiTDQJKEYSBJwjCQJGEYSJKA/w/vzek2lXuZhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(simulation_output,gp_output.mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation lengths = [7.55009206e-01 8.80164092e+02]\n",
      "Sigma = 10243.599711730027\n"
     ]
    }
   ],
   "source": [
    "print(\"Correlation lengths = {}\".format(np.sqrt(np.exp(-gp.theta[:2]))))\n",
    "print(\"Sigma = {}\".format(np.sqrt(np.exp(gp.theta[2]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_valid = 10\n",
    "validation_points = lhd.sample(n_valid)\n",
    "validation_output = np.array([simulator(p) for p in validation_points])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = gp.predict(validation_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWYElEQVR4nO3de4zd5Z3f8fdnjSFWLmsTJsjYTiEbi8hpWkNOgVVWFWUVbLyr2htFKfljcVO03jahTZSKBm+kZnNZhSzaXFBJVuyGjamycShhwaJQ102o0v6BYRw7GMO6nkAiPDh4NsaQNBYX99s/zmP2eOLLmfF4Lp73S/ppfuf7u5zn98D4c87vec6cVBWSpNnt16a6AZKkqWcYSJIMA0mSYSBJwjCQJAFnTXUDxuu8886rCy+8cKqbIUkzyrZt2/6uqgZG12dsGFx44YUMDg5OdTMkaUZJ8pNj1b1NJEkyDCRJhoEkCcNAkoRhIEliBs8mkqTZ5N7tw9yyeTfPHjzEBfPnceOKi1lzyaIJO79hIEnT3L3bh1l/z04OvXIYgOGDh1h/z06ACQsEbxNJ0jR3y+bdrwXBEYdeOcwtm3dP2HMYBpI0zT178NCY6uNhGEjSNHfB/Hljqo+HYSBJ09yNKy5m3tw5R9XmzZ3DjSsunrDncABZkqa5I4PEziaSpFluzSWLJvQf/9G8TSRJMgwkSYaBJAnDQJKEYSBJwjCQJNFHGCR5XZJHkvwwya4kn271byR5OsmOtixv9SS5NclQkseSXNpzrrVJ9rRlbU/93Ul2tmNuTZLTcK2SpOPo53MGLwFXVdUvkswF/neSB9u2G6vq7lH7XwMsbcvlwNeAy5OcC3wK6AAFbEuyqaqeb/v8AbAVeABYCTyIJGlSnPSdQXX9oj2c25Y6wSGrgTvbcQ8D85MsBFYAW6rqQAuALcDKtu1NVfVwVRVwJ7Bm/JckSRqrvsYMksxJsgPYT/cf9K1t05+0W0FfSnJOqy0Cnuk5fG+rnai+9xj1Y7VjXZLBJIMjIyP9NF2S1Ie+wqCqDlfVcmAxcFmSfwisB94B/BPgXOATp6uRPe24vao6VdUZGBg43U8nSbPGmGYTVdVB4CFgZVXta7eCXgL+Cris7TYMLOk5bHGrnai++Bh1SdIk6Wc20UCS+W19HvBe4G/bvX7azJ81wOPtkE3AdW1W0RXAC1W1D9gMXJ1kQZIFwNXA5rbtxSRXtHNdB9w3kRcpSTqxfmYTLQQ2JJlDNzzuqqr7k3wvyQAQYAfwr9v+DwCrgCHgl8CHAKrqQJLPAo+2/T5TVQfa+oeBbwDz6M4iciaRJE2idCfwzDydTqcGBwenuhmSNKMk2VZVndF1P4EsSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIk+vsO5NcleSTJD5PsSvLpVr8oydYkQ0m+neTsVj+nPR5q2y/sOdf6Vt+dZEVPfWWrDSW56TRcpyTpBPp5Z/AScFVV/WNgObCyfdH9F4AvVdXbgeeB69v+1wPPt/qX2n4kWQZcC7wTWAl8Ncmc9t3KtwHXAMuAD7Z9JUmT5KRhUF2/aA/ntqWAq4C7W30DsKatr26Padt/O0lafWNVvVRVTwNDwGVtGaqqp6rqZWBj21eSNEn6GjNor+B3APuBLcCPgINV9WrbZS+wqK0vAp4BaNtfAN7cWx91zPHqx2rHuiSDSQZHRkb6abokqQ99hUFVHa6q5cBiuq/k33E6G3WCdtxeVZ2q6gwMDExFEyTpjDSm2URVdRB4CPhNYH6Ss9qmxcBwWx8GlgC07b8O/Ky3PuqY49UlSZOkn9lEA0nmt/V5wHuBJ+mGwvvbbmuB+9r6pvaYtv17VVWtfm2bbXQRsBR4BHgUWNpmJ51Nd5B50wRcmySpT2edfBcWAhvarJ9fA+6qqvuTPAFsTPI5YDvw9bb/14H/nGQIOED3H3eqaleSu4AngFeBj1TVYYAkNwCbgTnAHVW1a8KuUJJ0Uum+aJ95Op1ODQ4OTnUzJGlGSbKtqjqj634CWZJkGEiSDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSfT3HchLkjyU5Ikku5J8tNX/OMlwkh1tWdVzzPokQ0l2J1nRU1/ZakNJbuqpX5Rka6t/u30XsiRpkvTzzuBV4N9X1TLgCuAjSZa1bV+qquVteQCgbbsWeCewEvhqkjntO5RvA64BlgEf7DnPF9q53g48D1w/QdcnSerDScOgqvZV1Q/a+s+BJ4FFJzhkNbCxql6qqqeBIeCytgxV1VNV9TKwEVidJMBVwN3t+A3AmnFejyRpHMY0ZpDkQuASYGsr3ZDksSR3JFnQaouAZ3oO29tqx6u/GThYVa+Oqh/r+dclGUwyODIyMpamS5JOoO8wSPIG4DvAx6rqReBrwG8Ay4F9wJ+djgb2qqrbq6pTVZ2BgYHT/XSSNGuc1c9OSebSDYJvVtU9AFX1XM/2vwDubw+HgSU9hy9uNY5T/xkwP8lZ7d1B7/6SpEnQz2yiAF8HnqyqL/bUF/bs9nvA4219E3BtknOSXAQsBR4BHgWWtplDZ9MdZN5UVQU8BLy/Hb8WuO/ULkuSNBb9vDN4D/D7wM4kO1rtj+jOBloOFPBj4A8BqmpXkruAJ+jORPpIVR0GSHIDsBmYA9xRVbva+T4BbEzyOWA73fCRJE2SdF+YzzydTqcGBwenuhmSNKMk2VZVndF1P4EsSTIMJEmGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkkSf32cg6eTu3T7MLZt38+zBQ1wwfx43rriYNZec6BtipenDMJAmwL3bh1l/z04OvXIYgOGDh1h/z04AA0EzgreJpAlwy+bdrwXBEYdeOcwtm3dPUYuksTEMpAnw7MFDY6pL041hIE2AC+bPG1Ndmm4MA2kC3LjiYubNnXNUbd7cOdy44uIpapE0NicNgyRLkjyU5Ikku5J8tNXPTbIlyZ72c0GrJ8mtSYaSPJbk0p5zrW3770mytqf+7iQ72zG3JsnpuFjpdFlzySI+/753sWj+PAIsmj+Pz7/vXQ4ea8Y46XcgJ1kILKyqHyR5I7ANWAP8S+BAVd2c5CZgQVV9Iskq4N8Cq4DLga9U1eVJzgUGgQ5Q7TzvrqrnkzwC/DtgK/AAcGtVPXiidvkdyJI0duP+DuSq2ldVP2jrPweeBBYBq4ENbbcNdAOCVr+zuh4G5rdAWQFsqaoDVfU8sAVY2ba9qaoerm4y3dlzLknSJBjTmEGSC4FL6L6CP7+q9rVNPwXOb+uLgGd6Dtvbaieq7z1G/VjPvy7JYJLBkZGRsTRdknQCfYdBkjcA3wE+VlUv9m5rr+hPfL9pAlTV7VXVqarOwMDA6X46SZo1+gqDJHPpBsE3q+qeVn6u3eI5Mq6wv9WHgSU9hy9utRPVFx+jLkmaJP3MJgrwdeDJqvpiz6ZNwJEZQWuB+3rq17VZRVcAL7TbSZuBq5MsaDOPrgY2t20vJrmiPdd1PeeSJE2Cfv420XuA3wd2JtnRan8E3AzcleR64CfAB9q2B+jOJBoCfgl8CKCqDiT5LPBo2+8zVXWgrX8Y+AYwD3iwLZKkSXLSqaXTlVNLJWnsxj21VJJ05jMMJEmGgSTJMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSSJ/r4D+Y4k+5M83lP74yTDSXa0ZVXPtvVJhpLsTrKip76y1YaS3NRTvyjJ1lb/dpKzJ/ICJUkn1887g28AK49R/1JVLW/LAwBJlgHXAu9sx3w1yZwkc4DbgGuAZcAH274AX2jnejvwPHD9qVyQJGnsThoGVfV94MDJ9mtWAxur6qWqehoYAi5ry1BVPVVVLwMbgdVJAlwF3N2O3wCsGdslSJJO1amMGdyQ5LF2G2lBqy0CnunZZ2+rHa/+ZuBgVb06qn5MSdYlGUwyODIycgpNlyT1Gm8YfA34DWA5sA/4s4lq0IlU1e1V1amqzsDAwGQ8pSTNCmeN56Cqeu7IepK/AO5vD4eBJT27Lm41jlP/GTA/yVnt3UHv/pKkSTKudwZJFvY8/D3gyEyjTcC1Sc5JchGwFHgEeBRY2mYOnU13kHlTVRXwEPD+dvxa4L7xtEmSNH4nfWeQ5FvAlcB5SfYCnwKuTLIcKODHwB8CVNWuJHcBTwCvAh+pqsPtPDcAm4E5wB1Vtas9xSeAjUk+B2wHvj5RFydJ6k+6L85nnk6nU4ODg1PdDEmaUZJsq6rO6LqfQJYkGQaSJMNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRJ9hEGSO5LsT/J4T+3cJFuS7Gk/F7R6ktyaZCjJY0ku7Tlmbdt/T5K1PfV3J9nZjrk1SSb6IiVJJ9bPO4NvACtH1W4CvltVS4HvtscA1wBL27IO+Bp0w4PudydfDlwGfOpIgLR9/qDnuNHPJUk6zU4aBlX1feDAqPJqYENb3wCs6anfWV0PA/OTLARWAFuq6kBVPQ9sAVa2bW+qqoer+2XMd/acS5I0ScY7ZnB+Ve1r6z8Fzm/ri4Bnevbb22onqu89Rv2YkqxLMphkcGRkZJxNlySNdsoDyO0VfU1AW/p5rturqlNVnYGBgcl4SkmaFcYbBs+1Wzy0n/tbfRhY0rPf4lY7UX3xMeqSpEk03jDYBByZEbQWuK+nfl2bVXQF8EK7nbQZuDrJgjZwfDWwuW17MckVbRbRdT3nkiRNkrNOtkOSbwFXAucl2Ut3VtDNwF1Jrgd+Anyg7f4AsAoYAn4JfAigqg4k+SzwaNvvM1V1ZFD6w3RnLM0DHmyLJGkSpXvLf+bpdDo1ODg41c2QpBklybaq6oyu+wlkSZJhIEkyDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkcYphkOTHSXYm2ZFksNXOTbIlyZ72c0GrJ8mtSYaSPJbk0p7zrG3770my9njPJ0k6PSbincE/q6rlPV+jdhPw3apaCny3PQa4BljalnXA16AbHnS/V/ly4DLgU0cCRJI0OU7HbaLVwIa2vgFY01O/s7oeBuYnWQisALZU1YGqeh7YAqw8De2SJB3HqYZBAf89ybYk61rt/Kra19Z/Cpzf1hcBz/Qcu7fVjleXJE2Ss07x+N+qquEkbwG2JPnb3o1VVUnqFJ/jNS1w1gG89a1vnajTStKsd0rvDKpquP3cD/wN3Xv+z7XbP7Sf+9vuw8CSnsMXt9rx6sd6vturqlNVnYGBgVNpuiSpx7jDIMnrk7zxyDpwNfA4sAk4MiNoLXBfW98EXNdmFV0BvNBuJ20Grk6yoA0cX91qkqRJciq3ic4H/ibJkfP8dVX9tySPAncluR74CfCBtv8DwCpgCPgl8CGAqjqQ5LPAo22/z1TVgVNolyRpjFI1Ybf0J1Wn06nBwcGpboYkzShJtvV8FOA1fgJZkmQYSJJOfWqpgHu3D3PL5t08e/AQF8yfx40rLmbNJX5UQtLMYRiconu3D7P+np0ceuUwAMMHD7H+np0ABoKkGcPbRKfols27XwuCIw69cphbNu+eohZJ0tgZBqfo2YOHxlSXpOnIMDhFF8yfN6a6JE1HhkGf7t0+zHtu/h4X3fRfec/N3+Pe7d2/mHHjiouZN3fOUfvOmzuHG1dcPBXNlKRxcQC5D/0MEjubSNJMZhj04USDxGsuWfTaIkkzlbeJ+uAgsaQznWHQBweJJZ3pDIM+OEgs6UznmEEfHCSWdKYzDPrkILGkM9msDgP/wJwkdc3aMPAPzEnS35t1YXDk3cDwMaaF9n52QJJmk2kTBklWAl8B5gB/WVU3T/RzXP4nW3ju5y+fcB8/OyBpNpoWU0uTzAFuA64BlgEfTLJsIp/jvV/8nycNAvCzA5Jmp2kRBsBlwFBVPVVVLwMbgdUT+QR79v/fk+7jZwckzVbTJQwWAc/0PN7bakdJsi7JYJLBkZGRiW3A/Hl8/n3vcrxA0qw0bcYM+lFVtwO3A3Q6nZqo8375Xyw3BCTNatPlncEwsKTn8eJWmzBL3/L6Y9bPf+PZBoGkWW+6hMGjwNIkFyU5G7gW2DSRT7Dl41f+SiAsfcvr2frJ907k00jSjDQtbhNV1atJbgA2051aekdV7Zro59ny8Ssn+pSSdEaYFmEAUFUPAA9MdTskaTaaLreJJElTyDCQJBkGkiTDQJIEpGrCPrs1qZKMAD8Z5+HnAX83gc05U9lP/bOv+mM/9ed09tM/qKqB0cUZGwanIslgVXWmuh3Tnf3UP/uqP/ZTf6ain7xNJEkyDCRJszcMbp/qBswQ9lP/7Kv+2E/9mfR+mpVjBpKko83WdwaSpB6GgSRpdoVBkpVJdicZSnLTVLdnsiS5I8n+JI/31M5NsiXJnvZzQasnya2tjx5LcmnPMWvb/nuSrO2pvzvJznbMrUkyuVc4MZIsSfJQkieS7Ery0Va3r3okeV2SR5L8sPXTp1v9oiRb27V9u/05epKc0x4Pte0X9pxrfavvTrKip37G/K4mmZNke5L72+Pp2U9VNSsWun8a+0fA24CzgR8Cy6a6XZN07f8UuBR4vKf2p8BNbf0m4AttfRXwIBDgCmBrq58LPNV+LmjrC9q2R9q+acdeM9XXPM5+Wghc2tbfCPwfYJl99Sv9FOANbX0usLVd013Ata3+58C/aesfBv68rV8LfLutL2u/h+cAF7Xfzzln2u8q8HHgr4H72+Np2U+z6Z3BZcBQVT1VVS8DG4HVU9ymSVFV3wcOjCqvBja09Q3Amp76ndX1MDA/yUJgBbClqg5U1fPAFmBl2/amqnq4uv/n3tlzrhmlqvZV1Q/a+s+BJ+l+F7d91aNd7y/aw7ltKeAq4O5WH91PR/rvbuC32zui1cDGqnqpqp4Ghuj+np4xv6tJFgO/A/xlexymaT/NpjBYBDzT83hvq81W51fVvrb+U+D8tn68fjpRfe8x6jNae4t+Cd1XvfbVKO3Wxw5gP92w+xFwsKpebbv0Xttr/dG2vwC8mbH330z0ZeA/AP+vPX4z07SfZlMY6Djaq1TnGDdJ3gB8B/hYVb3Yu82+6qqqw1W1nO73lV8GvGNqWzT9JPldYH9VbZvqtvRjNoXBMLCk5/HiVputnmu3LWg/97f68frpRPXFx6jPSEnm0g2Cb1bVPa1sXx1HVR0EHgJ+k+5tsiPfnth7ba/1R9v+68DPGHv/zTTvAf55kh/TvYVzFfAVpms/TfXgymQtdL/i8ym6AzBHBlveOdXtmsTrv5CjB5Bv4ehB0T9t67/D0YOij7T6ucDTdAdEF7T1c9u20YOiq6b6esfZR6F7H//Lo+r21dH9MQDMb+vzgP8F/C7wXzh6YPTDbf0jHD0weldbfydHD4w+RXdQ9Iz7XQWu5O8HkKdlP015J03yf5BVdGeI/Aj45FS3ZxKv+1vAPuAVuvcVr6d7L/K7wB7gf/T8YxXgttZHO4FOz3n+Fd3BqyHgQz31DvB4O+Y/0T7ZPtMW4Lfo3gJ6DNjRllX21a/00z8Ctrd+ehz4j63+NrphN9T+wTun1V/XHg+17W/rOdcnW1/spmdm1Zn2uzoqDKZlP/nnKCRJs2rMQJJ0HIaBJMkwkCQZBpIkDANJEoaBJAnDQJIE/H8eteTbz4LaegAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(validation_output, predictions.mean)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Target Point                   Predicted mean            Actual Value\n",
      "------------------------------------------------------------------------------\n",
      "[ -3.00247231 500.24522638]      2013.7262830499285       2019.7659045735888\n",
      "[3.48210286e-01 7.02039887e+02]      11.96257048327243       2.3268476049186186\n",
      "[ -1.70204227 280.04888571]      101.41500365988031       139.13238287755163\n",
      "[ -4.15967925 974.84025616]      21035.43852115173       20449.836403470068\n",
      "[-1.06180863 68.07306775]      54.262952719436726       25.79549210931367\n",
      "[  0.77394004 674.20696712]      6.264385850925464       0.9267904241747875\n",
      "[ -0.70247942 153.61449563]      21.0349623519578       16.157892684661302\n",
      "[ -3.2677743 442.3562146]      3003.3950544493814       3021.710989325538\n",
      "[ -4.73866385 897.2140366 ]      37588.357348681166       40915.61586507107\n",
      "[ -2.37660086 374.32960132]      556.6393130187503       564.1552108234388\n",
      "Mean squared error: 337.8814234940752\n"
     ]
    }
   ],
   "source": [
    "print_results(validation_points, predictions.mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_predict = 10000\n",
    "prediction_points = lhd.sample(n_predict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "hm = mogp_emulator.HistoryMatching(gp=gp, coords=prediction_points, obs=[2000., 400.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "nroy_points = hm.get_NROY()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ruled out 9884 of 10000 points\n"
     ]
    }
   ],
   "source": [
    "print(\"Ruled out {} of {} points\".format(n_predict - len(nroy_points), n_predict))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyMElEQVR4nO3de3hU5bX48e8iREgECYI3AggqonJHvFSqgvYIrYqIt3o5itViL2pFTyyc+lM81Uobj1hqL4cWi22pFRUjWlsKCloV0ISES1CUQ1WI2CIaVAjHGNbvj70nJGEueyazZ/aeWZ/n4Ulmz86ed4BkZb9rvesVVcUYY4yJp0O2B2CMMSb4LFgYY4xJyIKFMcaYhCxYGGOMSciChTHGmIQ6ZnsAfujZs6f269cv28MwxphQqaqq+lBVD4n2XE4Gi379+lFZWZntYRhjTKiIyLuxnrNpKGOMMQlZsDDGGJOQBQtjjDEJWbAwxhiTkG/BQkQeFpF/icj6FscOFpElIvK2+7G7e1xEZLaIbBKRtSIyssXXXOOe/7aIXOPXeI0xxsTm553FPGB8m2PTgOdVdQDwvPsY4KvAAPfPFOCX4AQX4C7gFOBk4K5IgDEmGRXVdYye+QL9p/2Z0TNfoKK6LttDMiZUfCudVdWXRKRfm8MXAGPczx8BlgPfd4//Tp0WuCtFpEREjnDPXaKqHwGIyBKcAPSoX+M24VVRXUf54o28X99Ar5IiysYNZOKIUiqq65i+cB0NjU0A1NU3MH3hOgAmjijN5pCNCY1Mr7M4TFW3uZ9/ABzmfl4KbGlx3lb3WKzjxjQHh7r6BgRo2Wy/rr6BsifWMGNRLfUNjft9bUNjE+WLN1qwMMajrCW43buItG2mISJTRKRSRCq3b9+ersuaaNYugFmDYUaJ83HtgowPIXK3UFffAET/j9TYpFEDRURdfYNNSRnjUaaDxT/d6SXcj/9yj9cBfVqc19s9Fuv4flR1jqqOUtVRhxwSdbW6SYe1C+CZm2HnFkCdj8/cnPGAUb54Y/O0UntEpqQsYBgTX6aDxSIgUtF0DfB0i+NXu1VRpwI73emqxcA5ItLdTWyf4x4z2fL8f0FjQ+tjjQ3O8Qx6v74h8UkeRaakjDGx+ZazEJFHcRLUPUVkK05V00xggYhcB7wLXOqe/hzwNWATsBu4FkBVPxKRHwKvu+f9VyTZbbJk59bkjvukV0lR8xRUOqQz+BiTi/yshro8xlNnRzlXge/GuM7DwMNpHJppj2693SmoKMczaOxxhzB/5XtpS3r1KilK05Xy2NoFzh3mzq3O/4ez74Shlyb+OhMKtoI7H7UnQX32nVDY5gdrYZFzPEMqqut4sqoufdURQNm4gWm8Wh4KSC7L+MeCRb5p7zf10Evh/NnQrQ8gzsfzZ2f0N8h0JbcjRGy9RbsFJJdl/JOT+1mYOOJ9U3v9gT/00qxOL6Q7v3DlKX3Ter28FJBclvGP3Vnkmxz4pk5nfuGqU/tyz8Qhabte3oqVs8pwLsv4x4JFvsmBb+qycQMpKixIy7UsUKRJAHJZxl8WLPJNjnxTd+rY/v+6pVYBlT4ByGUZf1nOIt9EvnlDWuLYtilgqooKC4JTAZUrJadZzmUZf1mwyEch/qZORyVUgQj3TRoSjAqoSHVapOggUp0Gof03MrnJpqFMqKSjEmqvajACBVjJqQkNCxYmVNJRCRWo1do5UJ1m8oMFCxMqY49rX0fhQOUqICeq00x+sGBhQmXZm6nvVVJaUhScXEVEjlSnmdxnCe58FsIqnFRzFqUlRbwy7aw0jyYNQl6dZvKHBYt8FdIqnFRbkwe6BXmIq9NM/rBpqHwV0iqcsnEDKewgSX9doJLaJv8EYCvi9rI7i3wV5iqcJGNFYQcJVlLb5JeQ3sW3ZXcW+SqkVTjlizfS2JTcThZdOncMVlLb5JeQ3sW3ZcEiX4W0CieV3EP97kYfRmJayYFpFt+E+S6+BQsW+Sqkjd9SyT1YvsJntktefCG9i2/LgkU+G3opTF0PM+qdjwEPFJD89qeSwteYJOXINItvQnoX35YFCxMqE0eUUlJU6Pl8xbZM9V2OTLP4JqR38W1ZNZQJnRkTBnluU55Xe1Zka5Flt97uFFSU48aRA2tp7M7ChM7EEaVcdGLiu4XA9YHyUzbzBjkyzWLis2BhQqeiuo4nq+rinhOoPSsyIZt5gxyZZjHx2TSUCZ1EGyAVFgjlFw/Ln0AB2c8b5MA0i4nP7ixMqFRU1yXsDXXZSX2yGyiyseYgR8ozTXBZsDChEdl/O5Enq+qoqI4/TeWbbOUOLG9gfGbBwoSG1/23GxqbKF+8MQMjiiJbuQPLGxifWc7ChEYyrT6y1pI8m7kDyxsYH9mdhQmNZNp2ZK3FRzZyB9aXyWSABQsTGmXjBnrqTp7V9RWZzh1YXyaTIRYsTGhMHFFKoubkWd9nO9O5A+vLZDIkKzkLEZkKXI/TumcdcC1wBPAnoAdQBfy7qn4uIp2A3wEnAjuAy1T1nWyM22RfSVEh9Q3RW44HZp/tTOYOsr2+wuSNjN9ZiEgpcDMwSlUHAwXA14EfA7NU9RjgY+A690uuAz52j89yzzN5qrFpb8znAr3Ptl9sfYXJkGxNQ3UEikSkI1AMbAPOAp5wn38EmOh+foH7GPf5s0Uk+U2YTehVVNex6/PYpbN5uW+Fra8wGZLxYKGqdcD9wHs4QWInzrRTvap+4Z62FYhMOpcCW9yv/cI9v0fb64rIFBGpFJHK7du3+/smvLAKlbRLtHYib5oGtmTrK0yGZDxnISLdce4W+gP1wOPA+PZeV1XnAHMARo0aldwmzemWIxu0B02iaabKdz/Kr35QEba+wmRANqahvgL8Q1W3q2ojsBAYDZS401IAvYFIv4Y6oA+A+3w3nER3cFmFii9KiuNvejR/1XsZGokJJbvbb5dsBIv3gFNFpNjNPZwNbACWARe751wDPO1+vsh9jPv8C6qa3TuHRKxCJe0qquv4bM8Xcc8J+P8Kk022HqXdspGzWIWTqF6NUzbbAWf66PvArSKyCScnMdf9krlAD/f4rcC0TI85aVahknblizfSuNeigUmR3e23W1bWWajqXcBdbQ5vBk6Ocu4e4JJMjCttzr6zdc4CrEKlnbyUxRYX2hpTE4Pd7bebfXf5wSpU0q5bUfx8BcCPJg3NwEhMKNndfrtZ11m/5EKFytoFzm36zq3ON9XZd2btPXlZWZOXlVDGG7vbb7e4wUJEOgPnAacDvYAGYD3wZ1Wt9X94JmsCVv5bvzt6i4+I0nxckGe8i/yfDcgvP2EUM1iIyN04gWI5sAr4F9AZOBaY6QaS21R1bQbGaTItXkIwC99gJcWFfBwjYGS1y6wJj1y428+ieHcWr7mJ6GgeEJFDgb4+jMkEQcASgrHKYgWy22XWmDwRM8Gtqn9ue0xEOojIQe7z/1LVSj8HZ7IoYAnBnTE6zVoxrTGZkbAaSkT+KCIHiciBOPmKDSJS5v/QTFYFrEFdvCaB0xeuo6K6LubzaWMrgE0e81I6e4KqfoLTBfYvOD2d/t3PQZkACFj5b9m4gRQVFkR9rqGxKWGTwXazFcAmz3kpnS0UkUKcYPGQqjaKiN3954MAJQQjOYlbHquJ+rzve1kELOFvTKZ5ubP4H+Ad4EDgJRE5EvjEz0EZE83EEaUxS2R938siYAl/YzItZrAQkS+JiKjqbFUtVdWvuQ383gPGZm6IxuwTbTrKc+lse3IOAUv4G5Np8e4srgaqRORPIjJZRA4HUEf89p/G+GTiiFLumzSE0pIiBGcxnqfS2fbmHAKW8Dcm02LmLFT12wAichzwVWCeiHTDaSX+V+AVVY29x6UxaVRRXcfdz9Q2L8wrKSpk1mXDva+vaG/OwVYAmzyXMMGtqm8CbwKzRKQIZwrqEuABYJS/wzPGCRRlT6yhsWlfXUV9QyNlj68BPPaESkfOIUAJf2MyzVPXWRHpLiJDgeOBD4DfqqoFCpMR5Ys3tgoUEY171XvJrOUcTBgEeC1PwjsLEfkhMBlnv4m97mEFzvJvWMbsE68sts5ryax1HTVBF7DmnW15WWdxKXC0qn7u92CMiaZXSVHMoFDgpXc5WM7BBF/A1/J4CRbrgRKcrrPGZFzZuIExF+M1JbPxtuUcTJAFfC2Pl2BxH1AtIuuB/4scVNUJvo3KGI9sHwuTM7r1dku7oxwPAC/B4hHgx8A69uUsjMmYWElsAdvHIl8EaNdG3wQ8r+YlWOxW1dm+j8SYGGIluBXbSjUvBDzx2y5tg+CwK+DtvwUyKHoJFn8XkfuARbSehlrt26iMaSFWgtumoPJEwBO/KYsWBNf8MavdnePxEixGuB9PbXHMSmdNxpSNG8j0hetoaNzXMMC2Us0jAU/8pixkQdDLCm5rGmiyKjLVVL54I+/XN9CrpIiycQNtCipfBDzxm7KQBcGYwUJErgL+qKpRk9oicjRwhKq+7NfgjImYOKLUgkO+CnjiN2UhC4Lx7ix64JTMVgFVwHagM3AMcCbwITDN9xEaY/Jbri6oDFkQFI2zqElECnByE6OBI4AG4A3gL6r6XkZGmIJRo0ZpZWVltodhjDHxBawkWESqYvX9i5uzcFuQL3H/GGMSCdg3vwm4EHUV8FINZUzmhPmHbS6vBzB5z4KFCY4YP2xff+djbtkwIPiVUCErhTQmGQn3s3DzFsb4L8YP215VP6GuvgHFaUk+feE6KqrrsjLEuEJWCmlMMrxsfvS2iJSLyAnpelERKRGRJ0TkTRF5Q0S+JCIHi8gSEXnb/djdPVdEZLaIbBKRtSIyMl3jMAET44fqEexo9bihscn7pkeZZBssmRzmJVgMA94CfiMiK0Vkiogc1M7X/SnwV1U9zr3+GzhluM+r6gDgefaV5X4VGOD+mQL8sp2vbYIqxg/V97XH/se8bnqUSWff6ZQ+thTgUkhjkpEwWKjqp6r6a1U9Dfg+cBewTUQeEZFjkn1BEekGnAHMda//uarWAxfgdLjF/TjR/fwC4HfqWAmUiMgRyb6uCYEoP2wb6MRPvth/vr9XEPtCDb3U6evTrQ8gzseA9vkxJlletlUtAM4FrgX6Af8NzAdOB54Djk3yNfvjLPD7rYgMw1nw9z3gMFXd5p7zAXCY+3kp0HKZ41b32LYWxxCRKTh3HvTt2zfJIQVYmKuDkhVl8dX6o29iyetHwt6Q9IUKUSmkMcnwUg31NrAMKFfVV1scf0JEzkjxNUcCN6nqKhH5KW1WgquqikgSW6CBqs4B5oCzKC+FcQVPPpZitvlhexJwX5866wtlTJZ5yVlcrarXtQwUIjIaQFVvTuE1twJbVXWV+/gJnODxz8j0kvsxso1rHdCnxdf3do/lvnilmHlk4ohSysYNpFdJEe/XN1C+eGMwq6GMyWFegkW0jY9+luoLquoHwBYRicwjnA1swNkv4xr32DXA0+7ni4Cr3aqoU4GdLaarsmvtApg1GGaUOB/XLkjv9a0UE4A7KtYx9bGacJTPGpOj4nWd/RJwGnCIiNza4qmDgPauvbgJmC8iBwCbcfIhHYAFInId8C4QmYt4DvgasAnY7Z6bfZmYIgpZV0o/VFTX8YeV+7chi5TP2nSUMZkRL2dxANDFPadri+OfABe350VVtQaI1qzq7CjnKvDd9ryeLzKxWjdkXSn9cPcztTGfC2T5rDE5KmawUNUXgRdFZJ6qvpvBMYVDJqaIcrU1cxI+3t0Y87lAls8ak6PiTUM9qKq3AA9Fq0xS1Ql+DizwMjVFZKWYMQW2fNbknnwqYY8h3jTU792P92diIKFjU0S+i5fALirsYPkKkxn5WMIeRbxpqCr300qgIbK9qrtIr1MGxhZsNkXku3j9n+6bNDSDIzF5zboJA94W5T0PfAX4zH1cBPwNp1Iqv9kUka/iJbDtrsJkjJWwA97WWXRW1UigwP282L8hGeOIlcAu9TOx7ffaGRM+1k0Y8BYsdrVsCy4iJ+LsxW2Mr8rGDaSosPWSHl/7QkXmpnduAXTf3LQFjPxm3YQBb9NQtwCPi8j7gACHA5f5OShjYN9UU8b6QtnctInG8pMAiLPmLcFJIoVA5Ne5jaoau/g9AEaNGqWVlZXZHoYJmxklQLTvB4EZ9ZkdizFZICJVqhptwbSnFuWFwLdx9qAAWC4i/xP0gGFM0qy9ijExeclZ/BI4EfiF++dEbLc6k4tsbjo8rBAh47zkLE5S1WEtHr8gImv8GpAxWWNz0+Fgi+SywkuwaBKRo1X1fwFE5CigKcHXGBNOtnYm+KwQISu8BIsyYJmIbMaphjqSoLQJNzntjop1PLpqC02qFIhw+Sl9uGfikPS/kPX9CRdbJJcVCYOFqj4vIgNoXQ31f/4Oy+S7OyrWtdrHokm1+XFaA4ZNaYSPFSJkRcwEt4hMivwBzgWOcf+c6x4zxjfzV+2/4VG84ymzrWvDxwoRsiLencX5cZ5TYGGax2JMs1jLfzwsC0qOTWmEjxUiZEW8rrOWlzC5z6Y0wskKETIu4ToLETlMROaKyF/cxye4+2Qbk14taudf6XQzEzq8vN8pxYVelgYlwaY0jPHEy3fePGAx0Mt9/BZOvyhj0qdNE79S+ZCZhb9pFTA6CPwo3ftYDL0Uzp8N3foA4nw8f7b91mpMG15KZ3uq6gIRmQ6gql+IiK2zMOkVJdFcLJ9ze8cFLPr8ywhwxSl9/WkiaFMaxiTktUV5D9wOayJyKrDT11HlC2tZsE+MhHIv2QE4//mWvbk9gwMyxrTk5c7iNmARcLSIvAIcAlzs66jygdX3txYj0fy+9tj3eZyd84wx/kp4Z+HuxX0mzjaqNwCDVHWt3wPLeVbf31qURPNuPYCffLEvcMbaOc8Y4z8v1VBrgduBPaq63lqTp4nV97fWItGsCHXak2mN17No75cBn3fIM8Yk5GUa6nycnfEWiMhe4DFggaqmeSltnrH6/v25iWYBXq+uo2rxRiQTO+QZYxLytFNe88lOj6j/B1ypqgWJzs+WUOyU1zZnAc40jJVtmkyw5okminbtlOde4Eicu4vLcNqT356+4eUpa1lgssWKK0wKvGyrugooBBYAl6jqZt9HlS+svj+qiuo6yhdv5H2bgvKH7QdhUuDlzuJqVd3o+0iMwQkU0xeuo6HRWfdZV9/A9IXrACxgpIsVV5gUeCmdtUBhMqZ88cbmQBHR0NhE+WL7b5g2sYoo8rm4wiSU5q5s3olIgYhUi8iz7uP+IrJKRDaJyGMicoB7vJP7eJP7fL9sjdn4L9bCO88L8mxVfGLWPNGkIGvBAvge8EaLxz8GZqnqMcDHQKSz7XXAx+7xWe55Jkd1KyqMetzTgrw2zQibE7cWMFqz5okmBV6roU4D+rU8X1V/l+qLikhvnN337gVuFREBzgKucE95BJgB/BK4wP0c4AngIRERTabm14RCRXUduz7/Yr/jhR3E24I8S9x6Z8UVJkleqqF+DxwN1OCUzYLT1y3lYAE8iFN+29V93AOoV9XIT4qtQCSbWQpsgeaOtzvd8z9sM84pwBSAvn37tmNoJlvKF2+ksWn/3wG6dO7oLbkd9MStrW0wIeblzmIUcEK6fpMXkfOAf6lqlYiMScc1AVR1DjAHnEV56bquyZxYeYn63R47zAR5VbytbTAh5yVnsR44PI2vORqYICLvAH/CmX76KVAiIpHg1Ruocz+vA/oAuM93A3akcTz7syRpVsTKS3huIBjkxK01jjQhFzNYiMgzIrII6AlsEJHFIrIo8ifVF1TV6araW1X7AV8HXlDVK4Fl7Gt9fg3wtPv5Ivcx7vMv+JqvsCRp1pSNG0hRYesuMkk1EAxy4jboU2TGJBBvGur+jI3C8X3gTyJyD1ANzHWPzwV+LyKbgI9wAox/LEmaNZG8RLtWbwc1cRvkKTJjPIgZLFT1RXDWPwDbVHWP+7gIOCwdL66qy4Hl7uebgZOjnLMHuCQdr+eJ/QaYNXdUrOPRVVtoUqVAhLHHHZI7q7bPvjN648ggTJEZ44GXnMXjwN4Wj5vcY7nJVrdmxR0V6/jDyvdocmcYm1T5w8r3uKNiXZZHliZBniIzxgMv1VAdVfXzyANV/Tyyujon2W+AWTF/ZfTtUeavfI97Jg7J8Gh8ksoUmZXbmoDwcmexXUQmRB6IyAW0WeOQU+w3wKyIVbGQ1zXQVmxhAsTLncW3gPki8hAgOAvkrvZ1VNkW1CSpyS9WbGECJGGwUNX/BU4VkS7u4898H5XJKxXVdYlPykdWbGECxEu7j07ARbi9oZw2TqCqtprIpEW89uPdi6M3FswLVm5rAsRLzuJpnGZ+XwC7WvwxJi3itR+/6/xBGRxJwAR5RbrJO15yFr1VdbzvIzF5q6S4kI+j9H8qKSrMnXUWqcjVfdqtwiuUvASLV0VkiKrmSMG7CZI7KtZFDRSFBcKMCXl8VxGRa8UW1lAxtLxMQ30ZqBKRjSKyVkTWichavwdmcl9kIV40Bx7gsS25CRdrqBhaXu4svur7KEzeqaiui7kQD2Bng8e25CZcrMIrtLzcWWiMP8akrHzxxrj/iTy3JTfhYu10khOg7RK83Fn8GSc4CNAZ6A9sBGxC2aSsLk4FlID3tuTGH34loa2djncBy+8kvLNQ1SGqOtT9OACnM+wK/4dmclVFdR0S5/krT+1r+Yps8rPNiLXT8S5g+R0vdxatqOpqETnFj8GY3FdRXcdtC9bEnIK66tS+udM4MKz8bjOSaxVefglYfsfLCu5bWzzsAIwE3vdtRCZnVVTXUfb4muY25NFYoAiAgP2QylsBW8HvJcHdtcWfTjg5jAv8HJTJTf+5cC2Ne2MHigKJNzllMsaS0MEQsBX8XhoJ3p2JgZjcVlFdx+7GvXHPiXfHYTLIktDBELAV/F6moQ4BbsepfuocOa6qZ/k4LpNj4jULjCi1ctlgCNgPqbwWoPyOlwT3fOAx4DycvS2uAbb7OahAsn427RKvVDbCymUDJEA/pEwweMlZ9FDVuUCjqr6oqt8A8uuuwnYsaxcv+1V0L87zpoHGBJyXYBHpu7BNRM4VkRHAwT6OKXgCVu8cNommoIoKC/K7FbkxIeBlGuoeEekG3Ab8DDgIuMXPQQWOlRK2S7z9KgDumzTE7ipM+9g0se+8rOB+VlV3qup6VR2rqicCR2dgbMFhpYTtEq/PU2lJkQUK0z42TZwRXqahork18Sk5JGD1zmFTNm4ghQX7r6Eo7CCW1DbtZ9PEGZFqsMiv1VPWzyYpFdV1jJ75Av2n/ZnRM18AoPziYa320y4pKqT8kmF2V2Haz6aJMyLp3lCu/Fs9ZaWEnkRaekRWatfVN1D2+BrKLxlG9Z3nZHl0JicFrC1Grop5ZyEin4rIJ1H+fAr0yuAYTXtloCd+5G7ilsdq9mvp0bhXmbGoNu2vaQxg08QZEvPOQlW7ZnIgxicZ6IlfUV3H9IXraGhsinlOve18Z/xiK84zItVpKBMWfrebxllHES9QGOM7myb2XaoJbhMWGUj+JVpHAbRKbqcsQFtMGpNvLFjkugysEfGyX3a7V2hbLb0xWZXxYCEifURkmYhsEJFaEfmee/xgEVkiIm+7H7u7x0VEZovIJhFZKyIjMz3mjPHjN+cMJP/Kxg2kqLAg7jntLpG1WnpjsiobOYsvgNvc7Vm7AlUisgSYDDyvqjNFZBowDfg+8FVggPvnFOCX7sfc4lciOs3Jv4rqOsoXb+T9+gZ6lRRRNm5gcyC4bUH0XfDS0nrcaumNyaqMBwtV3QZscz//VETeAEpxdt8b4572CLAcJ1hcAPxOVRVYKSIlInKEe520iveD0Hd+JqLTlPxrW/VUV9/A9IXrgH13Dm2roooKC9KzSttq6Y3JqqxWQ4lIP2AEsAo4rEUA+AA4zP28FGj5U2Kre6xVsBCRKcAUgL59+yY9Fi8/CH0V0N+cWwbQDiL73Tk0NDZRvngjE0eUNv89+RJwbfc2Y7Iqa8FCRLoATwK3qOon0mL/ZVVVEUlqlbiqzgHmAIwaNSrpFebRyj8bGpuYsag2M3cbAfzNuaK6jrIn1tDY5Px1xtr2tGU1VMugkVZWS5+7rGNsKGQlWIhIIU6gmK+qC93D/4xML4nIEcC/3ON1QJ8WX97bPZZWsco/6xsamxeU1dU3MPWxGirf/Yh7Jg5J7wAC+Jvz3c/UNgeKeLxUQ6WF1dLnngwsGjXpkfFgIc4txFzgDVV9oMVTi3C2bJ3pfny6xfEbReRPOIntnX7kK7oVFXpaZazA/JXvMepIZ/+ntN11BOQ358i0k5dtUCGNOQmTnzKwaNSkRzbuLEYD/w6sE5Ea99h/4gSJBSJyHfAuEPmf8hzwNWATsBu41o9BSRJ9dBXnt+49jXuTynEkTKBn+TdnL207IgQyXwRgck9Ac3Vmf9mohnqZ2C3Oz45yvgLf9XVQQP3u5HoXfRzl/JbJ3raynkBPoKK6Lmbpa1slRYXU3GUdZE0aBDBXZ6Kz3lCuXiVFnqde4omV+4iVQI/sT521kl32BTIvgQJgxgTbL9ukSQBzdSnJgyS9tftweVmFHFFUWEBJUfReR7GSvbGCSOQOo66+AW3xuKI67Tn8mJJpBNi9uDAQd0ImR+TCxmJ50orG7ixcbdcIdCsqRMSZniopLkQVdjY0Nv/mD8ktQIt151IgEvOOI1M/lL00AgQoLJD293gypq2wV7nlSZLegkULqawR8Dp9VDZuYNTgEus3+sgP8HSuKo91LS9TcAceUMC9Fw5JewBrbGxk69at7NmzJ63XNcHRuXNnevfuTWFhGjoPp0O6p4zyJElvwaIdkgkusVY3xypT7VVSlNakeLxrRQtkglP1VepzDmXr1q107dqVfv36IcmUpJlQUFV27NjB1q1b6d+/f7aH48+6jjxJ0luwyKBYwSXWdFa8pHgqd0CxrvXKtLOaz8l0kn3Pnj0WKHKYiNCjRw+2b9+e7aE4/JgyypUkfQIWLLIsXj+lqY/VRP0arzkGL18TOe5bmw4PLFDktkD9+/oxZRSQBbV+s2ARALF+UMfKJaTSXiOd1zImtPyaMgp7kt4DK50NsGjlvKm210jntbKporqO0TNfoP+0PzN65gtpKTEWEW677bbmx/fffz8zZswAYMaMGZSWljJ8+HBOOOEEHn300ebzVJV77rmHAQMGcOyxxzJ27Fhqa2upra3l2GOPpaFhX3A+99xzW31tqk477bSE5zz44IPs3r273a+VkzKwGViusmARYBNHlHLfpCGUlhQhOMnm+yalVpGUzmslLU07AEaS9Olek9KpUycWLlzIhx9+GPX5qVOnUlNTw9NPP80NN9xAY6Ozev/nP/85r776KmvWrOGtt95i+vTpTJgwgaOPPppJkyZx7733OuOuqKCxsZHLL7+8XeMEePXVVxOeY8EijlxY15ElNg0VcOnMJWQlL5HG6pN0Jvxb6tixI1OmTGHWrFnNP+CjGTBgAMXFxXz88ccceuih/PjHP+bFF1+kuLgYgHPOOYfTTjuN+fPnc+eddzJixAguvvhipk2bxjPPPLPf9ebNm8dTTz3Fzp07qaur46qrruKuu+4C4IEHHuDhhx8G4Prrr+eWW24BoEuXLnz22WcsX76cGTNm0LNnT9avX8+JJ57IH/7wB372s5/x/vvvM3bsWHr27MnSpUu57rrrqKysRET4xje+wdSpU1P+u8oJeTBl5AcLFia92tawf74rbdUniZL07fHd736XoUOHcvvtt8c8Z/Xq1QwYMIBDDz2UTz75hF27dnHUUUe1OmfUqFHU1tZSXFzM/fffzxlnnMGtt97KgAEDol7ztddeY/369RQXF3PSSSdx7rnnIiL89re/ZdWqVagqp5xyCmeeeSYjRoxo9bXV1dXU1tbSq1cvRo8ezSuvvMLNN9/MAw88wLJly+jZsydVVVXU1dWxfv16AOrr69v3F2Xylk1DmfSJ1vag4aPo56ZQfRIrGZ+OJP1BBx3E1VdfzezZs/d7btasWQwaNIhTTjmFH/zgB56vef7551NSUsJ3vvOdmOf827/9Gz169KCoqIhJkybx8ssv8/LLL3PhhRdy4IEH0qVLFyZNmsTf//73/b725JNPpnfv3nTo0IHhw4fzzjvv7HfOUUcdxebNm7npppv461//ykEHHeR5/Ma0ZMHCpE+0GvZYUqg+8TtJf8sttzB37lx27drV6vjUqVOpra3lySef5LrrrmPPnj0cdNBBHHjggWzevLnVuVVVVQwatK8lSocOHejQIfa3Wduy0mTKTDt16tT8eUFBAV988cV+53Tv3p01a9YwZswYfvWrX3H99dd7vr4xLVmwMOnj9W4hxeoTv5P0Bx98MJdeeilz586N+vyECRMYNWoUjzzyCABlZWXcfPPNzVVPS5cu5eWXX+aKK67w/JpLlizho48+oqGhgYqKCkaPHs3pp59ORUUFu3fvZteuXTz11FOcfvrpnq/ZtWtXPv30UwA+/PBD9u7dy0UXXcQ999zD6tWrPV/HmJYsZ2HSJ1YNe9HBcMCBaVmw5HeS/rbbbuOhhx6K+fydd97JFVdcwTe/+U1uuukmPv74Y4YMGUJBQQGHH344Tz/9NEVF3qfFTj75ZC666CK2bt3KVVddxahRowCYPHkyJ598MuAkuNvmK+KZMmUK48ePp1evXjz44INce+217N27F4D77rvP83WMaUnU4x4GYTJq1CitrKzM9jDyT9vKJ3DuIuKUJr7xxhscf/zxGRpgsMybN4/Kysq4wSljdn8En26Dps+h4ADoegQUH5y2y+fzv3OYiEiVqo6K9pxNQ5n0sRr2cNr9kXNH2PS587jpc+fx7hjFCSYv2TSUSS+rYfds8uTJTJ48OdvDcO4odG/rY7rXOZ7GuwsTbnZnYUy+i9xReD1u8pIFizBLUxsNk+cKDkjuuMlLFizCKk/2/TUZ0PUIkDY/CqSDc9wYlwWLsIq3iYsxySg+2ClGiNxJFBzgPLZ8hWnBgkVY5cm+v/vxYert3nvvZdCgQQwdOpThw4ezatUqwFnfsGHDhnZfH6Bfv34xu9pG/OhHP2r12Es7ci8mT55M//79GT58OCNHjmTFihX7n1R8MBw2CHqN4LQLv5kwUFhn2/xjwSKsYrXLyLF9f1vxYeptxYoVPPvss6xevZq1a9eydOlS+vTpA8BvfvMbTjjhhDQNPrG2wcJLO3KvysvLqampYebMmdxwww1xz7U26CHlcw7TgkVY5eMmLj5MvW3bto2ePXs291nq2bMnvXr1AmDMmDFEFnd26dKFsrIyBg0axFe+8hVee+01xowZw1FHHcWiRYsAZ5HdjTfe2Hzt8847j+XLl+/3mhMnTuTEE09k0KBBzJkzB4Bp06bR0NDA8OHDufLKK5tfE5xNlsrKyhg8eDBDhgzhscceA2D58uWMGTOGiy++mOOOO44rr7ySRItszzjjDDZt2gQ4bdAHDx7M4MGDefDBB5vPibxurOvPnj27uQ362LFjaWpqYvLkyc3jmzVrlvd/AJMeGchhWrAIq3xcAOfD1Ns555zDli1bOPbYY/nOd77Diy++GPW8Xbt2cdZZZ1FbW0vXrl254447WLJkCU899RR33plcgH744YepqqqisrKS2bNns2PHDmbOnElRURE1NTXMnz+/1fkLFy6kpqaGNWvWsHTpUsrKyti2bRvgtCl/8MEH2bBhA5s3b+aVV16J+9rPPPMMQ4YMoaqqqrkN+sqVK/n1r39NdXX1fudHu/7NN99Mr169WLZsGcuWLaOmpqa5Dfq6deu49tprk/r7MGmQgRymBYswG3opTF0PM+qdj7kcKMCXqbcuXbpQVVXFnDlzOOSQQ7jsssuYN2/efucdcMABjB8/HoAhQ4Zw5plnUlhYyJAhQ6K2Bo9n9uzZDBs2jFNPPZUtW7bw9ttvxz3/5Zdf5vLLL6egoIDDDjuMM888k9dffx3w1qYcnKaHw4cPZ86cOcydO9faoOeaDOQwLViY8PBp6q2goIAxY8Zw991389BDD/Hkk0/ud05hYWFz+/AOHTo0T1t16NChuTV4x44dmxv2AezZs2e/6yxfvpylS5eyYsUK1qxZw4gRI6Ke55WXNuWwL2exZMkSBg8enNbr52Qb9LCtYcpADtOChQkPH6beNm7c2Oo3+5qaGo488siUrtWvXz9qamrYu3cvW7Zs4bXXXtvvnJ07d9K9e3eKi4t58803WblyZfNzhYWFzft7t3T66afz2GOP0dTUxPbt23nppZeaO9KmytqgxxHGNUwZyGFabygTLmnuPfXZZ59x0003UV9fT8eOHTnmmGOak87JGj16NP379+eEE07g+OOPZ+TIkfudM378eH71q19x/PHHM3DgQE499dTm56ZMmcLQoUMZOXJkq7zFhRdeyIoVKxg2bBgiwk9+8hMOP/xw3nzzzZTGCTBy5Ehrgx5LvPn/oE71RsbVckvjdmwFEE1oWpSLyHjgp0AB8BtVnRnrXGtRHh7Wujo/hOrfeUYJEO3nojj5wRwW+hblIlIA/Bz4KnACcLmIZK4A3hiTP/JxDZMHoQgWwMnAJlXdrKqfA38CLsjymIwxuSgf1zB5EJZgUQq03K9zq3usmYhMEZFKEancvn17Rgdn2icsU6EmNaH7983HNUwe5EyCW1XnAHPAyVlkeTjGo86dO7Njxw569OjRXJpqcoeqsmPHDjp37pztoSTHNvHaT1iCRR3Qp8Xj3u4xE3K9e/dm69at2N1g7urcuTO9e+f3fH8uCEuweB0YICL9cYLE14Ersjskkw6FhYX0798/28MwxiQQimChql+IyI3AYpzS2YdVtTbLwzLGmLwRimABoKrPAc9lexzGGJOPwlINZYwxJotCs4I7GSKyHXg3hS/tCcTfzixccun92HsJJnsvwZTqezlSVQ+J9kROBotUiUhlrKXuYZRL78feSzDZewkmP96LTUMZY4xJyIKFMcaYhCxYtJZab+rgyqX3Y+8lmOy9BFPa34vlLIwxxiRkdxbGGGMSsmBhjDEmIQsWbYjIDBGpE5Ea98/Xsj2m9hKR20RERaRntseSKhH5oYisdf9N/iYivbI9plSJSLmIvOm+n6dEpCTbY2oPEblERGpFZK+IhK70VETGi8hGEdkkItOyPZ72EJGHReRfIrI+3de2YBHdLFUd7v4JdYsREekDnAO8l+2xtFO5qg5V1eHAs0CYd6JZAgxW1aHAW8D0LI+nvdYDk4CXsj2QZOXgLpzzgPF+XNiCRe6bBdxO9E2FQ0NVP2nx8EBC/H5U9W+q+oX7cCVOy/3QUtU3VHVjtseRopzahVNVXwI+8uPaFiyiu9GdInhYRLpnezCpEpELgDpVXZPtsaSDiNwrIluAKwn3nUVL3wD+ku1B5LGEu3AaR2i6zqaTiCwFDo/y1A+AXwI/xPnN9YfAf+N8QwdSgvfynzhTUKEQ772o6tOq+gPgByIyHbgRuCujA0xCovfinvMD4AtgfibHlgov78fktrwMFqr6FS/nicivcebHAyvWexGRIUB/YI27XWlvYLWInKyqH2RwiJ55/XfB+eH6HAEOFonei4hMBs4DztYQLHZK4t8mbGwXTo9sGqoNETmixcMLcZJ3oaOq61T1UFXtp6r9cG6vRwY1UCQiIgNaPLwAeDNbY2kvERmPk0eaoKq7sz2ePNe8C6eIHICzC+eiLI8pkGwFdxsi8ntgOM401DvADaq6LZtjSgcReQcYpaqhbMEsIk8CA4G9OO3nv6WqofwNUEQ2AZ2AHe6hlar6rSwOqV1E5ELgZ8AhQD1Qo6rjsjqoJLjl8Q+ybxfOe7M7otSJyKPAGJwW5f8E7lLVuWm5tgULY4wxidg0lDHGmIQsWBhjjEnIgoUxxpiELFgYY4xJyIKFMcaYhCxYmFARkc/SdJ15IvIPEVkjIm+JyO9ExLceTSJyrIg8JyJvi8hqEVkgIoeleK1yt8truYgcIiKrRKRaRE53X6Mkztd+S0SuTvF1+4nIFal8rQk/K501oSIin6lqlzRcZx7wrKo+Ic4S91uAb+N0g/28zbkFqtrUjtfqDKwDblXVZ9xjY4APVTXpRZ8ishM4WFWbROTrwFdU9fpUx5fE644B/kNVz/P7tUzw2J2FCSVxlIvIehFZJyKXucc7iMgv3P0ilri/aV8c71rqmAV8gNOqGhH5TET+W0TWAF8SkTtF5HX39ea4AQYROanFPhvlMfYRuAJYEQkU7msuV9X1ItJZRH7rvodqERnrXrfAvd7r7vVvcI8vAroAVSLyfeAnwAXu6xeJyDvi7lsiIle7X7vGXWwa2a/lP9zPjxaRv4pIlYj8XUSOc4/PE5HZIvKqiGxu8fc3Ezjdfa2pSf+jmVDLy95QJidMwllpPwxnterrIvISMBroh7M3waHAG8DDHq+5GjgOeBqnDfoqVb0NQEQ2qOp/uZ//Hqev0zPAb4FvquoKEZkZ47qDgaoYz30XJ14NcX9Y/01EjgWuBnaq6kki0gl4RUT+pqoT3Lur4e5Y/omzMv9G9zHux0HAHcBpqvqhiBwc5bXn4KyEf1tETgF+AZzlPncE8GX372MR8AQwDbuzyFsWLExYfRl41J0e+qeIvAic5B5/XFX3Ah+IyLIkriktPm8CnmzxeKyI3A4UAwcDtSLyd6Crqq5wz/kjThBJ9n38DEBV3xSRd4FjcboFD23xW303YADwD4/XPQvn7+FD99qt9jgQkS7AacDjkQCD04IkosL9O9yQam7F5BYLFsbsMwJ43v18TyRP4eYcfoHzG/wWEZkBdE7iurXAmUmORYCbVHVxkl/nVQegPnKHEsX/tRmLyXOWszBh9XfgMndu/xDgDOA14BXgIjd3cRhOU7W43PzHzThTL3+NckokMHzo/kZ+MYCq1gOfulM44HQsjeaPwGkicm6L1zxDRAa77+NK99ixQF9gI7AY+LaIFEaeE5EDE72XFl4ALhGRHu7Xt5qGcnce/IeIXOI+LyIyLME1PwW6JjEGk0MsWJiwegpYC6zB+cF4u9t+/UmcduwbgD/g5CF2xrhGuZvAfgtnCmts20ooaA4Kv8ZpV78Yp611xHXAr0WkBifPsd9rqWoDzvTUTeKUzm4AvgNsx7lj6SAi64DHgMmq+n/Ab9z3sNpNmv8PScwEqGotcC/wovseH4hy2pXAde7ztSTeTnQt0OQmzC3BnWesdNbkHBHpoqqfub9VvwaM9msfj8hruZ9PA45Q1e/58VrGZJPlLEwuelachWkHAD/0ecOnc8XZ5rUjzj4bk318LWOyxu4sjDHGJGQ5C2OMMQlZsDDGGJOQBQtjjDEJWbAwxhiTkAULY4wxCf1/00JqhfkTHc4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(prediction_points[nroy_points,0], prediction_points[nroy_points,1], \"o\", label=\"NROY points\")\n",
    "plt.plot(simulation_points[:,0], simulation_points[:,1],\"o\", label=\"Simulation Points\")\n",
    "plt.xlabel(\"log Drag Coefficient\")\n",
    "plt.ylabel(\"Launch velocity (m/s)\")\n",
    "plt.legend()\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
}