ricardoV94/PyMC_optimization.ipynb

## PyMC_optimization.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1hBg4YmsqvSz"
   },
   "source": [
    "# PyMC Optimization workflow\n",
    "\n",
    "Notebook create for the following discourse thread: https://discourse.pymc.io/t/pymc-for-bayesian-optimization/11293"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "xhmWGgXSe5aC"
   },
   "outputs": [],
   "source": [
    "import pymc as pm\n",
    "import pytensor\n",
    "import pytensor.tensor as pt\n",
    "from pytensor.graph.rewriting.utils import rewrite_graph\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "s0M4UuGngCLg"
   },
   "outputs": [],
   "source": [
    "seed = sum(map(ord, \"PyMC Optimization\"))\n",
    "rng = np.random.default_rng(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "ARMR6-3De97W"
   },
   "outputs": [],
   "source": [
    "with pm.Model() as generative_model:\n",
    "    x = pm.MutableData(\"x\", np.random.normal(size=(100, 5)), dims=[\"batch\", \"features\"])\n",
    "    betas = pm.Normal(\"betas\", shape=5, dims=\"features\")\n",
    "\n",
    "    mu = pm.Deterministic(\"mu\", x @ betas, dims=[\"batch\"])\n",
    "    sigma = pm.HalfNormal(\"sigma\")\n",
    "    y = pm.Normal(\"y\", mu, sigma, shape=mu.shape, dims=[\"batch\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "blf77lysfpmr"
   },
   "outputs": [],
   "source": [
    "# Simulate data for specific parameter values\n",
    "fixed_parameters = {\n",
    "    \"betas\": [-2, -1, 0, 1, 2],\n",
    "    \"sigma\": 0.5,\n",
    "}\n",
    "with pm.do(generative_model, fixed_parameters) as synthetic_model:\n",
    "    synthetic_y = pm.draw(synthetic_model[\"y\"], random_seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 332
    },
    "id": "gfKm-lwpgQtw",
    "outputId": "91261faa-9f44-4290-b636-e73ca8c9ada0"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [betas, sigma]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.\n"
     ]
    },
    {
     "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>betas[0]</th>\n",
       "      <td>-1.997</td>\n",
       "      <td>0.054</td>\n",
       "      <td>-2.095</td>\n",
       "      <td>-1.890</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4625.0</td>\n",
       "      <td>3487.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>betas[1]</th>\n",
       "      <td>-0.990</td>\n",
       "      <td>0.064</td>\n",
       "      <td>-1.108</td>\n",
       "      <td>-0.870</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>6048.0</td>\n",
       "      <td>3008.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>betas[2]</th>\n",
       "      <td>0.038</td>\n",
       "      <td>0.060</td>\n",
       "      <td>-0.071</td>\n",
       "      <td>0.156</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>6143.0</td>\n",
       "      <td>3044.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>betas[3]</th>\n",
       "      <td>0.958</td>\n",
       "      <td>0.052</td>\n",
       "      <td>0.857</td>\n",
       "      <td>1.054</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4936.0</td>\n",
       "      <td>3266.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>betas[4]</th>\n",
       "      <td>2.022</td>\n",
       "      <td>0.055</td>\n",
       "      <td>1.925</td>\n",
       "      <td>2.131</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.000</td>\n",
       "      <td>6202.0</td>\n",
       "      <td>3688.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.536</td>\n",
       "      <td>0.039</td>\n",
       "      <td>0.465</td>\n",
       "      <td>0.609</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.000</td>\n",
       "      <td>5813.0</td>\n",
       "      <td>3339.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "betas[0] -1.997  0.054  -2.095   -1.890      0.001    0.001    4625.0   \n",
       "betas[1] -0.990  0.064  -1.108   -0.870      0.001    0.001    6048.0   \n",
       "betas[2]  0.038  0.060  -0.071    0.156      0.001    0.001    6143.0   \n",
       "betas[3]  0.958  0.052   0.857    1.054      0.001    0.001    4936.0   \n",
       "betas[4]  2.022  0.055   1.925    2.131      0.001    0.000    6202.0   \n",
       "sigma     0.536  0.039   0.465    0.609      0.001    0.000    5813.0   \n",
       "\n",
       "          ess_tail  r_hat  \n",
       "betas[0]    3487.0    1.0  \n",
       "betas[1]    3008.0    1.0  \n",
       "betas[2]    3044.0    1.0  \n",
       "betas[3]    3266.0    1.0  \n",
       "betas[4]    3688.0    1.0  \n",
       "sigma       3339.0    1.0  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Infer parameters\n",
    "with pm.observe(generative_model, {\"y\": synthetic_y}) as inference_model:\n",
    "    idata = pm.sample(random_seed=rng)\n",
    "\n",
    "pm.stats.summary(idata, var_names=[\"betas\", \"sigma\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "fwA9BXCfgc4d"
   },
   "outputs": [],
   "source": [
    "from typing import Sequence\n",
    "from pymc import Model\n",
    "from arviz import InferenceData\n",
    "from pytensor.tensor import TensorVariable\n",
    "\n",
    "def rewrite(*vars: TensorVariable) -> Sequence[TensorVariable]:\n",
    "    return rewrite_graph(vars, include=(\"ShapeOpt\", \"canonicalize\", \"stabilize\", \"specialize\"))\n",
    "\n",
    "def posterior_predictive_fn(\n",
    "    model: Model,\n",
    "    idata: InferenceData,\n",
    "    var_names: Sequence[str],\n",
    "    replace: dict[TensorVariable, TensorVariable],\n",
    ") -> Sequence[TensorVariable]:\n",
    "\n",
    "    # Replace every unobserved RV by a dummy with the same shape (excluding chain, draws)\n",
    "    dummy_posterior_vars = {\n",
    "        model_rv: pt.tensor(model_rv.name, shape=idata.posterior[model_rv.name].shape[2:])\n",
    "        for model_rv in model.free_RVs\n",
    "        if model_rv.name not in var_names\n",
    "    }\n",
    "\n",
    "    predicted_vars = pytensor.clone_replace(\n",
    "        [model[var_name] for var_name in var_names],\n",
    "        replace=dummy_posterior_vars | replace,\n",
    "        rebuild_strict=False,\n",
    "    )\n",
    "\n",
    "    # Clean up shape graph (or else blockwise fails)\n",
    "    predicted_vars = rewrite(*predicted_vars)\n",
    "\n",
    "    # Replace dummy_vars by posterior draws.\n",
    "    # Vectorize graph handles the new batch dims\n",
    "    posterior_vars = {\n",
    "        var: pt.constant(idata.posterior[var.name], name=var.name)\n",
    "        for var in dummy_posterior_vars.values()\n",
    "    }\n",
    "    batch_predicted_vars = pytensor.graph.replace.vectorize_graph(\n",
    "        predicted_vars,\n",
    "        replace=posterior_vars\n",
    "    )\n",
    "\n",
    "    return batch_predicted_vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "-TW5EsJ3TlP5",
    "outputId": "58b44ca5-3e8c-4825-9472-5b1ecad9613d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DropDims{axes=[2, 3]} [id A] <Matrix(float64, shape=(4, 1000))>\n",
      " └─ SpecifyShape [id B] <Tensor4(float64, shape=(4, 1000, 1, 1))>\n",
      "    ├─ Transpose{axes=[1, 2, 0, 3]} [id C] <Tensor4(float64, shape=(?, ?, 1, 1))>\n",
      "    │  └─ Reshape{4} [id D] <Tensor4(float64, shape=(1, ?, ?, 1))>\n",
      "    │     ├─ dot [id E] <Matrix(float64, shape=(1, 4000))>\n",
      "    │     │  ├─ ExpandDims{axis=0} [id F] <Matrix(float64, shape=(1, 5))>\n",
      "    │     │  │  └─ x [id G] <Vector(float64, shape=(5,))>\n",
      "    │     │  └─ [[-2.00337 ... 6598e+00]] [id H] <Matrix(float64, shape=(5, 4000))>\n",
      "    │     └─ [   1    4 ... 1000    1] [id I] <Vector(int64, shape=(4,))>\n",
      "    ├─ 4 [id J] <Scalar(int8, shape=())>\n",
      "    ├─ 1000 [id K] <Scalar(int16, shape=())>\n",
      "    ├─ 1 [id L] <Scalar(int8, shape=())>\n",
      "    └─ 1 [id M] <Scalar(int8, shape=())>\n",
      "DropDims{axis=2} [id N] <Matrix(float64, shape=(4, 1000))>\n",
      " └─ normal_rv{0, (0, 0), floatX, False}.1 [id O] <Tensor3(float64, shape=(4, 1000, 1))>\n",
      "    ├─ RandomGeneratorSharedVariable(<Generator(PCG64) at 0x7F2FEBC302E0>) [id P] <RandomGeneratorType>\n",
      "    ├─ [   4 1000    1] [id Q] <Vector(int64, shape=(3,))>\n",
      "    ├─ 11 [id R] <Scalar(int64, shape=())>\n",
      "    ├─ DropDims{axis=3} [id S] <Tensor3(float64, shape=(4, 1000, 1))>\n",
      "    │  └─ SpecifyShape [id B] <Tensor4(float64, shape=(4, 1000, 1, 1))>\n",
      "    │     └─ ···\n",
      "    └─ [[[0.51133 ... 4335869]]] [id T] <Tensor3(float64, shape=(4, 1000, 1))>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7f30760e1840>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create graphs from x to (y, mu), for every draw from the posterior\n",
    "x = pt.tensor(\"x\", shape=(5,))\n",
    "\n",
    "# Add the batch dim\n",
    "replace_x = x[None]\n",
    "# replace_x = pt.broadcast_to(x, idata.constant_data[\"x\"].shape)\n",
    "\n",
    "predicted_y, predicted_mu = posterior_predictive_fn(\n",
    "    model=generative_model,\n",
    "    idata=idata,\n",
    "    var_names=[\"y\", \"mu\"], \n",
    "    replace={generative_model[\"x\"]: replace_x},\n",
    ")\n",
    "\n",
    "# Drop the batch dim\n",
    "predicted_y, predicted_mu = predicted_y.squeeze(-1), predicted_mu.squeeze(-1)\n",
    "# predicted_y, predicted_mu = predicted_y[..., 0], predicted_mu[..., 0]\n",
    "\n",
    "predicted_y, predicted_mu = rewrite(predicted_y, predicted_mu)\n",
    "pytensor.dprint([predicted_mu, predicted_y], print_type=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Hco5BGG3UjFD",
    "outputId": "de1e46d0-e96c-45e3-fbb6-aebb04f9c490"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.5119952383345927, (4, 1000))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For debugging, we check the graph looks correct and that evaluating returns what we expect\n",
    "x_test = np.array([0, 0.25, 0.5, 0.75, 1.0])\n",
    "res = predicted_mu.eval({x: x_test})\n",
    "res.mean(), res.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "XqoDLzlYjVV_",
    "outputId": "6fdaa625-fec2-4ca3-a443-9d9e9e1b9433"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(24786.39509671)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a cost function\n",
    "# target_y could be a SharedVariable so we can optimize for different values\n",
    "# without having to recompile the cost and grad functions\n",
    "target_y = 5.0\n",
    "\n",
    "cost = pt.sum((target_y - predicted_mu) ** 2)\n",
    "cost_fn = pytensor.function([x], cost)\n",
    "cost_fn(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "8inHy0kZkyiS",
    "outputId": "0df87060-9f1b-4c9e-db61-b333ecd7ae49"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 39743.84888322,  19721.50155962,   -740.94660199, -19045.16624961,\n",
       "       -40223.42767489])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Can also get gradient function for scipy\n",
    "grad_fn = pytensor.function([x], pt.grad(cost, wrt=x))\n",
    "grad_fn(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "lR82IxHukE9q",
    "outputId": "75ef9aaa-03a1-43ff-9cef-65c6b1fb2c65"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 27.73603988253558\n",
       " hess_inv: array([[ 0.02245723, -0.00485788, -0.0006937 ,  0.00593183,  0.01698518],\n",
       "       [-0.00485788,  0.02792761, -0.00248036,  0.00198485,  0.00798208],\n",
       "       [-0.0006937 , -0.00248036,  0.03571901,  0.00078844, -0.00294736],\n",
       "       [ 0.00593183,  0.00198485,  0.00078844,  0.0409628 , -0.01258027],\n",
       "       [ 0.01698518,  0.00798208, -0.00294736, -0.01258027,  0.02670575]])\n",
       "      jac: array([-2.58878825e-06, -1.25715929e-06, -8.10551371e-08,  1.27604829e-06,\n",
       "        2.66765740e-06])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 14\n",
       "      nit: 10\n",
       "     njev: 14\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([-1.07362961, -0.42703831, -0.07130571,  0.55570669,  0.94079541])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Optimize x\n",
    "res = minimize(\n",
    "    fun=cost_fn,\n",
    "    x0=x_test,\n",
    "    jac=grad_fn,\n",
    ")\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_pred_fn = pm.pytensorf.compile_pymc([x], predicted_y, random_seed=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAG7CAYAAABeqpbSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbF0lEQVR4nO3dd3hUddrG8e9MeiMhpFETQif03pQmYEHBjugKFtTXteuqu7oqu7qu3V3XtqsudgVdUEBQpCrSS6ihhx7Se5857x8nhVADJDnJzP25rlxmzpk584TEmXt+1WYYhoGIiIi4LbvVBYiIiIi1FAZERETcnMKAiIiIm1MYEBERcXMKAyIiIm5OYUBERMTNKQyIiIi4OYUBERERN6cwICIi4uYUBkRERNycwoCIiIibUxgQERFxcwoDIiIibk5hQERExM0pDIi4mTvvvBObzcaoUaM41Q7mzzzzDDabja5du1JUVGRBhSJS12zGqV4NRMRl5ebm0r17d/bu3csbb7zBQw89VHFu1apVDB48GA8PD1avXk337t2tK1RE6oxaBkTcTGBgIJ9++ikeHh788Y9/ZOvWrQDk5+fzu9/9DofDwV//+lcFARE3ojAg4oYGDRrE448/TmFhIbfccgvFxcU88sgj7Nq1i4svvpjHHnvM6hJFpA6pm0DETZWUlNC/f382bNjAqFGjWLBgAY0aNWLTpk1ER0dbXZ6I1CGFARE3tm3bNnr37k1hYSEA06ZNY9KkSRZXJSJ1TWFAxI0VFxfTtWtXdu7cSXBwMIcOHSIwMNDqskSkjmnMgIgbe+qpp9i5cyd2u52srCwefvhhq0sSEQsoDIi4qWXLlvH666/j7+/PggULCAkJ4YMPPmD27NlWlyYidUxhQMQNZWdnM2nSJJxOJ6+88gojRozg7bffBsxFiVJSUiyuUETqksKAiBt64IEHSExMZPTo0dx7770ATJw4kRtvvJHk5GTuuusuiysUkbqkAYQibmbmzJlcc801NG7cmC1bttCsWbOKcxkZGXTp0oUjR47w0Ucfcdttt1lYqYjUFYUBETdy7NgxunTpQmpqKl9++SUTJkw46T4//fQTl156KYGBgWzatImYmJi6L1RE6pTCgIiIiJvTmAERERE3pzAgIiLi5hQGRERE3JzCgIiIiJtTGBAREXFzCgMiIiJuTmFARETEzSkMiIiIuDmFARERETenMCAiIuLmFAZERETcnMKAiIiIm1MYEBERcXMKAyIiIm5OYUBERMTNKQyIiIi4OYUBERERN6cwICIi4uYUBkRERNycwoCIiIibUxgQERFxcwoDIiIibs7T6gJEpP45klnAh7/uY/W+dPKKSokND+S63i0YExeJzWazujwRqWE2wzAMq4sQkfrj23WHeHrWFgpKHCedG9U5ktdv6E6Qr5cFlYlIbVEYEJEKX6w6wJ9mbgagb0xjJg9qTeMAL5btTOWjX/dR7HDSuWkjvrxrAMF+CgQirkJhQEQAmB1/hAe+2oBhwJSLWvPHyzpht1d2CWw6lMnt09aSmlvEwNgmfHZnfzzs6jIQcQUaQCgi7DqWwx++iccw4HcDovnT5VWDAEC3FiF8cns/Arw9WLE3jfeW7rGoWhGpaQoDIm7O4TR4dEY8hSVOLmoXxtSr4k47SLBzs0Y8e1UcAG8s2MmWw1l1WaqI1BKFARE39/mq/Ww6lEWQryevXd/9pBaBE11fNqug1GnwxLebcDjV0yjS0CkMiLix5JxCXpm/A4A/jOlARCPfsz7GZrPxwtVdCfLxZOuRbP63/lBtlykitUxhQMSN/ePnXeQUldKtRTA394+u9uPCAn24b0RbAF75cQf5xaW1VaKI1AGFARE3dTA9n+lrDwLw1OWdznlmwKRBMbQM9SM5p4gPf9lXGyWKSB1RGBBxU28v3k2Jw2BI2zD6xzY558f7ennw2OgOAHzw6z5yi9Q6INJQKQyIuKGD6fnMWGf29T88qt15X2dst2bEhgeQVVDCJysSa6g6EalrCgMibuiDX/bicBpc1C6M3tGh530dD7uN+8vGDvxn2V7y1Dog0iApDIi4mYy8YqavNVsF7hna5oKvd2W3ZsQ08Scjv4RvNbNApEFSGBBxM5+u3E9BiYO4Zo0Y1ObcxwqcyNPDzu1DWgMwbXkiTq07INLgKAyIuJGiUkdF3/5dF8fW2HbE1/RqQZCPJ3tT81i2K6VGrikidUdhQMSNzNucRGpuMU2Dfbm8a9Mau26gjyc39G0JwLTfEmvsuiJSNxQGRNzIpyv3A3BTv1Z4edTs//63DozGZoMlO1LYk5Jbo9cWkdqlMCDiJrYdyWbd/gw87TYmlH2Kr0nRTQIY2TECgE/UOiDSoCgMiLiJz1aZrQJjukRVaw+C8zF5kDmQ8Jt1hzTNUKQBURgQcQPZhSXM2nAYgN8NqP4eBOdqcNsmtA4LIK/YwffxR2rteUSkZikMiLiBmesPk1/soF1EIP1bn/8iQ2djs9m4qZ/ZBfHl6gO19jwiUrMUBkRcnGEYFQMHfzcwusamE57Odb1b4u1hZ9OhLLYczqrV5xKRmqEwIOLiVu5NZ3dyLv7eHlzds3mtP19ogDdjukQB8IVaB0QaBIUBERf3WVmrwPiezQny9aqT57ypbLbC7PgjFJY46uQ5ReT8KQyIuLDk7EJ+3JoEwC39a2/g4IkGxDahRWM/cgpLK55fROovhQERF/bl6oOUOg16Rzemc7NGdfa8druNa3u1AGDGWm1eJFLfKQyIuKjiUiefl60tcOvAumsVKHddbzMMLN+TyuHMgjp/fhGpPoUBERf149YkknOKCA/y4bIuNbcPQXW1DPVnQGwohgH/W6fWAZH6TGFAxEWV7044sV8rvD2t+V/9+t7mQMJv1h/CMLS1sUh9pTAg4oK2HsliTaK5D8HE/q0sq+OyrlEE+niyPy2f1fvSLatDRM5MYUDEBX3ymzlW4LKuTYmspX0IqsPf25MryrZK/kZdBSL1lsKAiIvJzC9m1kZzH4JJFgwcPNF1fcyBhHM3H9XmRSL1lMKAiIv5es1BikqddG7aiN7Rja0uhz7RjWkdFkB+sYMfNh+1uhwROQWFAREXUupwVuxDMHlQTK3vQ1AdNputYprhDHUViNRLCgMiLmTu5qMcyiggNMCbq3o0s7qcCtf0ao7NBqv3pbM/Lc/qckTkBAoDIi7CMAzeW7oXMFsFfL08LK6oUtNgP4a0DQPgW7UOiNQ7CgMiLmLpzhS2H83G39vDkhUHz+b6PuaaA9+uP4zTqTUHROoThQERF/He0j0A3NSvFSH+3hZXc7LRnSMJ8vXkcGYBK/amWV2OiBxHYUDEBWw4kMHKvel42m3cMaS11eWckq+XB1d1N8cxzFh70OJqROR4CgMiLuCtRbsBGNejOc1C/Cyu5vTKuwrmb00iu7DkrPefNm0aNpvttF8TJkyo7ZJF3IKn1QWIyIVZk5jOooRkPOw2fj+8jdXlnFH3FsG0jQhkd3Iuczcd5aZ+1VsquXv37vTo0eOk4/3796/hCkXck8KASANmGAYvz08A4IY+LYgND7S4ojOz2Wxc37sFL85LYMbag9UOA+PHj+e5556r3eJE3Ji6CUQasCU7UliTmIGPp50HRrazupxqubpnczzsNtYfyGRPSq7V5YgICgMidSoxMRGbzcawYcPIy8vjkUceoWXLlvj5+dGrVy9mz55dcd8ZM2bQr18/AgICiIyM5IEHHqCgoKDivNNp8PKPO3AWF9Ay8QdGD+mHv78/jRo1YujQocyaNeuUNcydO5fbb7+dTp060ahRIwICAujevTt/+9vfKCoqOun+5f32zz33HAcOHGDixImEh4fj5+dHnz59qtRcHRGNfBnaPhzQ5kUi9YXN0CbjInUmMTGR1q1bM3DgQJxOJ3v27GHAgAHk5uaybNkybDYb8+fPZ/PmzTz++OP07duXyMhIfvnlF9LS0pg4cSKff/45YL6RPjxtCSlfP0VRygGaN29O7969yc/PZ8WKFeTl5fHiiy/y5JNPVqkhKiqKvLw84uLiaNWqFdnZ2axevZqMjAxGjBjBTz/9hIdH5YJF06ZN47bbbmPSpEnMmzcPX19fevXqxbFjx1ixYgV2u5158+YxevToav87zNt8lP/7fD2RjXz47cmReNhPvWxy+XOPHTuWjh07kp2dTVRUFCNGjGDo0KHn8RsQkVMyRKTO7Nu3zwAMwBg2bJiRnp5ece6///2vARht27Y1QkNDjWXLllWcO3z4sBEREWEAxp49e4ysgmKj918XGL6xvQ3AePzxx43i4uKK++/Zs8do06aN4eHhYcTHx1epYebMmUZubm6VY9nZ2cbYsWMNwPj444+rnCuvCzDuv/9+o6SkpOLcm2++aQDGRRdddNLPGh0dXfG46n7t27fvtM994tfQoUONpKSk6v/ji8hpqWVApA6Vtwx4eHiQkJBA27ZtK845nU6ioqJISUnhmWeeYerUqVUe+8gjj/DGG2/w3//+l8Ph/Xj7m4UcnfYAAwcOYvnyX0/alOi7775j/Pjx3H///fzzn/88a227d++mXbt2XHPNNXz77bcVx8s/ncfGxpKQkICXl1fFudLSUiIiIsjNzSU3Nxdv78rFjh577DFSU1NP+3wbD2ayOzmXFo39GBDbBIBXX32VsLCwivv8+OOPrFy5knHjxhEbG0tBQQGrV6/m8ccfJyEhgd69e7Nq1aoqLRkicu40m0DEAjExMVWCAIDdbic6OpqUlBRGjRp10mPatDGnDW7csY85OyMoSNwAwPjx4065O+GQIUMAWLNmzUnndu3axQ8//MDu3bvJy8vD6XRS/rlg165dp6x52LBhVYIAgKenJ7Gxsaxbt460tDSaNm1ace7VV1897c8PsOVwFmPf+hWnh503nxp5ylUTx4wZw5gxYypuN2rUiCuvvJLhw4fTu3dv1q1bx9dff83EiRPP+FwicmYKAyIWaN68+SmPBwQEnPZ8+bk5Gw5Q2qMfzT1zyQSeeOIJnnjiidM+1/Gfzg3D4LHHHuONN97gdI2COTk5pzzeokWLUx4PDDSnM55q8OGZxDVrRMeoIBKScpgdf4TfDYyp9mMDAwN54IEHuO+++/jxxx8VBkQukMKAiAVO9Um+uudTc4toE+BN11bBbAUuuugiYmNjT3v/45vdv/76a15//XVatGjBm2++ycCBAwkPD8fLy4vi4mJ8fHxOGxLOVvOJztZNAJBzLJfUQ5n8abEXCztFntRNcCbt2plTKY8ePXpOdYnIyRQGRBqIg+n5Fd8/P74L675bBcB1113HAw88UK1rzJw5E4B3332XsWPHVjm3d+/eGqrU9M0337B///5q3TcP+Hg1PPfcc9UOAxkZGUBly4SInD+tMyDSABSVOvhq9QEA2kcGcVnXplxyySUAp11P4FTK30Bbtmx50rnp06dfeKHHSUxMxDCMs35N+XgN0U/M4fk5W4mJian29csHOfbu3btG6xZxRwoDIg3Amz/v4liO2Sc/vEMEAAMGDGDkyJEsXryYhx9+mNzcqqv5OZ1OfvrpJ3799deKY+3btwfg3//+d5XugF9++YVXXnmltn+MU7qutzkWYeaGw5Q4nFXO/fOf/zzp5yopKWHq1KnMmDEDPz8/Jk+eXFelirgshQGRem7jwUzeX7qn4rafd+U0us8//5xu3brx5ptvEh0dzciRI5kwYQIXXXQRUVFRjBkzhrVr11bc/4EHHiAgIIB33nmHLl26cNNNN3HxxRczdOhQ7rnnnjr9ucoN7xhBWKA3qbnFLEpIrnLuwQcfJCIigj59+nDttddyxRVXEBMTw3PPPYevry+fffbZaQdjikj1KQyI1GOFJQ4enb4RpwG9WoWcdD4yMpKVK1fy+uuv065dO9asWcOsWbM4dOgQPXv25O233+aWW26puH/79u1Zs2YNV155JampqXz//ffk5uby/vvvW9Yy4OVh59qy1oEvy7pCyj3zzDMMGTKE5ORk5s2bx6JFi/D39+fuu+9m48aNXHPNNVaULOJytOiQSD324g/beX/ZXsKDfFjw8MWnnIvvCval5jH81SXYbPDrEyNoHuJndUkibkUtAyL11Lr9GfznF3OE/9+u7uqyQQCgdVgAA2ObYBgwfc1Bq8sRcTsKAyL1UGGJgz/MiMdpwDW9mjOqc6TVJdW6Cf3MGQ7T1x7E4VSDpUhdUhgQqYde+2kHe1PziGzkw7Nj46wup06MiYsixN+Lo1mFLN2ZfPYHiEiNURgQqWfWH8jgg1/3AfDiNV0J9vc6yyNcg6+XB9f2Kh9IqK4CkbqkMCBSjxSWOHj8m00YZd0DIzq6fvfA8W4q6ypYlJDMsexCi6sRcR8KAyL1yFuLdrE7OZfwIB+eGdvZ6nLqXNuIIPrGNMbhNJixVq0DInVFYUCknthyOIv3lpqzB/46rotLzx44kwl9WwHw9dqDODWQUKROKAyI1ANOp8FTMzfjcBpc0a0pl3aJsroky1zetSlBvp4cTC9g+Z4z73ooIjVDYUCkHvh67UHiD2UR5OPJs1e6X/fA8fy8Pbimp7nE8FcaSChSJxQGRCyWkVfMS/MTAHh4VHsignwtrsh6E/qZXQU/bUsiNbfI4mpEXJ/CgIjFXvlpB5n5JXSMCuLWgdFWl1MvdGraiB4tQyhxGMxYe8jqckRcnsKAiIW2HM6q2JznL+O64Omh/yXL3dzfbB34bOV+rUgoUsv0yiNioZfmJ2AYcFX3ZvRrHWp1OfXKld2bEeLvxeHMAhYnaEVCkdqkMCBikeW7U/llVypeHjb+MKaD1eXUO75eHtzYx1yE6OMVidYWI+LiFAZELOB0Gvx9njlo8Ob+0bQM9be4ovrplgHR2Gzwy65U9qbkWl2OiMtSGBCxwA9bjrL5cBYB3h7cN6Kt1eXUWy1D/RnRIQKAT1fut7gaEdelMCBSxxxOg9d/2gnAlItjCQv0sbii+u13ZTMsvll3iPziUourEXFNCgMidWzu5qPsTc0jxN+LOy+Ktbqceu/iduHENPEnp7CUWRuOWF2OiEtSGBCpQ06nwTuLdwNw26DWBPp4WlzRCQwDSovML6N+TOez223cMsBsHfhkRSJGPalLxJUoDIjUoYUJySQk5RDo48nkQTFWl3MyRzHMutf8chRbXU2F63u3xM/Lg4SkHFbsSbO6HBGXozAgUkcMw+BfZa0CvxsYTbC/l8UVNRzB/l5c17sFAB/+us/iakRcj8KASB1ZvjuN+IOZ+HrZuWNIa6vLaXBuGxwDmK0rmmYoUrMUBkTqyL9/2QvAhL6tNIPgPMSGBzKyoznN8L/LE60tRsTFKAyI1IHdybks25mCzQa3D1arwPkqb1H5Zt0hMvPrz5gGkYZOYUCkDnz8WyIAl3SKpFUTrTZ4vga2aULHqCAKShx8ufqg1eWIuAyFAZFallVQwrfrzW14b6uPMwgaEJvNVtE68PFviZQ4nBZXJOIaFAZEatn0NQfJL3bQITKIgW2aWF1Og3dVj2aEBfqQlF3ID5uPWl2OiEtQGBCpRQ6nUbHj3m2DY7DZbNYW5AJ8PD34XdkiRB/9uk+LEInUAIUBkVq0OCGZQxkFhPh7Ma5Hc6vLcRk3D2iFt6ed+ENZrN2fYXU5Ig2ewoBILfpy9QEAru/dAj9vD4urcR1hgT5c28sMV+8t2WNxNSINn8KASC05klnA4h3JAEzo18rialzPlItisdnKl3jOtrockQZNYUCklkxfexCnAf1bh9ImPNDqclxObHggl3dpCqh1QORCKQyI1AKH0+DrNeY8+In91SpQW/5vWBsAZm86ysH0fIurEWm4FAZEasHSnckczSokxN+LMXFRVpfjsro0D+aidmE4nAbvL1PrgMj5UhgQqQVfrDJbBa7t1QJfLw0crE33DmsLwPS1h0jJKbK4GpGGSWFApIal5BRVDBy8qV9Li6txfQNiQ+nRMoTiUicfLdf2xiLnQ2FApIZ9t/EwDqdB95YhtI0Isrocl2ez2bi3bOzApyv2k5VfYnFFIg2PwoBIDft2/WEAruulRYbqyiWdImkfGUhuUWnFio8iUn0KAyI1aNuRbLYfzcbLw8bYbs2sLsdt2O02fj/cHDvw0fJ95BaVWlyRSMOiMCBSg/5XtjvhyI6RNA7wtrga9zK2WzNahwWQmV/C5yv3W12OSIOiMCBSQ0odTmZtPALAtb1bWFyN+/Gw2yrWHfjPL/soLHFYXJFIw6EwIFJDlu1KITW3iCYB3gzrEG51OW7p6p7NaR7iR2puUcWiTyJydgoDIjXk23XmwMGrejTDy0P/a1nBy8POPWWtA+8t3UNxqdPiikQaBr1iidSArPwSFmw7BpgLDYl1ru/dgoggH45mFVaM4RCRM1MYEKkBP2w5SrHDSYfIIOKaNbK6HLfm6+XBXRfHAvDOkj2UOtQ6IHI2CgMiNWB2vDlwcFzPZthsNourkYn9WxEa4M2B9HxmbzpidTki9Z7CgMgFSskpYuXeNADGdtXaAvWBv7cndwxpDcC/Fu3G6TQsrkikflMYELlA87YcxWlA9xbBtGrib3U5UubWgdE08vVkT0oe87cmWV2OSL2mMCBygebEHwXQioP1TJCvF5MHm60Dby3ajWGodUDkdBQGRC5AUlYha/anA3BFt6YWVyMnum1QDAHeHmw/ms2ihGSryxGptxQGRC7A3M1HMQzoHd2YZiF+VpcjJ2gc4M0tA6MBtQ6InInCgMgFmFs2Un2sWgXqrTuHxOLjaWfjwUyW706zuhyReklhQOQ8Hc4sYP2BTGw2uLyrwkB9FR7kw039WgHwr8W7LK5GpH6yGWo3k/pk0wz4351Vj417B3refPJ9DQO2fAsbPoOj8VCUAwFh0GoADPg9tOx77s//3ytg/6/VuuuRkN4MSnqUAbGhfHXXQNi7FJa+DEc2gOGAiM4w5CHoPO7kB5cWw7sDIW03jH8Xekw891prSnV/5j8eBp/A2q/neItfhKV/P/v9bvoaOlx68vHCbFjxL0q2zqYoZS8AHk1a49f1Khh4H/iexwJRW2fC7p/Nv7ncZCjIAJsHBIRD027Q9TroPB5OXG8ifR8s/AvsXWL+rYa0Mn/vQx4Gu8fJzzPzHoj/EnrcDOPfOfc6Rc6Bp9UFiFTIT4f5T1bvvo4SmDEZEuZUPZ5z1Hyx3joLRv0FBj9Q01VWSMstBspmEexZDJ9dA4YTvAPBwxeOrIfpt8L496DHTVUfvOItMwi06AvdbzrF1eWCpe2Bj6+C7EN4AV7l783p22HpdtjwOUz6Hpq0Obfr/vYvOLz25ONZB8yvhDnQ6Uq44dPKQJBzDD4cBXkpYPcE/yaQvgcW/dUMCePfrnqtg6sh/ivwCYZLnjvHH1zk3CkMSP0x/4+Qn1q9+y6cWjUIBEZBZJz5aS0/FTBgwZ8hohO0G1X9GmIGg3/oqc+l7oKU7RU3Vxa2xG6DS7tEwfS7zSDQuDXcvRQ8feHjK+HgKlj8QtUwkHUYlr0GNjtc/srJnyCt1LSH+YYF0KynWWM5u8UvF8GtoFmPU58Liqx6u7QIvrgRssv3JrCRG9WXrUey6Wvbgd1mmOe+uBH+7zfw9D7HYmwQGgvBLcx/l6TNkHfcbIXts2Hr/6DLtebtNf8x/11tHnDXEojqCj8+BSv+BRs/h4sfg1BzGiROJ8x9FDBg2JMQGHGOtYmcO4UBqR/2LIJNX5nfB7eErDNsP1uUA6v+XXk7NBbuWmo2+ealwTv9K9/QfvrzuYWB4X86/bmPLqv41omdTx2j6BsTSligDxzZaJ7ocDn4Bpvfx11jhoGsg5CXanZhAPz0FJTkQe/bzDfcuuAohT0L4dAaGPH06e/X53bYt8z8fvw74OlTt89/JjFD4Op3q3ffddMg7bjxAZe9TGD/u3jtvRV0PPglf/H62Dyetsu8b/+7ql/H0MehaXcIiqo85iiBmXeb3Vbl9q+oDAPlfx+RcWYQALNFaMW/AAOObqwMA+s+gqRNEN4J+p1DXSIXQAMIxXrF+TD7IfP7kGgY/OCZ739oDTiKKm/HXVPZ9xvQBDpeUXkuZbvZWnChjm2FA79V3Nzg05cDRiRj4qLO8KBTDMfZt8zsxvBrDCOfufC6ziZpM8z/E7zeCb64AXbMr/3nrA/PH/9V5fc+jaDPbQDccVFrvnCMJJfjpoHGf3lu124/pmoQAPDwgi7XVT3m5XuWC53i7yM/HRY9b35/+cvgoc9rUjf0lybWW/wCZO43v7/yTcg+eub7550wPcwvpOpt3xNuH1pjfpK7EKver3LzHznDARjTpexNoVlPMyzs+AGGPQGefuabPpgtHQFh5qfjH/5gHhvx9Om7Iy5UbjJsmm6+yR3bUvWc91mWS971IyRtNQdA/vgUtOwPncaC1zmsoXAhz38mydvM5vOCDPMNPjLO7Js/8Y25pND8ZF0uorP5Zg1c0imSZqGNSMhtSR/7TvN80ibzMWd98z4DRyls+abqsdhhld836wm7F5ihMmmz2TpQEVhsZvcMwM/PmT9f3NXQ+uLzr0fkHCkMiLWObICVZU2/3SdCmxHmwK4zOXEEePres9zed2E1FmTC5hkVN7MDYlhW2JVuLYJpXr7Q0NDHzQGEGfvg9c5mP3JhpnmuvFl81XuQkgBR3aD37RdW04lKCs0gEv8l7F5ovpmX8wqAjpdD1+uhzcgzXydhbuX3a/5jfgWEm4Mg211S+89/Jkc3ml/H+/FPMOyPcNEjlccy94OztPL2ceMJPOw2Jg2KIeXH4MrzzlLzMeEdzq2e+X8yu4CKcsw3+IrxLjZz4Grb4/69+t4J6/5rdl/9e5g5gDD3mHmu581mF8GRDbDhU/Pfa/QL51aLyAVSGBDrOErh+/vNN46AcBhTzRfAFn3Bwxsc5mh+Nk2HjmPNKYV7FsOOeVXvX5RzYXVu+AxK8itufud9BWCr2kXQZjj8blbl1MLSQmjWq3JqYc4xWPoSYIPLXwW7HRJ+MAeZ5SRBYKT5abDT2HOr7cAqiP/CbIUozKo8bvcyg1XX681ukwv5RJ6XAl/dBLfNhxa96/75z8RRbA4m9fKDAf9nHju+DjBbaY5zXa8WLPnphLEQhdnn/tx7FlUZUAqYAy6H/wkGndDVFRQJdywwa927xOwOCI01pw0OfsicJjv3MXMQ6sWPQnBzM2Cs/cicdeIdaLY09J5cc+M4RI6jMCDW+e2f5gsewGUvVb/Z3D8U+t8Nv71l3i7ONT+Vn86FvHgaBqz9sPKmTxCvJZtviCeNF4gdan6dyoI/Q1G2OWisVX9zYONv/6x6ny3fmHPfqxuKEn+FaceNj8BmBqKu15njKKr77xnWFqIHmQMtGzUzP/EWZJgtG+XN7Y5icxrcrbNq/vnPJCDMHNTY8QoI62DeTt8Hy14xg1S5xX8z3yhP2Z1RtW8+2N+LFo39IOsUd71QhtPs89+1ACZOr9qFFdoarp926set/8ScrhgaCwPvN1topk8CZ0nlfXb8YA5QnDRbgUBqnAYQijXy081P0QDtL6scdV1dI58zuxVOJfCEaWb+Tc65vAq7FlTpdtjbfByZDl/aRgTSNqKaC/AcWAmbvjb7uS+ZCofXVwaBPnfAE4nmGx6Yo8sPr6vedU9cL6zbDTD2DbNJ+lzeiK/8B4x4Clr2M1toPLzM6WwTvqj6qTrxV3OxpJp+/jPpN8W8ZttLIKSl+WYf2Rmu/bDqTIyibHNsCJj/zscrKTjpsrEhJyzycz6LD/1+JTyXBY/vg1u/M1usyh1cBUuqsVgSmN1QP081v7/0JTNQfP+AGQSa9YTHdsPV71de94TxKyI1QWFArFGUA6VlL9KJv8BLrSu/5j1e9b7zHjeP//pm5TEPT3Oa2Z2LYMgjZnN0v7th3Ntw1b+qPr58Ktf5WH3cFEZsfFJqTlO89IyzCI7jdMAPj5nfD33CbC4+vhtjyEPmzILBD1Ueq+6I+0bNqr4hbvoa3hkA7w6GX9+AzDNMz6yOwAgIb19521kCBel19/xnYrdD9OCqx3LL5vk3jqm6JkJO0kkPD3FUDkJ12jzNWSznyz/UbMK/eYa5jkC548dfnMniF8zxBu0vg/ajzTf88vEH/e6CwHDoPsFcZwFO7gYTqQHqJhDrFedW7/wpPuHRovfJ/djf3Vf5vaffyW8a1ZW+11x2towjdgTTd/kCjrNMKTzO2o/MrpDwjtD/HvNY7nFvTkFlexo0alZ5rHxg2dk0aWMuYJOcYPbbb5pursB4bIv59fPUymb7zleb0y5P5Cg9/fQ1wwnZR6oe8z6uNaQmnv9MnA6zD/50izJlHqh62yfI/K+XrxkAj2wwbydvN1s0yhcWKi3GllzZ17/T1poOnj5c8NJPPsHg5Q/FZWNUjl+E6HSStsCaD81Fqi590Tx2/O//+JkSjZqaKxxW9+9D5ByoZUAaptTdleMNyjmd5syEDZ9WHusx8eSph290heeCza//XsFprf6A4/ub45vdSEGJg+YhfnRpXo1m5by0yjnjl71U+abre9xI9vyyT9r5x02XPP58dUR0NJdefngr3PKtOd/d08+s/cAKczrea+3hs+vMlfGOt3k6fHWzuUDO8Qwn/Pxs5eJNAJFdTr03wYU8/5lkHYL3L4Yt/6vaPQHmseNXoLR5QPM+lbe7Taj8vjjHHMlfbu1HVQLoV0UDWbs/o/L8zP+r/Pt47oTfRcJcc5GiEwcpOh3w62uVQQCq19rwwx/MAbSDHqhcdOhUfx9Q+Tdyrn8fItWglgGxRuNos7/1VDZ8Dt/dW3n7VBsVHV4HM+8y5/CHtjY/QSYnVP3UHdzy/Bf2Kc6HjZ9V3g6N5auMDsBhRnWOxFadJYR/ftYchNd5XNU55zEXVQ5+3PgZXPQobPziuPNDzq9mu4fZt972EnN0/NaZ5lS/AyvM6XO7F5hN5p2urHyMYZhvqglzICDCXL45JcF8systrHr946fv1dTzn03SJvjmNnO6XdPuZt9+2p6qqwsC9J5UteWh92Sziyd9j3l73hPmfhVg1lMmxbsFXxSOJHPlfvrGVGOcQ8Z++PGP5sj/iE7QqLnZ3XXi3155TWcS/7W5NkVwS3OzonIt+oKHj7mwVvxX5mDMw2vNWQVgLpktUsMUBqRhyzp46qWLm7Qz+3BPbBWork1fV/n05+w7hcWLzH7cSzpFnu5RlQ6vM6ckevmfPGe87SgzECT+Yu5it/K9yiblmIug3ejzq/l4vo3MN6Pek8zujvivqq7KV+74UJOXDPtO0bRts8Pwp85tkGd1n/9Mjq+tJK/KCpBVdBwLl54wWM/LFyZ+bW5UlHMEs5XihMcHNSN9zCcUf5bED1uSmJpfQrC/V/Vqc5aYQeX4xY2O13sy9P+/0z++KAcWlAXVMS9UnXrpH2qGg6V/NwPUK23MAZJgDo4deH/1ahQ5BwoD0jC16GMOrtr/m9lPXZhtvgGFd4K48dBr0nlsPnOcNR9Ufu8dyLbIsaTkbCbA24N+rc/yCbJ8zjiG+Wk6pGXV83a7+Ua1+G/mp+fcY+YnzLirzTfdmt64KDTWnPs+7I/mp/7jdb3enG2xcz4c3WR+cs85YtbQONZspeh7x4UNwjzT859JSCuYsgi2zzFnCmQfNtdrcJaagxub9zKnana47NSPD2tnjvj/7S2zeT+jbJXLxtHmVMVB99PBN5iOUctISMph7uajTOzf6sw1dbjMXHPiwArzk3p+OhTnmWMpQlqZ22b3uNn8+zyTJX83WxJih516i+vhfzR/xtX/qVxnoM1wcwfDEzdlEqkBNsM4cX6QiJzojQU7+cfCXVwaF8V7v+t99gc0VKVFMKusi6YmNyqqx95fuocX5yXQN6YxM+4ZZHU5IpbQAEKRaliUYDafj+ik7WRdzbgezbHZYE1iBgfT88/+ABEXpDAgchbJ2YVsPmyOHxjeQWHA1UQF+zK4jbm99MwNhy2uRsQaCgMiZ7F4h9kq0L1lCOFBrt9s7o6u7tkcgP+tP4R6TsUdKQyInMXC7WVdBGoVcFmXdonCz8uDxLR8NhzMtLockTqnMCByBoUlDn7dbU4pHKnxAi4rwMeTMXHmKP2Z69VVIO5HYUDkDFbtSye/2EFkIx/imp3HZjbSYFzdqwUAszcdobjUaXE1InVLYUDkDBaXzyLoGFG9VQelwRrcpgnhQT5k5pewZEc19hUQcSEKAyKnYRgGCxPMTWFGdNRCL67O08POuO7mhlHfxx85y71FXIvCgMhp7E7O5WB6Ad6edga3Pccd96RBurIsDCzcnkx+canF1YjUHYUBkdNYWNZFMDC2Cf7eWrnbHXRrEUzLUD8KShwsTkg5+wNEXITCgMhpLCqbUqhZBO7DZrNxRVezdWDuZnUViPtQGBA5hcz8YtYdMPe416qD7mVst6aAuQR1XpG6CsQ9KAyInMLSnSk4nAYdIoNoGep/9geIy4hr1ojoJv4UljgruopEXJ3CgMgplG9MNLyjWgXcjc1mq2gdmLtJXQXiHhQGRE5Q6nCyZIc5eEzjBdxT+biBxTtSyFVXgbgBhQGRE2w4mElWQQkh/l70bBlidTligU5Ng4gNC6C41MnC7cesLkek1ikMiJygfGOiYe3D8fTQ/yLu6PiugtnxRy2uRqT26ZVO5ASLylYd1HgB93ZFN7OrYNnOFLILSyyuRqR2KQyIHOdgej47j+XiYbcxtH241eWIhdpHBtI2IpBih5Oft6mrQFybwoDIccpnEfSObkyIv7fF1YiVqs4qUFeBuDaFAZHjlIeBkeoiECoXIFq2K4WsAnUViOtSGBApk1dUyoo9aYCmFIqpbUQQ7SMDKXEYLFBXgbgwhQGRMst3p1LscNIy1I824YFWlyP1xOVdzdaBHzarq0Bcl8KASJnKLoJIbDabxdVIfXFFWRj4RV0F4sIUBkQAwzAqwsAIjReQ47SLrOwq0KwCcVUKAyLA1iPZJOcU4e/tQf/YUKvLkXpGXQXi6hQGRKhcdfCidmH4eHpYXI3UN5VdBalagEhcksKACJWrDqqLQE6lXWQQ7bQAkbgwhQFxe8k5hcQfygJgeAeFATm18q4CLUAkrkhhQNxe+XbF3VoEE9HI1+JqpL66opu6CsR1KQyI21u0XbMI5OzaRwZprwJxWQoD4taKSh38sstsGRjZMdLiaqS+06wCcVUKA+LWVu1NJ6/YQXiQD3HNGlldjtRz5bMKlu1UV4G4FoUBcWsLt5vNvSM7RmC3a9VBObP2kYG0CQ+g2OGs+NsRcQUKA+K2DMPg57LxAiM7qYtAzs5ms1W0DszdlGRxNSI1R2FA3NaOYzkczizAx9POkLZhVpcjDcQV3ZoBsGxnCjnqKhAXoTAgbqt81cHBbcPw89aqg1I9x3cV/KyuAnERCgPitspfyEd20pRCqT51FYgrUhgQt5SaW8TGg5mAphTKubu8bAGiZbvUVSCuQWFA3NLihGQMA+KaNSIqWKsOyrnpEBlEbHgAxaXOiu4mkYZMYUDc0kLNIpALUKWrQAsQiQtQGBC3c/yqg5dovICcp/LVCJfuTNECRNLgKQyI21lZtupgRJAPXZoFW12ONFAdo8r2Kih18tNWzSqQhk1hQNzOwuNmEWjVQTlfNpuNcd3NNQe+23jY4mpELozCgLgVwzAqxwtoFoFcoCvLwsDy3amk5BRZXI3I+VMYELdy/KqDg7XqoFygmLAAurcMwWloJ0Np2BQGxK1o1UGpaeoqEFegMCBuRasOSk0b260pdhusP5DJgbR8q8sROS8KA+I2tOqg1IaIRr4MbNMEgNmbjlhcjcj5URgQt7Fou7nqYJfmWnVQata47s0B+H6jwoA0TAoD4jbmbzU3lRndOcriSsTVjOkShbeHnR3HckhIyra6HJFzpjAgbiG3qJRfd6UCMCZOYUBqVrCfF8M6hAPwnVoHpAFSGBC3sGRHMsUOJzFN/GkfGWh1OeKCxvWo7CowDMPiakTOjcKAuIUfy5aLHRMXhc2mVQel5o3sFEGAtweHMwtYtz/D6nJEzonCgLi8olIHixPM9QXGdFEXgdQOXy+Pii4odRVIQ6MwIC7vtz1p5BaVEhHkQ48WIVaXIy7sqh7mAkRzNx+lxOG0uBqR6lMYEJf345ayWQRxkdqYSGrVkLZhhAV6k55XXLFNtkhDoDAgLs3hNFiwrXK8gEht8vSwM7ab2Towa4O6CqThUBgQl7ZufwZpecU08vVkQGwTq8sRNzCurKtgwbZj5BWVWlyNSPUoDIhL+7FsoaGRnSLx8tCfu9S+Hi1DiG7iT0GJg5+2JVldjki16NVRXJZhGMwvGy8wJk57EUjdsNlsFWsOqKtAGgqFAXFZGw9mcjizAH9vD4a21y6FUnfGl3UV/Lo7ldTcIourETk7hQFxWXM2HQXgkk6R+Hl7WFyNuJPY8EC6tQjG4TSYE6/WAan/FAbEJTmdBnPLwsAV3ZpaXI24o4quAi1AJA2AwoC4pPUHMkjKLiTIx5Oh7cOtLkfc0JXdm2K3md1Vial5VpcjckYKA+KSyrsIRnWOxNdLXQRS9yKCfBncNgzQ8sRS/ykMiMtxOA3mbjbDwNju6iIQ65R3FXy38bB2MpR6TWFAXM7qfemk5BTRyNeTIW3VRSDWGRMXiY+nnb2peWw5nG11OSKnpTAgLmfuZrNJdkxcFN6e+hMX6wT5enFJZ3ONi1kbD1tcjcjp6ZVSXEqpw8m8zeZCQ2O7N7O4GhEYX9ZVMDv+CA6nugqkflIYEJeycm86aXnFNPb3YlAb7UUg1hvaPpwQfy+Sc4pYsSfN6nJETklhQFxKeVPspV2aai8CqRe8Pe1c0dUcyKquAqmv9GopLqOg2MG8slkE1/RqbnE1IpXG9zT/HudvSaKwxGFxNSInUxgQl7Fg+zHyih20aOxH71aNrS5HpELvVo1pHuJHblEpC7cnW12OyEkUBsRlzNpgNsFe3bM5drvN4mpEKtntNsaVbV6krgKpjxQGxCWk5haxdGcKULnQi0h9Ut5VsGRHMpn5xRZXI1KVwoC4hDll07a6tQimbUSg1eWInKR9ZBCdmjaixGEwb0uS1eWIVKEwIC5h5nFdBCL11ZVly2OX76gpUl8oDEiDtycll/hDWXjYbVyphYakHhvb1fz7/G1PKqm5RRZXI1JJYUAavPKBgxe3CyMs0MfiakROr1UTf7q1CMZpmNMMReoLhQFp0JxOo6KLYLy6CKQBKF+ASF0FUp8oDEiDtnJfGocyCgjy8WR05yiryxE5qyu6mWFg1b40knMKLa5GxKQwIA3ajLWHALiyRzP8vD0srkbk7Fo09qdHyxB1FUi9ojAgDVZ2YQk/lC0/fEOflhZXI1J9Y8taB+bEq6tA6geFAWmwZscfoajUSfvIQLq3CLa6HJFqu7xs3MCa/ekkZamrQKynMCAN1vSyLoIb+rTEZtPyw9JwNAvxo090YwwD5m1R64BYT2FAGqQdSTnEH8zE027TLAJpkMoHEs7RrAKpBxQGpEGasfYgACM7RWhtAWmQLu/aFJsN1u3P4EhmgdXliJtTGJAGp8ThrFhbQAMHpaGKbORL35hQgIqBsCJWURiQBmfh9mTS8ooJD/JhaPtwq8uplpUrVzJu3DjCwsLw9fWlffv2PP300+Tn51f7Gpdccgk2mw2bzUZS0slT0goLC/n9739PWFgYAQEBXHXVVezfv/+U18rKyiIqKoqbbrrpnH+WxMREbDYbMTExZ7zf5MmTsdlsTJs27ZTHy7/sdjvBwcHExMRw5ZVX8vLLL3Ps2LFzvm5DNFZdBVJPKAxIg/PF6gMAXNurBZ4e9f9P+PPPP2fIkCF8//33xMTEcPnll1NYWMgLL7zAoEGDyMnJOes1pk2bxsKFC884UPLBBx/knXfeITo6mosuuog5c+Zw+eWX43A4TrrvM888Q15eHq+++uoF/WwXYvDgwUyaNIlbb72V0aNH06JFCxYuXMgTTzxBq1ateOmllzAMw7L66sKlXaKw22DjwUwOplc/GIrUtPr/SipynANp+fyyKwWAif1aWVzN2R06dIg777wTh8PBRx99xNq1a/nf//7Hrl27uP7664mPj+fxxx8/4zVSUlJ47LHHGD16NK1anfpnPnr0KB999BGXXXYZa9euZf78+fz1r39l27ZtzJw5s8p9t2zZwjvvvMOf//xnmje3bvDlnXfeybRp05g2bRozZszg119/JS0tjX/+8594enry5JNP8tRTT1lWX12ICPKlf+smgLoKxFoKA9KgfLnmAIYBF7ULo1UTf6vLOatp06ZRWFjIqFGjuO222yqO+/j48Pbbb+Pv78+HH35IWlraaa/x0EMPkZeXxzvvvHPa+2zZsoXS0lJuvfXWitaD22+/HYCNGzdWue99991HmzZtePjhhy/gJ6sdfn5+3H///cydOxcPDw9efPFF4uPjrS6rVpXPKlAYECspDEiDUVzqrJhFcHP/aIurqZ5169YBMGzYsJPOhYeH07lzZ0pKSvjhhx9O+fgff/yRL774gqeeeoo2bdqc9nkyMjIAaNy4ccWx8u/T09Mrjn3xxRcsXbqUt956Cy8vr3P+eerKsGHDKsYzvPXWWxZXU7vKuwriD2Wpq0AsozAgDcZP25JIzS0mIsiHkZ0irC6nWvLy8oCqb9LHCw01R5Of6tNvfn4+99xzDx07djxrV0J598GuXbsqju3cuROA6GgzOOXm5vKHP/yBa6+9llGjRp3jT1L3JkyYAMDixYstrqR2hQX6qKtALKcwIA3GF6vMgYMT+rbEqwEMHATz0z9w2lH95ccTExNPOvfnP/+ZxMRE3n33Xby9vc/4PD169KBp06a8/vrrbNmyhWPHjvH4449js9m47LLLAPjLX/5CZmYmr7/++gX8RHWnR48eAOzdu5fi4mJri6lll6urQCzWMF5Rxe3tTcnltz1p2G1wYwMYOFhu6NChAHz55ZcnvaGtXLmSHTt2AJw0o2D9+vX84x//YNKkSafsYjiRr68vr7zyComJiXTt2pWoqCh+/PFH7rnnHrp168aOHTt48803+dOf/lRlEGJBQcF5j9jfv39/lSmCJ359/PHH53XdcmFhYRXfl3eDuKpL49RVINbytLoAker4smw64fAOETQP8bO4muq7+eabeeGFFzhw4ADjxo3j1VdfpVWrVixfvpwpU6bg6elJaWkpdntlLnc4HEyZMoWQkJBzmvp38803Exsby4wZMygsLGTEiBFce+21ANx///20atWKxx57DICvvvqKJ598kv379xMcHMx9993HX/7yl3P6dBAQEMB111132vO//vore/bsOYcrVnV8SHH1vSfCg3zo1zqUlXvTmbflKHddfPrxISK1QWFA6r3CEgcz1pmbEk3s33BaBcB8w5wzZw5jx45l/vz5zJ8/v+Jcq1ateOSRR3j55ZerjCl48803Wb9+PR9++GGVT8fVMXDgQAYOHFjl2LfffsuCBQuYM2cOPj4+rFu3jokTJzJmzBj+8Y9/sHTpUl544QUiIiJ44N67q/1cYWFhZ1z4Z/LkyRcUBlJTUyu+P92YC1dyRdemrNybztzNSQoDUucUBqTem7XhMJn5JbRo7MewDg1j4ODxunbtSkJCAjNmzGDt2rWUlpbSvXt3Jk6cyPPPPw9AXFxcxf1nz55d0cz+ySefVLlW+cqD11xzDd7e3jz//PMMGTLktM9dUFDAo48+ypVXXskVV1wBwGuvvUZgYCDTp08nKCiIcePGsX79el555ZVzCgO1rXxKZLt27er1zIeaMqZLFM98v5X4g5kcysinReP6P3VWXIfCgNRrhmEw7bdEACYNjMHD3jCbi/38/Lj11lu59dZbqxz/+eefgZOnHhqGwbJly057vRUrVgBVPz2fyt/+9jeOHTvGm2++WXEsISGBjh07EhQUVHGsX79+LF26lOzsbBpV5weqA1999RUAw4cPt7iSuhER5Eu/mFBW7Utn3uYkplwca3VJ4kY0gFDqtZV700lIysHPy8PlNiVaunQp69evJy4ujsGDB1ccX7JkCYZhnPKrfJrg0aNHMQyD8ePHn/b6e/bs4ZVXXuHxxx8nNrbqG8uJeyKUT4GsL33zS5Ys4auvvsJms3H//fdbXU6dKV+AaK5mFUgdUxiQeu2/y/cBcG3v5gT7N8ym4o0bN1JaWlrl2Pr165k4cSI2m63WFtV58MEHadq0KU8++WSV43FxcWzbto0NGzYA5kyG2bNn06pVqyqtBVYoLCzkX//6F1dccQUOh4M///nPdOnSxdKa6tKlXaKwle1VcFjbGksdUjeB1FsH0/P5ebu5e93kQTHWFnMBHnroIbZt20aPHj0ICwsjMTGRVatWYbfbef/992ulGXzu3LnMnTuXmTNn4udXdfbFH/7wB7744guGDx/OiBEj2LBhAwcPHuS9996r8TrO5IMPPmDJkiWA2VKRlJTEunXryM/Px8fHh5dffrli9oO7iAgytzVevS+deZuPcudF6iqQuqEwIPXWJysScZbtQ9A2wtpPrBfilltu4bPPPmPjxo1kZmYSHh7OhAkT+MMf/lCxsE5NKioq4sEHH2TMmDGn7Ebo1q0bs2bN4umnn2bOnDlERUXx97//nbvvvhtKi2q8ntNZvnw5y5cvx2azERgYSGhoKMOHD2fo0KFMmjSJiIiGN1i0JlzRtSmr96UzV2FA6pDNcPU9QqVByisqZcCLC8kpLOWjyX0Y0THS6pLcQ2kRzLrX/H78O+DpY209big5u5D+Ly7EMOC3J0fQrAGtqyENl8YMSL30vw2HySksJaaJP8Pau+cnRHFPEY186Rtt7lmh5YmlrigMSL3jcBp8+MtewBwrYG+g0wlFztflXaMAhQGpOwoDUu8s2JZEYlo+If5e3NDXtaYTilTHZV2bYrPB+gOZHNGsAqkDCgNSrxiGwfvLzFaB3w2Ixt9bY1zF/UQ28qVPtLkE87wtSRZXI+5AYUDqlXX7M9hwIBNvTzu3DoyxuhwRy1zeVdsaS91RGJB6pbxV4NpezQkP0kh2cV+XdTHDwLr9GRzNUleB1C6FAak3difnVCwypPnV4u6igo/rKtisrgKpXQoDUm+8s3gPhgGjO0fSJjzQ6nJELKeuAqkrCgNSL+xPy+O7+CMA3D+incXViNQPl5VNMVy7P4OkrEKLqxFXpjAg9cK7S/bgcBoM6xBO1xbBVpcjUi80Dfajd8WsArUOSO1RGBDLHc4s4Nv1hwC4f0Rbi6sRqV+uKOsq+L6s5UykNigMiOX+vXQPJQ6DgbFN6F22DKuImMZ2b4rdBhsOZLI/Lc/qcsRFKQyIpZJzCvlyzUEA7h+pVgGRE0UE+TK4bRgAszaodUBqh8KAWOo/y/ZSXOqkd3RjBsY2sbockXppfI/mAHy38TDaaFZqg8KAWOZYdiGfrNgPwH0j2mKzaUMikVMZ0yUKXy87e1Pz2HQoy+pyxAUpDIhl3l68m6KyVoFh7cOtLkek3gr08WRUZ3Oa4ayNhy2uRlyRwoBY4lBGPl+uPgDAo6Pbq1VA5CzG92gGwOz4o5Q6nBZXI65GYUAs8dbC3ZQ4DAa3bcKgNmFWlyNS713cPpzG/l6k5haxfE+a1eWIi1EYkDq3LzWPb8rWFXhkVAeLqxFpGLw87IztZrYOfLdBXQVSsxQGpM794+edOJwGIzpGVKyuJiJnN76nGQbmb00iv7jU4mrElSgMSJ1KSMqu2IPgkVHtLa5GpGHp1aox0U38yS928IN2MpQapDAgdepvPyRgGOYSq12aaw8CkXNhs9m4oU9LAL5ec8DiasSVKAxInVm2M4VlO1Pw8rDx+KUaKyByPq7r3QK7DdYkZrA7OdfqcsRFKAxInXA4Df72w3YAfjcghugmARZXJNIwRTbyZUTHCACmrz1ocTXiKhQGpE58u/4QCUk5NPL11M6EIhfoxr6tAPh23SGKS7XmgFw4hQGpdQXFDl77aQdgLjvcOMDb4opEGrbhHcKJCPIhLa+YhduPWV2OuACFAal1H/yyl2PZRbRo7MetA2OsLkekwfP0sHNd7xYAfLFaAwnlwikMSK1KySnivaV7AHj80o74enlYXJGIa5jQtxU2G/yyK5U9KRpIKBdGYUBq1esLdpJX7KB7i2Cu7NbU6nJEXEarJv6MLBtI+MlvidYWIw2ewoDUms2HsviqbC7002M7azMikRo2eVBrAL5Zd4icwhKLq5GGTGFAaoVhGDz7/RYMA8b1aEbfmFCrSxJxOYPbNqFtRCB5xQ6+WXfI6nKkAVMYkFoxc8Nh1h/IxN/bgz9e1snqckRcks1mY9KgGAA+/i0Rp9OwtiBpsBQGpMblFJbw4rwEwJxKGBXsa3FFIq7rmp7NCfL1JDEtnyU7k60uRxoohQGpcf9atJuUnCJimvhzx5DWVpcj4tICfDyZ0Nfcr+C9JXstrkYaKoUBqVF7UnL5aPk+AJ65sjM+nppKKFLb7hgSi7eHndWJ6axJTLe6HGmAFAakxhiGwdTZ2yhxGIzoGMGIjpFWlyTiFqKCfbm2d3MA3lm82+JqpCFSGJAa8/P2ZJbtTMHbw84zYztbXY6IW7n74jbYbbB4Rwpbj2RZXY40MAoDUiMKSxz8dc42AO64qDUxYdqVUKQuxYQFMLZbMwDeWbLH4mqkoVEYkBrxwS97OZCeT2QjH+4brl0JRazwf8PaAPDD5qPsSMqxuBppSBQG5IIdySzg7cXmJ5E/Xd6JAB9PiysScU+dmjbi0rgoDANe+XGH1eVIA6IwIBfshR+2U1DioF9MKFd1b2Z1OSJu7bExHbDb4Oftx1irmQVSTQoDckFW7Elj7qaj2G3w7FXaf0DEam0jArmhj7nuwEvzEzAMrUooZ6cwIOet1OFk6uytAEzs34q4ZsEWVyQiAA9e0g4fTztrEjNYvEOrEsrZKQzIeft81QESknII8ffi0VEdrC5HRMo0DfZjctmeBS/M3U5xqdPagqTeUxiQ85KWW8RrP5kDlB4d3YHGAd4WVyQix7t3eFuaBHizJyWvYlVQkdNRGJDz8upPO8guLKVz00ZM7NfK6nJE5ATBfl788XJzx9B//LyLI5kFFlck9ZnCgJyzzYey+GrNQQCmjovDw65BgyL10bW9mtM3pjEFJQ6en7vN6nKkHlMYkHPidBo8+/0WDAPG9WhG35hQq0sSkdOw2Wz8ZVwXPOw2fticxOIEDSaUU1MYkHMyc8Nh1h/IxN/bgz9e1snqckTkLDo1bcTtg2MAePJ/m8jKL7G2IKmXFAak2nIKS/j7/AQA7hvRlqhgX4srEpHqeHR0B2LDAjiWXcRzZdOBRY6nMCDV9tai3aTkFNE6LIA7hrS2uhwRqSZfLw9evaE7dpvZujd/S5LVJUk9ozAg1bIjKYePfjWnJ/15bCd8PD0srkhEzkWvVo25e6i5kdGfZm4mKavQ4oqkPlEYkLMyDIM/z9pCqdNgVOdIRnSMtLokETkPD13Sjk5NG5GeV8x9X6ynxKHFiMSkMCBn9c26Q6xOTMfPy4PnroqzuhwROU8+nh68e3Mvgnw8Wbs/g5fLxgCJKAzIGWXkFfPiPPMF48FL2tE8xM/iikTkQsSEBfDK9d0B+M8v+5i3+ajFFUl9oDAgZ/TyjztIzyumfWSgBg2KuIhLu0Rx18WxADwyPZ7Nh7IsrkispjAgp7X+QAZfrj4AwPPju+LloT8XEVfx+JgODG0fTkGJgzs+XqPlit2cXt3llEodTp6auQWAa3u1oF9rrTQo4ko8Pez8a2JPOkQGkZxTxO3T1pBTqAWJ3JXCgJzSxyv2s/1oNsF+Xvzp8o5WlyMitSDI14sPJ/chLNCHhKQc7vh4LYUlDqvLEgsoDMhJkrIKeb1se+InLu1Ik0AfiysSkdrSorE/027rS5CPJ6v3pXPv55py6I4UBuQkf52zjbxiBz1bhTChb0uryxGRWtaleTAfTOqDj6edRQnJPDo9HofTsLosqUMKA1LFz9uOMXfzUew2eH58F+zanljELfSPbcJ7t/TG027j+/gjPD1rM04FArehMCAVsgtLeHqWOWhwykWxxDULtrgiEalLwztG8MaNPbDZ4MvVB/nDN5vUQuAmFAakwos/JJCUXUhME38eHtXe6nJExAJXdm/Gmzf2wMNu49v1h3jwqw0aQ+AGFAbcWHp6OhEREdhsNqLbtKtYU+Cla7vh66WNiETcVWBaApEr/8mht27h7Vv7ExIWxaWXXcb3339vdWlSSxQG3NgjjzxCamoqACk5RQDcMqAV/WObWFmWiFjoySef5JJLLmHd8sV07RpHYIfBOAIjWLh4Kf+b9Z3V5UktURhwUwsXLuTjjz9mypQpAJQ4DJoF+/LEpVpTQMRdvfvuu7z00kv07duX3bt3s2Hlr8z/7htaT36Vpr//lKSWw8ktKrW6TKkFCgNuqKCggHvuuYfOnTsz/PrbK47/7ZquBPl6WViZiFglMzOTJ554gqCgIL777jtatjSnFQ9uG8Ynd/SjUWAgCUWNufXDVWQVaKVCV6Mw4IamTp3Knj17eOWNt3jpp90ABPt5MqxDhMWViYhVvvzyS3Jycrjpppto2rRplXN9Y0L5/M7+BPt5sf5AJjd/sJKMvGKLKpXa4Gl1AVK3Nm3axGuvvcZtt93Gj+mNSclJBKBJgFYZFHFnCxcuBGDUqFEcO3aMzz//nJ07dxIUFMTgwYO58sor+XLKAH734Sq2HM5mwr9X8umd/YgI8rW4cqkJahlwI06nkylTphASEsJFEx/iu41HsNvMRYVsWltIxK1t3boVgP3799OhQwceffRR3n//fV599VWuvvpq+vbtS7CRw9d3DyAiyIcdx3K48f2VHNZuhy5BYcCNvPXWW6xevZrH/vxXXlpyGIDJg2KsLUpE6oWMjAzAnE3Qtm1bVq5cSXZ2NitWrKBnz55s2LCB6667jjbhgcy4ZyDNQ/zYl5rH9e/+xt6UXIurlwulMOAmDh48yNNPP81FF1/MIkdnCkocDGrThEkKAyICOBzmboV+fn7Mnz+f/v37ExQUxIABA5g/fz4BAQGsXLmShQsXEt0kgG/+byCx4QEcySrkhvdXsP1otsU/gVwIhQE3ce+991JcXEy7qx9ix7EcwgJ9eHOCucqYiEhQUBAAV111FWFhYVXORUREcMUVVwCwZMkSAJoG+zH97oF0btqI1Nxibnx/BesPZNRpzVJzNIDQTcyZMwf/wEZ8/vozADRu2ogbfvaisLAQgAMHDjBs2LCK+wYGBlpVqohYICYmhn379hEdHX3a8wDJyckVx8ICffjyrgHcPm0N6/ZncMsHq/jg1j4Maht2ymtI/aUw4Ebyc7Mh19yIaOPBqucKCgpYunQpAKWlWlRExN307NmTxYsXk56efsrzaWlpACd9UAj28+LTO/px96fr+GVXKpOnreHtib0Y1Tmy1muWmqNuAjewNyWXrs/OJ/qJOTz89QacTieGYWAYBvv27QOgQ4cOFcdCQkKsLVhE6txVV10FwNKlS3E6q25M5HA4+OWXXwDo1avXSY/19/bkg0l9GBMXSXGpk3s+W8d3Gw/XftFSYxQGXNyx7EJu/Wg12YWl9GwVwt+u7opN8whF5ARDhw5l4MCBbN++neeff77KualTp7Jz504iIiK4+uqrT/l4H08P3p7Yi2t6NsfhNHjo6418vmp/XZQuNUDdBC4sK7+EWz9czaGMAlqHBfCfW/toN0IROa1PP/2UQYMG8eyzz/LVV1/RuXNntm7dSkJCAn5+fnz++ecEBASc9vGeHnZevb47gb6efLJiP0/N3EJOYSn3DG1Thz+FnA+1DLionMISbv94DTuO5RDZyIdPbu9HWKBWGRSR02vTpg3x8fHcfffd5OTk8P3335ORkcFNN93EmjVruOSSS856DbvdxtSr4rh3mBkA/j4vgVd/3IFhGLVdvlwAm6HfkMtJzytm0ker2Xw4i0a+nsy4ZxAdooKsLksagtIimHWv+f34d8BTAVLO37tL9vDS/ATAXODsmbGdsWs6c72klgEXk1S2AMjmw1mEBnjzxZQBCgIiYon/G9aGv47vgs0G035L5LFv4ilxOM/+QKlzCgMuZMOBDMa9/Su7k3NpGuzL9LsH0qV5sNVliYgb+92AaN64wVzg7H/rD3Pbf9eQXagtkOsbhQEX8fWaA9z4/kqOZRfRNsJcO7xthBYOEhHrje/ZnA9u7YO/twe/7k7lhvdWcEQbHNUrCgMNXHpeMfd/uYEnvt1MscPJmLhIZv1+MC0a+1tdmohIheEdI5h+90DCg3xISMrh6neWs/lQltVlSRmFgQbKMAzmbT7K6DeWMjv+CB52G4+Nbs+7N/cm0EczRkWk/unSPJiZ9w6iXUQgx7KLuPa93/h6zQGryxI0m6BBij+Yyd/nJbBir7k8aLuIQF69vjvdW4ZYW5g0fJpNIHUgq6CER6dv5Oft5j4HE/q25Lmr4rQOioUUBhoIp9Ng6a4U/r10b0UI8Pa0c9dFsdw3oq3+J5KaoTAgdcTpNHhnyW5eW7ATw4C4Zo1448YetI/U7CcrKAzUY4ZhsONYDrPjjzBrwxEOlw248bTbuKp7Mx4e1Z6WoRobIDVIYUDq2LKdKTz41QYy8kvw9rDzyOj2TLkoVtur1zGFgXrE4TTYfjSbVfvSWb0vjTWJGaTnFVecD/Lx5Ia+Lbl9SGuah/hZWKm4LIUBscCx7EKe/HYTi3ekAGYrwZ/HdmZAbBOLK3MfCgMWKnE42Xw4i9X70lm9L501ienkFFbdPtjXy87gNmFc3as5l3SKVHeA1C6FAbGIYRjMWHeIv87ZVvE6eGlcFA+NakfHqEYWV+f6FAbqUFGpgw0HMive/Nftz6CgxFHlPoE+nvSJaUy/1qH0b92Ers2D8fbUpA+pIwoDYrG03CJeX7CTL1cfwFn27jSiYwRTLoqlf+tQLWdcSxQGapHTabDtaDa/7k5l+e5U1iSmU1hSdSnOEH8v+sWEVrz5d2oahKeH3vzFIgoDUk8kJGXz1qLd/LD5KOXvUi1D/bimZwsu79qU9pGB2o69BikM1LDCEge/7kplwbZjLEw4RmpucZXzYYE+DIgNpX/rUPq1bkK7iEAlXak/FAakntmXmsd/ftnL9xuPkFtU2Y3aPMSPoR3C6RPdmJ6tGhPTxF/h4AIoDNSAzPxift6ezE9bk/hlV2qVpv9AH08GxIYyqE0YQ9qF0S6i4afZlJQUq0uQ2uIoImjhkwDkjPw7eCgMuIrw8HCrS7ggBcUOftqWxKwNh/ltTxpFpVVbWRv7e9GzVWN6tAyhW4tgurcIoXGAt0XVNjwKA+cpKauQn7Yl8ePWJFbuTcfhrPxnbBbsy6jOkYzqHEX/2FC8XKzZv6GHGTk9bw/47zhzpspt3xVQ7DjLA6TBcKWX+oJiByv2prJ8dxobDmSw5XA2xafYDbFlqB/dW4TQvYUZELo0DyZAK7SeksJANRmGwb7UPH7ceowftyax8WBmlfMdo4IYHRfF6M6RxDVr5NJvmK78s7k7hQHX5cov9UWlDrYfzWH9/gw2Hcpk06Es9qbmnXQ/uw3aRgTSrUUI3VuG0LNlCB2jNE4LFAbOKL+4lBV70liyI4WlO1M4kJ5f5Xzv6MaMiYtkTFwU0U0CLKqy7ikMuC6FAdflbi/1WQUlbD6URfyhzIqAcDSr8KT7+Xl50K1FMD1bNaZnqxB6tgohIsjXgoqtpTBwnKz8EtYdSGdNYgZrE9OJP5hVpenJy8PGgNgmjClrAYho5H5/MKAw4MoUBlyXXuohObuQ+ENZbDqUycaD5teJa7uAOTixZ6sQepUFhM7NGuHj6dprvLhtGMgqKCHhaDbbj2az/WgOGw9msuNYzkn3a9HYj2EdwhnaPoJBbZqovwkNIHRpGkDoshr6AMLa4HQa7EnJZcOBTDYczGDDAfN94MR3RW8PO3HNG9GjZQjtI4NoHRZAbFgA4UE+LvPhyKXDgMNpcCy7kP1p+RxMz2d/eh47knLZfjS7Yp3/E8WGBdAnpjF9YkLpFxNKtKariDvR1EJxczmFJWw6lMWGAxllISGzyrLwxwv08aR5iB/hQT6EBXoTHuRDaIAPgT4e+Hl7EuDtgZ+3B14edmw28LDZsNtt2G1mC6sNMKAifDQN9qWZRUvNWxoGvl13iG/XH8LDbsPTbsOj7MvTbsfTw4aXhx2viv+ax7yP+76k1KCw1EFhiYOiUic5haVk5BWTlldMel4R6XnFlDhO/+M1D/GjU9MgOjVtRFyzYPrENCYsUC9+4sYUBkSqMAyDA+n5bDiQSfyhTPam5LEvNY9DGfk4a/jd86FL2vHQJe1r9qLVZGmb94H0fH7bk1arz+Fpt9GisR8tQ/2JbuJPm/BAOjVtRKeoRgT7e9Xqc4uISMNms9mIbhJAdJMAxvdsXnG8qNTBwfR8jmYVkpJTREpOEam5RaTlFpNf7CC/xEFBcSn5xQ5KHQZOo/wLnIZRZTq6zQY2bAT7WfeeZGnLQEJSNjuScnA4DUqdBs6y/zqcBiUOJ6VOg5JSJyVlt0tKzWPFZd97edrx9fTAx8v8b4CPB6EB3oQGeNMkwIfQQG8ig3w0bUSkugwDHGVNoh7e5quUiLg8lx4zICIiImenj8wiIiJuTmFARETEzSkMiIiIuDmFARERETenMCAiIuLmFAZERETcnMKAiIiIm1MYEBERcXMKAyIiIm5OYUBERMTNKQyIiIi4OYUBERERN6cwICIi4uYUBkRERNycwoCIiIib86zOnQzDoLi4uLZrERERkVrg7e2NzWY77flqhYHi4mL+/ve/11hRIiIiUneefPJJfHx8TnveZhiGcbaL1HTLQFJSEtOmTWPy5MlERUXV2HXFGvp9uhb9Pl2Pfqeu5Xx+nzXSMmCz2c6YKM6Vt7d3xX9r8rpiDf0+XYt+n65Hv1PXUhu/Tw0gFBERcXOWhIHAwECGDh1KYGCgFU8vNUy/T9ei36fr0e/UtdTG77NaYwZERETEdambQERExM0pDIiIiLg5hQERERE3pzAgIiLi5hQGRERE3FydhIHMzEweeOABBg4cSFRUFD4+PjRv3pwRI0bw7bffogkNDd/LL7+MzWbDZrOxcuVKq8uR8xATE1PxOzzx65577rG6PDlPM2fOZNSoUTRp0gQ/Pz9at27NTTfdxMGDB60uTc7BtGnTTvv/Z/nXyJEjz/v61VqB8EKlpqby0UcfMWDAAMaPH09oaCjJycnMnj2b6667jilTpvDvf/+7LkqRWrB9+3aeeeYZAgICyMvLs7ocuQDBwcE89NBDJx3v06dP3RcjF8QwDO655x7+/e9/06ZNGyZMmEBQUBBHjhxh6dKl7N+/n5YtW1pdplRTjx49ePbZZ0957ptvvmHr1q2MGTPmvK9fJ+sMOBwODMPA07Nq9sjJyWHAgAFs27aNLVu2EBcXV9ulSA1zOBwMHDgQm81G+/bt+eyzz1ixYgUDBgywujQ5RzExMQAkJiZaWofUjH/+8588+OCD/P73v+cf//gHHh4eVc6Xlpae9JosDU9xcTHNmjUjKyuLQ4cOERkZeV7XqZNuAg8Pj1P+0QUFBVUkmd27d9dFKVLDXnrpJeLj4/noo49OerEREWsUFBQwdepUYmNjefPNN0/5/6aCgGuYOXMmaWlpjB079ryDANRRN8HpFBYWsmjRImw2G507d7ayFDkPW7ZsYerUqTz99NNq1XERRUVFfPzxxxw+fJjGjRszaNAgunfvbnVZco4WLFhAeno6kydPxuFw8P3337Nz505CQkK45JJLaNu2rdUlSg358MMPAbjzzjsv6Dp1GgYyMzN58803cTqdJCcn88MPP3Dw4EGeffZZ2rVrV5elyAUqLS1l8uTJdOrUiSeffNLqcqSGJCUlMXny5CrHLr30Uj799FPCwsKsKUrO2dq1awHz03/37t3ZsWNHxTm73c7DDz/Mq6++alV5UkP279/PwoULad68OZdeeukFXavOw8DUqVMrbnt5efHKK6/w6KOP1mUZUgP+9re/ER8fz6pVq/Dy8rK6HKkBt99+O0OHDiUuLg4fHx+2bdvG1KlTmTdvHldddRXLly8/437oUn8kJycD8Nprr9GrVy9Wr15Np06d2LBhA3fddRevvfYabdq04f/+7/8srlQuxH//+1+cTie33XbbhXfTGhYoLS019u3bZ7z44ouGt7e3cfXVVxslJSVWlCLnYePGjYaXl5fx5JNPVjk+adIkAzBWrFhhUWVS0xwOhzFkyBADMObMmWN1OVJNU6ZMMQDDz8/POHz4cJVzW7ZsMex2u9GmTRuLqpOa4HA4jFatWhk2m83Yu3fvBV/PkkWHPDw8iImJ4cknn+T5559n5syZ/Oc//7GiFDkPkyZNok2bNjz33HNWlyK1zG63c9tttwGwfPlyi6uR6goODgbMKaHNmjWrci4uLo7Y2Fj27NlDZmamBdVJTViwYAEHDhxgxIgRtG7d+oKvZ/kKhKNHjwZgyZIl1hYi1RYfH09CQgK+vr5VFrz4+OOPASqmGs6aNcvaQqVGlI8VyM/Pt7gSqa4OHToAEBIScsrz5ccLCgrqqCKpaTU1cLCc5XNLjhw5AmiaS0Nyxx13nPL4smXL2LVrF1dddRXh4eEV89alYVu1ahWAfp8NyPDhwwFzQbATlZSUsHv3bgICAggPD6/r0qQGpKWl8d133xEaGsrVV19dMxetga6Ls9qwYYORmZl50vG0tDSjR48eBmB8+umndVGK1CKNGWi4tm7damRkZJx0/JdffjF8fX0NHx8fY//+/XVfmJy30aNHG4Dxn//8p8rxv/zlLwZg3HLLLRZVJhfqjTfeMADjgQceqLFr1snH8WnTpvHBBx8wfPhwoqOjCQgIYP/+/cydO5fc3FyuvfZaJk6cWBeliMgpTJ8+nZdffpmRI0cSExODj48PW7Zs4aeffsJut/Pee+/RqlUrq8uUc/DOO+8waNAgpkyZwqxZs+jYsSMbNmxg0aJFREdH88orr1hdopynmu4igDrqJrjuuuvIyspi5cqVLFu2jPz8fEJDQxkyZAi33norEyZM0JQlEQsNHz6c7du3s379epYuXUphYSGRkZHceOONPPzww/Tr18/qEuUctWnThrVr1/LMM88wf/58fvrpJ6Kiovj973/PM888Q0REhNUlynlYvXo1W7ZsoV+/fnTt2rXGrlsnexOIiIhI/WX5bAIRERGxlsKAiIiIm1MYEBERcXMKAyIiIm5OYUBERMTNKQyIiIi4OYUBERERN6cwICIi4uYUBkRERNycwoCIiIibUxgQERFxcwoDIiIibu7/ARJ9yb29N7jvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.plots.plot_posterior(\n",
    "    posterior_pred_fn(res.x),\n",
    "    ref_val=target_y,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RQ87MQE2o13F"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "pymc",
   "language": "python",
   "name": "pymc"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}