A iPython notebook showing how to use SVI for logistic regression in Pyro

pyro_logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.special as ssp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.distributions.constraints as constraints\n",
    "\n",
    "from torch.utils.data import DataLoader\n",
    "from torch.utils.data.sampler import SubsetRandomSampler\n",
    "\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "\n",
    "from pyro.infer import SVI, Trace_ELBO\n",
    "from pyro.optim import Adam, SGD\n",
    "\n",
    "pyro.enable_validation(True)\n",
    "torch.set_default_dtype(torch.double) # this was necessary on the CPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_logistic_dataset(N, p=1, noise_std=0.01):\n",
    "    X = np.random.randn(N, p)\n",
    "    \n",
    "    w = np.random.randn(p)\n",
    "    w += 2 * np.sign(w)\n",
    "\n",
    "    y = np.round(ssp.expit(np.matmul(X, w) \n",
    "                           + np.repeat(1, N) \n",
    "                           + np.random.normal(0, noise_std, size=N)))\n",
    "    y = y.reshape(N, 1)\n",
    "    return X, y, w\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these were adapted from the Pyro VAE tutorial\n",
    "\n",
    "def train(svi, train_loader, n_train):\n",
    "    # initialize loss accumulator\n",
    "    epoch_loss = 0.\n",
    "    # do a training epoch over each mini-batch x returned\n",
    "    # by the data loader\n",
    "    for _, xs in enumerate(train_loader):\n",
    "        # do ELBO gradient and accumulate loss\n",
    "        epoch_loss += svi.step(*xs)\n",
    "\n",
    "    # return epoch loss\n",
    "    total_epoch_loss_train = epoch_loss / n_train\n",
    "    return total_epoch_loss_train\n",
    "\n",
    "\n",
    "def evaluate(svi, test_loader, n_test):\n",
    "    # initialize loss accumulator\n",
    "    test_loss = 0.\n",
    "    # compute the loss over the entire test set\n",
    "    for _, xs in enumerate(test_loader):\n",
    "        # compute ELBO estimate and accumulate loss\n",
    "        test_loss += svi.evaluate_loss(*xs)\n",
    "\n",
    "    total_epoch_loss_test = test_loss / n_test\n",
    "    return total_epoch_loss_test\n",
    "\n",
    "\n",
    "def plot_llk(train_elbo, test_elbo, test_int):\n",
    "    plt.figure(figsize=(8, 6))\n",
    "\n",
    "    x = np.arange(len(train_elbo))\n",
    "\n",
    "    plt.plot(x, train_elbo, marker='o', label='Train ELBO')\n",
    "    plt.plot(x[::test_int], test_elbo, marker='o', label='Test ELBO')\n",
    "    plt.xlabel('Training Epoch')\n",
    "    plt.legend()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogRegressionModel(nn.Module):\n",
    "    def __init__(self, p):\n",
    "        super(LogRegressionModel, self).__init__()\n",
    "        \n",
    "        self.p = p\n",
    "\n",
    "        # hyperparameters for normal priors\n",
    "        self.alpha_h_loc = torch.zeros(1, p)\n",
    "        self.alpha_h_scale = 10.0 * torch.ones(1, p)\n",
    "        self.beta_h_loc = torch.zeros(1)\n",
    "        self.beta_h_scale = 10.0 * torch.ones(1)\n",
    "        \n",
    "        # initial values of variational parameters\n",
    "        self.alpha_0 = np.zeros((1, p))\n",
    "        self.alpha_0_scale = np.ones((1, p))\n",
    "        self.beta_0 = np.zeros((1,))\n",
    "        self.beta_0_scale = np.ones((1,))\n",
    "\n",
    "    def model(self, x, y):\n",
    "        # sample from prior\n",
    "        a = pyro.sample(\n",
    "            \"weight\", dist.Normal(self.alpha_h_loc, self.alpha_h_scale, validate_args=True).independent(1)\n",
    "        )\n",
    "        b = pyro.sample(\n",
    "            \"bias\", dist.Normal(self.beta_h_loc, self.beta_h_scale, validate_args=True).independent(1)\n",
    "        )\n",
    "\n",
    "        with pyro.iarange(\"data\", x.size(0)):\n",
    "            model_logits = (torch.matmul(x, a.permute(1, 0)) + b).squeeze()\n",
    "            \n",
    "            pyro.sample(\n",
    "                \"obs\", \n",
    "                dist.Bernoulli(logits=model_logits, validate_args=True),\n",
    "                obs=y.squeeze()\n",
    "            )\n",
    "            \n",
    "    def guide(self, x, y):\n",
    "        # register variational parameters with pyro\n",
    "        alpha_loc = pyro.param(\"alpha_loc\", torch.tensor(self.alpha_0))\n",
    "        alpha_scale = pyro.param(\"alpha_scale\", torch.tensor(self.alpha_0_scale),\n",
    "                                 constraint=constraints.positive)\n",
    "        beta_loc = pyro.param(\"beta_loc\", torch.tensor(self.beta_0))\n",
    "        beta_scale = pyro.param(\"beta_scale\", torch.tensor(self.beta_0_scale),\n",
    "                                constraint=constraints.positive)\n",
    "\n",
    "        pyro.sample(\n",
    "            \"weight\", dist.Normal(alpha_loc, alpha_scale, validate_args=True).independent(1)\n",
    "        )\n",
    "        pyro.sample(\n",
    "            \"bias\", dist.Normal(beta_loc, beta_scale, validate_args=True).independent(1)\n",
    "        )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyro.clear_param_store()\n",
    "\n",
    "optim = Adam({'lr': 0.01})\n",
    "\n",
    "num_epochs = 1000\n",
    "batch_size = 50\n",
    "\n",
    "N = 1000\n",
    "p = 3\n",
    "\n",
    "X, y, w = build_logistic_dataset(N, p)\n",
    "\n",
    "example_indices = np.random.permutation(N)\n",
    "n_train = int(0.9 * N)  # 90%/10% train/test split\n",
    "n_test = N - n_train\n",
    "test_iter = 50\n",
    "\n",
    "X = torch.from_numpy(X)\n",
    "y = torch.from_numpy(y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[epoch 000]  average training loss: 1.2700\n",
      "[epoch 050]  average training loss: 0.3476\n",
      "[epoch 100]  average training loss: 0.2844\n",
      "[epoch 150]  average training loss: 0.2586\n",
      "[epoch 200]  average training loss: 0.2638\n",
      "[epoch 250]  average training loss: 0.2504\n",
      "[epoch 300]  average training loss: 0.2299\n",
      "[epoch 350]  average training loss: 0.2368\n",
      "[epoch 400]  average training loss: 0.2376\n",
      "[epoch 450]  average training loss: 0.2421\n",
      "[epoch 500]  average training loss: 0.2412\n",
      "[epoch 550]  average training loss: 0.2241\n",
      "[epoch 600]  average training loss: 0.2325\n",
      "[epoch 650]  average training loss: 0.2162\n",
      "[epoch 700]  average training loss: 0.2396\n",
      "[epoch 750]  average training loss: 0.2344\n",
      "[epoch 800]  average training loss: 0.2515\n",
      "[epoch 850]  average training loss: 0.2385\n",
      "[epoch 900]  average training loss: 0.2319\n",
      "[epoch 950]  average training loss: 0.2378\n"
     ]
    }
   ],
   "source": [
    "lr_model = LogRegressionModel(p=p)\n",
    "\n",
    "svi = SVI(\n",
    "    lr_model.model, lr_model.guide, optim,\n",
    "    loss=Trace_ELBO(), use_cuda=False\n",
    ")\n",
    "\n",
    "\n",
    "lr_dataset = torch.utils.data.TensorDataset(X, y)\n",
    "\n",
    "data_loader_train = DataLoader(\n",
    "    dataset=lr_dataset, batch_size=batch_size, pin_memory=False,\n",
    "    sampler=SubsetRandomSampler(example_indices[:n_train]),\n",
    ")\n",
    "    \n",
    "data_loader_test = DataLoader(\n",
    "    dataset=lr_dataset, batch_size=batch_size, pin_memory=False,\n",
    "    sampler=SubsetRandomSampler(example_indices[n_train:]),\n",
    ")\n",
    "\n",
    "train_elbo = []\n",
    "test_elbo = []\n",
    "for epoch in range(num_epochs):\n",
    "    total_epoch_loss_train = train(svi, data_loader_train, n_train)\n",
    "    train_elbo.append(-total_epoch_loss_train)\n",
    "\n",
    "    if epoch % test_iter == 0:\n",
    "        print(\"[epoch %03d]  average training loss: %.4f\" % (epoch, total_epoch_loss_train))\n",
    "        # report test diagnostics\n",
    "        total_epoch_loss_test = evaluate(svi, data_loader_test, n_test)\n",
    "        test_elbo.append(-total_epoch_loss_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAF3CAYAAACMvMPjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8lOW9///XZyYDhDUQkCVAQYqyC5haFD0UN2rrQq3oqcuxm56etsdWW6p+26PW2lMqbbXa86vluNRWj0otjVqrqLjhikiQHUGtQBIWgZAAIcvM9fvjnoQsM0kmM0lmcr+fjwePzNxzz31fczMz77mu+76uy5xziIiISOYIdHYBREREJDEKbxERkQyj8BYREckwCm8REZEMo/AWERHJMApvERGRDKPwFhERyTAKbxERkQyj8BYREckwCm8REZEMk9XZBWjOwIED3ahRozq7GCIiIh3i3Xff/cQ5N6il9dI6vEeNGsXKlSs7uxgiIiIdwsw+bs16ajYXERHJMApvERGRDKPwFhERyTAKbxERkQyj8BYREckwCm8REZEMo/AWERHJMApvERGRDKPwFhERyTAKb5FMtmYx3DEJbsnx/q5Z3NklEpEOkNbDo4pIM9YshqeugeoK7/6B7d59gCkXd165RKTdqeYtkqmW3Xo0uGtVV3jLRaRLU3iLZKoDOxJbLlJLp1synprNRTJNJAwr74//eHZ/b51AsOPKJJlDp1u6BNW8RTLJrg1w/xz4xw9h0DjI6tFoBYOKffD7mbDpaXCuU4opaSze6ZbnfgLluyBck/g2VZPvcEnVvM1sAPAYMAr4J3Cxc25/o3WmAr8H+gJh4OfOuceS2a+I71QfgVcXwut3Qo9+8KVFXi1p7V+8L+MDO6DfcDj9vyAYgpd+Do9eCnn5cObNMPpfOvsVSGer2A8b/+7VtGM5uAt+fZx3O3sA9BoU/Tcwzu3o/S3PqSbfCcwl8cvczG4H9jnnFpjZDUB/59z1jdY5DnDOuS1mNgx4FxjvnCttafv5+flu5cqVbS6fSJfw0XJ46nuw7wM44Stw9s+hV27zzwnXwOqH4ZVfQlkRHDsbzrgJ8qZ3TJklPVQehPefhXV/hS3PQ6TaO50SCTddt2cuzP5/cOgTOLQHDu4+evvQHjjS4ld2Q/1GwLXrUvM6fMTM3nXO5be4XpLhvRn4nHOuxMyGAi87545v4TnvARc557a0tH2Ft3SqNYsb1mrPuKljaxKH98HzN0Hhn6H/KDj3ThgzO7FtVB+Bd+6F5b/2mtPHnw+n/wQGNfsxlUxWfQS2vuAF9vvPQvVh6DMMJl0Ik74Me7c2rCkDhLLhvLuaf3/XVMHhvUfDvDbYn/tx/Od87VkY8VkIZN4Z2oLCIhYu3UxxaQXDcrKZP+d45k7La/f9dlR4lzrncurd3++c69/M+icBDwITnXORlrav8JZO0/iiHmjdF1wqOOd98T57gxfgp/wnzLoeuvVs+zaPlMGb/wNv/s77Mj/hUvjc9ZAzMnXlzlSd/SMtFcLV8NErsG4JbHwKKsu8mvSEuTD5Ihgxo2GApvI13zEpTlO8AQ76jfTKMOUSOGZc2/bRwQoKi7hxyVoqqo+2UGSHgvziwsntHuApC28zewEYEuOhHwMPtja8a2vmwJXOubea2d/VwNUAI0eOPPHjjz9u6TWIpI5z3hfaolleLaOxfsPh2vXtt//SbfD362Dr8zBsOpx/FwyZnLrtH/oElv/Gq43jIP/rcNoPofeg1O2jnaW0RtSZP9ISEStsJ10E2970fuhtKPDer937wfjzvFr26FkQ7IAORfGO4Tm3Q7A7rF0MH7wILuK9l6dc4rUA9B3W/mVrXM5W/mCZueBFikormizPy8nm9RtOb9diplWzuZn1xQvuXzjn/tLa7avmLe3u8D4oWgXFq6DoXe/foT3NP2f8eXDs57zzyAOOBbPkyxGugRV/gBdvAwzO+C846er26+51YAe8vMA7L56VDSd/26vhv7+0xS+42vAsKq0gaEbYOfI6qFkx5TWieLXGdDpfGyscA1kQ6g2VpRDqCcef4wXip8+ErO6dU8bm3jcHd3utAmse8z5rmHcR5ZSLvc9Tj37tX74EfqSNvuFpYiWjAR8t+GK7FRM6LrwXAnvrXbA2wDn3o0brdAOeAZ5yzt2ZyPYV3pKQlr5Aqg5ByRovoGvDev8/ow8aDDwO8k70Lup69XbvC6exUC/oOeDoF36/kTDmc16Qj57V8oVksZS8B09eAyWrYewc+OKvIWdE4ttpRv3aar/sEGZQeriaGX338auBT5FXvNQLgXC1d1FTrUZfcLHCs1ZHNCumtEZUsR/3y1HE+unlAPuvvSmvucb7f2i2BeGOibEH3snKhgt+5wV3t14pLWd7en7562x7+UHOqH6ZUYFdHCHEsvB0Xu0xm1M+/69ccOLolDXrF6zaxl1PvsUjbj6DrekFd4ezh9Lz+k1Nlsd7nwXNiDjXrufAOyq8c4HFwEhgGzDPObfPzPKBbznnvmlmlwMPAPXbGr/qnFvd0vYV3mkk3c8LxvplndUdJs3zzvUVrYLdG7ymO4C+w72Qrg3roVOhR9/mt1cbZJPnwd4P4MOX4MOXvavBKw8ABkOnHK2Vj5wBoey6L+z8sue5sdtfGMwnWL/hMOsG+GSzdy66Zy6c80uY+KWYNfnmmooLCou45cn1lFZ4odu/Z4ibz5tY9/hPCtby8FvbYtYkwAvd/+/0ALNfvwJqjjRdofdg+I83oGcuM3/5UswvtVpBM3598QnMnZZHQWERq59exDerHmJYYC9HsofQ85xbY75vWtsU3pYaUf1tj+jXjQXTPuGU8qVeP/hwVdzXQo8ctufO5P7dYykoH0/PnGPqytX4mAcMIo4GLRCNX9PscYP467tFMX/4QL0fP1OOgeJC+OjV6L9XYq4fwRhz5OEGx6tNPw6aEes1vLRpD0WlFbVntIGm77l42+iXHeJQVQ3VYQc4ptoHzA2+xrnBtxhoZZS63lQNmsgx+1dDuLJuGzXBHtxm3+LBgyfRr0cWfeww3Y58wrg+FXx1Sk8+M7Da6+p2cDcc3EXpnh1UH9jJAHeAoMXPuIgznpy7Pu5nqTnt9WO1Q8K7vSm800QHnBdM+jxm3Itm8EYcG1YvqIdNhz6DW95ma3+whGu8L9sPX/YCffsKr/aa1YPd/afxp12jiYSr+W7WE/S0o2HhMAwH0/8NzrqVgk2HYx6D5pqKV368j4fe2hb3JeRkh1r1RZSXk81rRy70yhNHJSF2RvpTQi4lbgA73QCKXS473YDo/Vz20gcXHftpbtbr/Hfwfxu85krrzi3u33nkyAzgaOjVD4JaPUMBuoeClB6urgui/Ydjv5Z4Ne/aYzesZhsXBV/lS8HXGGL7qQr1o9u0f2Xh6/v4TqP/lwrXjUfCs/l03zATDq1goJURdsYqN5blTGf3kFk8uq0vxKyze/8300f24/UP9sU9lvUFiDDePuaUwHpmd9/MKVmboeogABsjIxlhu+ltTX9U7YgM5NSqu+r2+eUT85r9cRAKGgsvOqHB5ypeOBeXVtAjFKCiusXrimNuP5EQBMiihlMDa5kbfJ3zA28QiHFoK10Wu1x/BtkBsi3Gj65AFvQezP5Af1bv78bOcD/20I89LofvZ/2VXCtv8pQdkYGcVnUX2aEAhxN4rbXa4xy4wltSJ04w7mQQb13wStK/PNt6HrOgsIiFz25ieFkhj3b/WZzmT+OJ89ex8Ln3Y35B1YYkEPdLLJEfEwWFRfzu2dWMKC/k89mbmB5ezVji/KgAjnTPpceNH8atHffqFuRQVbyaWmJfri15rfs1DLdPmizf6/rwu/CXGMw+hto+htg+hrGXwbaPbtawbLVfsCXkMtk+oqdVNtle/dBJhdrgr/8DIGDQ2x3kvOBbXBR8lWmBrdS4AC9HTuDx8CxejExjUE5fDlfVcNqRl/hR1mKG2V6KXS6311zMk5FTo9uOcIJ9yOxgIWcECpkU+CcARS6XF8PTeDEyjTciE6mkW115zg+8Ft3eJxS7gQ2253GMtSJOCaznlMB6PhvYSI4dAuCDyFDWd5/KM4eO4+3IePbRl/MDr7EgdG+DHxiHXTduqP5mo+0mdsxCQaMqnNrv/9prIJLxYfdLY4a3c/C3yKnscV4g73H92ENO3e0D9CKnp3e+v/GPvPY4hrUunzGS2+am7qJShbekRiQCt8bu/eccnB75H7534eykArw15zEbN725mkpmh9/gG1n/YHLgn4SdxWweS1VQGHDKmAH8c29FzJpx7QVcsQxiPyu6fyfmdW0RZ1wx/JlW19LaU6JfcEaEXMoZYnsZZnsZYvXC3fZykm2K+Zqdgz+Hz2KjG8nGyKfY7IZTQeNhXtsmQITTAmu5KPgKZwfepbtVszkynL+EZ/FEeCZ7yGl5I80YzD5mB1dzeqCQmYF19LJKKlw33ohM5MXINLJcDdeHHm1yDH9dcxEVZHNyYD0zAhsYZGUAbI8M4o3IRN6ITODNyER2E/uzdvQHQdMfGF3Na92uYXig6Y/IZD/L7XkMUxngCm9JXsl78PQPYceKmA875513ezN4IqdeMh/GnhX36ujmmsXjnccEL8DrnyvMoZxLg8v4t6znGWL72RoZxv3hc6h0QX4W+mO7/LJOhfb6Qkq1VH7BxXvNlS6LKkL0Me/HTsQZ/3SD2eRGsjEyko3uU2xyI9nhBhKraTpWzXa9G82Xg8u5MLicIbaf/a43T4RP4fHwv7DOjY65nWR1p4rPBjYyO7CaMwKrGBlooZcCsNP1583IBN6ITOTNyAR2uGNSXq5M15615PaSyqvQFd7Sdof3eWNjr7zfG+N43Be9vpr1znkfdt24veZiBlo5Fwdf4Rgr5XD2EP6v+nPce/BUgjl5rTpnC/CDxe+12NQ2xor4evBZLgwuJ9uqeDU8mfvD5/BKZErdOdZ0rp1k4hdSspp/zTMZbnuYYB8zzrYzPvAx420bowK76tYtc9kNAz0ykk/bdm4NPdhgm7WtLjUuwCuRE/hLtFm8ilAHvlrHGCvmhW7z47Y2nFH1Kz50Q2mPHxJdTTp/luO585KpKbl4TeEtiYtEKHzybkat/hV9XTmP2ef5vV3CjopuXNl7BVfXPMQQ1/TDFLIazrBVXBZcxmnBtYSd8WJkOv8XPp3hnzmXFzfva/YK5fgcpwXW8o3gM3wu+B6VLsSS8Kk8EP4877vUdqXqCJn4hZSsRF9zT44wzrYxLrCd8fYx4wPbGGfb6i7Yci52t/pS15OzKheyJ06zc0fJlBYWSb2c7BCrbz476e0ovCUxRavY95drGFC6lhWR47m5+qtsdJ9qsEoo4F2MEmnmLTPCdvGV4EvMC77MICtjhxvIozWzWRz+XNzzeY11p4q5wdf5evAZjg/sYI/rx59qzuLh8Jnso2/LG5AuxYhEa+nbuCd0R9xrB46tfLjjC9eIH1tY5KhU1L4V3tJA3HPOh/bCsp/Cqj+xh378vOpSCiIzSbZpL0QNZwVWcmlwGacG11PjArwQOZH/C5/O8shkzgu80eS85ZuRSVyR9RyXBZeRa+VsiHyK+2rO4anIyR3cBCrpKhNqtn5sYRFPKrqOKbylTqxzzr1CxsPTNnLc+jvpVnOQP4bncGfNlzlIEpNfxDHKSvjX4EvMC75CrpXzSaQP/ewwoXrdjGpcAKJnr5dFpnNf+BzeioxH5welPtVsJZ2l4sK11oZ3B4xaL+0hVk0aqOuyFGvQi1rTbAu32gNMXvNP3oqM56bqr7brOeR/uqEsqLmU39TMY07gHRaG7mkQ3ABZFuGg68G5VT/nn25ou5VFMtuTkVOhGtVsJS0Ny8nusH2p5p2Bmhtfujm5HOD6rEe5OOsVdrr+/Hf1ZTwZOZmOrt3GG4QhXc5biogkyoA7OvCct2reGeinT61vMbjr94Utcbm8EZnA2cF36Ukl99Scy901X+IQHfcrsb5iNzDmSF7Frg2Tekidscf04sM9h5Me4aoryms0al6/7BBVNeE2DYkpqdEtaNREmr8AtqMZkBW06NjriT3vshkj231WvfoU3hmmoLAo7vjOtRqfF8yzvcwLLGdTJI/vVH+fD1zHvcFiub3m4pjnLW+vSaOJTtJcdihAVY0j7BxBM77y2RENRniadutzLb5POlu3oNGre1bd2OVlR6pT/kXemmF2mxskSNqm9rRdvOFSa7tVxWpFDAWM3j2y6iZUafyjq3aylZbGu2/u1GFzFl50Qt0pyUAz5e/VPSs1c8q3kcI7w9zy5PoW17k+67EGwVirN5WdHtzQNc5b1v7SfuTt7XUBOuPY/qwvLm/1ZAxt1ZpAuvm8iS2eWunfM8QXpwzl7++VtKnMsYaMHZWb3eJQr7F+bIAXoons+7IYQ1K2ZYKbYTnZccchyGkwC1bscrRH8NeGXiLj19dO8tLS+OKtGX88FDQu+cyImO+N2glIgBZnMIs3QNMt508EqPu/uXHJmrrX2btHVtwZymKJt48vn5gXcwa05gzLyWbutLwGs4zFK39Hh3VjOuedRmLNxNO/Z4gJQ/vw1of7m/3AhahhZmAtXwy8zUXBV9O6L2xX0FyXkHgzKoUCRnUrqpbNfdHkJThJSmunh4w3rSi0/AXd2E8K1tb9qEmk3M2NcT9/zvHJzTrXjJYmxmnpODY3rv3lM0aS/6kBcabEbF79meNqj6cBPaOT1dSGcONj25qJfgoKi5j/l/divh97dQvy8y81XDeZY9/S8wsKi5j/+HsNjkmi0222Zh8tzXIWb59Jz3iYIHUVyzDNfZji8QJ7HV8IvM2c4Dv0s8OUuWyCROjVAbM5JSsA9GjjVHxHp0DcEbdmUn+O5ZbmUo4nFDAw2vTFEutDv/LjfU1mD4vVTNi4rO01d3A6aevscqnad1u/oGOVO17LQON9tXWq01S8npbmge8orZmYKFVSPd95e1B4Z5CCwqJWje8N3ry3MwPr+WLgLc4OriTHDlHmsnk+ciJPh2fwWmQynw+83Sl9YQMG/bJDrTrX2praQv1Qa+6D1tov3sbrHa6qiVnWoBkR5+JOF5rsh7y1X6wd+Ws/XWTq60623PHOu6dywot05efXHovCO0209KH+ScFaHnprW7PbyKKGUwLr+WLgbc4OrqS/HaS8LrA/y/LIlCYjkLX3KE+Nr95t6TwXtPzLvqO/uDuzpidSX0fWPtONn197LArvNNBSiK38eB8PvbUt5hSH/4jM4OTABr4YeIs59QL7hch0ng7PYHlkMpV0a7eyZ4eC9AgFYtZMW/OhypQaVKaUU7o2P/+Q9PNrj0XhnQbi/aKsL9ZwjzUuwBFC9LZKDroePB85kX+EP8urkSkpCezmghmONmkD+lCJdBA//5D082tvTIO0pIHiVkyD+aOsxU26dWVZhICDq6uu5ZXICSmtYde/iri1wawPlUj7q99FyW/8/NrbSuHdjnJ6xr94y4jwL4G15MUYaQygB1U8F/lM0mWo7XYUr6tOS8GsD5WISPpReLeTnxSsjRncfTnEvOArXB58ntGBXYSdEYxxrWUqhgptqW+tgllEJDMpvFMs3mAA42wb/xZ8jrnB1+lplayMHMcdVReRRQ23hR5I+VChfr1SU0TEDxTeKVB7sUXji9OyqOHzgXe4Iut5PhvYxBEXoiA8kz+Hz2a9G1W3Xk11sM3dunKyQ1TWRJqcu6694ExERLoehXeSYnVzOIb9XJq1jK8EX2SwlfJx5Bhuq76Mv4RncYDeTbbxZORUnqxKvA92/TGCdVGZiIh/KLzbYs1iWHYrHNjBDAZyVngeTzKTz9hmrsx6jjmBdwhZmJfCJ3B9+CpeiZyAI5DQLi6fMbLJEJmNh9GsH9IKaxER/1B4J2rNYnjqGqj2msiHsIdfhf7A9e5R8gL7OOB68sfwHB4Kn8nHbkibdpGXk81tcyc3mNBANWoREaml8E7UslvrgrtWNwsziANcX30VT4RP4Qjd27z5+uerdTW4iIjEovBOkDuwgxizbZJFmMfCsxPeXu1kHuk2s42IiKQvhXeCdjGQIexpsrzYDUx4W501BZ+IiGS2xK6iEv67ah6HXcPhStvSLzsvJ5vCm85WcIuISMJU805AQWERT0VOhWq4I/R7AjiKorOAJTLdpvphi4hIMlTzTsDCpZtxwLORzxI0xx01F3Fq1V2tCm4zb5zxvJxszcolIiJJUc27lQoKi+pGUDvG9gNQwoBWPVfTaIqISCqp5t0KtaOo1RrCPgB2uf4tPtcMBbeIiKSUwrsVFi7d3GCks6HmhXdJCzN/GXDHxVMV3CIiklIK71ZoPOHI4GizeXM1717dgtxxiYJbRERST+e8WyFoRtgdnXN7iO3jsOtOGT3rloWCxsKLTlBYi4hIu1PNuwUFhUUNghu88C5xAyA61pqBgltERDqMwrsZBYVFzP/Le02WD7H9TZrMFdwiItJR1GweQ0FhEbc8uZ7SiuqYjw+xfbztxtXdH5aT3VFFExERUXg3Vlvbro64mI8bEQazn13O6+MdCppGSxMRkQ6lZvNGFi7dHDe4AXIpJ2RhdkabzXWuW0REOprCu5HiRt3CGhtiewHY6QaQl5Ot4BYRkQ6n8G4kp2eo2ceHRPt473QD1FwuIiKdIunwNrMBZva8mW2J/o07comZ9TWzIjP7XbL7bS+HKmuafXxIdHS1w90Hq9YtIiKdIhU17xuAZc65scCy6P14fga8koJ9touCwiKqwvHPd4MX3jUuwHfPP7mDSiUiItJQKq42vwD4XPT2g8DLwPWNVzKzE4HBwLNAfgr2mzIFhUUsXLq5yTCosYzMKqWqxzHMnT6yA0omIiLSVCrCe7BzrgTAOVdiZsc0XsHMAsCvgSuAM1Kwz5SpnTGs/sQj8YSCxoyBlfTsMaIDSiYiIhJbq8LbzF4AhsR46Met3M+3gX8457abWUv7uhq4GmDkyPav3TaeMaw51WHH4b3bYdyJ7VwqERGR+FoV3s65M+M9Zma7zGxotNY9FNgdY7WTgdPM7NtAb6CbmR10zjU5P+6cWwQsAsjPz2/+BHQKtNQ1rLGBkb3QVxeqiYhI50nFBWtPAldGb18JPNF4BefcZc65kc65UcAPgT/FCu7OkMjQpr05TG87An2HtmOJREREmpeK8F4AnGVmW4Czovcxs3wzuzcF229X8+ccT3Yo2Kp1PxU64N3oM6wdSyQiItK8pC9Yc87tJcZFaM65lcA3Yyz/I/DHZPebKrV9tb//2Oq46xheDf3/TesHbwJ9Fd4iItJ5NMIaXoAH4lxH179niI8WfJHXbzidmcdUeQvVbC4iIp1I4Q38pGAt8eYicfWXlxV7f9VsLiIinUjhDTzy9va4jx2oP6d3eTFkD4BQjw4olYiISGwKbyDs4vdIa3A1elmJzneLiEin8314FxQWNft4g5nDyooU3iIi0ul8H94//tvauI/NHDOg4cxh5SXQRxeriYhI5/J1eBcUFnGoKv7QqA9fVW/msJpKOLRHo6uJiEin83V4L1y6ufUrl+/0/qqbmIiIdDJfh3dz45rnZIcaLigv8f6qm5iIiHQyX4d3Ts9Q3MduOX9iwwVl0QvbdMGaiIh0Mt+Gd0FhEQeP1MR87PIZIxteqAZeNzFQs7mIiHQ634b3T59aT3WMYdVyskPcNndy0yeUl0CoJ/TI6YDSiYiIxOfL8C4oLGL/4eqYjzUYUa2+siKvm5jFGQRdRESkg/gyvJu7yrxf4wvVaml0NRERSRO+DO/mrjIvr6yJPepaWbHCW0RE0oIvw7vBeOWNhCOuac08EtHoaiIikjZ8Gd4NxiuPoUnN/PBeiFSr5i0iImnBl+HdkiY1c/XxFhGRNOLL8G7ugrVQ0JrWzDW6moiIpBFfhndzF6xd8pkRMQZoUc1bRETShy/Du7kL1l7atKfpwrISsCD0PqYdSyUiItI6vgzv5i5Yi1krLy+B3oMhEGzHUomIiLSOL8N77rQ8+mVnxXwsZq28rEhN5iIikjZ8Gd4A35o1psmy7FAwdq28rEQTkoiISNrwbXifPGYgALm9umFAXk42v7hwctOL1SA6ulqM5SIiIp0gdtuxD1RUhQG4+9JpnBIN8pgqy6GqXKOriYhI2vBtzftIjRfe2aEWLkKrm8db57xFRCQ9+DK8CwqL+MFjqwH49z+/G3siklrq4y0iImnGd83mBYVF3LhkLRXVXs17d3klNy5ZCxD7fHfd6GpqNhcRkfTgu5r3wqWb64K7VkV1OP6Qqap5i4hImvFdeMcbGjXukKllJZDdH0LxR2UTERHpSL4L73hDo8YdMrW8RBOSiIhIWvFdeM+fczyhoDVYFnMmsVoaXU1ERNKM78IbANfC/fo0upqIiKQZ34X3T59aT3WkYVpXR1zsC9ZqquDQbo2uJiIiacVX4V1QWMT+w9UxH4t5wdrBnd5fdRMTEZE04qvwjtsdjHiziWl0NRERST++Cu+43cGA2eMGNV2oPt4iIpKGfBXecbuDAS9t2tN0oUZXExGRNOSr8I7bHYw4tfKyYsjK9gZpERERSRO+Cu+50/Lo3zMU87HY57yLvW5iZk0fExER6SS+Cm+Am8+bSKBRFmeHgrFr5RpdTURE0pDvwnvutDwmDO1LVsAwIC8nm19cODn2jGIaXU1ERNKQ76YEBRjSL5uIg39877T4K0UiUL5To6uJiEja8V3NG8A5R6ClV354L4SrNLqaiIikHV+Gd8Q5Ai1dhFZe7P1VNzEREUkzSYW3mQ0ws+fNbEv0b8w+VWY20syeM7ONZrbBzEYls99kRRxYS+Gt0dVERCRNJVvzvgFY5pwbCyyL3o/lT8BC59x44CRgd5L7TYpX825hJY2uJiIiaSrZ8L4AeDB6+0FgbuMVzGwCkOWcex7AOXfQOXc4yf0mxTla0WxeAhaAXsd0TKFERERaKdnwHuycKwGI/o2VdMcBpWa2xMwKzWyhmQXjbdDMrjazlWa2cs+eGEOWpkA40pqadzH0HgJBX16QLyIiaazFZDKzF4AhMR76cQL7OA2YBmwDHgO+CtwXa2Xn3CJgEUB+fr6LtU6yWnXBWu3oaiIiImmmxfB2zp0Z7zEz22VmQ51zJWbC6dm7AAAfF0lEQVQ2lNjnsncAhc65D6PPKQBmECe8O4JzEGyp6l1eArmf7pgCiYiIJCDZZvMngSujt68EnoixzjtAfzOrnXPzdGBDkvtNSqQ1/bzLitXHW0RE0lKy4b0AOMvMtgBnRe9jZvlmdi+Acy4M/BBYZmZrAQP+N8n9JqXFZvPKcqgsU7O5iIikpaSuxnLO7QXOiLF8JfDNevefB6Yks69UarGfd10fb9W8RUQk/fhyhDXXUj9vja4mIiJpzDf9oAoKi1i4dDPFpRVkBY2xx/SOv7JGVxMRkTTmi5p3QWERNy5ZS1FpBQ6oDjs27SynoLAo9hNqR1dTzVtERNKQL8J74dLNVFSHGyyLOG95TOUl0CMHuvXsgNKJiIgkxhfhXVxakdByykp0sZqIiKQtX4T3sJzshJZTVqRuYiIikrZ8Ed7z5xxPdqjhcOpB85bHVF6i890iIpK2fBHec6fl8YsLJ9O3x9GL66cMz2HutBhN4+FqOLhbzeYiIpK2fBHe4AX4NWeMrbs/MjfOxWjlOwGnZnMREUlbvglvgFDw6MuNOzxqebSPdx/18RYRkfTkq/DOCh4N7Lijo9b28dYALSIikqZ8Fd6helOJBeOlt0ZXExGRNOev8M46Gthxm83LiiCrB2T376BSiYiIJMZX4Z1Vr+Yddz7v2m5izc06JiIi0ol8Fd71a9txpwQtK1GTuYiIpDVfhXfEubrbcacELStSeIuISFrzVXiv+Ghv3e2/rSpqOquYcxpdTURE0p5vwrugsIhH39led/9QVZgbl6xtGOCH90K4SqOriYhIWvNNeC9cupnqsGuwrKI63HBa0LJi769GVxMRkTTmm/Bu1bSgGl1NREQygG/Cu1XTgmp0NRERyQC+Ce/Z4wa1vLysBCwAvQd3UKlEREQS55vwfmnTnpaXlxV7wR3MirmuiIhIOvBNeLfunHexuomJiEja8014t+6ct0ZXExGR9Oeb8J4/53hCwYbDqmWHgsyfc/zRBWXFCm8REUl7vgnvudPyOO+Eo8HcPSvALy6czNxp0QFZKg9C5QE1m4uISNrzTXgDTMnrV3d7cl6/o8ENR/t4a3Q1ERFJc74K70i9AdZWfryfmQtePDo8qkZXExGRDOGrPlFrdpQ2uF9UWsGNS9YCMDeg0dVERCQz+Krm/eKm3U2W1Y1vXje6mmreIiKS3nwV3mVHamIuLy6t8LqJ9egH3Xp1cKlEREQS46vw7tsj9lmCYTnZ0W5iulhNRETSn6/C+7SxA5ssq+vrrdHVREQkQ/gqvCcM87qKDevXAwPycrKP9vUuK9H5bhERyQi+uto8Eu0r9uqPZpMVrPe7JVwNB3ep2VxERDKCr2retf28A9ZwmFQO7gKcms1FRCQj+Cy8vfRunN1HB2hRzVtERNKfr8LbOYcZWOP01uhqIiKSQXwV3hEXo8kcjo5rrtHVREQkA/gsvB0xotsbXS3YHXoO6OgiiYiIJMxn4R2n5l3bTSzWYyIiImnGV+HtcLHzWaOriYhIBvFXeMc9563R1UREJHP4KrwjEUegcXY7p9HVREQkoyQd3mY2wMyeN7Mt0b/946x3u5mtN7ONZnaXNemv1f5invM+vA/ClWo2FxGRjJGKmvcNwDLn3FhgWfR+A2Z2CjATmAJMAj4DzErBvhMScTHOeZdH+3ir2VxERDJEKsL7AuDB6O0Hgbkx1nFAD6Ab0B0IAbtSsO+EOOcING431+hqIiKSYVIR3oOdcyUA0b/HNF7BOfcm8BJQEv231Dm3MQX7TkjE0bSft0ZXExGRDNOqWcXM7AVgSIyHftzK538aGA8Mjy563sz+xTn3aox1rwauBhg5cmRrNt9qDtf0nHd5CWDQe3BK9yUiItJeWhXezrkz4z1mZrvMbKhzrsTMhgK7Y6z2JeAt59zB6HOeAWYATcLbObcIWASQn5/vWlO+1oq4WOOaF3nBHQylclciIiLtJhXN5k8CV0ZvXwk8EWOdbcAsM8sysxDexWod3mzuXIyuYuomJiIiGSYV4b0AOMvMtgBnRe9jZvlmdm90nceBD4C1wHvAe865p1Kw74REIjG6ipWX6GI1ERHJKK1qNm+Oc24vcEaM5SuBb0Zvh4F/T3ZfyYrErHkXwadmdkp5RERE2sJfI6w1PudddQiOHFCzuYiIZBRfhbdrPEhLWXQebzWbi4hIBvFXeNPonLdGVxMRkQzkm/AuKCzi2XU72bbvMDMXvEhBYVG9mvewzi2ciIhIAnwR3gWFRdy4ZC0V1WEAikoruHHJWtZvivZWU81bREQyiC/Ce+HSzXXBXauiOszG9zdD937QvXcnlUxERCRxvgjv4tKKmMv7VO1Wk7mIiGQcX4T3sJzsmMtHZJWqm5iIiGQcX4T3/DnHkx0KNliWHQoyunsZ9FHNW0REMosvwnvutDx+ceFkumd5LzcvJ5sFc8eTXfmJms1FRCTjJD08aqaYOy2Ph9/+mFAwwP9dNQMOFIGLqNlcREQyji9q3rWqwo5QMPqSy6IDtGh0NRERyTC+Cu/qmsjR8NboaiIikqH8Fd7hCN2yosOjanQ1ERHJUL4L76PN5kUQ7AY9czu3UCIiIgnyWXjXO+ddXuI1mVvjCb5FRETSm6/Cu6pBzbtYF6uJiEhG8lV4V4cjdAvWnvMuVjcxERHJSP4K79qrzZ072mwuIiKSYfwV3mFHKCsAFfuh5oiazUVEJCP5Jrydc0fPedcN0KKat4iIZB7fDI+6ZNUOAO5atoXiFRv4FajmLSIiGckXNe+CwiJ+XLCu7n7W4Z0ALN2mbmIiIpJ5fBHeC5du5kh1pO7+ENtHxBk/f3V/J5ZKRESkbXwR3sWlFQ3uD2Efn9CP7QeqO6lEIiIibeeL8B6Wk93g/hDbz07Xv8lyERGRTOCL8J4/53i6Zx19qUNsH7stl/lzju/EUomIiLSNL8J77rQ8vjpzFAAGDAvs59jRY5k7TVebi4hI5vFFV7GCwiIWv7MdgOG9oW/NQfoe++lOLpWIiEjbdPmad0FhETcuWcv+w97FacFD3jze7+7v2ZnFEhERabMuH94Ll26mojpcd3+Ied3DHlxX1VlFEhERSUqXD+9Y3cQA1h/s1RnFERERSVqXD++m3cS88La+wzqjOCIiIknr8uE9f87xZIeCdfcH237KXTbf/fzUTiyViIhI23X5q81ru4P91xPrKD9Sw+hupbhew9RNTEREMlaXr3mDF+BfnzkagFlDa+g7aGQnl0hERKTtfBHeAC7618qKNRWoiIhkNN+EN84RJAwHd0HfoZ1dGhERkTbzTXg7YJAdABeBPgpvERHJXP4Jb+eNaQ6o2VxERDKab8I74lzdAC1qNhcRkUzmm/B2wNBAbXir5i0iIpnLP+HtYDD7INgNeuZ2dnFERETazD/hjWOw7Yc+Q8Css4sjIiLSZr4JbxwMsb1qMhcRkYyXVHib2TwzW29mETPLb2a9z5vZZjPbamY3JLPPtnLAYParm5iIiGS8ZGve64ALgVfjrWBmQeB/gHOACcBXzGxCkvtNWCQc8c55azYxERHJcElNTOKc2whgzZ9DPgnY6pz7MLruo8AFwIZk9p2o7uEysq1K4S0iIhmvI8555wHb693fEV3WofpU7YneULO5iIhkthZr3mb2AjAkxkM/ds490Yp9xKqWuxjLavd3NXA1wMiRqZv9q0/Vbu+GLlgTEZEM12J4O+fOTHIfO4AR9e4PB4qb2d8iYBFAfn5+3JBPVF3NW6OriYhIhuuIZvN3gLFmNtrMugH/CjzZAfttoE91NLx7x2pEEBERyRzJdhX7kpntAE4GnjazpdHlw8zsHwDOuRrgu8BSYCOw2Dm3PrliJ2jNYk7Z9X9eW/3d02HN4g7dvYiISCqZcylrmU65/Px8t3LlyuQ2smYxPHUNVFccXRbKhvPugikXJ7dtERGRFDKzd51zccdNqdX1R1hbdmvD4Abv/rJbO6c8IiIiSer64X1gR2LLRURE0lzXD+9+wxNbLiIikua6fnifcZN3jru+ULa3XEREJAN1/fCecjGcdxf7Q4OJYNBvhC5WExGRjJbU2OYZY8rF/HLLOF7avJu3r012zBkREZHO1fVr3lER57CYI7WKiIhkFt+Et3PQ/ORnIiIimcE/4U3sGVJEREQyjX/C27U477iIiEhG8E94x5+FVEREJKP4JrzROW8REekifBPeDggovUVEpAvwTXhHnFPNW0REugTfhLdzutpcRES6Bv+EN7raXEREugb/hLdzqnmLiEiX4J/wBrWbi4hIl+Cb8EbnvEVEpIvwTXg7nLqKiYhIl+Cb8I5ENEiLiIh0Db4Jb4emBBURka7BP+Gt4VFFRKSL8E94d3YBREREUsQ/4a0pQUVEpIvwTXiDBmkREZGuwTfhHXEQ8M2rFRGRrsw3ceYNj6q6t4iIZD7/hDe62lxERLoG/4S3hkcVEZEuwj/hDap6i4hIl+Cf8NaUoCIi0kX4JrwBAkpvERHpAnwT3hHnNEiLiIh0Cb4Jb12wJiIiXYW/wlvpLSIiXYB/wltTgoqISBfhn/B2qN1cRES6BP+EN8puERHpGnwT3jgI6KS3iIh0Ab4Jb6+rWGeXQkREJHm+CW9NTCIiIl2Ff8JbU4KKiEgX4Z/wRjVvERHpGvwT3q6zSyAiIpIaWck82czmAbcA44GTnHMrY6wzAvgTMASIAIucc79NZr9t4dW8VfUWEWlOdXU1O3bs4MiRI51dlC6tR48eDB8+nFAo1KbnJxXewDrgQuAPzaxTA/zAObfKzPoA75rZ8865DUnuOzHOaVYxEZEW7Nixgz59+jBq1ChVeNqJc469e/eyY8cORo8e3aZtJNVs7pzb6Jzb3MI6Jc65VdHb5cBGIC+Z/bZFRBOTiIi06MiRI+Tm5iq425GZkZubm1TrRoee8zazUcA04O2O3C9ExzbXm1FEpEX6rmx/yR7jFsPbzF4ws3Ux/l2QyI7MrDfwV+D7zrmyZta72sxWmtnKPXv2JLKLZmlKUBGR9Ld3716mTp3K1KlTGTJkCHl5eXX3q6qqWrWNr33ta2ze3GyjcAP33nsvgwYNqtvP1KlT2bx5M1u3bmXq1KlN1r/88ssZPXo0U6dOZdy4cdx22211j1VWVvKf//mfjBkzhrFjxzJ37lyKi4tbXZbWavGct3PuzGR3YmYhvOB+2Dm3pIX9LQIWAeTn56fsGnFNCSoiknoFhUUsXLqZ4tIKhuVkM3/O8cyd1vYzo7m5uaxevRqAW265hd69e/PDH/6wwTrOOZxzBAKx658PPPBAwvu97LLLuPPOOxss27p1a9z177jjDubOnUtFRQXjxo3jyiuvZMSIEVx//fVUVlby/vvvEwwG+d///V++/OUv8+abbyZcpua0e7O5eW0D9wEbnXO/ae/9xeP9ClB6i4ikSkFhETcuWUtRaQUOKCqt4MYlaykoLEr5vrZu3cqkSZP41re+xfTp0ykpKeHqq68mPz+fiRMncuutt9ate+qpp7J69WpqamrIycnhhhtu4IQTTuDkk09m9+7dKS1XRUUFZkbPnj0pLy/noYce4je/+Q3BYBCAq666CoBXXnklpftNtqvYl4C7gUHA02a22jk3x8yGAfc6574AzASuANaa2eroU/+fc+4fyew7UU5jm4uIJOSnT61nQ3Hcs5wUbiulKhxpsKyiOsyPHl/DIyu2xXzOhGF9ufm8iW0qz4YNG3jggQe45557AFiwYAEDBgygpqaG2bNnc9FFFzFhwoQGzzlw4ACzZs1iwYIFXHfdddx///3ccMMNTbb98MMP8/LLL9fdX7FiRbNlufbaa7nlllvYsmULP/jBD8jNzWXVqlWMHj2a3r17N1g3Pz+f9evXM2vWrDa97liSCm/n3N+Av8VYXgx8IXr7NdKkyquuYiIiqdM4uFtanqwxY8bwmc98pu7+I488wn333UdNTQ3FxcVs2LChSXhnZ2dzzjnnAHDiiSeyfPnymNuO1WzenNpm8/LycmbPns25555LMBiMeSGaV3lMbQAl2887Y0Q0trmISEJaqiHPXPAiRaUVTZbn5WTz2L+fnPLy9OrVq+72li1b+O1vf8uKFSvIycnh8ssvj9n1qlu3bnW3g8EgNTU1KS1Tnz59mDVrFq+99hpXXXUVH374IQcPHmxQ+161ahXz5s1L6X59NTyqms1FRFJn/pzjyQ4FGyzLDgWZP+f4dt93WVkZffr0oW/fvpSUlLB06dJ232cs1dXVrFixgjFjxtCnTx8uvfRS5s+fTyTitT7cf//9hMPhlDaZg49q3pqYREQktWqvKk/l1eatNX36dCZMmMCkSZM49thjmTlzZlLba3zO+w9/+AO5ubls2LCB4cOH1y2/++67gaPnvCsrK5kzZw7nn38+ALfffjs/+MEPGDt2LGbGhAkTWLKk2U5WbWIujWfsyM/PdytXNhkuvU3O+PXLjBvSl/+5bHpKtici0hVt3LiR8ePHd3YxfCHWsTazd51z+S091z/N5pAml82JiIgkxzfhjUZYExGRLsI34e2AgE56i4hIF+Cb8I5okBYREekifBPemphERES6Cv+Et6YEFRGRLsI/4a2at4hI2kvFlKDgDY6yc+fOmI/Vn9Jz6tSpnHbaaYA3Nej3v//9JusPHz6cyZMnM3XqVCZPnsxTTz1V99j27ds5//zzGTt2LGPGjOG6666juro6wVedOF+Ed0FhESWlR1hSWMTMBS+2y4w3IiK+tGYx3DEJbsnx/q5ZnNTmaqcEXb16Nd/61re49tpr6+7XH+q0Jc2FN3hjk9duN9545/UtX76c1atX8+ijj9YFvHOOCy64gHnz5rFlyxY2b97M3r17uemmm1pdzrbq8uFdO2VdODoYTXtOWSci4itrFsNT18CB7YDz/j51TdIBHs+DDz7ISSedxNSpU/n2t79NJBKhpqaGK664gsmTJzNp0iTuuusuHnvsMVavXs0ll1yScI29JWVlZfTv3x+A5557jpycHK644goAsrKy+O1vf8uiRYtijrOeSl1+eNSFSzdTUR1usKyiOszCpZs7ZAg/EZGM9cwNsHNt/Md3vAPhyobLqivgie/Cuw/Gfs6QyXDOgoSLsm7dOv72t7/xxhtvkJWVxdVXX82jjz7KmDFj+OSTT1i71itnaWkpOTk53H333fzud79j6tSpMbdXO7wpwJQpU/jTn/7U7P5PO+00IpEIH330Ud1wp+vXr+fEE09ssF5OTg7Dhg3jww8/bDLDWSp1+fAujjHjTXPLRUSklRoHd0vLk/DCCy/wzjvvkJ/vjRxaUVHBiBEjmDNnDps3b+Z73/seX/jCFzj77LNbtb3aKT1ba/ny5eTk5PD+++8zZ84c1q9fH3eqz/aYArSxLh/ew3KyY05ZNywnuxNKIyKSQVqqId8xKdpk3ki/EfC1p1NaFOccX//61/nZz37W5LE1a9bwzDPPcNddd/HXv/6VRYsWpXTf9R133HEMGDCATZs2MXHiRJ5+uuHrLC0tpbi4mNGjR7dbGcAH57w7c8o6EZEu7YybINSoIhTK9pan2JlnnsnixYv55JNPAO+q9G3btrFnzx6cc8ybN4+f/vSnrFq1CvDm2S4vL095OXbu3Mm2bdsYOXIkZ599Nvv37+fhhx8GoKamhuuuu46rrrqKHj16pHzf9XX5mndnTlknItKlTbnY+7vsVjiwA/oN94K7dnkKTZ48mZtvvpkzzzyTSCRCKBTinnvuIRgM8o1vfKOuqfqXv/wlAF/72tf45je/SXZ2NitWrGhypXr9c94A7777LgD33Xcfjz/+eN3y2pktTzvtNILBINXV1fzqV79i4MCBABQUFPCd73yHW265hUgkwrnnnhuzdSDVfDMlqIiItExTgnYcTQkqIiLiIwpvERGRDKPwFhERyTAKbxERaSCdr4XqKpI9xgpvERGp06NHD/bu3asAb0fOOfbu3ZtUd7Iu31VMRERab/jw4ezYsYM9e/Z0dlG6tB49ejB8+PA2P1/hLSIidUKhULuPDibJU7O5iIhIhlF4i4iIZBiFt4iISIZJ6+FRzWwP8HEKNzkQ+CSF2/MjHcPk6Rimho5j8nQMk5fqY/gp59ygllZK6/BONTNb2ZoxYyU+HcPk6Rimho5j8nQMk9dZx1DN5iIiIhlG4S0iIpJh/Bbeizq7AF2AjmHydAxTQ8cxeTqGyeuUY+irc94iIiJdgd9q3iIiIhnPF+FtZp83s81mttXMbujs8qQrMxthZi+Z2UYzW29m34suH2Bmz5vZlujf/tHlZmZ3RY/rGjOb3rmvIH2YWdDMCs3s79H7o83s7egxfMzMukWXd4/e3xp9fFRnljudmFmOmT1uZpui78mT9V5MjJldG/0srzOzR8ysh96LLTOz+81st5mtq7cs4feemV0ZXX+LmV2ZyjJ2+fA2syDwP8A5wATgK2Y2oXNLlbZqgB8458YDM4DvRI/VDcAy59xYYFn0PnjHdGz039XA7zu+yGnre8DGevd/CdwRPYb7gW9El38D2O+c+zRwR3Q98fwWeNY5Nw44Ae946r3YSmaWB1wD5DvnJgFB4F/Re7E1/gh8vtGyhN57ZjYAuBn4LHAScHNt4KeEc65L/wNOBpbWu38jcGNnlysT/gFPAGcBm4Gh0WVDgc3R238AvlJv/br1/PwPGB79cJ8O/B0wvEEcsqKP170ngaXAydHbWdH1rLNfQ2f/A/oCHzU+FnovJnQM84DtwIDoe+vvwBy9F1t9/EYB6+rdT+i9B3wF+EO95Q3WS/Zfl695c/QNXGtHdJk0I9pkNg14GxjsnCsBiP49Jrqajm1sdwI/AiLR+7lAqXOuJnq//nGqO4bRxw9E1/e7Y4E9wAPR0w/3mlkv9F5sNedcEfArYBtQgvfeehe9F9sq0fdeu74n/RDeFmOZLrFvhpn1Bv4KfN85V9bcqjGW+frYmtm5wG7n3Lv1F8dY1bXiMT/LAqYDv3fOTQMOcbSZMhYdx0aiTbQXAKOBYUAvvCbexvReTE6849aux9MP4b0DGFHv/nCguJPKkvbMLIQX3A8755ZEF+8ys6HRx4cCu6PLdWybmgmcb2b/BB7Fazq/E8gxs6zoOvWPU90xjD7eD9jXkQVOUzuAHc65t6P3H8cLc70XW+9M4CPn3B7nXDWwBDgFvRfbKtH3Xru+J/0Q3u8AY6NXWHbDu2DjyU4uU1oyMwPuAzY6535T76EngdorJa/EOxdeu/zfoldbzgAO1DYr+ZVz7kbn3HDn3Ci899qLzrnLgJeAi6KrNT6Gtcf2ouj6vq/tOOd2AtvN7PjoojOADei9mIhtwAwz6xn9bNceQ70X2ybR995S4Gwz6x9tBTk7uiw1OvuigA668OALwPvAB8CPO7s86foPOBWvWWcNsDr67wt4572WAVuifwdE1ze8K/k/ANbiXdXa6a8jXf4BnwP+Hr19LLAC2Ar8BegeXd4jen9r9PFjO7vc6fIPmAqsjL4fC4D+ei8mfAx/CmwC1gF/Brrrvdiq4/YI3nUC1Xg16G+05b0HfD16PLcCX0tlGTXCmoiISIbxQ7O5iIhIl6LwFhERyTAKbxERkQyj8BYREckwCm8REZEMo/AW6URmlmtmq6P/dppZUb373Vq5jQfq9YeOt853zOyyFJX5NfNm6ast52Op2G697e8ws5xUblOkq1FXMZE0YWa3AAedc79qtNzwPquRmE/sYGb2GvBd59zqdtr+DmCSc660PbYv0hWo5i2Shszs09E5mO8BVgFDzWyRma2Mzs98U711XzOzqWaWZWalZrbAzN4zszfN7JjoOreZ2ffrrb/AzFZEa9CnRJf3MrO/Rp/7SHRfUxMo80Nm9nszW25m75vZOdHl2Wb2oJmtNbNVZvYv0eVZZnZH9HWuMbNv19vc96MTkqwxs+OSPqAiXYzCWyR9TQDuc85Nc94MUTc45/Lx5rY+K8689P2AV5xzJwBv4o3wFIs5504C5gO1PwT+E9gZfe4CvFnl4nmsXrP5gnrLRwCzgPOARWbWHW9O6Srn3GTgCuDP0VMC/4E3YcYJzrkpeGPB19rlvAlJ7gWua6YcIr6U1fIqItJJPnDOvVPv/lfM7Bt4n9theOG+odFzKpxzz0RvvwucFmfbS+qtMyp6+1TglwDOuffMbH0zZbskTrP54mjz/mYz2w6MjW53YXS7682sGPg03sQZdzrnwtHH6k+CUb98X2imHCK+pPAWSV+Ham+Y2Vjge8BJzrlSM3sIbyzqxqrq3Q4T/zNeGWOdWFMYJqrxRTTxpkas3V+8i25ilU9EotRsLpIZ+gLlQFl0OsI57bCP14CLAcxsMl7NPlHzorMrHYfXhL4FeBW4LLrd8cBQvIkangP+w8yC0ccGJP0KRHxCv2hFMsMqvCbydcCHwOvtsI+7gT+Z2Zro/tYBB+Ks+5iZVURv73LO1f6Y2IoX1scAVzvnqszsbuAPZrYWb5amf4su/wNes/oaM6sBfg/c0w6vS6TLUVcxEQG8q7+BLOfckWgz/XPAWOdcTSuf/xDwuHOuoD3LKSKqeYvIUb2BZdEQN+DfWxvcItKxVPMWERHJMLpgTUREJMMovEVERDKMwltERCTDKLxFREQyjMJbREQkwyi8RUREMsz/D3O0UowHJ+iSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_llk(train_elbo, test_elbo, test_iter)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:dose_ensemble]",
   "language": "python",
   "name": "conda-env-dose_ensemble-py"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}