notebooks/solve_inverse_linear.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w-Tp5UzlOzso"
      },
      "source": [
        "# Solving Inverse Problems - Linear Example\n",
        "\n",
        "Consider a scenario where we have some observed output data, and our goal is to determine the underlying conditions (inputs) that led to these outcomes.\n",
        "\n",
        "This can be accomplished by performing experiments to obtain some sample inputs and ouputs, using twinlab to train a model on the data, and applying twinlab's `Emulator.calibrate` function to predict the missing inputs.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SFnNvQv_Ozsp"
      },
      "source": [
        "First, import the required libraries.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KTjc23BCOzsq",
        "outputId": "268287fe-4e72-48d0-880a-395d2d848b95"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "          ====== TwinLab Client Initialisation ======\n",
            "          Version     : 2.10.0\n",
            "          User        : alexander@digilab.co.uk\n",
            "          Server      : https://twinlab.digilab.co.uk/v3\n",
            "          Environment : /Users/mead/digiLab/twinLab/tutorials/.env\n",
            "\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.patches as patches\n",
        "import twinlab as tl"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y_bl4UnROzsq"
      },
      "source": [
        "## Establish Parameters and Create Datasets\n",
        "\n",
        "We want to establish data and model parameters, and also generate the noisy observed and experimental data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "_xWC4akAOzsq"
      },
      "outputs": [],
      "source": [
        "# Data and model parameters\n",
        "random_seed = 8675309\n",
        "err_sig = 2 / 3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HewUUBlXSZ9W"
      },
      "source": [
        "In order to solve an inverse problem (using MCMC behind the scenes), there are a few essential ingredients.\n",
        "\n",
        "**An observation**: The algorithm computes a posterior distribution of the inputs that is most likely to have resulted in an observation. To begin with, we need to have an observation that we are interested in.\n",
        "\n",
        "**The standard deviation of the observation**: one of the assumptions in solving a Bayesian inverse problem is that the measurement error/noise level of the observation is known. This may come from the known precision of a sensor module, or multiple measurements can be taken and the mean and standard deviation provided instead.\n",
        "\n",
        "In our example, we can simply choose a noise level to add to the observation.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 143
        },
        "id": "qmDICytPOzsr",
        "outputId": "ba05df58-2ba3-419f-ecf7-39009d04f515"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>10.25</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "       y\n",
              "0  10.25"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0.666667</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "          y\n",
              "0  0.666667"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Observed datapoint we would like to know the inputs for.\n",
        "obs = pd.DataFrame({\"y\": [10.25]})\n",
        "display(obs)\n",
        "obs_std = pd.DataFrame({\"y\": [err_sig]})\n",
        "display(obs_std)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jb_6C-30UToZ"
      },
      "source": [
        "We need a forward model before we can start solving calibration problems. Let's imagine we've sampled from our simulation, done some experiments, or collected some data:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 700
        },
        "id": "0t4iPTOlOzsr",
        "outputId": "902bd0b0-51cf-4f48-92f5-096e49be3af2"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>X</th>\n",
              "      <th>y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0.0</td>\n",
              "      <td>3.392682</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.5</td>\n",
              "      <td>4.488746</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>1.0</td>\n",
              "      <td>4.225207</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>1.5</td>\n",
              "      <td>5.628789</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>2.0</td>\n",
              "      <td>6.485008</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>2.5</td>\n",
              "      <td>7.887852</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>3.0</td>\n",
              "      <td>8.722331</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>3.5</td>\n",
              "      <td>9.081046</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>4.0</td>\n",
              "      <td>11.499504</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>4.5</td>\n",
              "      <td>12.119256</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>10</th>\n",
              "      <td>5.0</td>\n",
              "      <td>13.462674</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>11</th>\n",
              "      <td>5.5</td>\n",
              "      <td>12.681298</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>12</th>\n",
              "      <td>6.0</td>\n",
              "      <td>14.444124</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>13</th>\n",
              "      <td>6.5</td>\n",
              "      <td>16.376249</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>14</th>\n",
              "      <td>7.0</td>\n",
              "      <td>17.208009</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>15</th>\n",
              "      <td>7.5</td>\n",
              "      <td>17.417770</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>16</th>\n",
              "      <td>8.0</td>\n",
              "      <td>19.104009</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>17</th>\n",
              "      <td>8.5</td>\n",
              "      <td>19.741933</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>18</th>\n",
              "      <td>9.0</td>\n",
              "      <td>20.840528</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>19</th>\n",
              "      <td>9.5</td>\n",
              "      <td>22.592334</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>20</th>\n",
              "      <td>10.0</td>\n",
              "      <td>23.862569</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "       X          y\n",
              "0    0.0   3.392682\n",
              "1    0.5   4.488746\n",
              "2    1.0   4.225207\n",
              "3    1.5   5.628789\n",
              "4    2.0   6.485008\n",
              "5    2.5   7.887852\n",
              "6    3.0   8.722331\n",
              "7    3.5   9.081046\n",
              "8    4.0  11.499504\n",
              "9    4.5  12.119256\n",
              "10   5.0  13.462674\n",
              "11   5.5  12.681298\n",
              "12   6.0  14.444124\n",
              "13   6.5  16.376249\n",
              "14   7.0  17.208009\n",
              "15   7.5  17.417770\n",
              "16   8.0  19.104009\n",
              "17   8.5  19.741933\n",
              "18   9.0  20.840528\n",
              "19   9.5  22.592334\n",
              "20  10.0  23.862569"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Seed the random-number generator\n",
        "np.random.seed(random_seed)\n",
        "\n",
        "\n",
        "# Create a linear function to generate experimental data\n",
        "def f(x):\n",
        "    return 2 * x + 3\n",
        "\n",
        "\n",
        "# Add Gaussian noise to the linear ouput\n",
        "def model(X):\n",
        "    return np.random.normal(f(X), err_sig)\n",
        "\n",
        "\n",
        "# Set up experimental data dataframe\n",
        "X = np.linspace(0, 10, 21)\n",
        "y = model(X)\n",
        "df_exp = pd.DataFrame({\"X\": X, \"y\": y})\n",
        "display(df_exp)\n",
        "\n",
        "# Upload the dataset to the cloud\n",
        "dataset = tl.Dataset(\"inverse-methods-dataset\")\n",
        "dataset.upload(df_exp)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 449
        },
        "id": "_NEjPQ3GOzsr",
        "outputId": "e514e4c8-01a8-4f38-bb80-90702551f110"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAy0lEQVR4nO3de1yUZf7/8fcMJzkJaCkHAUnxVB7zkGmrlYbm+sgO+rUvKVbf2oeiieZu7pbpVkbp1+yo2WHVdi2sLa3Wb7pqqYt5StS0TBEpTFErBQQMcOb+/cHPyRFFUGDmnnk9H495PJq5r7nnwx06b6/ruq/LYhiGIQAAABOyuroAAACAy0WQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApuXr6gLqm91u15EjRxQaGiqLxeLqcgAAQA0YhqFTp04pOjpaVuvF+108PsgcOXJEsbGxri4DAABchkOHDqlFixYXPe7xQSY0NFRS5YVo3Lixi6sBAAA1UVRUpNjYWMf3+MV4fJA5O5zUuHFjggwAACZzqWkhTPYFAACmRZABAACmRZABAACm5fFzZGrKZrOpoqLC1WUAtebn5ycfHx9XlwEALuH1QcYwDB09elQFBQWuLgW4bOHh4YqMjGStJABex+uDzNkQ06xZMwUFBfFFAFMxDEOlpaU6fvy4JCkqKsrFFQFAw/LqIGOz2RwhpmnTpq4uB7gsgYGBkqTjx4+rWbNmDDMB8CpePdn37JyYoKAgF1cCXJmzv8PM8wLgbbw6yJzFcBLMjt9hAN7Kq4eWAADA5bHbDWXnn1FBqV3hQVYlRvnKam34f1QRZAAAQK1k5ZQrI7NEJ0sMx2sRwRaN7Busbq38G7QWhpZQ58aMGaNhw4a5uowr0rJlS7344ou1eo8n/NwAcClZOeWav6rYKcRI0skSQ/NXFSsrp7xB6yHI1AG73dC+wxXakl2mfYcrZLcbl37TFRgzZowsFkuVx6BBg+r1c2vqpZde0qJFi1xdhqTKuSPLly93dRkX9P3338tisWjnzp2uLgUAasRuN5SRWVJtm4yNJfX+PXguhpaukKu61wYNGqSFCxc6vRYQEFBvn1cTNptNFotFYWFhLq0DAFA/svPPVOmJOd/J4sq5M21j/BqkJnpkroAru9cCAgIUGRnp9IiIiJAkrVu3Tv7+/vrPf/7jaD9r1iw1a9ZMx44dkyT1799f48eP1/jx4xUWFqarrrpK06ZNk2H89rOUlZVpypQpiomJUXBwsHr16qV169Y5ji9atEjh4eH65JNP1KFDBwUEBCgvL6/KEEv//v01YcIEpaWlKSIiQs2bN9ebb76pkpIS3X///QoNDVXr1q312WefOf2Me/bs0eDBgxUSEqLmzZtr1KhR+vnnn53O+8gjj+hPf/qTmjRposjISM2YMcNxvGXLlpKkO++8UxaLxfE8JydHd9xxh5o3b66QkBD16NFDa9asqdX1t9lsmjx5ssLDw9W0aVP96U9/crp2krRy5Ur17dvX0eb3v/+9cnJyHMcTEhIkSV27dpXFYlH//v0lSdu2bdPAgQN11VVXKSwsTP369VNWVlat6gOA+lBQaq/TdnWBIHOZ3LF77az+/fsrLS1No0aNUmFhoXbs2KFp06bprbfeUvPmzR3tFi9eLF9fX23dulUvvfSSXnjhBb311luO4+PHj9emTZuUkZGhr7/+WsOHD9egQYOUnZ3taFNaWqrnn39eb731lr755hs1a9bsgjUtXrxYV111lbZu3aoJEyZo7NixGj58uG688UZlZWXptttu06hRo1RaWipJKigo0C233KKuXbvqq6++0sqVK3Xs2DGNGDGiynmDg4O1ZcsWzZo1S0899ZRWr14tqTIQSNLChQuVn5/veF5cXKzbb79da9eu1Y4dOzRo0CANHTpUeXl5Nb7Gc+bM0aJFi/S3v/1NmZmZOnHihJYtW+bUpqSkRJMnT9ZXX32ltWvXymq16s4775TdXvkHfOvWrZKkNWvWKD8/Xx999JEk6dSpU0pJSVFmZqY2b96sxMRE3X777Tp16lSN6wOA+hAeVLPYUNN2dcLwcIWFhYYko7CwsMqx06dPG99++61x+vTpWp/3ux/Ljf957ZdLPr77sbwufgwnKSkpho+PjxEcHOz0mDlzpqNNWVmZ0aVLF2PEiBFGhw4djIceesjpHP369TPat29v2O12x2uPPfaY0b59e8MwDOOHH34wfHx8jMOHDzu979ZbbzX+/Oc/G4ZhGAsXLjQkGTt37qxS3x133OH0WX379nU8P3PmjBEcHGyMGjXK8Vp+fr4hydi0aZNhGIbx9NNPG7fddpvTeQ8dOmRIMvbt23fB8xqGYfTo0cN47LHHHM8lGcuWLbvAVXR27bXXGq+88orjeXx8vDF37tyLto+KijJmzZrleF5RUWG0aNHC6ec+308//WRIMnbv3m0YhmHk5uYakowdO3ZUW5vNZjNCQ0ONTz/99KJtruR3GQBqymazG39cdKLa770/Lj5h2Gz2S5/sEqr7/j4Xc2Quk6u7126++WbNnz/f6bUmTZo4/tvf319LlixRp06dFB8fr7lz51Y5xw033OC0kFrv3r01Z84c2Ww27d69WzabTW3atHF6T1lZmdN2Dv7+/urUqdMl6z23jY+Pj5o2baqOHTs6XjvbU3R2z6Bdu3bpiy++UEhISJVz5eTkOOo6/7OjoqIc57iY4uJizZgxQytWrFB+fr7OnDmj06dP17hHprCwUPn5+erVq5fjNV9fX3Xv3t1peCk7O1tPPvmktmzZop9//tnRE5OXl6frrrvuouc/duyYnnjiCa1bt07Hjx+XzWZTaWlprXqMAKA+WK2Vc0Dnryq+aJuRfYIbdD0ZgsxlcnX3WnBwsFq3bl1tmy+//FKSdOLECZ04cULBwcE1Pn9xcbF8fHy0ffv2Knv3nBsuAgMDa7SqrJ+f86Qvi8Xi9NrZc5z9si8uLtbQoUP1/PPPVznXuRsjXui8Z89xMVOmTNHq1av1v//7v2rdurUCAwN1zz33qLy8buc0DR06VPHx8XrzzTcVHR0tu92u66677pKfk5KSol9++UUvvfSS4uPjFRAQoN69e9d5fQBwObq18tfYpJCqN7qEWDSyT8OvI0OQuUyJUb6KCLZUO3s7IsSixCjXXOKcnBxNmjRJb775ppYuXaqUlBStWbNGVutvwWrLli1O7zk7H8PHx0ddu3aVzWbT8ePHddNNNzV0+erWrZs+/PBDtWzZUr6+l38N/fz8ZLPZnF7buHGjxowZozvvvFNSZWj6/vvva3zOsLAwRUVFacuWLfrd734nSTpz5oy2b9+ubt26SZJ++eUX7du3T2+++abj+mVmZjqdx9+/8g/7heqbN2+ebr/9dknSoUOHnCY5A4CrdWvlry4Jfm6xsi+TfS/T2e616tRn91pZWZmOHj3q9Dj7ZWez2XTfffcpKSlJ999/vxYuXKivv/5ac+bMcTpHXl6eJk+erH379um9997TK6+8ookTJ0qS2rRpo+TkZI0ePVofffSRcnNztXXrVqWnp2vFihX18jOdKzU1VSdOnNC9996rbdu2KScnR6tWrdL9999f5Yu/Oi1bttTatWt19OhRnTx5UpKUmJiojz76SDt37tSuXbv03//935fsxTnfxIkT9dxzz2n58uX67rvvNG7cOBUUFDiOR0REqGnTpnrjjTd04MABff7555o8ebLTOZo1a6bAwEDHRObCwkJHfX//+9+1d+9ebdmyRcnJyY4drgHAXVitFrWN8VOvxAC1jfFzSYiRCDJX5Gz3WkSw8/+8iBCLxiaF1Gv32sqVKxUVFeX06Nu3ryRp5syZ+uGHH7RgwQJJlUMxb7zxhp544gnt2rXLcY7Ro0fr9OnT6tmzp1JTUzVx4kQ9/PDDjuMLFy7U6NGj9eijj6pt27YaNmyYtm3bpri4uHr7uc6Kjo7Wxo0bZbPZdNttt6ljx45KS0tTeHi4U6/SpcyZM0erV69WbGysunbtKkl64YUXFBERoRtvvFFDhw5VUlKSoyelph599FGNGjVKKSkp6t27t0JDQx09PJJktVqVkZGh7du367rrrtOkSZM0e/Zsp3P4+vrq5Zdf1oIFCxQdHa077rhDkvT222/r5MmT6tatm0aNGqVHHnnkoneDAYC3sxiG0fD3BzegoqIihYWFqbCwUI0bN3Y69uuvvyo3N1cJCQlq1KjRZX+Gu2ycVRv9+/dXly5dar0MP9xTXf0uA4C7qO77+1zMkakDZ7vXAABAw2JoCQAAmBY9Ml7q3K0GAAAwK3pkAACAaRFkAACAaRFkAACAaRFkAACAaTHZFwAAD2bGtc5qgyDj4Vq2bKm0tDSlpaW5upQ64Wk/DwDUp6yc8qqbOwZXbrHT0Js71heGlkzq0KFDeuCBBxQdHS1/f3/Fx8dr4sSJ+uWXX1xdmsudOHFCaWlpio+Pl7+/v6Kjo/XAAw8oLy+v1ueyWCxavnx53RepylDGysoA6ktWTrnmryqusrnxyRJD81cVKyun3EWV1S2CjAkdPHhQ3bt3V3Z2tt577z0dOHBAr7/+utauXavevXvrxIkTLqvNZrPVegPGunTixAndcMMNWrNmjV5//XUdOHBAGRkZOnDggHr06KGDBw+6rDYAaCh2u6GMzJJq22RsLJHdbv5diggyJpSamip/f3/9+9//Vr9+/RQXF6fBgwdrzZo1Onz4sB5//HGn9qdOndK9996r4OBgxcTE6LXXXnMcMwxDM2bMUFxcnAICAhQdHa1HHnnEcbysrExTpkxRTEyMgoOD1atXL6fF9BYtWqTw8HB98skn6tChgwICAvTWW2+pUaNGTrtBS5U7Rt9yyy2O55mZmbrpppsUGBio2NhYPfLIIyop+e0P3vHjxzV06FAFBgYqISFBS5YsueS1efzxx3XkyBGtWbNGgwcPVlxcnH73u99p1apV8vPzU2pqqqPthXpEunTpohkzZjiOS9Kdd94pi8XieD5jxgx16dJFCxYsUGxsrIKCgjRixAjH7tVS5V5W5w9/DRs2TGPGjHEc/+GHHzRp0iRZLBZZLJ4zXg3A9bLzz1TpiTnfyeLKuTNmR5A5l2FIZb+65lHDvTtPnDihVatWady4cQoMDHQ6FhkZqeTkZC1dulTn7gU6e/Zsde7cWTt27NDUqVM1ceJErV69WpL04Ycfau7cuVqwYIGys7O1fPlydezY0fHe8ePHa9OmTcrIyNDXX3+t4cOHa9CgQcrOzna0KS0t1fPPP6+33npL33zzjZKTkxUeHq4PP/zQ0cZms2np0qVKTk6WJOXk5GjQoEG6++679fXXX2vp0qXKzMzU+PHjHe8ZM2aMDh06pC+++EL//Oc/NW/ePB0/fvyi18ZutysjI0PJycmKjIx0OhYYGKhx48Zp1apVNe6x2rZtm6TKXcDz8/MdzyXpwIEDev/99/Xpp59q5cqV2rFjh8aNG1ej80rSRx99pBYtWuipp55Sfn6+8vPza/xeALiUgtKa9YzXtJ07Y7LvucrLpNRhrvns15ZLAZfetTg7O1uGYah9+/YXPN6+fXudPHlSP/30k5o1ayZJ6tOnj6ZOnSpJatOmjTZu3Ki5c+dq4MCBysvLU2RkpAYMGCA/Pz/FxcWpZ8+ekqS8vDwtXLhQeXl5io6OliRNmTJFK1eu1MKFC/Xss89KkioqKjRv3jx17tzZUcfIkSP17rvv6sEHH5QkrV27VgUFBbr77rslSenp6UpOTnb0WiQmJurll19Wv379NH/+fOXl5emzzz7T1q1b1aNHD0nS22+/fdGfW5J++uknFRQUVHttDMPQgQMHHD9jda6++mpJUnh4eJVg9Ouvv+qdd95RTEyMJOmVV17RkCFDNGfOnCptL6RJkyby8fFRaGhojdoDQG2EB9Wsn6Km7dyZ+X8CL2XUsAdHknr37l3l+d69eyVJw4cP1+nTp3XNNdfooYce0rJly3TmTGVX4+7du2Wz2dSmTRuFhIQ4HuvXr1dOTo7jfP7+/urUqZPTZyQnJ2vdunU6cuSIJGnJkiUaMmSIwsPDJUm7du3SokWLnM6blJQku92u3Nxc7d27V76+vrr++usd52zXrp3j/XV1bS5XXFycI8RIldfUbrdr37599f7ZAHApiVG+igiufsg6IsSixCjz92eY/yeoS/4BlT0jrvrsGmjdurUsFov27t2rO++8s8rxvXv3KiIiwtGbcCmxsbHat2+f1qxZo9WrV2vcuHGaPXu21q9fr+LiYvn4+Gj79u3y8fFxel9ISIjjvwMDA6vM8ejRo4datWqljIwMjR07VsuWLdOiRYscx4uLi/WHP/zBaT7OWXFxcdq/f3+N6j/X1VdfrfDwcEdIO9/evXtlsVjUunVrSZLVaq0SeioqKmr9uRdSn+cGgEuxWitvsZ6/qviibUb2CfaI9WQIMueyWGo0vONKTZs21cCBAzVv3jxNmjTJaZ7M0aNHtWTJEo0ePdopWGzevNnpHJs3b3YafgkMDNTQoUM1dOhQpaamql27dtq9e7e6du0qm82m48eP66abbqp1rcnJyVqyZIlatGghq9WqIUOGOI5169ZN3377rSNUnK9du3Y6c+aMtm/f7hha2rdvX5UJxOeyWq0aMWKElixZoqeeesppyOb06dOaN2+ekpKS1KRJE0mVwefcuSlFRUXKzc11Oqefn59sNluVz8rLy9ORI0ccQ26bN2+W1WpV27ZtL3hum82mPXv26Oabb3a85u/vf8FzA0Bd6NbKX2OTQqquIxNi0cg+rCMDF3r11VdVVlampKQkbdiwQYcOHdLKlSs1cOBAxcTEaObMmU7tN27cqFmzZmn//v167bXX9MEHH2jixImSKu86evvtt7Vnzx4dPHhQ//jHPxQYGKj4+Hi1adNGycnJGj16tD766CPl5uZq69atSk9P14oVKy5ZZ3JysrKysjRz5kzdc889Cgj4rdfpscce05dffqnx48dr586dys7O1scff+yY7Nu2bVsNGjRIf/jDH7RlyxZt375d//M//1NlgvP5nn32WUVGRmrgwIH67LPPdOjQIW3YsEFJSUmqqKhwumPrlltu0d///nf95z//0e7du5WSklKl56lly5Zau3atjh49qpMnTzpeb9SokVJSUrRr1y795z//0SOPPKIRI0Y4wtMtt9yiFStWaMWKFfruu+80duzYKiGsZcuW2rBhgw4fPqyff/75ktcTAGqrWyt/PTcqXFPuCNX/DAzWlDtC9dx94R4TYiRJhocrLCw0JBmFhYVVjp0+fdr49ttvjdOnT7ugsivz/fffGykpKUbz5s0NPz8/IzY21pgwYYLx888/O7WLj483/vrXvxrDhw83goKCjMjISOOll15yHF+2bJnRq1cvo3HjxkZwcLBxww03GGvWrHEcLy8vN5588kmjZcuWhp+fnxEVFWXceeedxtdff20YhmEsXLjQCAsLu2idPXv2NCQZn3/+eZVjW7duNQYOHGiEhIQYwcHBRqdOnYyZM2c6jufn5xtDhgwxAgICjLi4OOOdd94x4uPjjblz51Z7bX766SdjwoQJRmxsrOHn52c0b97cGDNmjPHDDz84tSssLDT+67/+y2jcuLERGxtrLFq0yOjcubMxffp0R5tPPvnEaN26teHr62vEx8cbhmEY06dPNzp37mzMmzfPiI6ONho1amTcc889xokTJ5yu29ixY40mTZoYzZo1M9LT04077rjDSElJcbTZtGmT0alTJyMgIMC40j+KZv5dBoALqe77+1wWw2iAmZEuVFRUpLCwMBUWFqpx48ZOx3799Vfl5uYqISFBjRq595AS3MeMGTO0fPly7dy509WlOPC7DMDTVPf9fS6GlgAAgGkRZAAAgGkRZIBamjFjhlsNKwGANyPIAAAA0yLIAAAA0yLIAAAA0yLIAAAA0yLIAAAA0yLIAAAA02LTyAsp+1U604A7Ffv61dtmlS1btlRaWprS0tLq5fwNbd26dbr55pt18uRJhYeHu7ocAICLEWTOV/ar9NUGqeRUw31mcKjU/Xe1CjOHDh3S9OnTtXLlSv3888+KiorSsGHD9OSTT6pp06b1WGzD6d+/v7p06aIXX3zR8dqNN96o/Px8hYWFua4wAIDbIMic70xFZYjx85f8G2B30PLyys87U1HjIHPw4EH17t1bbdq00XvvvaeEhAR98803+uMf/6jPPvtMmzdvVpMmTeq58Auz2WyyWCyyWutn1NLf39+xw7Q7MQxDNptNvr7Of6TKy8vlfxm/R5f7PgDwNsyRuRh/fykgsP4fl/FllZqaKn9/f/373/9Wv379FBcXp8GDB2vNmjU6fPiwHn/8caf2p06d0r333qvg4GDFxMTotddecxwzDEMzZsxQXFycAgICFB0drUceecRxvKysTFOmTFFMTIyCg4PVq1cvrVu3znF80aJFCg8P1yeffKIOHTooICBAb731lho1aqSCggKnOiZOnKhbbrlFkvTLL7/o3nvvVUxMjIKCgtSxY0e99957jrZjxozR+vXr9dJLL8lischisej777/XunXrZLFYnM794Ycf6tprr1VAQIBatmypOXPmOH1uy5Yt9eyzz+qBBx5QaGio4uLi9MYbb1R7je12u9LT05WQkKDAwEB17txZ//znPx3Hz9bx2Wef6frrr1dAQIAyMzPVv39/jR8/XmlpabrqqquUlJQkSVq/fr169uypgIAARUVFaerUqTpz5ozjfBd7HwCgegQZkzlx4oRWrVqlcePGKTAw0OlYZGSkkpOTtXTpUp27qfns2bPVuXNn7dixQ1OnTtXEiRO1evVqSZUhYO7cuVqwYIGys7O1fPlydezY0fHe8ePHa9OmTcrIyNDXX3+t4cOHa9CgQcrOzna0KS0t1fPPP6+33npL33zzjZKTkxUeHq4PP/zQ0cZms2np0qVKTk6WVLlb8/XXX68VK1Zoz549evjhhzVq1Cht3bpVkvTSSy+pd+/eeuihh5Sfn6/8/HzFxsZWuR7bt2/XiBEjNHLkSO3evVszZszQtGnTtGjRIqd2c+bMUffu3bVjxw6NGzdOY8eO1b59+y56ndPT0/XOO+/o9ddf1zfffKNJkybpvvvu0/r1653aTZ06Vc8995z27t2rTp06SZIWL14sf39/bdy4Ua+//roOHz6s22+/XT169NCuXbs0f/58vf3223rmmWecznX++wAAl8bQkslkZ2fLMAy1b9/+gsfbt2+vkydP6qefflKzZs0kSX369NHUqVMlSW3atNHGjRs1d+5cDRw4UHl5eYqMjNSAAQPk5+enuLg49ezZU5KUl5enhQsXKi8vT9HR0ZKkKVOmaOXKlVq4cKGeffZZSVJFRYXmzZunzp07O+oYOXKk3n33XT344IOSpLVr16qgoEB33323JCkmJkZTpkxxtJ8wYYJWrVql999/Xz179lRYWJj8/f0VFBRU7VDSCy+8oFtvvVXTpk1z/HzffvutZs+erTFjxjja3X777Ro3bpwk6bHHHtPcuXP1xRdfqG3btlXOWVZWpmeffVZr1qxR7969JUnXXHONMjMztWDBAvXr18/R9qmnntLAgQOd3p+YmKhZs2Y5nj/++OOKjY3Vq6++KovFonbt2unIkSN67LHH9OSTTzqG4c5/HwDg0uiRMalze1wu5eyX8bnP9+7dK0kaPny4Tp8+rWuuuUYPPfSQli1b5hjy2L17t2w2m9q0aaOQkBDHY/369crJyXGcz9/f39EbcVZycrLWrVunI0eOSJKWLFmiIUOGOO40stlsevrpp9WxY0c1adJEISEhWrVqlfLy8mp1Hfbu3as+ffo4vdanTx9lZ2fLZrM5Xju3PovFosjISB0/fvyC5zxw4IBKS0s1cOBAp5/7nXfecfq5Jal79+5V3n/99ddXqbF3796yWCxONRYXF+vHH3+86PsAAJdGj4zJtG7dWhaLRXv37tWdd95Z5fjevXsVERGhq6++ukbni42N1b59+7RmzRqtXr1a48aN0+zZs7V+/XoVFxfLx8dH27dvl4+Pj9P7QkJCHP8dGBjo9CUtST169FCrVq2UkZGhsWPHatmyZU7DPbNnz9ZLL72kF198UR07dlRwcLDS0tJUXl5ei6tRc35+fk7PLRaL7Hb7BdsWFxdLklasWKGYmBinYwEBAU7Pg4ODq7z/Qq/VxOW+D4D52e2GsvPPqKDUrvAgqxKjfGW1Wi79RhBkzKZp06YaOHCg5s2bp0mTJjnNkzl69KiWLFmi0aNHOwWLzZs3O51j8+bNTkNTgYGBGjp0qIYOHarU1FS1a9dOu3fvVteuXWWz2XT8+HHddNNNta41OTlZS5YsUYsWLWS1WjVkyBDHsY0bN+qOO+7QfffdJ6lycu3+/fvVoUMHRxt/f3+nXpULad++vTZu3Oj02saNG9WmTZsq4aumzk5azsvLcxpGulzt27fXhx9+KMMwHP9fNm7cqNDQULVo0eKKzw/A3LJyypWRWaKTJb/1tEcEWzSyb7C6teLuxUthaMmEXn31VZWVlSkpKUkbNmzQoUOHtHLlSg0cOFAxMTGaOXOmU/uNGzdq1qxZ2r9/v1577TV98MEHmjhxoqTKu47efvtt7dmzRwcPHtQ//vEPBQYGKj4+Xm3atFFycrJGjx6tjz76SLm5udq6davS09O1YsWKS9aZnJysrKwszZw5U/fcc49Tb0ZiYqJWr16tL7/8Unv37tUf/vAHHTt2zOn9LVu21JYtW/T999/r559/vmAPyqOPPqq1a9fq6aef1v79+7V48WK9+uqrTvNvais0NFRTpkzRpEmTtHjxYuXk5CgrK0uvvPKKFi9eXOvzjRs3TocOHdKECRP03Xff6eOPP9b06dM1efLkertNHYA5ZOWUa/6qYqcQI0knSwzNX1WsrJz66aX2JPTIXEw9DXHUxeckJibqq6++0vTp0zVixAidOHFCkZGRGjZsmKZPn15lDZlHH31UX331lf7617+qcePGeuGFFxy394aHh+u5557T5MmTZbPZ1LFjR3366aeORfUWLlyoZ555Ro8++qgOHz6sq666SjfccIN+//vfX7LO1q1bq2fPntq6davTonaS9MQTT+jgwYNKSkpSUFCQHn74YQ0bNkyFhYWONlOmTFFKSoo6dOig06dPKzc3t8pndOvWTe+//76efPJJPf3004qKitJTTz3lNNH3cjz99NO6+uqrlZ6eroMHDyo8PFzdunXTX/7yl1qfKyYmRv/3f/+nP/7xj+rcubOaNGmiBx98UE888cQV1QjA3Ox2QxmZJdW2ydhYoi4JfgwzVcNi1GbWqAkVFRUpLCxMhYWFaty4sdOxX3/9Vbm5uUpISFCjRv9/MTqTrOwLnOuCv8sA3Nq+wxX6348v/V0z5Y5QtY3xu2Q7T1Pd9/e56JE5X0CjylDhIXstAQDcU0HphW84uNx23oogcyEBjQgWAIB6FR5UszlyNW3nrbg6AAC4QGKUryKCq5/7EhFiUWIUfQ7VIcgAAOACVmvlLdbVGdknmIm+l0CQUe1WyQXcEb/DgDl1a+WvsUkhVXpmIkIsGpsUwjoyNeDV/VVnV3stLS2tsgEjYCalpaWSqq5gDMD9dWvlry4Jfqzse5m8Osj4+PgoPDzcsedOUFBQlaX2AXdmGIZKS0t1/PhxhYeHX/ZqxgBcy2q1eOUt1nXBq4OMJMfOyhfbQBAwg/Dw8Gp3CQcAT+X1QcZisSgqKkrNmjVTRUUDrh0D1BE/Pz96YgB4La8PMmf5+PjwZQAAgMm49K6l9PR09ejRQ6GhoWrWrJmGDRumffv2ObX59ddflZqaqqZNmyokJER33313lc0FAQCAd3JpkFm/fr1SU1O1efNmrV69WhUVFbrttttUUvLbJlqTJk3Sp59+qg8++EDr16/XkSNHdNddd7mwagAA4C7catPIn376Sc2aNdP69ev1u9/9ToWFhbr66qv17rvv6p577pEkfffdd2rfvr02bdqkG264oco5ysrKVFZW5nheVFSk2NjYS246BQAA3EdNN410qwXxCgsLJUlNmjSRJG3fvl0VFRUaMGCAo027du0UFxenTZs2XfAc6enpCgsLczxiY2Prv3AAgFew2w3tO1yhLdll2ne4Qna72/QFeC23mexrt9uVlpamPn366LrrrpMkHT16VP7+/goPD3dq27x5cx09evSC5/nzn/+syZMnO56f7ZEBAOBKZOWUKyOzRCdLfgsvEcGV2wywAq/ruE2QSU1N1Z49e5SZmXlF5wkICFBAQEAdVQUAQGWImb+quMrrJ0sMzV9VzHYCLuQWQ0vjx4/Xv/71L33xxRdq0aKF4/XIyEiVl5eroKDAqf2xY8dY/AsA0CDsdkMZmSXVtsnYWMIwk4u4NMgYhqHx48dr2bJl+vzzz5WQkOB0/Prrr5efn5/Wrl3reG3fvn3Ky8tT7969G7pcAIAXys4/4zScdCEniw1l559poIpwLpcOLaWmpurdd9/Vxx9/rNDQUMe8l7CwMAUGBiosLEwPPvigJk+erCZNmqhx48aaMGGCevfufcE7lgAAqGsFpfY6bYe65dIgM3/+fElS//79nV5fuHChxowZI0maO3eurFar7r77bpWVlSkpKUnz5s1r4EoBAN4qPKhmgxc1bYe65VbryNSHmt6HDgDAhdjthqb+vaDa4aWIEIueuy9cVqulASvzbKZcRwYAAHdjtVbeYl2dkX2CCTEuQpABAOASurXy19ikEEUEO4eViBALt167mNusIwMAQF2w2yvvICootSs8yKrEKN866S3p1spfXRL86uXcuHwEGQCAx6jv1XetVovaxvhd8XlQdxhaAgB4hLOr754/Kffs6rtZOeUuqgz1iSADADA9Vt/1XgQZAIDpsfqu9yLIAABMj9V3vRdBBgBgeqy+6734PwoAML3EKN8qa7ycLyLEosQobtb1NAQZAIDpsfqu9yLIAAA8Aqvveif62AAAHoPVd70PQQYA4FFYfde7MLQEAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMiyADAABMy9fVBQAAvJPdbig7/4wKSu0KD7IqMcpXVqvF1WXBZAgyAIAGl5VTrozMEp0sMRyvRQRbNLJvsLq18ndhZTAbhpYAAA0qK6dc81cVO4UYSTpZYmj+qmJl5ZS7qDKYEUEGAHBRdruhfYcrtCW7TPsOV8huNy79pkucLyOzpNo2GRtLrvhz4D0YWgIAXFB9DP9k55+p0hNzvpPFlXNn2sb4XdZnwLvQIwMAqKK+hn8KSu112g4gyAAAnNTn8E94UM2+dmraDuA3BQDgpDbDP7WVGOWriODqb7GOCLEoMYqZD6gZggwAwEl9Dv9YrZVzbKozsk8w68mgxggyAAAn9T38062Vv8YmhVTpmYkIsWhsUgjryKBW6LsDADg5O/xT3fDSlQ7/dGvlry4JfqzsiytGjwwAwElDDf9YrRa1jfFTr8QAtY3xI8TgshBkAABVMPwDs2BoCQBwQQz/wAwIMgCAizo7/AO4K4aWAACAaRFkAACAaTG0BAAewG43mMsCr0SQAQCTq49dqgGzYGgJAEysvnapBsyCIAMAJlWfu1QDZkGQAQCTqs9dqgGzIMgAgEnV5y7VgFkQZADApOp7l2rADPjtBgCTOrtLdXWudJdqwN0RZADApBpql2rAnRFkAMDE2KUa3o7+RgAwOXaphjcjyACAB2CXangrhpYAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpuTTIbNiwQUOHDlV0dLQsFouWL1/udHzMmDGyWCxOj0GDBrmmWAAA4HZcGmRKSkrUuXNnvfbaaxdtM2jQIOXn5zse7733XgNWCAAA3JmvKz988ODBGjx4cLVtAgICFBkZWeNzlpWVqayszPG8qKjosusDAADuze3nyKxbt07NmjVT27ZtNXbsWP3yyy/Vtk9PT1dYWJjjERsb20CVAgCAhmYxDMNwdRGSZLFYtGzZMg0bNszxWkZGhoKCgpSQkKCcnBz95S9/UUhIiDZt2iQfH58LnudCPTKxsbEqLCxU48aN6/vHAAAAdaCoqEhhYWGX/P526dDSpYwcOdLx3x07dlSnTp3UqlUrrVu3TrfeeusF3xMQEKCAgICGKhEAALiQ2w8tneuaa67RVVddpQMHDri6FAAA4AZMFWR+/PFH/fLLL4qKinJ1KQAAwA24dGipuLjYqXclNzdXO3fuVJMmTdSkSRP99a9/1d13363IyEjl5OToT3/6k1q3bq2kpCQXVg0Al8duN5Sdf0YFpXaFB1mVGOUrq9Xi6rIAU3NpkPnqq6908803O55PnjxZkpSSkqL58+fr66+/1uLFi1VQUKDo6Gjddtttevrpp5kDA8B0snLKlZFZopMlv91fERFs0ci+werWyt+FlQHm5jZ3LdWXms56BoD6kpVTrvmrii96fGxSCGEGOE9Nv79NNUcGAMzGbjeUkVlSbZuMjSWy2z3635RAvSHIAEA9ys4/4zScdCEniyvnzgCoPYIMANSjglJ7nbYD4IwgAwD1KDyoZn/N1rQdAGf8yQGAepQY5auI4OpvsY4IsSgxyq0XWgfcFkEGAOqR1Vp5i3V1RvYJZj0Z4DIRZACgnnVr5a+xSSFVemYiQizceg1cIfoyAaABdGvlry4JfqzsC9QxggwANBCr1aK2MX6uLgPwKAwtAQAA0yLIAAAA0yLIAAAA0yLIAAAA06p1kElJSdGGDRvqoxYAqDG73dC+wxXakl2mfYcr2HQR8FK1vmupsLBQAwYMUHx8vO6//36lpKQoJiamPmoDgAvKyilXRmaJ02aMEcGVC8+xJgvgXWrdI7N8+XIdPnxYY8eO1dKlS9WyZUsNHjxY//znP1VRUVEfNQKAQ1ZOueavKq6yo/TJEkPzVxUrK6fcRZUBcIXLmiNz9dVXa/Lkydq1a5e2bNmi1q1ba9SoUYqOjtakSZOUnZ1d13UCgOx2QxmZJdW2ydhYwjAT4EWuaLJvfn6+Vq9erdWrV8vHx0e33367du/erQ4dOmju3Ll1VSMASJKy889U6Yk538liQ9n5ZxqoIgCuVusgU1FRoQ8//FC///3vFR8frw8++EBpaWk6cuSIFi9erDVr1uj999/XU089VR/1AvBiBaX2Om0HwPxqPdk3KipKdrtd9957r7Zu3aouXbpUaXPzzTcrPDy8DsoDgN+EB9Xs3141bQfA/GodZObOnavhw4erUaNGF20THh6u3NzcKyoMAM6XGOWriGBLtcNLESEWJUaxjRzgLWr9z5ZRo0ZVG2IAoL5YrZW3WFdnZJ9gdpQGvAj9rwBMpVsrf41NClFEsHNYiQixaGxSCOvIAF6G/lcAptOtlb+6JPgpO/+MCkrtCg+yKjHKl54YwAsRZACYktVqUdsYP1eXAcDFCDIAcA673aCnBzARggwA/H/s4QSYD5N9AUDs4QSYFUEGgNdjDyfAvAgyALweezgB5kWQAeD12MMJMC+CDACvxx5OgHnxpxKA1zu7h1N12MMJcE8EGQBejz2cAPMiyACA2MMJMCv6SQHUG7OtksseToD5EGQA1AuzrpLLHk6AuTC0BKDOsUougIZCkAFQp1glF0BDIsgAqFOskgugIRFkANQpVskF0JAIMgDqFKvkAmhI/E0CoE6xSi6AhkSQAVCnWCUXQEMiyACoc6ySC6Ch0LcLoF5W4GWVXAANgSADeLn6XIGXVXIB1DeGlgAvxgq8AMyOIAN4KVbgBeAJCDKAl2IFXgCegCADeClW4AXgCQgygJdiBV4AnoC/oQAvxQq8ADwBQQbwUqzAC8ATEGQAL8YKvADMjj5jwMuxAi8AMyPIAGAFXgCmxdASAAAwLe/pkSk5JfnQVQ4AgCmUnKpRM+8JMps/l4KDXF0FAACoiZLSGjXzniDj5y8Fh7i6CgAAUBPlNdsexXuCjL+fFBDo6ioAAEBN+Fe/qe1ZTPYFAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACm5T23X1eckSrKXV0FAACoiQrWkXG29A3J33t+XAAATK2GC+IxtAQAAEzLpV0UGzZs0OzZs7V9+3bl5+dr2bJlGjZsmOO4YRiaPn263nzzTRUUFKhPnz6aP3++EhMTa/9h//Ww1KRJ3RUPAADqz4kT0t/XXrKZS4NMSUmJOnfurAceeEB33XVXleOzZs3Syy+/rMWLFyshIUHTpk1TUlKSvv32WzVq1Kh2H+bnW7nfEgAAcH9+NYsoLg0ygwcP1uDBgy94zDAMvfjii3riiSd0xx13SJLeeecdNW/eXMuXL9fIkSMv+L6ysjKVlZU5nhcVFdV94QAAwC247RyZ3NxcHT16VAMGDHC8FhYWpl69emnTpk0XfV96errCwsIcj9jY2IYoFwAAuIDbBpmjR49Kkpo3b+70evPmzR3HLuTPf/6zCgsLHY9Dhw7Va50AAMB1PO5+5ICAAAUEBLi6DAAA0ADctkcmMjJSknTs2DGn148dO+Y4BgAAvJvbBpmEhARFRkZq7drfbr0qKirSli1b1Lt3bxdWBgAA3IVLh5aKi4t14MABx/Pc3Fzt3LlTTZo0UVxcnNLS0vTMM88oMTHRcft1dHS001ozAADAe7k0yHz11Ve6+eabHc8nT54sSUpJSdGiRYv0pz/9SSUlJXr44YdVUFCgvn37auXKlbVfQ0aSyiukstN1VToAAKhP5RU1amYxDMOo51JcqqioSGFhYSr86B01Dg5ydTkAAKAGikpKFXbXaBUWFqpx48YXbedxdy1d1A23SNVcCAAA4EZquKCt9wSZ4NDKBwAAcH+2mg0Yue1dSwAAAJfiPT0ygMnZ7Yay88+ooNSu8CCrEqN8ZbVaXF0WALgUQQYwgayccmVkluhkyW9drRHBFo3sG6xurdjVHYD3YmgJcHNZOeWav6rYKcRI0skSQ/NXFSsrp9xFlQGA6xFkADdmtxvKyCyptk3GxhLZ7R69igIAXBRBBnBj2flnqvTEnO9kceXcGQDwRsyRAepYXU7KLSi112k7APA0BBmgDtX1pNzwoJp1mta0HQB4Gv72A+pIfUzKTYzyVURw9b05ESEWJUbxbxIA3okgA9SB+pqUa7VW9uZUZ2SfYNaTAeC1CDJAHajPSbndWvlrbFJIlZ6ZiBCLxiaFsI4MAK9GfzRQB+p7Um63Vv7qkuDHyr4AcB6CDFAHGmJSrtVqUdsYv8t+PwB4IoaWgDrApFwAcA2CDFAHmJQLAK5BkAHqCJNyAaDh0c8N1CEm5QJAwyLIAHWMSbkA0HAYWgIAAKZFkAEAAKbF0BK8Ul3uUA0AcB2CDLxOXe9QDQBwHYaW4FXqY4dqAIDrEGTgNeprh2oAgOsQZOA16nOHagCAaxBk4DXqe4dqAEDDI8jAazTEDtUAgIbF39jwGuxQDQCehyADr8EO1QDgeQgy8CrsUA0AnoU+dHgddqgGAM9BkIFXYodqAPAMDC0BAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTIsgAAADTYq8luDW73WBzRwDARRFk4LaycsqVkVmikyWG47WIYItG9g1Wt1b+LqwMAOAuGFqCW8rKKdf8VcVOIUaSTpYYmr+qWFk55S6qDADgTggycDt2u6GMzJJq22RsLJHdblTbBgDg+QgycDvZ+Weq9MSc72Rx5dwZAIB3I8jA7RSU2uu0HQDAcxFk4HbCg2r2a1nTdgAAz8U3AdxOYpSvIoKrv8U6IsSixChuugMAb0eQgduxWitvsa7OyD7BrCcDACDIwD11a+WvsUkhVXpmIkIsGpsUwjoyAABJLIgHN9atlb+6JPixsi8A4KIIMnBrVqtFbWP8XF0GAMBNMbQEAABMiyADAABMi6ElXDF2qAYAuApBBleEHaoBAK7E0BIuGztUAwBcjSCDy8IO1QAAd0CQwWVhh2oAgDsgyOCysEM1AMAdEGRwWdihGgDgDviWwWVhh2oAgDsgyOCysEM1AMAdEGRw2dihGgDgavT744qwQzUAwJUIMrhi7FANAHAVtx5amjFjhiwWi9OjXbt2ri4LAAC4Cbfvkbn22mu1Zs0ax3NfX7cvGQAANBC3TwW+vr6KjIx0dRkAAMANufXQkiRlZ2crOjpa11xzjZKTk5WXl1dt+7KyMhUVFTk9AACAZ3LrINOrVy8tWrRIK1eu1Pz585Wbm6ubbrpJp06duuh70tPTFRYW5njExsY2YMUAAKAhWQzDMM32xAUFBYqPj9cLL7ygBx988IJtysrKVFZW5nheVFSk2NhYFRYWqnHjxg1VKgAAuAJFRUUKCwu75Pe328+ROVd4eLjatGmjAwcOXLRNQECAAgICGrAqAADgKm49tHS+4uJi5eTkKCoqytWl1Bu73dC+wxXakl2mfYcrZLebpsMMAIAG59Y9MlOmTNHQoUMVHx+vI0eOaPr06fLx8dG9997r6tLqRVZOuTIyS3Sy5LfwEhFcuacRy/0DAFCVW/fI/Pjjj7r33nvVtm1bjRgxQk2bNtXmzZt19dVXu7q0OpeVU675q4qdQowknSwxNH9VsbJyyl1UGQAA7sute2QyMjJcXUKDsNsNZWSWVNsmY2OJuiT4sYcRAADncOseGW+RnX+mSk/M+U4WG8rOP3NFn8P8GwCAp3HrHhlvUVBqr9N2F8L8GwCAJ6JHxg2EB9Xsf0NN252P+TcAAE9FkHEDiVG+igiufu5LRIhFiVG170Cr6fwbhpkAAGZEkHEDVmvlEE91RvYJvqyJvg01/wYAAFcgyLiJbq38NTYppErPTESIRWOTQi57HktDzL8BAMBVmOzrRrq18leXBD9l559RQald4UFWJUb5XtEt1/U9/wYAAFciyLgZq9WitjF+dXa+s/Nvqhteutz5NwAAuBr/DPdw9Tn/BgAAVyPIeIH6mn8DAICrMZ7gJepj/g0AAK5GkPEidT3/BgAAV2NoCQAAmBZBBgAAmBZDS5fBbjeYawIAgBsgyNQSu0gDAOA+GFqqBXaRBgDAvRBkaohdpAEAcD8EmRpiF2kAANwPQaaG2EUaAAD3Q5CpIXaRBgDA/fCtW0Nnd5GuDrtIAwDQsAgyNcQu0gAAuB+CTC2wizQAAO6FcZBaYhdpAADcB0HmMrCLNAAA7oGhJQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoev7KvYRiSpKKiIhdXAgAAaurs9/bZ7/GL8fggc+rUKUlSbGysiysBAAC1derUKYWFhV30uMW4VNQxObvdriNHjig0NFQWS91t7FhUVKTY2FgdOnRIjRs3rrPzoiqudcPgOjcMrnPD4Do3jPq8zoZh6NSpU4qOjpbVevGZMB7fI2O1WtWiRYt6O3/jxo35Q9JAuNYNg+vcMLjODYPr3DDq6zpX1xNzFpN9AQCAaRFkAACAaRFkLlNAQICmT5+ugIAAV5fi8bjWDYPr3DC4zg2D69ww3OE6e/xkXwAA4LnokQEAAKZFkAEAAKZFkAEAAKZFkAEAAKZFkLlMr732mlq2bKlGjRqpV69e2rp1q6tL8ijp6enq0aOHQkND1axZMw0bNkz79u1zdVke77nnnpPFYlFaWpqrS/FIhw8f1n333aemTZsqMDBQHTt21FdffeXqsjyKzWbTtGnTlJCQoMDAQLVq1UpPP/30JffrQfU2bNigoUOHKjo6WhaLRcuXL3c6bhiGnnzySUVFRSkwMFADBgxQdnZ2g9RGkLkMS5cu1eTJkzV9+nRlZWWpc+fOSkpK0vHjx11dmsdYv369UlNTtXnzZq1evVoVFRW67bbbVFJS4urSPNa2bdu0YMECderUydWleKSTJ0+qT58+8vPz02effaZvv/1Wc+bMUUREhKtL8yjPP/+85s+fr1dffVV79+7V888/r1mzZumVV15xdWmmVlJSos6dO+u111674PFZs2bp5Zdf1uuvv64tW7YoODhYSUlJ+vXXX+u/OAO11rNnTyM1NdXx3GazGdHR0UZ6eroLq/Jsx48fNyQZ69evd3UpHunUqVNGYmKisXr1aqNfv37GxIkTXV2Sx3nssceMvn37uroMjzdkyBDjgQcecHrtrrvuMpKTk11UkeeRZCxbtszx3G63G5GRkcbs2bMdrxUUFBgBAQHGe++9V+/10CNTS+Xl5dq+fbsGDBjgeM1qtWrAgAHatGmTCyvzbIWFhZKkJk2auLgSz5SamqohQ4Y4/V6jbn3yySfq3r27hg8frmbNmqlr16568803XV2Wx7nxxhu1du1a7d+/X5K0a9cuZWZmavDgwS6uzHPl5ubq6NGjTn9/hIWFqVevXg3yvejxm0bWtZ9//lk2m03Nmzd3er158+b67rvvXFSVZ7Pb7UpLS1OfPn103XXXubocj5ORkaGsrCxt27bN1aV4tIMHD2r+/PmaPHmy/vKXv2jbtm165JFH5O/vr5SUFFeX5zGmTp2qoqIitWvXTj4+PrLZbJo5c6aSk5NdXZrHOnr0qCRd8Hvx7LH6RJCB20tNTdWePXuUmZnp6lI8zqFDhzRx4kStXr1ajRo1cnU5Hs1ut6t79+569tlnJUldu3bVnj179PrrrxNk6tD777+vJUuW6N1339W1116rnTt3Ki0tTdHR0VxnD8XQUi1dddVV8vHx0bFjx5xeP3bsmCIjI11UlecaP368/vWvf+mLL75QixYtXF2Ox9m+fbuOHz+ubt26ydfXV76+vlq/fr1efvll+fr6ymazubpEjxEVFaUOHTo4vda+fXvl5eW5qCLP9Mc//lFTp07VyJEj1bFjR40aNUqTJk1Senq6q0vzWGe/+1z1vUiQqSV/f39df/31Wrt2reM1u92utWvXqnfv3i6szLMYhqHx48dr2bJl+vzzz5WQkODqkjzSrbfeqt27d2vnzp2OR/fu3ZWcnKydO3fKx8fH1SV6jD59+lRZQmD//v2Kj493UUWeqbS0VFar81ebj4+P7Ha7iyryfAkJCYqMjHT6XiwqKtKWLVsa5HuRoaXLMHnyZKWkpKh79+7q2bOnXnzxRZWUlOj+++93dWkeIzU1Ve+++64+/vhjhYaGOsZZw8LCFBgY6OLqPEdoaGiVeUfBwcFq2rQp85Hq2KRJk3TjjTfq2Wef1YgRI7R161a98cYbeuONN1xdmkcZOnSoZs6cqbi4OF177bXasWOHXnjhBT3wwAOuLs3UiouLdeDAAcfz3Nxc7dy5U02aNFFcXJzS0tL0zDPPKDExUQkJCZo2bZqio6M1bNiw+i+u3u+L8lCvvPKKERcXZ/j7+xs9e/Y0Nm/e7OqSPIqkCz4WLlzo6tI8Hrdf159PP/3UuO6664yAgACjXbt2xhtvvOHqkjxOUVGRMXHiRCMuLs5o1KiRcc011xiPP/64UVZW5urSTO2LL7644N/JKSkphmFU3oI9bdo0o3nz5kZAQIBx6623Gvv27WuQ2iyGwXKHAADAnJgjAwAATIsgAwAATIsgAwAATIsgAwAATIsgAwAATIsgAwAATIsgAwAATIsgAwAATIsgAwAATIsgA8BUbDabbrzxRt11111OrxcWFio2NlaPP/64iyoD4ApsUQDAdPbv368uXbrozTffVHJysiRp9OjR2rVrl7Zt2yZ/f38XVwigoRBkAJjSyy+/rBkzZuibb77R1q1bNXz4cG3btk2dO3d2dWkAGhBBBoApGYahW265RT4+Ptq9e7cmTJigJ554wtVlAWhgBBkApvXdd9+pffv26tixo7KysuTr6+vqkgA0MCb7AjCtv/3tbwoKClJubq5+/PFHV5cDwAXokQFgSl9++aX69eunf//733rmmWckSWvWrJHFYnFxZQAaEj0yAEyntLRUY8aM0dixY3XzzTfr7bff1tatW/X666+7ujQADYweGQCmM3HiRP3f//2fdu3apaCgIEnSggULNGXKFO3evVstW7Z0bYEAGgxBBoCprF+/XrfeeqvWrVunvn37Oh1LSkrSmTNnGGICvAhBBgAAmBZzZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGkRZAAAgGn9P+nA9koVbUz7AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Plotting routine\n",
        "def plot_data(X, y, obs=None, obs_std=None, inverse_df=None):\n",
        "    fig, ax = plt.subplots()\n",
        "    ax.scatter(X, y, color=\"cornflowerblue\", label=\"Experimental data\")\n",
        "    if (obs is not None) and (inverse_df is None):\n",
        "        y_mean = obs[\"y\"][0]\n",
        "        ax.axhline(y=y_mean, color=\"tomato\", label=\"Observed Output\")\n",
        "        if obs_std is not None:\n",
        "            y_std = obs_std[\"y\"][0]\n",
        "            ax.axhspan(\n",
        "                ymin=y_mean - y_std,\n",
        "                ymax=y_mean + y_std,\n",
        "                color=\"tomato\",\n",
        "                alpha=1 / 3,\n",
        "                label=\"Observation error\",\n",
        "            )\n",
        "    elif (obs is not None) and (obs_std is not None) and (inverse_df is not None):\n",
        "        y_mean = obs[\"y\"][0]\n",
        "        y_std = obs_std[\"y\"][0]\n",
        "        X_mean = inverse_df[\"mean\"][0]\n",
        "        X_std = inverse_df[\"sd\"][0]\n",
        "        ax.scatter(X_mean, y_mean, color=\"tomato\", label=\"Observed data\")\n",
        "        ax.add_patch(\n",
        "            patches.Ellipse(\n",
        "                (X_mean, y_mean),\n",
        "                X_std * 2,\n",
        "                y_std * 2,\n",
        "                color=\"tomato\",\n",
        "                alpha=1 / 3,\n",
        "                label=\"Uncertainty\",\n",
        "            )\n",
        "        )\n",
        "    ax.set_xlabel(\"X\")\n",
        "    ax.set_ylabel(\"y\")\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "    plt.close()\n",
        "\n",
        "\n",
        "# Plot the data\n",
        "plot_data(X=X, y=y, obs=obs, obs_std=obs_std)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eJQfXzF1Ozss"
      },
      "source": [
        "## Training the Model\n",
        "\n",
        "In this section, we'll be establishing the model training parameters, training our model on our experimental data, and producing the output of that training and prediction.\n",
        "\n",
        "First we need to establish what criteria our model training has. We do this via a `TrainParams` and `Emulator` object, which will include what we want to name our model, specify what the model inputs and outputs are, and fine-tune modelling parameters, such as the test_train_ratio.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "KoavhpYbOzss"
      },
      "outputs": [],
      "source": [
        "# Training parameters\n",
        "train_params = tl.TrainParams(train_test_ratio=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F9tF5d4wOzss"
      },
      "source": [
        "Next, we can actually train our model in the Cloud.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "KVd_nwglOzss"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Emulator inverse-methods-model has begun training.\n",
            "0:00:00: Job status: processing\n",
            "0:00:01: Job status: success\n",
            "Training of emulator inverse-methods-model is complete!\n"
          ]
        }
      ],
      "source": [
        "# Train model\n",
        "emulator = tl.Emulator(\"inverse-methods-model\")\n",
        "emulator.train(dataset, [\"X\"], [\"y\"], train_params)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ShX5pWQAOzss"
      },
      "source": [
        "## Inverse Problem\n",
        "\n",
        "Solve the inverse problem to obtain information about the input value that created the observed output value.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "pColRYXQOzst"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Job calibrate process ID: calibrate-lovely-ambitious-veterinarian\n",
            "0:00:00: Job status: processing\n",
            "0:00:01: Job status: processing\n",
            "0:00:03: Job status: processing\n",
            "0:00:05: Job status: processing\n",
            "0:00:07: Job status: processing\n",
            "0:00:09: Job status: processing\n",
            "0:00:12: Job status: processing\n",
            "0:00:15: Job status: processing\n",
            "0:00:18: Job status: processing\n",
            "0:00:21: Job status: processing\n",
            "0:00:25: Job status: processing\n",
            "0:00:30: Job status: processing\n",
            "0:00:34: Job status: processing\n",
            "0:00:39: Job status: processing\n",
            "0:00:45: Job status: success\n",
            "Calibration summary:\n",
            "    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
            "X  3.636  0.367   2.951    4.323      0.008    0.006    1939.0    2270.0   \n",
            "\n",
            "   r_hat  \n",
            "X    1.0  \n"
          ]
        }
      ],
      "source": [
        "# Solve the inverse problem\n",
        "inverse_df = emulator.calibrate(obs, obs_std)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 449
        },
        "id": "e0oxL44NOzst",
        "outputId": "78cce172-3d45-456b-b5d4-de9d41fb8445"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/var/folders/ml/754v6mn14dx2yfqknbvqmlmw0000gn/T/ipykernel_28132/1142193926.py:20: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
            "  X_mean = inverse_df[\"mean\"][0]\n",
            "/var/folders/ml/754v6mn14dx2yfqknbvqmlmw0000gn/T/ipykernel_28132/1142193926.py:21: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
            "  X_std = inverse_df[\"sd\"][0]\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAH0lEQVR4nO3deXiU1f3//9dMNrKQhbAkgRACREDZElGKYAFFAmguQcULTCFoq/0hKAGpy6dqbV2oKKi4REUF26pRq1Ctn4KBAn5D2SOLiiGmWBACKGQhiSRh5v79MZ9MHQJhEpKZuSfPx3XNBXPPmXveGULmlXPOfY7FMAxDAAAAJmT1dgEAAADNRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmFejtAlqb3W7X4cOH1b59e1ksFm+XAwAA3GAYhk6ePKmEhARZrefud/H7IHP48GElJiZ6uwwAANAMBw8eVLdu3c75uN8Hmfbt20tyvBGRkZFergYAALijoqJCiYmJzs/xc/H7IFM/nBQZGUmQAQDAZM43LYTJvgAAwLQIMgAAwLQIMgAAwLT8fo6Mu2w2m+rq6rxdBtBkQUFBCggI8HYZAOAVbT7IGIahI0eOqKyszNulAM0WHR2tuLg41koC0Oa0+SBTH2I6d+6ssLAwPghgKoZhqLq6WseOHZMkxcfHe7kiAPCsNh1kbDabM8TExsZ6uxygWUJDQyVJx44dU+fOnRlmAtCmtOnJvvVzYsLCwrxcCXBh6r+HmecFoK1p00GmHsNJMDu+hwG0VW16aAkAADSP3W6oqOS0yqrtig6zKiU+UFar53+pIsgAAIAmKSiuVW5+lUqrDOexmHCLpowIV1qvYI/WwtASWtyMGTM0ceJEb5dxQXr06KFnn322Sc/xh68bAM6noLhWOasrXUKMJJVWGcpZXamC4lqP1kOQaQF2u6HCQ3XaUlSjwkN1stuN8z/pAsyYMUMWi6XBbdy4ca36uu567rnntHz5cm+XIckxd2TlypXeLuOsvv32W1ksFu3cudPbpQCAW+x2Q7n5VY22yd1Y1eqfgz/F0NIF8lb32rhx47Rs2TKXYyEhIa32eu6w2WyyWCyKioryah0AgNZRVHK6QU/MmUorHXNn+nQN8khN9MhcAG92r4WEhCguLs7lFhMTI0lav369goOD9f/+3/9ztl+4cKE6d+6so0ePSpJGjRql2bNna/bs2YqKilLHjh310EMPyTD++7XU1NRo/vz56tq1q8LDwzV06FCtX7/e+fjy5csVHR2tjz76SBdffLFCQkJ04MCBBkMso0aN0l133aXs7GzFxMSoS5cuWrp0qaqqqnTrrbeqffv26t27t/7xj3+4fI1ffPGFxo8fr4iICHXp0kXTpk3TDz/84HLeu+++W/fee686dOiguLg4PfLII87He/ToIUmaNGmSLBaL835xcbGuv/56denSRREREbrsssu0Zs2aJr3/NptN8+bNU3R0tGJjY3Xvvfe6vHeStGrVKo0YMcLZ5rrrrlNxcbHz8eTkZElSamqqLBaLRo0aJUnatm2brrnmGnXs2FFRUVEaOXKkCgoKmlQfALSGsmp7i7ZrCQSZZvLF7rV6o0aNUnZ2tqZNm6by8nJ9/vnneuihh/Taa6+pS5cuznZvvvmmAgMDtXXrVj333HNavHixXnvtNefjs2fP1qZNm5Sbm6vdu3dr8uTJGjdunIqKipxtqqur9eSTT+q1117Tl19+qc6dO5+1pjfffFMdO3bU1q1bddddd2nmzJmaPHmyrrjiChUUFGjs2LGaNm2aqqurJUllZWW66qqrlJqaqu3bt2vVqlU6evSobr755gbnDQ8P15YtW7Rw4UL94Q9/UF5eniRHIJCkZcuWqaSkxHm/srJSEyZM0Nq1a/X5559r3LhxysjI0IEDB9x+jxctWqTly5frjTfeUH5+vk6cOKEVK1a4tKmqqtK8efO0fft2rV27VlarVZMmTZLd7vgPvnXrVknSmjVrVFJSog8//FCSdPLkSWVlZSk/P1+bN29WSkqKJkyYoJMnT7pdHwC0hugw92KDu+1ahOHnysvLDUlGeXl5g8d+/PFH46uvvjJ+/PHHJp/36+9qjV+9ePy8t6+/q22JL8NFVlaWERAQYISHh7vcHn/8cWebmpoaY/DgwcbNN99sXHzxxcbtt9/uco6RI0ca/fr1M+x2u/PYfffdZ/Tr188wDMP4z3/+YwQEBBiHDh1yed7VV19tPPDAA4ZhGMayZcsMScbOnTsb1Hf99de7vNaIESOc90+fPm2Eh4cb06ZNcx4rKSkxJBmbNm0yDMMwHn30UWPs2LEu5z148KAhySgsLDzreQ3DMC677DLjvvvuc96XZKxYseIs76KrSy65xHj++eed95OSkoxnnnnmnO3j4+ONhQsXOu/X1dUZ3bp1c/m6z/T9998bkow9e/YYhmEY+/fvNyQZn3/+eaO12Ww2o3379sbHH398zjYX8r0MAO6y2ezGb5afaPRz7zdvnjBsNvv5T3YejX1+/xRzZJrJ291ro0ePVk5OjsuxDh06OP8eHByst956SwMHDlRSUpKeeeaZBuf42c9+5rKQ2rBhw7Ro0SLZbDbt2bNHNptNF110kctzampqXLZzCA4O1sCBA89b70/bBAQEKDY2VgMGDHAeq+8pqt8zaNeuXVq3bp0iIiIanKu4uNhZ15mvHR8f7zzHuVRWVuqRRx7RJ598opKSEp0+fVo//vij2z0y5eXlKikp0dChQ53HAgMDNWTIEJfhpaKiIj388MPasmWLfvjhB2dPzIEDB9S/f/9znv/o0aN68MEHtX79eh07dkw2m03V1dVN6jECgNZgtTrmgOasrjxnmynDwz26ngxBppm83b0WHh6u3r17N9rmX//6lyTpxIkTOnHihMLDw90+f2VlpQICArRjx44Ge/f8NFyEhoa6tapsUJDrpC+LxeJyrP4c9R/2lZWVysjI0JNPPtngXD/dGPFs560/x7nMnz9feXl5evrpp9W7d2+FhobqpptuUm1ty85pysjIUFJSkpYuXaqEhATZ7Xb179//vK+TlZWl48eP67nnnlNSUpJCQkI0bNiwFq8PAJojrVewZqZHNLzQJcKiKcM9v44MQaaZUuIDFRNuaXT2dkyERSnx3nmLi4uLNXfuXC1dulTvvvuusrKytGbNGlmt/w1WW7ZscXlO/XyMgIAApaamymaz6dixY7ryyis9Xb7S0tL0wQcfqEePHgoMbP57GBQUJJvN5nJs48aNmjFjhiZNmiTJEZq+/fZbt88ZFRWl+Ph4bdmyRT//+c8lSadPn9aOHTuUlpYmSTp+/LgKCwu1dOlS5/uXn5/vcp7gYMd/9rPV99JLL2nChAmSpIMHD7pMcgYAb0vrFazByUE+sbIvk32bqb57rTGt2b1WU1OjI0eOuNzqP+xsNpt+8YtfKD09XbfeequWLVum3bt3a9GiRS7nOHDggObNm6fCwkK98847ev755zVnzhxJ0kUXXaTMzExNnz5dH374ofbv36+tW7dqwYIF+uSTT1rla/qpWbNm6cSJE5o6daq2bdum4uJirV69WrfeemuDD/7G9OjRQ2vXrtWRI0dUWloqSUpJSdGHH36onTt3ateuXbrlllvO24tzpjlz5uiPf/yjVq5cqa+//lp33nmnysrKnI/HxMQoNjZWr776qr755hv985//1Lx581zO0blzZ4WGhjonMpeXlzvr+/Of/6y9e/dqy5YtyszMdO5wDQC+wmq1qE/XIA1NCVGfrkFeCTESQeaC1HevxYS7/uPFRFg0Mz2iVbvXVq1apfj4eJfbiBEjJEmPP/64/vOf/+iVV16R5BiKefXVV/Xggw9q165dznNMnz5dP/74oy6//HLNmjVLc+bM0R133OF8fNmyZZo+fbruuece9enTRxMnTtS2bdvUvXv3Vvu66iUkJGjjxo2y2WwaO3asBgwYoOzsbEVHR7v0Kp3PokWLlJeXp8TERKWmpkqSFi9erJiYGF1xxRXKyMhQenq6syfFXffcc4+mTZumrKwsDRs2TO3bt3f28EiS1WpVbm6uduzYof79+2vu3Ll66qmnXM4RGBioJUuW6JVXXlFCQoKuv/56SdLrr7+u0tJSpaWladq0abr77rvPeTUYALR1FsMwPH99sAdVVFQoKipK5eXlioyMdHns1KlT2r9/v5KTk9WuXbtmv4avbJzVFKNGjdLgwYObvAw/fFNLfS8DgK9o7PP7p5gj0wLqu9cAAIBnMbQEAABMix6ZNuqnWw0AAGBW9MgAAADTIsgAAADTIsgAAADTIsgAAADTYrIvAAB+zIxrnTUFPTJ+rkePHn616F1zvp4ZM2Zo4sSJrVIPAPiyguJa3f/nMj39t5N6La9KT//tpO7/c5kKiv1nE1qCjEkdPHhQt912mxISEhQcHKykpCTNmTNHx48f93Zppvftt9/KYrFo586d3i4FAJqtoLhWOasrG2xuXFplKGd1pd+EGYJMS7DbpK93SVvWOf60u7+pYXP8+9//1pAhQ1RUVKR33nlH33zzjV5++WWtXbtWw4YN04kTJ1r19Rtjs9mavAEjAKBl2e2GcvOrGm2Tu7FKdrv5dykiyFyoHfnSfVnS0/dJS590/HlfluN4K5k1a5aCg4P16aefauTIkerevbvGjx+vNWvW6NChQ/rtb3/r0v7kyZOaOnWqwsPD1bVrV7344ovOxwzD0COPPKLu3bsrJCRECQkJuvvuu52P19TUaP78+eratavCw8M1dOhQl8X0li9frujoaH300Ue6+OKLFRISotdee03t2rVz2Q1acuwYfdVVVznv5+fn68orr1RoaKgSExN19913q6rqv//xjh07poyMDIWGhio5OVlvvfXWed8bm82mefPmKTo6WrGxsbr33nt15nZiq1at0ogRI5xtrrvuOhUXFzsfT05OliSlpqbKYrFo1KhRkqRt27bpmmuuUceOHRUVFaWRI0eqoKDgvDUBgKcVlZxu0BNzptJKx9wZsyPIXIgd+VLOY1LpD67HS39wHG+FMHPixAmtXr1ad955p0JDQ10ei4uLU2Zmpt59912XD++nnnpKgwYN0ueff677779fc+bMUV5eniTpgw8+0DPPPKNXXnlFRUVFWrlypQYMGOB87uzZs7Vp0ybl5uZq9+7dmjx5ssaNG6eioiJnm+rqaj355JN67bXX9OWXXyozM1PR0dH64IMPnG1sNpveffddZWZmSpKKi4s1btw43Xjjjdq9e7feffdd5efna/bs2c7nzJgxQwcPHtS6dev017/+VS+99JKOHTvW6PuzaNEiLV++XG+88Yby8/N14sQJrVixwqVNVVWV5s2bp+3bt2vt2rWyWq2aNGmSsydp69atkqQ1a9aopKREH374oSRHIMzKylJ+fr42b96slJQUTZgwQSdPnjzPvxoAeFZZtXs94+6282VctdRcdpuU+3LjbXJfkVKHSdaAFnvZoqIiGYahfv36nfXxfv36qbS0VN9//706d+4sSRo+fLjuv/9+SdJFF12kjRs36plnntE111yjAwcOKC4uTmPGjFFQUJC6d++uyy+/XJJ04MABLVu2TAcOHFBCQoIkaf78+Vq1apWWLVumJ554QpJUV1enl156SYMGDXLWMWXKFL399tv65S9/KUlau3atysrKdOONN0qSFixYoMzMTGVnZ0uSUlJStGTJEo0cOVI5OTk6cOCA/vGPf2jr1q267LLLJEmvv/76Ob/ues8++6weeOAB3XDDDZKkl19+WatXr3ZpU19DvTfeeEOdOnXSV199pf79+6tTp06SpNjYWMXFxTnb/bQ3SZJeffVVRUdHa8OGDbruuusarQsAPCk6zL1+Cnfb+TLzfwXesu+Lhj0xZyr93tGuFZw5XNKYYcOGNbi/d+9eSdLkyZP1448/qmfPnrr99tu1YsUKnT7t6Grcs2ePbDabLrroIkVERDhvGzZscBmKCQ4O1sCBA11eIzMzU+vXr9fhw4clSW+99ZauvfZaRUdHS5J27dql5cuXu5w3PT1ddrtd+/fv1969exUYGKhLL73Uec6+ffs6n3825eXlKikp0dChQ53HAgMDNWTIEJd2RUVFmjp1qnr27KnIyEj16NFDkiO4Nebo0aO6/fbblZKSoqioKEVGRqqysvK8zwMAT0uJD1RMeOOXWMdEWJQSb/7+DPN/Bd5S7uaEWnfbual3796yWCzau3evJk2a1ODxvXv3KiYmxtmrcD6JiYkqLCzUmjVrlJeXpzvvvFNPPfWUNmzYoMrKSgUEBGjHjh0KCHDtVYqIiHD+PTQ0VBaL63+Yyy67TL169VJubq5mzpypFStWaPny5c7HKysr9etf/9plPk697t27a9++fW7V3xwZGRlKSkrS0qVLlZCQILvdrv79+6u2tvEZ/FlZWTp+/Liee+45JSUlKSQkRMOGDTvv8wDA06xWi6aMCFfO6spztpkyPNwv1pMhyDRXVIeWbeem2NhYXXPNNXrppZc0d+5cl3kyR44c0VtvvaXp06e7BIvNmze7nGPz5s0uQzShoaHKyMhQRkaGZs2apb59+2rPnj1KTU2VzWbTsWPHdOWVVza51szMTL311lvq1q2brFarrr32WudjaWlp+uqrr9S7d++zPrdv3746ffq0duzY4RxaKiwsbDCB+KeioqIUHx+vLVu26Oc//7kkOc+RlpYmSTp+/LgKCwu1dOlS59eUn+86lyk4OFiSY17PT23cuFEvvfSSJkyYIMlxCfwPP5ynVw4AvCStV7BmpkcoN7/KZeJvTIRFU4aHK61XsBerazkEmea6qL8U07Hx4aWYTo52LeyFF17QFVdcofT0dD322GNKTk7Wl19+qd/85jfq2rWrHn/8cZf2Gzdu1MKFCzVx4kTl5eXp/fff1yeffCLJcdWRzWbT0KFDFRYWpr/85S8KDQ1VUlKSYmNjlZmZqenTp2vRokVKTU3V999/r7Vr12rgwIEuweRsMjMz9cgjj+jxxx/XTTfdpJCQEOdj9913n372s59p9uzZ+tWvfqXw8HB99dVXysvL0wsvvKA+ffpo3Lhx+vWvf62cnBwFBgYqOzu7wQTnM82ZM0d//OMflZKSor59+2rx4sUu4ScmJkaxsbF69dVXFR8frwMHDjjnD9Xr3LmzQkNDtWrVKnXr1k3t2rVTVFSUUlJS9Oc//1lDhgxRRUWFfvOb35y3HgDwprRewRqcHMTKvjgLa4A05f9rvM2UX7foRN96KSkp2r59u3r27Kmbb75ZvXr10h133KHRo0dr06ZN6tDBtRfonnvu0fbt25WamqrHHntMixcvVnp6uiQpOjpaS5cu1fDhwzVw4ECtWbNGH3/8sWJjYyVJy5Yt0/Tp03XPPfeoT58+mjhxorZt26bu3buft87evXvr8ssv1+7du51XK9UbOHCgNmzYoH379unKK69UamqqHn74Yeek4vrXTkhI0MiRI3XDDTfojjvucE5gPpd77rlH06ZNU1ZWloYNG6b27du7DMFZrVbl5uZqx44d6t+/v+bOnaunnnrK5RyBgYFasmSJXnnlFSUkJOj666+X5JhsXFpaqrS0NE2bNk133333eesBAG+zWi3q0zVIQ1NC1KdrkF+FGEmyGE2ZNWpCFRUVioqKUnl5uSIjI10eO3XqlPbv36/k5GS1a9eueS+wI99x9dJPe2ZiOjlCzKUjLqBywH0t8r0MAD6ksc/vn2Jo6UJdOsJxifW+LxwTe6M6OIaTWqEnBgAAuCLItARrgNR30PnbAQCAFsUcGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFqsI3M2Naek03Wee73AICmE1VglqUePHsrOzlZ2dra3SwEAmABB5kw1p6Ttn0lVJz33muHtpSE/b1KYGTVqlAYPHqxnn33W5fjy5cuVnZ3d6C7RnmCxWLRixQpNnDixSc/btm2bwsPD3W6/fv16jR49WqWlpYqOjm5akQAA0yPInOl0nSPEBAVLwR7Y4ry21vF6p+v8olemtrZWwRfwvnXq1KkFqwEA+DvmyJxLcLAUEtr6t1YMSzNmzNDEiRP19NNPKz4+XrGxsZo1a5bq6v47bFZTU6P77rtPiYmJCgkJUe/evfX66687H//iiy80fvx4RUREqEuXLpo2bZp++OG/G2SOGjVKs2fPVnZ2tjp27Kj09HT16NFDkjRp0iRZLBbn/eLiYl1//fXq0qWLIiIidNlll2nNmjUuNffo0cOll8lisei1117TpEmTFBYWppSUFH300UeSpG+//VajR4+WJMXExMhisWjGjBn605/+pNjYWNXU1Lice+LEiZo2bdoFv68AAN9BkPFz69atU3FxsdatW6c333xTy5cv1/Lly52PT58+Xe+8846WLFmivXv36pVXXlFERIQkqaysTFdddZVSU1O1fft2rVq1SkePHtXNN9/s8hpvvvmmgoODtXHjRr388svatm2bJGnZsmUqKSlx3q+srNSECRO0du1aff755xo3bpwyMjJ04MCBRr+G3//+97r55pu1e/duTZgwQZmZmTpx4oQSExP1wQcfSJIKCwtVUlKi5557TpMnT5bNZnMGHkk6duyYPvnkE912220X/J4CAHwHQ0t+LiYmRi+88IICAgLUt29fXXvttVq7dq1uv/127du3T++9957y8vI0ZswYSVLPnj2dz33hhReUmpqqJ554wnnsjTfeUGJiovbt26eLLrpIkpSSkqKFCxc2eO3o6GjFxcU57w8aNEiDBv13c81HH31UK1as0EcffaTZs2ef82uYMWOGpk6dKkl64okntGTJEm3dulXjxo1Thw4dJEmdO3d2mSNzyy23aNmyZZo8ebIk6S9/+Yu6d++uUaNGufvWAQBMgB4ZP3fJJZcoICDAeT8+Pl7Hjh2TJO3cuVMBAQEaOXLkWZ+7a9curVu3ThEREc5b3759JTmGiepdeumlbtVSWVmp+fPnq1+/foqOjlZERIT27t173h6ZgQMHOv8eHh6uyMhI59dwLrfffrs+/fRTHTp0SJJjEvSMGTNksVjcqhUAYA70yJhUZGSkysvLGxwvKytTVFSU835QUJDL4xaLRXa7XZIUGhra6GtUVlYqIyNDTz75ZIPH4uPjnX939yqj+fPnKy8vT08//bR69+6t0NBQ3XTTTaqtrW30eY19DeeSmpqqQYMG6U9/+pPGjh2rL7/8Up988olbdQKAp9nthopKTqus2q7oMKtS4gNltfKLlzsIMibVp08fffrppw2OFxQUOId8zmfAgAGy2+3asGGDc2jpp9LS0vTBBx+oR48eCgxs2rdKUFCQbDaby7GNGzdqxowZmjRpkiRHUPr222+bdN4z1V8hdeZrSdKvfvUrPfvsszp06JDGjBmjxMTEC3otAGgNBcW1ys2vUmmV4TwWE27RlBHhSuvlgatnTY6hJZOaOXOm9u3bp7vvvlu7d+9WYWGhFi9erHfeeUf33HOPW+fo0aOHsrKydNttt2nlypXav3+/1q9fr/fee0+SNGvWLJ04cUJTp07Vtm3bVFxcrNWrV+vWW289a3A489xr167VkSNHVFpaKskxl+bDDz/Uzp07tWvXLt1yyy3n7Vk5n6SkJFksFv3973/X999/r8rKSudjt9xyi7777jstXbqUSb4AfFJBca1yVle6hBhJKq0ylLO6UgXFjfdYgyBzbrW1Us2PrX87z7DKufTs2VOfffaZvv76a40ZM0ZDhw7Ve++9p/fff1/jxo1z+zw5OTm66aabdOedd6pv3766/fbbVVVVJUlKSEjQxo0bZbPZNHbsWA0YMEDZ2dmKjo6W1dr4t86iRYuUl5enxMREpaamSpIWL16smJgYXXHFFcrIyFB6errS0tKa9fXX69q1q37/+9/r/vvvV5cuXVwmDUdFRenGG29UREREkxfmA4DWZrcbys2varRN7sYq2e1Go23aOothGH79DlVUVCgqKkrl5eWKjIx0eezUqVPav3+/kpOT1a7d/y1GZ5KVfeGeq6++WpdccomWLFni7VJa1Vm/lwH4tMJDdXr6b+f/rJl/fXv16Rp03nb+prHP759ijsyZQto5QgV7LZlaaWmp1q9fr/Xr1+ull17ydjkA0EBZtXtD6+62a6sIMmcT0o5gYXKpqakqLS3Vk08+qT59+ni7HABoIDrMvdkd7rZrqwgy8EsXejUUALS2lPhAxYRbGkz0/amYCItS4vmobgwxDwAAL7BaHZdYN2bK8HDWkzkPgowkP5/vjDaA72HAnNJ6BWtmeoRiwl3DSkyERTPTI1hHxg1tur+qfsXY6urq865yC/iy6upqSQ1XQQbg+9J6BWtwchAr+zZTmw4yAQEBio6Odu7bExYWxl48MBXDMFRdXa1jx44pOjraZV8tAOZhtVra5CXWLaFNBxlJzt2Zz7cJIeDLztxpHADaijYfZCwWi+Lj49W5c2fV1Xlw7RighQQFBdETA6DNavNBpl5AQAAfBgAAmIxXr1pasGCBLrvsMrVv316dO3fWxIkTVVhY6NLm1KlTmjVrlmJjYxUREaEbb7xRR48e9VLFAADAl3g1yGzYsEGzZs3S5s2blZeXp7q6Oo0dO9a5aaEkzZ07Vx9//LHef/99bdiwQYcPH9YNN9zgxaoBAICv8KlNI7///nt17txZGzZs0M9//nOVl5erU6dOevvtt3XTTTdJkr7++mv169dPmzZt0s9+9rMG56ipqVFNTY3zfkVFhRITE8+76RQAAPAd7m4a6VML4pWXl0uSOnToIEnasWOH6urqNGbMGGebvn37qnv37tq0adNZz7FgwQJFRUU5b4mJia1fOACgTbDbDRUeqtOWohoVHqqT3e4zfQFtls9M9rXb7crOztbw4cPVv39/SdKRI0cUHBys6Ohol7ZdunTRkSNHznqeBx54QPPmzXPer++RAQDgQhQU1yo3v8plb6SYcMc2A6zA6z0+E2RmzZqlL774Qvn5+Rd0npCQEIWEhLRQVQAAOEJMzurKBsdLqwzlrK5kOwEv8omhpdmzZ+vvf/+71q1bp27dujmPx8XFqba2VmVlZS7tjx49yuJfAACPsNsN5eZXNdomd2MVw0xe4tUgYxiGZs+erRUrVuif//ynkpOTXR6/9NJLFRQUpLVr1zqPFRYW6sCBAxo2bJinywUAtEFFJaddhpPOprTSUFHJaQ9VhJ/y6tDSrFmz9Pbbb+tvf/ub2rdv75z3EhUVpdDQUEVFRemXv/yl5s2bpw4dOigyMlJ33XWXhg0bdtYrlgAAaGll1fYWbYeW5dUgk5OTI0kaNWqUy/Fly5ZpxowZkqRnnnlGVqtVN954o2pqapSenq6XXnrJw5UCANqq6DD3Bi/cbYeW5VPryLQGd69DBwDgbOx2Q/f/uazR4aWYCIv++ItoWa0WD1bm30y5jgwAAL7GanVcYt2YKcPDCTFeQpABAOA80noFa2Z6hGLCXcNKTISFS6+9zGfWkQEAoCXY7Y4riMqq7YoOsyolPrBFekvSegVrcHJQq5wbzUeQAQD4jdZefddqtahP16ALPg9aDkNLAAC/UL/67pmTcutX3y0orvVSZWhNBBkAgOmx+m7bRZABAJgeq++2XQQZAIDpsfpu20WQAQCYHqvvtl38iwIATC8lPrDBGi9niomwKCWei3X9DUEGAGB6rL7bdhFkAAB+gdV32yb62AAAfoPVd9seggwAwK+w+m7bwtASAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwLYIMAAAwrUBvFwAAaJvsdkNFJadVVm1XdJhVKfGBslot3i4LJkOQAQB4XEFxrXLzq1RaZTiPxYRbNGVEuNJ6BXuxMpgNQ0sAAI8qKK5VzupKlxAjSaVVhnJWV6qguNZLlcGMCDIAgHOy2w0VHqrTlqIaFR6qk91unP9J5zlfbn5Vo21yN1Zd8Oug7WBoCQBwVq0x/FNUcrpBT8yZSisdc2f6dA1q1mugbaFHBgDQQGsN/5RV21u0HUCQAQC4aM3hn+gw9z523G0H8J0CAHDRlOGfpkqJD1RMeOOXWMdEWJQSz8wHuIcgAwBw0ZrDP1arY45NY6YMD2c9GbiNIAMAcNHawz9pvYI1Mz2iQc9MTIRFM9MjWEcGTULfHQDARf3wT2PDSxc6/JPWK1iDk4NY2RcXjB4ZAIALTw3/WK0W9ekapKEpIerTNYgQg2YhyAAAGmD4B2bB0BIA4KwY/oEZEGQAAOdUP/wD+CqGlgAAgGkRZAAAgGkxtAQAfsBuN5jLgjaJIAMAJtcau1QDZsHQEgCYWGvtUg2YBUEGAEyqNXepBsyCIAMAJtWau1QDZkGQAQCTas1dqgGzIMgAgEm19i7VgBnw3Q0AJlW/S3VjLnSXasDXEWQAwKQ8tUs14MsIMgBgYuxSjbaO/kYAMDl2qUZbRpABAD/ALtVoqxhaAgAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApkWQAQAApuXVIPPZZ58pIyNDCQkJslgsWrlypcvjM2bMkMVicbmNGzfOO8UCAACf49UgU1VVpUGDBunFF188Z5tx48appKTEeXvnnXc8WCEAAPBlgd588fHjx2v8+PGNtgkJCVFcXJzb56ypqVFNTY3zfkVFRbPrAwAAvs3n58isX79enTt3Vp8+fTRz5kwdP3680fYLFixQVFSU85aYmOihSgEAgKdZDMMwvF2EJFksFq1YsUITJ050HsvNzVVYWJiSk5NVXFys//mf/1FERIQ2bdqkgICAs57nbD0yiYmJKi8vV2RkZGt/GQAAoAVUVFQoKirqvJ/fXh1aOp8pU6Y4/z5gwAANHDhQvXr10vr163X11Vef9TkhISEKCQnxVIkAAMCLfH5o6ad69uypjh076ptvvvF2KQAAwAeYKsh89913On78uOLj471dCgAA8AFeHVqqrKx06V3Zv3+/du7cqQ4dOqhDhw76/e9/rxtvvFFxcXEqLi7Wvffeq969eys9Pd2LVQNA89jthopKTqus2q7oMKtS4gNltVq8XRZgal4NMtu3b9fo0aOd9+fNmydJysrKUk5Ojnbv3q0333xTZWVlSkhI0NixY/Xoo48yBwaA6RQU1yo3v0qlVf+9viIm3KIpI8KV1ivYi5UB5uYzVy21FndnPQNAaykorlXO6spzPj4zPYIwA5zB3c9vU82RAQCzsdsN5eZXNdomd2OV7Ha//p0SaDUEGQBoRUUlp12Gk86mtNIxdwZA0xFkAKAVlVXbW7QdAFcEGQBoRdFh7v2YdbcdAFf8zwGAVpQSH6iY8MYvsY6JsCgl3qcXWgd8FkEGAFqR1eq4xLoxU4aHs54M0EwEGQBoZWm9gjUzPaJBz0xMhIVLr4ELRF8mAHhAWq9gDU4OYmVfoIURZADAQ6xWi/p0DfJ2GYBfYWgJAACYFkEGAACYFkEGAACYFkEGAACYVpODTFZWlj777LPWqAUA3Ga3Gyo8VKctRTUqPFTHpotAG9Xkq5bKy8s1ZswYJSUl6dZbb1VWVpa6du3aGrUBwFkVFNcqN7/KZTPGmHDHwnOsyQK0LU3ukVm5cqUOHTqkmTNn6t1331WPHj00fvx4/fWvf1VdXV1r1AgATgXFtcpZXdlgR+nSKkM5qytVUFzrpcoAeEOz5sh06tRJ8+bN065du7Rlyxb17t1b06ZNU0JCgubOnauioqKWrhMAZLcbys2varRN7sYqhpmANuSCJvuWlJQoLy9PeXl5CggI0IQJE7Rnzx5dfPHFeuaZZ1qqRgCQJBWVnG7QE3Om0kpDRSWnPVQRAG9rcpCpq6vTBx98oOuuu05JSUl6//33lZ2drcOHD+vNN9/UmjVr9N577+kPf/hDa9QLoA0rq7a3aDsA5tfkyb7x8fGy2+2aOnWqtm7dqsGDBzdoM3r0aEVHR7dAeQDwX9Fh7v3u5W47AObX5CDzzDPPaPLkyWrXrt0520RHR2v//v0XVBgAnCklPlAx4ZZGh5diIixKiWcbOaCtaPKvLdOmTWs0xABAa7FaHZdYN2bK8HB2lAbaEPpfAZhKWq9gzUyPUEy4a1iJibBoZnoE68gAbQz9rwBMJ61XsAYnB6mo5LTKqu2KDrMqJT6QnhigDSLIADAlq9WiPl2DvF0GAC8jyADAT9jtBj09gIkQZADg/7CHE2A+TPYFALGHE2BWBBkAbR57OAHmRZAB0OaxhxNgXgQZAG0eezgB5kWQAdDmsYcTYF78rwTQ5tXv4dQY9nACfBNBBkCbxx5OgHkRZABA7OEEmBX9pABajdlWyWUPJ8B8CDIAWoVZV8llDyfAXBhaAtDiWCUXgKcQZAC0KFbJBeBJBBkALYpVcgF4EkEGQItilVwAnkSQAdCiWCUXgCfxkwRAi2KVXACeRJAB0KJYJReAJxFkALQ4VskF4Cn07QJolRV4WSUXgCcQZIA2rjVX4GWVXACtjaEloA1jBV4AZkeQAdooVuAF4A8IMkAbxQq8APwBQQZoo1iBF4A/IMgAbRQr8ALwB/yEAtooVuAF4A8IMkAbxQq8APwBQQZow1iBF4DZ0WcMtHGswAvAzAgyAFiBF4BpMbQEAABMix4ZwMzqaqWaGslukwxDslikwEApuJ3jTwDwc/ykA8zCbpNOVkgVpVLpD9LJMqnmlHS6TrLbJRmSrJLVKgUFS+1CpcgYKSbW8Wd4e0fQAQA/QpABfN2paunoIem7/Y4gU3vq/3pdQqTAYCm8nSO81LPZHD01leXS8WNSsV1q106KipUSk6VO8Y6gAwB+gCAD+KqaU9KBYulgsVRZ8X89LFFScOfGnxckR9t6hiHV/CidOCYdOyRFdZB6XCQlJDH8BMD0+CkG+KIfjkiFe6TjR6WISCk+sfnDQhaL1C7McbPZHENTuzZLxw5LfQZI7aNbtHQA8CSCDOBrDn0rfVngmPvSpasUENBy5w4IkGI6OoaeSg5IVSelgZc7jgGACXH5NeBLDv9H+nKHZJHUOb5lQ8xPBQVLcd0cQ1a7tjgmDwOACRFkAF9xslz6epfj703pIbHbpUP/kYq+dPxpt7v3PIvFEZYqK6S9uxy9NABgMgwtAb7Abpf+vdcx1BPXzf3n/ftraWOe43n1wttLw6+RevY9//MtFqlTnHT0sPSfb6TeFze9dgDwInpkAF9wskwq+U7q0NH9Sb3//lr69EPXECM57n/6oeNxdwQESOER0sF/0ysDwHQIMoAvKDsu1dZIIaHnbys5enA25jXeZuMa94eZIto7AlDZCffaA4CPIMgAvqC8tGlrupQcbNgTc6aqCkc7dwQESob9/OcEAB9DkAF8gc3mujrv+VRXtmw7SZLF/R4cAPARXg0yn332mTIyMpSQkCCLxaKVK1e6PG4Yhh5++GHFx8crNDRUY8aMUVFRkXeKBVpTULAjzLgrLKJl20mOFYBZ6ReAyXg1yFRVVWnQoEF68cUXz/r4woULtWTJEr388svasmWLwsPDlZ6erlOnTnm4UqCVxcQ6goxhuNc+PtFxdVJjwv9vRWB31NZIwUFSZLR77QHAR3j116/x48dr/PjxZ33MMAw9++yzevDBB3X99ddLkv70pz+pS5cuWrlypaZMmXLW59XU1KimpsZ5v6KiouULB1padKwUFu6YoxIRef72VqvjEutPPzx3m+Fj3B+uOlnu2KqAIAPAZHx2jsz+/ft15MgRjRkzxnksKipKQ4cO1aZNm875vAULFigqKsp5S0x08zdSwJvC20vde0kVZe4PMfXsK429oWHPTHik47g768hIjs0pT9dJyX0layutJAwArcRnB8SPHDkiSerSpYvL8S5dujgfO5sHHnhA8+bNc96vqKggzMAcklKkH45Kx0qkLgnu9ab07OvYybrkoGNib1iEYzjJ3Z6Y03XS8WNSUu+mLcQHAD7CZ4NMc4WEhCgkJMTbZQBNF9JOuuRSafdWx0q77oYZq1XqmtT016urk74vkRKSpD6DmnbVFAD4CJ/9yRUXFydJOnr0qMvxo0ePOh8D/E5UjDRoqGOF3yPfSaeqW+d1TpY7QkzXHtKAy6R2bi7EBwA+xmeDTHJysuLi4rR27VrnsYqKCm3ZskXDhg3zYmVAK4uMli4d4Rg2qqhw9M6c+vHCz2sYjsnEJQel06eli9McoSmk3YWfGwC8xKtDS5WVlfrmm2+c9/fv36+dO3eqQ4cO6t69u7Kzs/XYY48pJSVFycnJeuihh5SQkKCJEyd6r2jAE9qFOYaZunR17IF0rEQq/UEKDXPMgwl2c/jUMKSaH6XqKkcYCgt3TOpNTJaiOrTu1wAAHuDVILN9+3aNHj3aeb9+km5WVpaWL1+ue++9V1VVVbrjjjtUVlamESNGaNWqVWrXjt8g0QZYLFKneKljnFR+Qjr6nSPQVJRLtaccjwcGSoFBjquNLJIMSbbTjkm8NpvjfrtQKSJK6nWxIxi5c3k3AJiExTDcXYHLnCoqKhQVFaXy8nJFRvIDHCZnt0knKxy7Zdf86BgqqjzpCC6G4ZiwGxTsCCthEY6enchoxyXa7u6qDQA+wN3Pb7+7agnwa9YAx4TgqBhvVwIAPsFnJ/sCAACcD0EGAACYFkNLgEnY7YaKSk6rrNqu6DCrUuIDZbUy7wVA20aQAUygoLhWuflVKq3679z8mHCLpowIV1qvYC9WBgDexdAS4OMKimuVs7rSJcRIUmmVoZzVlSoorvVSZQDgfQQZwIfZ7YZy86sabZO7sUp2u1+vogAA50SQAXxYUcnpBj0xZyqtdMydAYC2iDkyQAtryUm5ZdX2Fm0HAP6GIAO0oJaelBsd5l6nqbvtAMDf8NMPaCGtMSk3JT5QMeGN9+bERFiUEs/vJADaJoIM0AJaa1Ku1erozWnMlOHhrCcDoM0iyAAtoDUn5ab1CtbM9IgGPTMxERbNTI9gHRkAbRr90UALaO1JuWm9gjU4OYiVfQHgDAQZoAV4YlKu1WpRn65BzX4+APgjhpaAFsCkXADwDoIM0AKYlAsA3kGQAVoIk3IBwPPo5wZaEJNyAcCzCDJAC2NSLgB4DkNLAADAtAgyAADAtBhaQpvUkjtUAwC8hyCDNqeld6gGAHgPQ0toU1pjh2oAgPcQZNBmtNYO1QAA7yHIoM1ozR2qAQDeQZBBm9HaO1QDADyPIIM2wxM7VAMAPIuf2Ggz2KEaAPwPQQZtBjtUA4D/IcigTWGHagDwL/Sho81hh2oA8B8EGbRJ7FANAP6BoSUAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBaBBkAAGBa7LUEn2a3G2zuCAA4J4IMfFZBca1y86tUWmU4j8WEWzRlRLjSegV7sTIAgK9gaAk+qaC4VjmrK11CjCSVVhnKWV2pguJaL1UGAPAlBBn4HLvdUG5+VaNtcjdWyW43Gm0DAPB/BBn4nKKS0w16Ys5UWumYOwMAaNsIMvA5ZdX2Fm0HAPBfBBn4nOgw974t3W0HAPBffBLA56TEByomvPFLrGMiLEqJ56I7AGjrCDLwOVar4xLrxkwZHs56MgAAggx8U1qvYM1Mj2jQMxMTYdHM9AjWkQEASGJBPPiwtF7BGpwcxMq+AIBzIsjAp1mtFvXpGuTtMgAAPoqhJQAAYFoEGQAAYFoMLeGCsUM1AMBbCDK4IOxQDQDwJoaW0GzsUA0A8DaCDJqFHaoBAL6AIINmYYdqAIAvIMigWdihGgDgCwgyaBZ2qAYA+AI+ZdAs7FANAPAFBBk0CztUAwB8AUEGzcYO1QAAb6PfHxeEHaoBAN5EkMEFY4dqAIC3+PTQ0iOPPCKLxeJy69u3r7fLAgAAPsLne2QuueQSrVmzxnk/MNDnSwYAAB7i86kgMDBQcXFx3i4DAAD4IJ8eWpKkoqIiJSQkqGfPnsrMzNSBAwcabV9TU6OKigqXGwAA8E8+HWSGDh2q5cuXa9WqVcrJydH+/ft15ZVX6uTJk+d8zoIFCxQVFeW8JSYmerBiAADgSRbDMEyzPXFZWZmSkpK0ePFi/fKXvzxrm5qaGtXU1DjvV1RUKDExUeXl5YqMjPRUqQAA4AJUVFQoKirqvJ/fPj9H5qeio6N10UUX6Ztvvjlnm5CQEIWEhHiwKgAA4C0+PbR0psrKShUXFys+Pt7bpbQau91Q4aE6bSmqUeGhOtntpukwAwDA43y6R2b+/PnKyMhQUlKSDh8+rN/97ncKCAjQ1KlTvV1aqygorlVufpVKq/4bXmLCHXsasdw/AAAN+XSPzHfffaepU6eqT58+uvnmmxUbG6vNmzerU6dO3i6txRUU1ypndaVLiJGk0ipDOasrVVBc66XKAADwXT7dI5Obm+vtEjzCbjeUm1/VaJvcjVUanBzEHkYAAPyET/fItBVFJacb9MScqbTSUFHJ6Qt6HebfAAD8jU/3yLQVZdX2Fm13Nsy/AQD4I3pkfEB0mHv/DO62OxPzbwAA/oog4wNS4gMVE9743JeYCItS4pvegebu/BuGmQAAZkSQ8QFWq2OIpzFThoc3a6Kvp+bfAADgDQQZH5HWK1gz0yMa9MzERFg0Mz2i2fNYPDH/BgAAb2Gyrw9J6xWswclBKio5rbJqu6LDrEqJD7ygS65be/4NAADeRJDxMVarRX26BrXY+ern3zQ2vNTc+TcAAHgbv4b7udacfwMAgLcRZNqA1pp/AwCAtzGe0Ea0xvwbAAC8jSDThrT0/BsAALyNoSUAAGBaBBkAAGBaDC01g91uMNcEAAAfQJBpInaRBgDAdzC01ATsIg0AgG8hyLiJXaQBAPA9BBk3sYs0AAC+hyDjJnaRBgDA9xBk3MQu0gAA+B4+dd1Uv4t0Y9hFGgAAzyLIuIldpAEA8D0EmSZgF2kAAHwL4yBNxC7SAAD4DoJMM7CLNAAAvoGhJQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFoEGQAAYFp+v7KvYRiSpIqKCi9XAgAA3FX/uV3/OX4ufh9kTp48KUlKTEz0ciUAAKCpTp48qaioqHM+bjHOF3VMzm636/Dhw2rfvr0slpbb2LGiokKJiYk6ePCgIiMjW+y8aIj32jN4nz2D99kzeJ89ozXfZ8MwdPLkSSUkJMhqPfdMGL/vkbFarerWrVurnT8yMpL/JB7Ce+0ZvM+ewfvsGbzPntFa73NjPTH1mOwLAABMiyADAABMiyDTTCEhIfrd736nkJAQb5fi93ivPYP32TN4nz2D99kzfOF99vvJvgAAwH/RIwMAAEyLIAMAAEyLIAMAAEyLIAMAAEyLINNML774onr06KF27dpp6NCh2rp1q7dL8isLFizQZZddpvbt26tz586aOHGiCgsLvV2W3/vjH/8oi8Wi7Oxsb5filw4dOqRf/OIXio2NVWhoqAYMGKDt27d7uyy/YrPZ9NBDDyk5OVmhoaHq1auXHn300fPu14PGffbZZ8rIyFBCQoIsFotWrlzp8rhhGHr44YcVHx+v0NBQjRkzRkVFRR6pjSDTDO+++67mzZun3/3udyooKNCgQYOUnp6uY8eOebs0v7FhwwbNmjVLmzdvVl5enurq6jR27FhVVVV5uzS/tW3bNr3yyisaOHCgt0vxS6WlpRo+fLiCgoL0j3/8Q1999ZUWLVqkmJgYb5fmV5588knl5OTohRde0N69e/Xkk09q4cKFev75571dmqlVVVVp0KBBevHFF8/6+MKFC7VkyRK9/PLL2rJli8LDw5Wenq5Tp061fnEGmuzyyy83Zs2a5bxvs9mMhIQEY8GCBV6syr8dO3bMkGRs2LDB26X4pZMnTxopKSlGXl6eMXLkSGPOnDneLsnv3HfffcaIESO8XYbfu/baa43bbrvN5dgNN9xgZGZmeqki/yPJWLFihfO+3W434uLijKeeesp5rKyszAgJCTHeeeedVq+HHpkmqq2t1Y4dOzRmzBjnMavVqjFjxmjTpk1erMy/lZeXS5I6dOjg5Ur806xZs3Tttde6fF+jZX300UcaMmSIJk+erM6dOys1NVVLly71dll+54orrtDatWu1b98+SdKuXbuUn5+v8ePHe7ky/7V//34dOXLE5edHVFSUhg4d6pHPRb/fNLKl/fDDD7LZbOrSpYvL8S5duujrr7/2UlX+zW63Kzs7W8OHD1f//v29XY7fyc3NVUFBgbZt2+btUvzav//9b+Xk5GjevHn6n//5H23btk133323goODlZWV5e3y/Mb999+viooK9e3bVwEBAbLZbHr88ceVmZnp7dL81pEjRyTprJ+L9Y+1JoIMfN6sWbP0xRdfKD8/39ul+J2DBw9qzpw5ysvLU7t27bxdjl+z2+0aMmSInnjiCUlSamqqvvjiC7388ssEmRb03nvv6a233tLbb7+tSy65RDt37lR2drYSEhJ4n/0UQ0tN1LFjRwUEBOjo0aMux48ePaq4uDgvVeW/Zs+erb///e9at26dunXr5u1y/M6OHTt07NgxpaWlKTAwUIGBgdqwYYOWLFmiwMBA2Ww2b5foN+Lj43XxxRe7HOvXr58OHDjgpYr8029+8xvdf//9mjJligYMGKBp06Zp7ty5WrBggbdL81v1n33e+lwkyDRRcHCwLr30Uq1du9Z5zG63a+3atRo2bJgXK/MvhmFo9uzZWrFihf75z38qOTnZ2yX5pauvvlp79uzRzp07nbchQ4YoMzNTO3fuVEBAgLdL9BvDhw9vsITAvn37lJSU5KWK/FN1dbWsVtePtoCAANntdi9V5P+Sk5MVFxfn8rlYUVGhLVu2eORzkaGlZpg3b56ysrI0ZMgQXX755Xr22WdVVVWlW2+91dul+Y1Zs2bp7bff1t/+9je1b9/eOc4aFRWl0NBQL1fnP9q3b99g3lF4eLhiY2OZj9TC5s6dqyuuuEJPPPGEbr75Zm3dulWvvvqqXn31VW+X5lcyMjL0+OOPq3v37rrkkkv0+eefa/Hixbrtttu8XZqpVVZW6ptvvnHe379/v3bu3KkOHTqoe/fuys7O1mOPPaaUlBQlJyfroYceUkJCgiZOnNj6xbX6dVF+6vnnnze6d+9uBAcHG5dffrmxefNmb5fkVySd9bZs2TJvl+b3uPy69Xz88cdG//79jZCQEKNv377Gq6++6u2S/E5FRYUxZ84co3v37ka7du2Mnj17Gr/97W+Nmpoab5dmauvWrTvrz+SsrCzDMByXYD/00ENGly5djJCQEOPqq682CgsLPVKbxTBY7hAAAJgTc2QAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAAIBpEWQAmIrNZtMVV1yhG264weV4eXm5EhMT9dvf/tZLlQHwBrYoAGA6+/bt0+DBg7V06VJlZmZKkqZPn65du3Zp27ZtCg4O9nKFADyFIAPAlJYsWaJHHnlEX375pbZu3arJkydr27ZtGjRokLdLA+BBBBkApmQYhq666ioFBARoz549uuuuu/Tggw96uywAHkaQAWBaX3/9tfr166cBAwaooKBAgYGB3i4JgIcx2ReAab3xxhsKCwvT/v379d1333m7HABeQI8MAFP617/+pZEjR+rTTz/VY489Jklas2aNLBaLlysD4En0yAAwnerqas2YMUMzZ87U6NGj9frrr2vr1q16+eWXvV0aAA+jRwaA6cyZM0f/+7//q127diksLEyS9Morr2j+/Pnas2ePevTo4d0CAXgMQQaAqWzYsEFXX3211q9frxEjRrg8lp6ertOnTzPEBLQhBBkAAGBazJEBAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACmRZABAACm9f8DoXSsRFiFShcAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Plot the data\n",
        "plot_data(X=X, y=y, obs=obs, obs_std=obs_std, inverse_df=inverse_df)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mf6LMlDSOzst"
      },
      "source": [
        "## Validate Expected Result\n",
        "\n",
        "Since we fabricated the experimental data, we can manually validate our computed result. Note this is a sanity check and would not be possible with real-world data of an unknown model.\n",
        "\n",
        "$$\n",
        "\\begin{align}\n",
        "    y = 2X+3 \\;\\;\\;\\rightarrow \\;\\;\\; \\frac{y-3}{2} = X\n",
        "\\end{align}\n",
        "$$\n",
        "\n",
        "The difference between the computed inverse model mean and X in above equation (1) should be less than the determined standard deviation.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "K2PfO5LVOzst",
        "outputId": "ef2e48d9-ac74-4f24-c1a0-490f1920883a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Success! The difference between the inverse model mean X (3.636) and the algebraically computed X (3.625) is 0.011 which is less than the inverse model standard deviation of 0.367.\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/var/folders/ml/754v6mn14dx2yfqknbvqmlmw0000gn/T/ipykernel_28132/3142967819.py:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
            "  inverse_mean = inverse_df[\"mean\"][0]\n",
            "/var/folders/ml/754v6mn14dx2yfqknbvqmlmw0000gn/T/ipykernel_28132/3142967819.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
            "  inverse_std = inverse_df[\"sd\"][0]\n"
          ]
        }
      ],
      "source": [
        "inverse_mean = inverse_df[\"mean\"][0]\n",
        "inverse_std = inverse_df[\"sd\"][0]\n",
        "y_obs = obs[\"y\"][0]\n",
        "\n",
        "actual_X = (y_obs - 3) / 2\n",
        "difference = abs(actual_X - inverse_mean)\n",
        "\n",
        "if difference <= inverse_std:\n",
        "    print(\n",
        "        f\"Success! The difference between the inverse model mean X ({inverse_mean}) and the algebraically computed X ({actual_X}) is {round(difference, 5)} which is less than the inverse model standard deviation of {inverse_std}.\"\n",
        "    )\n",
        "else:\n",
        "    print(\n",
        "        f\"Something went wrong! The difference between the inverse model mean X ({inverse_mean}) and the algebraically computed X ({actual_X}) is {round(difference, 5)} which is greater than the inverse model standard deviation {inverse_std}.\"\n",
        "    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q0HFQm7cOzst"
      },
      "source": [
        "## Finishing up\n",
        "\n",
        "This section covers how to delete your trained model and dataset from the cloud. Note that you don't need to delete your model and data to rerun a campaign or dataset of the same name--if, say, you rerun the functions above as-is, with no name changes, will simply overwrite your existing model and dataset on the Cloud.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "id": "EFwOyDBNOzst"
      },
      "outputs": [],
      "source": [
        "# Delete campaign and dataset if necessary\n",
        "dataset.delete()\n",
        "emulator.delete()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.14"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}