the scVI VAE model, but using Pyro-PPL for inference

pyro_scvi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "J0QZYD_HuDJF"
   },
   "source": [
    "# scVI with Pyro\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "\n",
    "from torch.distributions import constraints\n",
    "from torch.distributions.distribution import Distribution\n",
    "from torch.distributions.utils import broadcast_all, probs_to_logits, lazy_property, logits_to_probs\n",
    "\n",
    "import pyro\n",
    "import pyro.distributions as dist\n",
    "from pyro.infer import SVI, Trace_ELBO\n",
    "from pyro.optim import Adam, SGD\n",
    "\n",
    "import collections\n",
    "from numbers import Number\n",
    "from typing import Iterable\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import anndata\n",
    "\n",
    "from scvi.dataset import GeneExpressionDataset\n",
    "\n",
    "from torch.utils.data import DataLoader\n",
    "from torch.utils.data.sampler import SubsetRandomSampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NegativeBinomial distribution is not in the current version of PyTorch so I reimplemented it\n",
    "\n",
    "class _GreaterThanEq(constraints.Constraint):\n",
    "    \"\"\"\n",
    "    Constrain to a real half line `[lower_bound, inf)`.\n",
    "    \"\"\"\n",
    "    def __init__(self, lower_bound):\n",
    "        self.lower_bound = lower_bound\n",
    "\n",
    "    def check(self, value):\n",
    "        return self.lower_bound <= value\n",
    "\n",
    "class _HalfOpenInterval(constraints.Constraint):\n",
    "    \"\"\"\n",
    "    Constrain to a real interval `[lower_bound, upper_bound)`.\n",
    "    \"\"\"\n",
    "    def __init__(self, lower_bound, upper_bound):\n",
    "        self.lower_bound = lower_bound\n",
    "        self.upper_bound = upper_bound\n",
    "\n",
    "    def check(self, value):\n",
    "        return (self.lower_bound <= value) & (value < self.upper_bound)\n",
    "\n",
    "\n",
    "class NegativeBinomial(torch.distributions.Distribution,\n",
    "                       pyro.distributions.torch.TorchDistributionMixin):\n",
    "    r\"\"\"\n",
    "    Creates a Negative Binomial distribution parameterized by `total_count` and\n",
    "    either `probs` or `logits` (but not both). `total_count` must be\n",
    "    broadcastable with `probs`/`logits`.\n",
    "    :param (Tensor) total_count: number of Bernoulli trials\n",
    "    :param (Tensor) probs: Event probabilities\n",
    "    :param (Tensor) logits: Event log-odds\n",
    "    \"\"\"\n",
    "    arg_constraints = {'total_count': _GreaterThanEq(0),\n",
    "                       'probs': _HalfOpenInterval(0., 1.)}\n",
    "    support = constraints.nonnegative_integer\n",
    "\n",
    "    def __init__(self, total_count, probs=None, logits=None, validate_args=None):\n",
    "        if (probs is None) == (logits is None):\n",
    "            raise ValueError(\"Either `probs` or `logits` must be specified, but not both.\")\n",
    "        if probs is not None:\n",
    "            self.total_count, self.probs, = broadcast_all(total_count, probs)\n",
    "            self.total_count = self.total_count.type_as(self.probs)\n",
    "            is_scalar = isinstance(self.probs, Number)\n",
    "        else:\n",
    "            self.total_count, self.logits, = broadcast_all(total_count, logits)\n",
    "            self.total_count = self.total_count.type_as(self.logits)\n",
    "            is_scalar = isinstance(self.logits, Number)\n",
    "\n",
    "        self._param = self.probs if probs is not None else self.logits\n",
    "        if is_scalar:\n",
    "            batch_shape = torch.Size()\n",
    "        else:\n",
    "            batch_shape = self._param.shape\n",
    "        super(NegativeBinomial, self).__init__(batch_shape, validate_args=validate_args)\n",
    "\n",
    "    def _new(self, *args, **kwargs):\n",
    "        return self._param.new(*args, **kwargs)\n",
    "\n",
    "    @property\n",
    "    def mean(self):\n",
    "        return self.total_count * self.probs\n",
    "\n",
    "    @property\n",
    "    def variance(self):\n",
    "        return self.total_count * self.probs * (1 - self.probs)\n",
    "\n",
    "    @lazy_property\n",
    "    def logits(self):\n",
    "        return probs_to_logits(self.probs, is_binary=True)\n",
    "\n",
    "    @lazy_property\n",
    "    def probs(self):\n",
    "        return logits_to_probs(self.logits, is_binary=True)\n",
    "\n",
    "    @property\n",
    "    def param_shape(self):\n",
    "        return self._param.shape\n",
    "\n",
    "    @lazy_property\n",
    "    def _gamma(self):\n",
    "        return torch.distributions.Gamma(concentration=self.total_count,\n",
    "                                         rate=torch.exp(-self.logits))\n",
    "\n",
    "    def sample(self, sample_shape=torch.Size()):\n",
    "        with torch.no_grad():\n",
    "            rate = self._gamma.sample(sample_shape=sample_shape)\n",
    "            return torch.poisson(rate)\n",
    "\n",
    "    def log_prob(self, value):\n",
    "        if self._validate_args:\n",
    "            self._validate_sample(value)\n",
    "\n",
    "        log_unnormalized_prob = (self.total_count * F.logsigmoid(-self.logits)\n",
    "                                 + value * F.logsigmoid(self.logits))\n",
    "\n",
    "        log_normalization = (-torch.lgamma(self.total_count + value)\n",
    "                             + torch.lgamma(self.total_count)\n",
    "                             + torch.lgamma(1. + value))\n",
    "\n",
    "        return log_unnormalized_prob - log_normalization\n",
    "\n",
    "    def expand(self, batch_shape):\n",
    "        try:\n",
    "            return super(NegativeBinomial, self).expand(batch_shape)\n",
    "        except NotImplementedError:\n",
    "            validate_args = self.__dict__.get('_validate_args')\n",
    "            total_count = self.total_count.expand(batch_shape)\n",
    "            if 'probs' in self.__dict__:\n",
    "                probs = self.probs.expand(batch_shape)\n",
    "                return type(self)(total_count, probs=probs, validate_args=validate_args)\n",
    "            else:\n",
    "                logits = self.logits.expand(batch_shape)\n",
    "                return type(self)(total_count, logits=logits, validate_args=validate_args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# class for creating a bunch of fully connected layers\n",
    "# (taken from https://github.com/YosefLab/scVI/blob/master/scvi/models/modules.py)\n",
    "\n",
    "class FCLayers(nn.Module):\n",
    "    \"\"\"\n",
    "    n_in : int \n",
    "        shape of input data\n",
    "    n_out : int\n",
    "        shape of output data\n",
    "    n_layers : int\n",
    "        number of hidden layers\n",
    "    n_hidden : int\n",
    "        nodes in the hidden layers\n",
    "    dropout_rate : float\n",
    "        dropout rate for edges (when training)\n",
    "    \"\"\"\n",
    "    def __init__(self, n_in, n_out, n_layers=1, n_hidden=128, dropout_rate=0.1):\n",
    "        super(FCLayers, self).__init__()\n",
    "        layers_dim = [n_in] + (n_layers - 1) * [n_hidden] + [n_out]\n",
    "        self.fc_layers = nn.Sequential(collections.OrderedDict(\n",
    "            [('Layer {}'.format(i), nn.Sequential(\n",
    "                nn.Linear(n_in, n_out),\n",
    "                nn.BatchNorm1d(n_out, eps=1e-3, momentum=0.99),\n",
    "                nn.ReLU(),\n",
    "                nn.Dropout(p=dropout_rate)\n",
    "            )) for i, (n_in, n_out) in enumerate(zip(layers_dim[:-1], layers_dim[1:]))]))\n",
    "\n",
    "    def forward(self, x):\n",
    "        for layers in self.fc_layers:\n",
    "            for layer in layers:\n",
    "                if isinstance(layer, nn.BatchNorm1d) and x.dim() == 3:\n",
    "                    x = torch.cat([(layer(slice_x)).unsqueeze(0) for slice_x in x], dim=0)\n",
    "                else:\n",
    "                    x = layer(x)\n",
    "        return x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is scvi.models.modules.DecoderSCVI minus some stuff for batch effects\n",
    "# which I took out for simplicity\n",
    "# taken from https://github.com/YosefLab/scVI/blob/master/scvi/models/modules.py\n",
    "\n",
    "\n",
    "class DecoderSCVI(nn.Module):\n",
    "    r\"\"\"Decodes data from latent space of ``n_input`` dimensions ``n_output``\n",
    "    dimensions using a fully-connected neural network of ``n_hidden`` layers.\n",
    "\n",
    "    :param n_input: The dimensionality of the input (latent space)\n",
    "    :param n_output: The dimensionality of the output (data space)\n",
    "    :param n_layers: The number of fully-connected hidden layers\n",
    "    :param n_hidden: The number of nodes per hidden layer\n",
    "    :param dropout_rate: Dropout rate to apply to each of the hidden layers\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, n_input: int, n_output: int, n_layers: int = 1,\n",
    "                 n_hidden: int = 128, dropout_rate: float = 0.1):\n",
    "        super(DecoderSCVI, self).__init__()\n",
    "        self.px_decoder = FCLayers(n_in=n_input, n_out=n_hidden, n_layers=n_layers,\n",
    "                                   n_hidden=n_hidden, dropout_rate=dropout_rate)\n",
    "\n",
    "        # mean gamma\n",
    "        self.px_scale_decoder = nn.Sequential(nn.Linear(n_hidden, n_output), \n",
    "                                              nn.Softmax(dim=-1))\n",
    "\n",
    "        # dispersion: here we only deal with gene-cell dispersion case\n",
    "        self.px_r_decoder = nn.Linear(n_hidden, n_output)\n",
    "\n",
    "    def forward(self, z: torch.Tensor, library: torch.Tensor):\n",
    "        r\"\"\"The forward computation for a single sample.\n",
    "\n",
    "         #. Decodes the data from the latent space using the decoder network\n",
    "         #. Returns parameters for the ZINB distribution of expression\n",
    "         #. If ``dispersion != 'gene-cell'`` then value for that param will be ``None``\n",
    "\n",
    "        :param z: tensor with shape ``(n_input,)``\n",
    "        :param library: library size\n",
    "        :return: parameters for the ZINB distribution of expression\n",
    "        :rtype: 4-tuple of :py:class:`torch.Tensor`\n",
    "        \"\"\"\n",
    "\n",
    "        # The decoder returns values for the parameters of the NB distribution\n",
    "        px = self.px_decoder(z)\n",
    "        px_scale = self.px_scale_decoder(px)\n",
    "        # Clamp to high value: exp(12) ~ 160000 to avoid nans (computational stability)\n",
    "        px_rate = torch.exp(torch.clamp(library, max=12)) * px_scale\n",
    "        px_r = self.px_r_decoder(px)\n",
    "        return px_r, px_rate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is scvi.models.modules.Encoder, again minus batch effects\n",
    "\n",
    "class Encoder(nn.Module):\n",
    "    def __init__(self, n_input, n_output, n_layers=1, n_hidden=128, dropout_rate=0.1):\n",
    "        super(Encoder, self).__init__()\n",
    "        self.encoder = FCLayers(n_in=n_input, n_out=n_hidden, n_layers=n_layers,\n",
    "                                n_hidden=n_hidden, dropout_rate=dropout_rate)\n",
    "        self.mean_encoder = nn.Linear(n_hidden, n_output)\n",
    "        self.var_encoder = nn.Linear(n_hidden, n_output)\n",
    "\n",
    "    def forward(self, x):\n",
    "        # Parameters for latent distribution\n",
    "        q = self.encoder(x)\n",
    "        q_m = self.mean_encoder(q)\n",
    "        q_v = torch.exp(torch.clamp(self.var_encoder(q), -5, 5))  # (computational stability safeguard)\n",
    "        return q_m, q_v\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adapted from the Pyro VAE tutorial (http://pyro.ai/examples/vae.html)\n",
    "# but using the encoder and decoder modules from scVI as in\n",
    "# https://github.com/YosefLab/scVI/blob/master/scvi/models/vae.py\n",
    "\n",
    "class SCVAE(nn.Module):\n",
    "    def __init__(self, n_input, z_dim=8, n_layers=3, hidden_dim=64, use_cuda=False):\n",
    "        super(SCVAE, self).__init__()\n",
    "        # create the encoder and decoder networks\n",
    "        self.z_encoder = Encoder(n_input, z_dim, n_layers=n_layers, n_hidden=hidden_dim)\n",
    "        self.l_encoder = Encoder(n_input, 1, n_layers=n_layers, n_hidden=hidden_dim)\n",
    "        self.decoder = DecoderSCVI(z_dim, n_input, n_layers=n_layers, n_hidden=hidden_dim)\n",
    "\n",
    "        if use_cuda:\n",
    "            # calling cuda() here will put all the parameters of\n",
    "            # the encoder and decoder networks into gpu memory\n",
    "            self.cuda()\n",
    "        self.use_cuda = use_cuda\n",
    "        self.z_dim = z_dim\n",
    "\n",
    "    # define the model p(x|z)p(z)\n",
    "    def model(self, x):\n",
    "        # register PyTorch module `decoder` with Pyro\n",
    "        pyro.module(\"decoder\", self.decoder)\n",
    "        with pyro.iarange(\"data\", x.size(0)):\n",
    "            # setup hyperparameters for prior p(z)\n",
    "            z_loc = x.new_zeros(torch.Size((x.size(0), self.z_dim)))\n",
    "            z_scale = x.new_ones(torch.Size((x.size(0), self.z_dim)))\n",
    "\n",
    "            l_loc = x.new_zeros(torch.Size((x.size(0), 1)))\n",
    "            l_scale = x.new_ones(torch.Size((x.size(0), 1)))\n",
    "\n",
    "            # sample from prior (value will be sampled by guide when computing the ELBO)\n",
    "            z = pyro.sample(\"latent\", dist.Normal(z_loc, z_scale, validate_args=True).independent(1))\n",
    "            l = pyro.sample(\"library\", dist.Normal(l_loc, l_scale, validate_args=True).independent(1))\n",
    " \n",
    "            # decode the latent code z\n",
    "            px_r, px_rate = self.decoder.forward(z, l)\n",
    "\n",
    "            # the scVI model uses a different parameterization of the NB distribution\n",
    "            # so here I am converting to something the PyTorch version can use.\n",
    "        \n",
    "            # exp(px_r) = theta, px_rate = mu\n",
    "            # p = mu / (mu + theta)\n",
    "            # logit = log(p) - log(1 - p) \n",
    "            #       = log(mu) - log(mu + theta) - (log(theta) - log(mu + theta))\n",
    "            #       = log(mu) - log(theta)\n",
    "            px_logit = torch.log(px_rate) - px_r\n",
    "\n",
    "            # score against data\n",
    "            pyro.sample(\"obs\", NegativeBinomial(total_count=torch.exp(px_r), logits=px_logit,\n",
    "                                                validate_args=True).independent(1), obs=x)\n",
    "\n",
    "    # define the guide (i.e. variational distribution) q(z|x)\n",
    "    def guide(self, x):\n",
    "        # register PyTorch module `encoder` with Pyro\n",
    "        pyro.module(\"z_encoder\", self.z_encoder)\n",
    "        pyro.module(\"l_encoder\", self.l_encoder)\n",
    "        with pyro.iarange(\"data\", x.size(0)):\n",
    "            # use the encoder to get the parameters used to define q(z|x)\n",
    "            z_loc, z_scale = self.z_encoder.forward(x)\n",
    "            l_loc, l_scale = self.l_encoder.forward(x)\n",
    "            \n",
    "            # sample the latent code z and library size\n",
    "            pyro.sample(\"latent\", dist.Normal(z_loc, z_scale, validate_args=True).independent(1))\n",
    "            pyro.sample(\"library\", dist.Normal(l_loc, l_scale, validate_args=True).independent(1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing SCVAE on real data\n",
    "\n",
    "Data is available in a variety of places, e.g. s3://czbiohub-tabula-muris/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "droplet_metadata = pd.read_csv('TM_droplet_metadata.csv', index_col=0, dtype=str)\n",
    "tmd = anndata.read_h5ad('TM_droplet_mat.h5ad')\n",
    "\n",
    "def get_droplet_tissue(tissue):\n",
    "    tmd_scvi = GeneExpressionDataset(\n",
    "        *GeneExpressionDataset.get_attributes_from_matrix(tmd.X),\n",
    "        gene_names=np.array(tmd.var.values, dtype=str)\n",
    "    )\n",
    "    \n",
    "    tissue_metadata = droplet_metadata.loc[((droplet_metadata['tissue'] == tissue).values\n",
    "                                            & np.asarray(tmd_scvi.X.sum(1) > 1000).flatten()\n",
    "                                            & np.asarray((tmd_scvi.X > 0).sum(1) > 500).flatten())]\n",
    "\n",
    "    \n",
    "    # downsample to cells that passed QC and were (mostly) annotated\n",
    "    tmd_scvi.update_cells(np.where(droplet_metadata['tissue'] == tissue)[0])\n",
    "    tmd_scvi.update_cells(np.where(tmd_scvi.X.sum(1) > 1000)[0])\n",
    "    tmd_scvi.update_cells(np.where((tmd_scvi.X > 0).sum(1) > 500)[0])\n",
    "    \n",
    "    classes = sorted(lbl if lbl is not np.nan else 'na' \n",
    "                     for lbl in set(tissue_metadata['cell_ontology_class']))\n",
    "    \n",
    "    return tissue_metadata, tmd_scvi, classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downsampling from 70118 to 21384 cells\n",
      "Downsampling from 21384 to 11386 cells\n",
      "Downsampling from 11386 to 11269 cells\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "((11269, 23433),\n",
       " (11269, 8),\n",
       " ['blood cell',\n",
       "  'endothelial cell',\n",
       "  'epithelial cell',\n",
       "  'mesenchymal cell',\n",
       "  'neuroendocrine cell'])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "tissue_metadata, tmd_scvi, classes = get_droplet_tissue('Trachea')\n",
    "\n",
    "tmd_scvi.X.shape, tissue_metadata.shape, classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these were adapted from the Pyro VAE tutorial\n",
    "\n",
    "def train(svi, train_loader, n_train, use_cuda=False):\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",
    "        # if on GPU put mini-batch into CUDA memory\n",
    "        if use_cuda:\n",
    "            xs = (xs[0].cuda(),)\n",
    "        # do ELBO gradient and accumulate loss\n",
    "        epoch_loss += svi.step(xs[0])\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, use_cuda=False):\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",
    "        # if on GPU put mini-batch into CUDA memory\n",
    "        if use_cuda:\n",
    "            xs = (xs[0].cuda(),)\n",
    "        # compute ELBO estimate and accumulate loss\n",
    "        test_loss += svi.evaluate_loss(xs[0])\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": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to stop if I hit NaNs\n",
    "import warnings\n",
    "warnings.simplefilter(\"error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyro.validation_enabled(True)\n",
    "pyro.clear_param_store() # need to do this before every run in a notebook\n",
    "\n",
    "example_indices = np.random.permutation(len(tmd_scvi))\n",
    "n_train = int(0.9 * len(tmd_scvi))  # 90%/10% train/test split\n",
    "n_test = len(tmd_scvi) - n_train\n",
    "\n",
    "data_loader_train = DataLoader(\n",
    "    dataset=tmd_scvi, batch_size=64, pin_memory=True,\n",
    "    sampler=SubsetRandomSampler(example_indices[:n_train]),\n",
    "    collate_fn=tmd_scvi.collate_fn\n",
    ")\n",
    "    \n",
    "data_loader_test = DataLoader(\n",
    "    dataset=tmd_scvi, batch_size=64, pin_memory=True,\n",
    "    sampler=SubsetRandomSampler(example_indices[n_train:]),\n",
    "    collate_fn=tmd_scvi.collate_fn\n",
    ")\n",
    "\n",
    "# setup the VAE\n",
    "vae = SCVAE(tmd_scvi.nb_genes, use_cuda=True, z_dim=10, n_layers=1, hidden_dim=128)\n",
    "\n",
    "# setup the optimizer\n",
    "optimizer = Adam({\"lr\": 1e-4})\n",
    "# optimizer = SGD({'lr': 1e-6, 'momentum': 0.9, 'nesterov': True})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[epoch 000]  average training loss: 14337.5087\n",
      "[epoch 005]  average training loss: 5347.3135\n",
      "[epoch 010]  average training loss: 4987.7193\n",
      "[epoch 015]  average training loss: 4880.3074\n",
      "[epoch 020]  average training loss: 4833.1097\n"
     ]
    }
   ],
   "source": [
    "# setup the inference algorithm\n",
    "\n",
    "# compare to https://github.com/YosefLab/scVI/tree/master/scvi/inference (!)\n",
    "\n",
    "svi = SVI(vae.model, vae.guide, optimizer, loss=Trace_ELBO())\n",
    "\n",
    "train_elbo = []\n",
    "test_elbo = []\n",
    "test_iter = 5\n",
    "\n",
    "# training loop\n",
    "for epoch in range(21):\n",
    "    total_epoch_loss_train = train(svi, data_loader_train, n_train, use_cuda=True)\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, use_cuda=True)\n",
    "        test_elbo.append(-total_epoch_loss_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAF3CAYAAACrEkILAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX9//HXZ2ayQYBAAJGtIsUFUBHR1lprrQtoVXDXulD1K61L61Jx+dmqtX5bl7burVrFavWrUAuoVYtLrVutCoIIIhLRStgyBBISSCaZmfP7Y24ghMk+k5nJvJ+Pxzwyc+65N+dmknnnnHvvueacQ0RERLKHL9UNEBERka6l8BcREckyCn8REZEso/AXERHJMgp/ERGRLKPwFxERyTIKfxERkSyj8BcREckyCn8REZEso/AXERHJMoFUNyBZ+vfv73bbbbdUN0NERKRLLFiwYINzbkBb6nbb8N9tt92YP39+qpshIiLSJczsv22tq2F/ERGRLKPwFxERyTIKfxERkSzTbY/5x1NfX09paSm1tbWpbkq3lp+fz9ChQ8nJyUl1U0REJI6sCv/S0lJ69erFbrvthpmlujndknOO8vJySktLGTFiRKqbIyIicWTVsH9tbS3FxcUK/iQyM4qLizW6IiKSxrIq/AEFfxfQz1hEJL1lXfinUnl5OePGjWPcuHEMGjSIIUOGbHtdV1fXpm2cd955LF++vM3f8+GHH2bAgAHbvs+4ceNYvnw5JSUljBs3bqf6Z599NiNGjGDcuHHstdde3HLLLduWhUIhfvKTnzBy5EhGjRrFlClTWLNmTZvbIiIi6SGrjvm319yFq7lj3nLWVNQwuKiA6RP3ZMr+Qzq8veLiYhYtWgTATTfdRGFhIVddddUOdZxzOOfw+eL/X/boo4+2+/ueddZZ3HXXXTuUlZSUNFv/zjvvZMqUKdTU1LDXXnsxdepUhg0bxjXXXEMoFOKzzz7D7/fzpz/9iZNPPpl333233W0SEZHUUc+/GXMXrua62R+zuqIGB6yuqOG62R8zd+HqhH+vkpISxo4dy49//GPGjx/P2rVrmTZtGhMmTGDMmDHcfPPN2+p++9vfZtGiRYTDYYqKirj22mvZb7/9OPjggykrK0tou2pqajAzevToQVVVFU888QS///3v8fv9AFx44YUAvPHGGwn9viIi3driWXDnWLipKPZ18awub0LW9vx/+fxSPlmzudnlC7+qoC4S3aGspj7C1c8s5qn3v4q7zujBvbnx+DEdas8nn3zCo48+ygMPPADArbfeSr9+/QiHwxx++OGccsopjB49eod1KisrOeyww7j11lu58sormTFjBtdee+1O237yySf517/+te31+++/32JbrrjiCm666SZWrFjBz372M4qLi/nwww8ZMWIEhYWFO9SdMGECS5cu5bDDDuvQfouIpLuEjgIvnkX42Z8QiHgnRVeuir0G2Pe0RDW5VVkb/q1pGvytlXfWyJEjOfDAA7e9fuqpp3jkkUcIh8OsWbOGTz75ZKfwLygo4JhjjgHggAMO4K233oq77XjD/i1pGPavqqri8MMP57jjjsPv98c9kc85pxP8RNLYB889yLAP72CgC1JmA1g1fjoHnvCjDm8v0YdDE729RG+zYRS4pj4CwOqKrfx89iJ8kRAn7DMQomGIRryvjR8RcJEdX0fD1L5wLfmRHa+GCkRq2frSDfRQ+Cdfaz30Q279J6sranYqH1JUwMwfHZzw9vTs2XPb8xUrVnD33Xfz/vvvU1RUxNlnnx330rnc3Nxtz/1+P+FwOKFt6tWrF4cddhhvv/02F154IStXrqS6unqH3v+HH37IqaeemtDvK5Ip0j0IP3juQcYu+DkFVgcGgwjSZ8HP+QA69A/A3A9X8Ys5H1FXX08PIlRVVHP77A3k1e7BMaP7NwnCnYOv6ev3Pl/P6++uZHwkzIG+KP7NUT6Y8xKDS3ZhzKCeuGgYF4mtE42GIRLGNf4aDeMabzsSZl3FFvzrK7jORQjkRPFviZA7J8rn/8yhT54P3Pbvb9EIuNhXc2HMRbBoBJ+Lvfa5CEeE61ngi+LPi+InQsC8DuDfvUc75TdXXrOu/RvrhKwN/9ZMn7jnDv/tARTk+Jk+cc+kf+/NmzfTq1cvevfuzdq1a5k3bx6TJk1K+vdtqr6+nvfff5+rrrqKXr168YMf/IDp06dz//334/P5mDFjBpFIREP+kjGS2yOMnRcEtG+bkTBEQryw8Et+//ePsHCI3a2evMp6npq9lOINwzl0RCFE6nDhEJG6GiL1ISL1tUTqa4nW1RKtD+HCIaL1tbhw7DnhOsZ++VIs+BspsDr2+/D/sWHFDCy6PfAawq9x8DU8jCh+F2EKUab4AX+TfZjnPdrpG8A34m1vqfeI9+NyRhg/UXyE8RNp8rXQ+RiDn4j5CeMjQuxrRaWfcnxEXOPynG3LY1/9RPHjzE/U58dZgK0Rt/P3cX4i+Cnu3YOoefW9h7MAznxNXm9/fl7ZrxlgOx9yXhMtZmj7f4QdpvBvRsMfb6KHo9pi/PjxjB49mrFjx7L77rtzyCGHdGp7TY/5P/jggxQXF/PJJ58wdOj2X7d7770X2H7MPxQKMXHiRE444QQAbr/9dn72s58xatQozIzRo0cze/bsTrVNupHFs+C1m6GyFPoMhSNu6PQxzOSEdZgAETZVbOK22f8mv2Ykk/bqC5E6CNdCuA4iIQjHHtFwiHBdLeFQDZH6WsJe2AY/+JyfulpyA2HyqCOXMLnU0+vZCKv/nYdF67BwCF+0Dl+kDl+0Dn80hD9aj9/V4Y/Wk+Pq8BHrSX4f+H68IHzHewBG7EM73gd3nfMTIpc6AtSRQ50LMMDr8TeV48IsrOixQ/BFaRxgfqJeiDn/9gALbonEDdwIfvbYtQhnfsznx/kCmC+AswD4/OAPgFdmvgD4Y1+fXrBm23a2bcsL5ysn7o35c8AXwOePbcd8AXw+P36fYQZ+n+E3w8zw+wyfwQWPxb+VuwEvX/EdAn4fAZ+RG4h9Dfh95Pp9BPxGwGc7HcZsaRT4nWu+1/wvXDNuumUdV9f/gR6N/inb6nJ5OPdsbmr31jrOnHNd+O26zoQJE9z8+Tv+Eixbtoy99947RS3KLvpZZ5mmJzEBYX8+gcn3xv4BcA4i9V6oNoTrzkG7bXm4lgUr1zN3/kosEtoWrIX+MN8bVcRuRQEidV7P1+v1Oq/XGwvxEL5ICPNC1xcNEakLkUs9udTjt8R87oVcgBA5sbAlQJ3L8V7HAjjkti8LWw71lkvEl0vEcoj484j6cmMPfy4rN4W31W1YL+Rt56h9hmOBfHyBXCw3H18gDwvk48/Nw59TgD83l9xADjl+IyfgI8/vIyfgY+ifD2RXNuzU7rUMoMc1y8jxGwFfLAR9vtbP3WkxCK9tfxAmenvJ2GbTER6IjQL/5qR9OvSP6NyFq3l7zh+4nKcZbOWsccXcxRl8+8SLO925NLMFzrkJbamrnr9IlupQr7q+BqrXQ3XQ+7oetgSpf+tucuKcxORmT8M9eykWqcNoX+AeABzgY6cLkkMlsYANNerhNu7xxoIzQB19YkFKLmHLoSYa2B7MXp2G9UbuWgz+PCyQC16wWk4evkAevpw8fDkFBHLz8eXkEcjNJ5BbwI0vrGDj1vqd2r1L7zzmXHwIuQEfeQEfuYFYz7K1E2NbCq3zzuxYEH4w/mqKGo75e2pcLqUHTOfAgvbfeCvRh0OTcXg10dtM9ChwbL2LOX3eEV0+qtyYwl8kQyTreHUOYSIVpTwxexn91/Rmv7711G9eR3TzOqgO4ttaRk5NkLzacnIj1XG3F3DEHV7GOR4KHR0L6m2Bm7NDaEf9uVggN9abzcnHn5OPPyePD1bXUOcC23vDsT47YNx4/GgKcvwU5PrJz/Fvfx7wU5Tr26EsLxAbJm4pXH9zafvDNeLLjRsy1x2zN4OLCtq9vWQE4YEn/IgPwDvbfwNl1p9VB3T8bP/kBGFiD68ma5uJDOdEb68jNOwvSaGfdWK1Z+gxHIlSuaWWqo3r2LpxDXWV66ivXAfVZfi2lJFTu4GtG9fQz1UwwCrpa/EDvdL1IOiK2EAfgq4PQVdE0BWx2V/EltxiavP6U5ffn2hBMbd8dS5DfTsPL5dG+7PktHfokeunR24sjHvkBra/zvET8Mefayzdh28btpnOZ/tLdtGwv0iKJepDvD4Spby6jv99cRm19fX0pZr+tpkBVsGASAWfzvk7L/yzloJQOYXhjfSJbKIfm+hHFcVxjmtXu3w2Wh9Crg8lbgj/iY6OBTtFbHCxkD/3qIMI9B5Izx6F9MoP0Ds/hyH5AXrlByjMC8QN65tuObv5k5jGDmr3fkP6D982bLO79QglOyj8RRKstUvAauoibKgOxR5VISoryqndtIb6yvVEq9fj21JGXu0GCurK6RPZxACr4DnbTP+8SnIsstP3q6vKYbO/H9W5xdTmDWNdwXhW9xiIFQ7E33sXcot2pWe/XSnsN5hevYso9BlnttCrPvl77Z/HYtz3p3HDnDCXuyYnMX1/Wru31SATwlokUyn8RRJoc209//viMqjfwnCrpD+VsV56tJKvZs9k1txK+kYr6G+x8tFUkmc7nzQWwc+WnL7U9OxPuMcw3t8QoLS+V2wYvmEInj4Eeg/i5WuPo78Z/dvRzuT0qhN/EpPCWiQ5FP5dqLy8nCOOOAKAdevW4ff7GTBgABCbb7/xjH0tmTFjBsceeyyDBu08nHr22Wfzzjvv0KdPHyA2S99bb73Fww8/zJIlS3aa5nfo0KH07dsXv99PJBLh17/+NccffzwAq1at4pJLLmHZsmVEo1EmT57MbbfdRk5O+88STndtGqYPh4hWrad8/So2ri+lasMaaivWEK0qI7A1SEHdBvpGK3jdKinM33lGxqgztuYUUZvXn3DBAFzhGCp7DSKvaBAF/QaT22cQ9BwIhbvgL+hLb5+P3t66HyxczX2zP6Ym3OR49aR9oAPTK6tXLZLdFP4tSfCkJW25pW9bzJgxg/Hjx8cNf9g+N39bvfXWWxQVFbF06VJOOOEEjj/+eJxzTJ48mSuuuIJzzjmHcDjMBRdcwA033MBvfvObdrc54RL13kTC/OO9j3n8pf8wKrKRb/krGFBVSeWczSx9s55+bhOBmg30qC+nZ7QaHzDAezTYTCFVgX7UFvYn3GM3ng0GWFXXyztRrmjbyXIFfQby5nVHUdhMU1qisBaRRFL4N2fxLHj+p7HrmgEqV8VeQ1LuvPTYY49x//33U1dXx7e+9S3uu+8+otEo5513HosWLcI5x7Rp09hll11YtGgRp59+OgUFBe0aMWjN5s2b6du3LwAvv/wyRUVFnHPOOQAEAgHuvvtuRo4cyY033kh+fnMzVHeB1t6baBRqNm6/Dr3JNemx52W46jLYWs4kHJOazKxW5QoIlvdhFX3YZIMI5e+D6zmAnN6xXnrRwCEMHDScgbsOpXdewbYeOsCShat5LM6Q+vWTdrwxU3sprEUkUbI3/F+6FtZ93Pzy0g9is401Vl8Dz14KCx6Lv86gfeCYW9vdlCVLljBnzhz+/e9/EwgEmDZtGk8//TQjR45kw4YNfPxxrJ0VFRUUFRVx7733ct999zFu3Li422uYnhdg33335fHHH2/x+x966KFEo1G++OKLbdP1Ll26lAMOOGCHekVFRQwePJiVK1fudIfBLvXazduDv0F9Dcy9GF7+RSzg3c4nxkX9+dTmFbPJ15c14SK+rB3C6nDj4+h9KPPOeq8lDwPm//xI+vXMbdedC1M5NbSISFtkb/i3pmnwt1beCa+++ioffPABEybELs+sqalh2LBhTJw4keXLl3PZZZdx7LHHcvTRR7dpex0d9v/ss8+YOHEiS5cubfZWvWlxC9/K0vjl0XoYdRSu50A2WhFfhAr5dHM+H27M5d0yP2u35MAWI8dv7DmoF/vs0YexQ/rw5KsrCFbt/L4OLiqguDCvQ01UL11E0ln2hn9rPfQ7x8aGk5vqMwzOeyGhTXHOcf755/OrX/1qp2WLFy/mpZde4p577uFvf/sbDz30UEK/d2N77LEH/fr149NPP2XMmDG88MKO+1lRUcGaNWsYMWJE0trQqkh9bM74yM6XqW3wD+TS9WexdM1mqmpjtzfO9fvYa9deHL5fH/YZEnuM2qWQvMD2Mf6euYGU3cFRRCQVsjf8W3PEDTseVwbIKYiVJ9iRRx7JKaecwmWXXUb//v0pLy9ny5YtFBQUkJ+fz6mnnsqIESP48Y9/DMTO4K+qqkp4O9atW8dXX33F8OHDGTduHNdddx1PPvkkZ511FuFwmCuvvJILL7wwdcf767bCX39IIFJDvfPvcM37VpfLzTWnUFMfZfK4wYwdHOvV77FLL3ID8WeQa6BhehHJNgr/5jSc1JfgW5TGs88++3DjjTdy5JFHEo1GycnJ4YEHHsDv93PBBRdsG2q/7bbbADjvvPP4n//5n2ZP+Gt8zB9gwYIFADzyyCM888wz28obpj8+9NBD8fv91NfX89vf/pb+/WNXjM+dO5dLLrmEm266iWg0ynHHHRd3dKIrlJaWkjPrDPpvXsov6i+g2uVxdWDWtgllbg+fxvPRb/PFJR27/bGG6UUkm2huf0mKzv6snXN8tr6afyxZx4KPF3PDpp8zzIL8ttd0ZlaPY7M3rN9YZ24DKiKS6TS3v2Qk5xwflVbyjyXrmLd0HV9s2MIevlKeyr+NXjkhKiY/zfX7HsmYZm7QomP0IiJto/CXLtHcDHqRqOODLzduC/y1lbUEfMbBI4u5bmwFRy78X3w5BXD2cwwcNBbQMXoRkc5S+EvSxbvRzdXPLObp979iRVk15VvqyAv4+M4eA7jq6D05cu9d6PPVK/DMj2LnWpw9G/p+bYdt6hi9iEjHZV34p8V16t1c0/NI7pi3fIcheoC6SJT3vtjI8fsNZtLYQRy2xwB65nm/jh8+Ds9fBruOg7P+Cj3bc8saERFpTVaFf35+PuXl5RQXF+sfgCRxzlFeXr7tckDnXNxbxza458z9G68Mb/0O/vkrGHkEnPY45HVkJnwREWlJVoX/0KFDKS0tJRgMprop3Vp+fj5Dhw5l/pcbufWlT5utN7ioYPuLaBT+cS28/yDsezqccB8EEnPPAhER2VFWhX9OTk5qZ6fLEivWV3HxUx/xyifrGdgrj1MnDOX5j9ZQWx/dVmeHs/PDIZjzI1g6Bw6+FI76FfhanphHREQ6LqvCX5JrbWUNd72ygr8uWEXP3ADTJ+7JeYfsRo/cAIeM7B//7PzazTDzLPjizVjoH/LTVO+GiEi3p/CXTqvcWs8f3ijhz+98iXNw/iEjuOTwr9O35/Zh+7hn51eXwRMnw/qlMOUBGHdmF7dcRCQ7Kfylw2rrIzz27y+5//USqkJhTtx/CFcetQdD+/ZofeWNK+EvJ0H1evjBTBh1VPIbLCIigMJfOiASdfztw1LufOUz1lbWcvieA7h60l7svWvvtm1gzSJ48hSIRmDq8zC0TbNRiohIgij8pc2cc7y6rIw75n3KZ+ur2W9YEb8/bRwHjyxu+0ZW/guePgsK+sYm7xmwR9LaKyIi8Sn8ZSfxpuId2reAW1/6lPn/3cTu/Xvyx7PGM2nsoPbNl7BkNsyeBv1Hwdl/g96Dk7cTIiLSLIW/7CDeVLxXzlpE1MGAXnn8+sR9OHXCUHL87bwU770H4aVrYPjBcOb/xXr+IiKSEkm9mNrMfmJmy81sqZnd3qj8OjMr8ZZNbFQ+ySsrMbNrG5WPMLP3zGyFmc00M83+kiTxpuKNOuidH+CN6d/lB98Y3r7gdw5euxleuhr2+j6cM1vBLyKSYkkLfzM7HJgM7OucGwP81isfDZwBjAEmAX8wM7+Z+YH7gWOA0cCZXl2A24A7nXOjgE3ABclqd7Zb08xUvFW1YXrktnOgKBKG5y6NTdk7fiqc+hjkFLS+noiIJFUye/4XAbc650IAzrkyr3wy8LRzLuSc+wIoAQ7yHiXOuZXOuTrgaWCyxQ4qfw94xlv/MWBKEtud1YoL4w+q7DAVb1vUbYWZZ8PCJ+A7V8Pxd4NfR5lERNJBMsN/D+BQb7j+DTM70CsfAqxqVK/UK2uuvBiocM6Fm5RLgv1jyTo2bamj6Sl8O0zF2xZbN8JfToTP/gHH/ha+dz3oRkoiImmjU10xM3sVGBRn0fXetvsC3wQOBGaZ2e6wU7YAOOL/I+JaqB+vPdOAaQDDhw9vrfnSyOPvfsmNzy1lv6FFnHLAUP74r893noq3LSpXwxMnxSbxOfXPMEaDNCIi6aZT4e+cO7K5ZWZ2ETDbxW7u/r6ZRYH+xHruwxpVHQqs8Z7HK98AFJlZwOv9N67ftD0PAQ8BTJgwIe4/CLKjaNRx+7zlPPDG5xy59y7ce+b+FOT6OfubX2v/xoLLY7P21VbGLuUb8Z3EN1hERDotmcP+c4kdq8fM9gByiQX5c8AZZpZnZiOAUcD7wAfAKO/M/lxiJwU+5/3z8DpwirfdqcCzSWx31qgLR7ly1iIeeONzzvrGcB44ezwFuf6ObWzV+zBjIkTq4LwXFPwiImksmWdgzQBmmNkSoA6Y6gX5UjObBXwChIFLnHMRADO7FJgH+IEZzrml3rauAZ42s1uAhcAjSWx3VthcW89FTyzgnZJypk/ck4u/O7J9E/Y09tk8mDUVeu8am7Wvn26bLCKSziyWx93PhAkT3Pz581PdjLS0rrKWHz76PiVl1dx28r6cfMDQjm9s4ZPw3E9g0D5w1jNQOCBxDRURkTYzswXOuTbdLEXXXmWZFeurmDrjfSpr6pnxwwP5zh4dDGvn4J274NWbYPfvwulPQF6vBLZURESSReGfRd5bWc6Fj88nL8fPzB8dzNghfTq2oWgU5v0/eO+PMPZkmPIABDTpoohIplD4Z4kXFq/lipmLGNqvgMfOO4hh/Xp0bEPhOph7ESx5Br5xEUz8NfiSOku0iIgkmMI/Czzy9hfc8sInjB/el4fPnUDfnh3spYeqYOY5sPJ1OPImOORyTd4jIpKBFP7dWDTq+PWLy3j47S+YOGYX7j5jf/JzOngpX3UQnjwF1n0Mk++H/c9ObGNFRKTLKPy7qVA4ws9mfcTfF69l6sFf44bjx+D3dbCXvvGL2Kx9m9fCGf8He05KbGNFRKRLKfy7ocqaeqY9Pp/3vtjItcfsxY++s3vHr+FfuzjW4w+HYOpzMOygxDZWRES6nMK/G5i7cDV3zFvOmooaBvbOw4DyLXXcfcY4Jo/rxD2QvngLnv5B7BK+85+DgXslrM0iIpI6Cv8MN3fhaq6b/TE19REA1m8OAXDxd0d2LviXzoXZF0K/3WPz9PfpxERAIiKSVnSNVoa7Y97ybcHf2LOL4t77qG3e/xP89YcweH847yUFv4hIN6Oef4ZbU1HTrvIWOQev/xrevB32mASnPAq5HZwPQERE0pZ6/hlu16L8uOWDiwrat6FIGP5+eSz49z8bTn9SwS8i0k0p/DPchOF9dyoryPEzfeKebd9IfQ38dSos+DMc+jM44T7wa1BIRKS70id8BluyupKXlq5j36G9Ka+uY01FLYOLCpg+cU+m7N/Gk/1qKuCpM+Grd+GY2+EbP0puo0VEJOUU/hmqtj7C5TMX0bdHLo+d942OTdm7eQ08cTJsWAGnPBK7SY+IiHR7Cv8MdetLn1JSVs1fLjioY8Ef/Cw2a1/NJjjrrzDy8MQ3UkRE0pLCPwP9a3kZf/73l5x/yAgOHTWg/RsonQ9Pngo+P/zwBRg8LvGNFBGRtKUT/jLMxi11TH9mMXvsUsjVk9pxUl+DFa/AY8dDfm84f56CX0QkCyn8M4hzjmv/tpjKrfXcdXoH7tC36Cl46gwoHgnnvxz7KiIiWUfhn0FmzV/Fy5+sZ/rEPRk9uHf7Vn7nHpj7Y/jat+CHL0KvXZLTSBERSXs65p8hvtywhV8+/wkH717MBd8e0fYVo1F45Rfw7n0w5kQ48UEI5CWvoSIikvYU/hkgHIly+cxFBHzG707bD5+vjbfnjdTDs5fA4plw0DSYdBv4NNgjIpLtFP4Z4L7XS1i0qoJ7z9y/7dP2hqph1rnw+WvwvZ/DoVeBtfGfBhER6dYU/mnuw682ce8/Szhx/yEcv9/gtq20ZUPsUr61i+D4e+CAqcltpIiIZBSFfxrbEgpzxcxFDOqdzy8nj2nbSpv+G5u8p7I0dnOevY5NbiNFRCTjKPzT2K/+/glfbdzKzGkH0zs/p/UV1i2JTdcbroFz5sLXDk5+I0VEJOPo7K809Y8l63j6g1VcdNhIDhrRr/UVvnwbHj0WzAfn/UPBLyIizVL4p6GyzbVcN3sxY4f05vIj92h9hU+eg7+cFLt2/4KXYZfRyW+kiIhkLIV/mnHOcdUzi6mpj3DX6fuTG2jlLZo/A/46FXbdNzZdb9GwrmmoiIhkLIV/mnn83f/y5mdBrj92b74+sLD5is7Bv26Dv18BXz8Szn0WerTh8ICIiGQ9nfCXRlasr+LXLy7j8D0HcPY3v9Z8xWgEXrwq1uvf7wdwwj3gb8MJgSIiIij800ZdOMplTy+iMC/A7afshzU3IU99Lcz+H1j2PBxyORx5kybvERGRdlH4p4nfvbKcT9Zu5k/nTmBAr2bm3q+thKd+AP99Gyb+Bg6+uGsbKSIi3YLCPw28+3k5D725kjMPGs5Ro5u5217VOnjiFAgug5Mehn1P7dpGiohIt6HwT7HKmnp+NmsRuxX35BfH7R2/0oYSeOJE2FIOP5gFXz+iaxspIiLdisI/xW54dgnrq0L87aJv0SM3ztuxekFsnn6AHz4PQw7o2gaKiEi3o/BPgbkLV3PHvOWsrqgB4Jixgxg3rGjniiWvwcxzoGcxnD0H+n+9i1sqIiLdka7z72JzF67mutkfbwt+gNeXlzF34eodKy7+K/zfadBvd7jgFQW/iIgkjMK/i90xbzk19ZEdymrro9wxb/n2gnfvj13ON+ybcN4L0GtQF7dSRES6Mw37d7E1jXr8O5U7B6/eCO/cDXufACf9CXLyu7iFIiLS3ann38UGFxXELR/WJwfmXhwL/gkXwKl/VvA2lLHOAAAc+0lEQVSLiEhSKPy72PSJe5Lr33FGvn459cwquhc++j84/Hr4/u/A509RC0VEpLtT+HexKfsP4YRxQwAwYO8+9bza/3cMKnsbjrsLDrta0/WKiEhS6Zh/CgwpKsAMPps+lpz/OwU2/RdOexz2Pj7VTRMRkSyg8E+BYHWIAwvWkfPoz6BuK5wzB3Y7JNXNEhGRLKHwT4E+wfk84n4Brhec/xLsMibVTRIRkSyiY/5d7dMXuXzN1VT7+8IFLyv4RUSkyyn8u9KCx2DmWZTY1/jDiPuh79dS3SIREclCCv+u4By8cQc8/1Pc7odzVt319OjbzK17RUREkkzhn2zRCLw4HV6/BfY9nc1TnqAiksuAXnmpbpmIiGQpnfCXTOEQzJ4Gn8yFgy+Fo35FcMMWAIW/iIikjMI/WWo3w8yz4Is34ahfwSE/BaCsKgQo/EVEJHUU/slQtR6ePBnKlsGJD8J+Z2xbFPTCf6DCX0REUkThn2jln8MTJ0F1GZw5E0YducPihvAfUKib9oiISGok7YQ/MxtnZv8xs0VmNt/MDvLKzczuMbMSM1tsZuMbrTPVzFZ4j6mNyg8ws4+9de4xS9PJ79csghkTY0P+U5/fKfghNrtfrt9H7wL93yUiIqmRzLP9bwd+6ZwbB9zgvQY4BhjlPaYBfwQws37AjcA3gIOAG82sr7fOH726DetNSmK7O+bz1+HP34dAfmzynqET4lYLVoUY0CuPdP3/RUREur9khr8DenvP+wBrvOeTgcddzH+AIjPbFZgIvOKc2+ic2wS8AkzylvV2zr3rnHPA48CUJLa7/T5+Bp48FYqGx4K//6hmqwarQvTX8X4REUmhZI49Xw7MM7PfEvsn41te+RBgVaN6pV5ZS+WlccrTw38egH9cA8O/BWc+BQVFLVYPVoUY2rdHFzVORERkZ50KfzN7FRgUZ9H1wBHAFc65v5nZacAjwJHEbmPflOtAebz2TCN2eIDhw4e32v5OcQ5euxne/j3sdRyc/DDkFLS62obqEPsP79tqPRERkWTpVPg753Y+o81jZo8Dl3kv/wo87D0vBYY1qjqU2CGBUuC7Tcr/5ZUPjVM/XnseAh4CmDBhQtx/EBIiEoa/XwYLn4ADfgjf/z34/K2uFo5EKd9Sp2v8RUQkpZJ5zH8NcJj3/HvACu/5c8C53ln/3wQqnXNrgXnA0WbW1zvR72hgnresysy+6Z3lfy7wbBLb3bK6rbHJexY+AYddA8fd1abgByjfUodzmuBHRERSK5nH/C8E7jazAFCLNxwPvAgcC5QAW4HzAJxzG83sV8AHXr2bnXMbvecXAX8GCoCXvEfXWDwrNrxfWQq9B4M/DzZ9Ad//HRz4P+3a1PZr/BX+IiKSOkkLf+fc28ABccodcEkz68wAZsQpnw+MTXQbW7V4Fjz/U6ivib3evDr29aAftTv4odHsfr0V/iIikjq6q19LXrt5e/A3tvzFDm1OPX8REUkHCv+WVJa2r7wVwWrd1EdERFJP4d+SPkPbV96KYFWIXvkB8nPadoKgiIhIMij8W3LEDTtfu59TECvvgIapfUVERFJJ4d+SfU+D4++BPsMAi309/p5YeQcEq0I63i8iIimnW8u1Zt/TOhz2TQWrQ4wZ3Lv1iiIiIkmknn8X0rC/iIikA4V/F9laF6Y6FFb4i4hIyin8u8iGqjpA1/iLiEjqKfy7SLC6FtA1/iIiknoK/y6ybXY/hb+IiKSYwr+LlCn8RUQkTSj8u0iwKoTPoLinwl9ERFJL4d9FglUhigvz8Pss1U0REZEsp/DvIprdT0RE0oXCv4sEqzXBj4iIpAeFfxfR7H4iIpIuFP5dIBp1bFDPX0RE0oTCvwtU1tRTH3E65i8iImlB4d8FgtW6xl9ERNKHwr8LaHY/ERFJJwr/LqDwFxGRdKLw7wIKfxERSScK/y4QrA6RF/DRKy+Q6qaIiIgo/LtC2eZaBvTKw0xT+4qISOop/LtAsDrEQA35i4hImlD4dwHN7iciIulE4d8FFP4iIpJOFP5JVheOsmlrPQMK81PdFBEREUDhn3TlW3SZn4iIpBeFf5LpGn8REUk3Cv8kU/iLiEi6UfgnmcJfRETSjcI/yRrCv39hbopbIiIiEqPwT7JgdYg+BTnkBfypboqIiAig8E86XeMvIiLpRuGfZMGqEAMKFf4iIpI+FP5JFqxWz19ERNKLwj+JnHOUbdZNfUREJL0o/JNoS12EmvqIev4iIpJWFP5JpGv8RUQkHSn8k0jhLyIi6Ujhn0QKfxERSUcK/yQKVtUC6FI/ERFJKwr/JApWh/D7jL49NLWviIikD4V/EgWrQvQvzMXns1Q3RUREZBuFfxJpal8REUlHCv8kClZral8REUk/Cv8kUs9fRETSkcI/SaJRx4bqOoW/iIikHYV/kmzaWkck6jTsLyIiaUfhnyRl3gQ/A3vnp7glIiIiO1L4J4lm9xMRkXTVqfA3s1PNbKmZRc1sQpNl15lZiZktN7OJjconeWUlZnZto/IRZvaema0ws5lmluuV53mvS7zlu3WmzV1lW/hr2F9ERNJMZ3v+S4CTgDcbF5rZaOAMYAwwCfiDmfnNzA/cDxwDjAbO9OoC3Abc6ZwbBWwCLvDKLwA2Oee+Dtzp1Ut7wWr1/EVEJD11Kvydc8ucc8vjLJoMPO2cCznnvgBKgIO8R4lzbqVzrg54GphsZgZ8D3jGW/8xYEqjbT3mPX8GOMKrn9aCVSF65PrpmRdIdVNERER2kKxj/kOAVY1el3plzZUXAxXOuXCT8h225S2v9OqnNV3jLyIi6arVbqmZvQoMirPoeufcs82tFqfMEf+fDddC/Za2tfM3NZsGTAMYPnx4M03rGsEqze4nIiLpqdXwd84d2YHtlgLDGr0eCqzxnscr3wAUmVnA6903rt+wrVIzCwB9gI3NtPUh4CGACRMmxP0HoasEq0OMGliYyiaIiIjElaxh/+eAM7wz9UcAo4D3gQ+AUd6Z/bnETgp8zjnngNeBU7z1pwLPNtrWVO/5KcA/vfppTcP+IiKSrjp7qd+JZlYKHAy8YGbzAJxzS4FZwCfAP4BLnHMRr1d/KTAPWAbM8uoCXANcaWYlxI7pP+KVPwIUe+VXAtsuD0xXoXCEypp6DfuLiEha6tSp6M65OcCcZpb9L/C/ccpfBF6MU76S2NUATctrgVM7086utqG6DtBlfiIikp40w18SaHY/ERFJZwr/JFD4i4hIOlP4J0FD+A/spZv6iIhI+lH4J0FZVS0AxYW5KW6JiIjIzhT+SRCsCtGvZy45fv14RUQk/SidkkCz+4mISDpT+CdBsFoT/IiISPpS+CeBZvcTEZF0pvBPMOecwl9ERNKawj/BqkJhQuGojvmLiEjaUvgnmCb4ERGRdKfwTzCFv4iIpDuFf4Ip/EVEJN0p/BNsW/jrmL+IiKQphX+CBatD5PiNoh45qW6KiIhIXAr/BGuY3c/MUt0UERGRuBT+CaZr/EVEJN0p/BOsTOEvIiJpTuGfYOr5i4hIulP4J1Ak6ti4RXf0ExGR9KbwT6DyLSGiTtf4i4hIelP4J5Am+BERkUyg8E8ghb+IiGQChX8CbZ/dLz/FLREREWmewj+BgtWx8O/fKzfFLREREWmewj+BglUhCvMC9MgNpLopIiIizVL4J5Cu8RcRkUyg8E+ghnn9RURE0pnCP4GC1SEG9Fb4i4hIelP4J5B6/iIikgkU/glSWx+hqjasY/4iIpL2FP4Jogl+REQkUyj8E6RM4S8iIhlC4Z8g22f3U/iLiEh6U/gnSMPsfgPV8xcRkTSn8E+QYFUIM+jXU1P7iohIelP4J0iwKkRxz1wCfv1IRUQkvSmpEiRYFaK/jveLiEgGUPgnSLBa8/qLiEhmUPgnyAbd1EdERDKEwj8BnHO6o5+IiGQMhX8CbK4JUxeJMrBXfqqbIiIi0iqFfwIEq2sBze4nIiKZQeGfAGWa3U9ERDKIwj8BdFMfERHJJAr/BFD4i4hIJlH4J0CwKkRuwEfv/ECqmyIiItIqhX8CBKtCDCjMw8xS3RQREZFWKfwTQLP7iYhIJlH4J4Am+BERkUyi8E8Ahb+IiGQShX8n1UeibNxap2v8RUQkYyj8O2njljqc02V+IiKSOToV/mZ2qpktNbOomU1oVH6UmS0ws4+9r99rtOwAr7zEzO4x7xR5M+tnZq+Y2Qrva1+v3Lx6JWa22MzGd6bNiaZr/EVEJNN0tue/BDgJeLNJ+QbgeOfcPsBU4C+Nlv0RmAaM8h6TvPJrgdecc6OA17zXAMc0qjvNWz9tNIT/QIW/iIhkiE6Fv3NumXNueZzyhc65Nd7LpUC+meWZ2a5Ab+fcu845BzwOTPHqTQYe854/1qT8cRfzH6DI205aUM9fREQyTVcc8z8ZWOicCwFDgNJGy0q9MoBdnHNrAbyvA73yIcCqZtZJuWB1LPz764Q/ERHJEK3OR2tmrwKD4iy63jn3bCvrjgFuA45uKIpTzbXWhLauY2bTiB0aYPjw4a1sNjHKNtfSOz9Afo6/S76fiIhIZ7Ua/s65IzuyYTMbCswBznXOfe4VlwJDG1UbCjQcHlhvZrs659Z6w/pljdYZ1sw6Tdv6EPAQwIQJE1r7pyIhNLufiIhkmqQM+5tZEfACcJ1z7p2Gcm84v8rMvumd5X8u0DB68ByxkwPxvjYuP9c76/+bQGXD4YF0oAl+REQk03T2Ur8TzawUOBh4wczmeYsuBb4O/MLMFnmPhmP4FwEPAyXA58BLXvmtwFFmtgI4ynsN8CKw0qv/J+DizrQ50WLhn5/qZoiIiLRZp+5B65ybQ2xov2n5LcAtzawzHxgbp7wcOCJOuQMu6Uw7k6nhjn4iIiKZQjP8dcKWUJgtdREN+4uISEZR+HfChmpd4y8iIplH4d8JmuBHREQykcK/E7aFv475i4hIBlH4d0JQw/4iIpKBFP6dEKwK4fcZ/XrmpropIiIibabw74RgVYjinrn4ffFmIBYREUlPCv9O0Ox+IiKSiRT+naB5/UVEJBMp/DuhbLNm9xMRkcyj8O+gaNSxQT1/ERHJQAr/DqqoqSccdQp/ERHJOAr/DtLsfiIikqkU/h2k2f1ERCRTKfw7KFhdC6jnLyIimUfh30Ea9hcRkUyl8O+gYFWI/BwfhXmBVDdFRESkXRT+HdQwu5+ZpvYVEZHMovDvoGB1iIG98lPdDBERkXZT+HdQsEqz+4mISGZS+HeQbuojIiKZSuHfAXXhKJu21iv8RUQkIyn8O6B8iy7zExGRzKXw74CyzZrdT0REMpfCvwM0wY+IiGQyhX8HBKsV/iIikrkU/h3Q0PMvLsxNcUtERETaT+HfAcGqEEU9csgL+FPdFBERkXZT+HeAJvgREZFMpvDvgGC1JvgREZHMpfDvAM3uJyIimUzh307OOYJVIQYq/EVEJEMp/NtpS12EmvqIev4iIpKxFP7tpAl+REQk0yn822lb+Bfmp7glIiIiHaPwbyf1/EVEJNMp/NuprKoWUPiLiEjmUvi3U7AqRMBnFBXkpLopIiIiHaLwb6dgVYj+hXn4fJbqpoiIiHSIwr+dNLufiIhkOoV/O2l2PxERyXQK/3bSTX1ERCTTKfzbIRJ1lG+pU89fREQymsK/HTZtrSMSdQp/ERHJaAr/dmiY4Ec39RERkUym8G8Hze4nIiLdgcK/HRT+IiLSHSj82yFYHQv//jrbX0REMpjCvx2CVSF65vrpmRdIdVNEREQ6TOHfDprgR0REugOFfzuUVdUq/EVEJOMp/NtBPX8REekOFP7toKl9RUSkO1D4t1FtfYTNtWH1/EVEJON1KvzN7FQzW2pmUTObEGf5cDOrNrOrGpVNMrPlZlZiZtc2Kh9hZu+Z2Qozm2lmuV55nve6xFu+W2fa3FEbqnWNv4iIdA+d7fkvAU4C3mxm+Z3ASw0vzMwP3A8cA4wGzjSz0d7i24A7nXOjgE3ABV75BcAm59zXve3d1sk2d4gm+BERke6iU+HvnFvmnFseb5mZTQFWAksbFR8ElDjnVjrn6oCngclmZsD3gGe8eo8BU7znk73XeMuP8Op3qW3hX5jf1d9aREQkoZJyzN/MegLXAL9ssmgIsKrR61KvrBiocM6Fm5TvsI63vNKrH+/7TjOz+WY2PxgMJmJXtmmY3W9gb/X8RUQks7Ua/mb2qpktifOY3MJqvyQ2hF/ddHNx6roWyltaZ+dC5x5yzk1wzk0YMGBAC81rv2BVCDPo1zM3odsVERHpaq3OU+ucO7ID2/0GcIqZ3Q4UAVEzqwUWAMMa1RsKrAE2AEVmFvB69w3lEBsFGAaUmlkA6ANs7ECbOiVYFaJfj1xy/LpAQkREMltSJql3zh3a8NzMbgKqnXP3eeE9ysxGAKuBM4AfOOecmb0OnELsPICpwLPeJp7zXr/rLf+ncy5uzz+ZNMGPiIh0F5291O9EMysFDgZeMLN5LdX3evWXAvOAZcAs51zDCYHXAFeaWQmxY/qPeOWPAMVe+ZXAtaRAsFrhLyIi3UOnev7OuTnAnFbq3NTk9YvAi3HqrSR2NUDT8lrg1M60MxGCVSFGFPdMdTNEREQ6TQew28A5R5mG/UVEpJtQ+LfB5towdeGowl9ERLoFhX8r5i5czVG/fwOA+18vYe7C1SlukYiISOck5Wz/7mLuwtVcN/tjauojAGzaWs91sz8GYMr+Q1paVUREJG2p59+CO+Yt3xb8DWrqI9wxL+6MxiIiIhlB4d+CNRU17SoXERHJBAr/FgwuKmhXuYiISCZQ+Ldg+sQ9Kcjx71BWkONn+sQ9U9QiERGRztMJfy1oOKnvjnnLWVNRw+CiAqZP3FMn+4mISEZT+Ldiyv5DFPYiItKtaNhfREQkyyj8RUREsozCX0REJMso/EVERLKMwl9ERCTLKPxFRESyjMJfREQkyyj8RUREsozCX0REJMso/EVERLKMOedS3YakMLMg8N8EbrI/sCGB20sl7Ut66i770l32A7Qv6ai77Ackfl++5pwb0JaK3Tb8E83M5jvnJqS6HYmgfUlP3WVfust+gPYlHXWX/YDU7ouG/UVERLKMwl9ERCTLKPzb7qFUNyCBtC/pqbvsS3fZD9C+pKPush+Qwn3RMX8REZEso56/iIhIllH4N2Fmk8xsuZmVmNm1cZbnmdlMb/l7ZrZb17eydWY2zMxeN7NlZrbUzC6LU+e7ZlZpZou8xw2paGtbmNmXZvax1875cZabmd3jvS+LzWx8KtrZEjPbs9HPepGZbTazy5vUSdv3xMxmmFmZmS1pVNbPzF4xsxXe177NrDvVq7PCzKZ2Xavja2Zf7jCzT73fnzlmVtTMui3+Lna1ZvblJjNb3ej36Nhm1m3x864rNbMfMxvtw5dmtqiZddPtPYn7+ZtWfy/OOT28B+AHPgd2B3KBj4DRTepcDDzgPT8DmJnqdjezL7sC473nvYDP4uzLd4G/p7qtbdyfL4H+LSw/FngJMOCbwHupbnMr++MH1hG7Ljcj3hPgO8B4YEmjstuBa73n1wK3xVmvH7DS+9rXe943DfflaCDgPb8t3r54y1r8XUyTfbkJuKqV9Vr9vEv1fjRZ/jvghgx5T+J+/qbT34t6/js6CChxzq10ztUBTwOTm9SZDDzmPX8GOMLMrAvb2CbOubXOuQ+951XAMmBIaluVVJOBx13Mf4AiM9s11Y1qwRHA5865RE5ElVTOuTeBjU2KG/89PAZMibPqROAV59xG59wm4BVgUtIa2gbx9sU597JzLuy9/A8wtMsb1gHNvC9t0ZbPuy7T0n54n7GnAU91aaM6qIXP37T5e1H472gIsKrR61J2DsxtdbwPikqguEta10HeoYn9gffiLD7YzD4ys5fMbEyXNqx9HPCymS0ws2lxlrflvUsnZ9D8B1mmvCcAuzjn1kLsAw8YGKdOpr03AOcTG0mKp7XfxXRxqXcIY0Yzw8uZ9L4cCqx3zq1oZnnavidNPn/T5u9F4b+jeD34ppdDtKVO2jCzQuBvwOXOuc1NFn9IbNh5P+BeYG5Xt68dDnHOjQeOAS4xs+80WZ4x74uZ5QInAH+NsziT3pO2ypj3BsDMrgfCwJPNVGntdzEd/BEYCYwD1hIbMm8qk96XM2m515+W70krn7/NrhanLOHvi8J/R6XAsEavhwJrmqtjZgGgDx0bcks6M8sh9ov3pHNudtPlzrnNzrlq7/mLQI6Z9e/iZraJc26N97UMmENsyLKxtrx36eIY4EPn3PqmCzLpPfGsbzi84n0ti1MnY94b7+Sq44CznHcAtqk2/C6mnHNuvXMu4pyLAn8ifhsz4n3xPmdPAmY2Vycd35NmPn/T5u9F4b+jD4BRZjbC652dATzXpM5zQMPZl6cA/2zuQyKVvGNkjwDLnHO/b6bOoIbzFczsIGK/D+Vd18q2MbOeZtar4TmxE7OWNKn2HHCuxXwTqGwYXktDzfZiMuU9aaTx38NU4Nk4deYBR5tZX2/4+WivLK2Y2STgGuAE59zWZuq05Xcx5Zqc73Ii8dvYls+7dHAk8KlzrjTewnR8T1r4/E2fv5dUnxWZbg9iZ41/Ruws2Ou9spuJfSAA5BMbri0B3gd2T3Wbm9mPbxMbKloMLPIexwI/Bn7s1bkUWErsLN//AN9Kdbub2ZfdvTZ+5LW34X1pvC8G3O+9bx8DE1Ld7mb2pQexMO/TqCwj3hNi/7CsBeqJ9U4uIHa+y2vACu9rP6/uBODhRuue7/3NlADnpem+lBA71trw99JwVc9g4MWWfhfTcF/+4v0dLCYWOLs23Rfv9U6fd+m0H175nxv+PhrVTff3pLnP37T5e9EMfyIiIllGw/4iIiJZRuEvIiKSZRT+IiIiWUbhLyIikmUU/iIiIllG4S+SwcysuNFdz9Y1uZNbbhu38aiZ7dlKnUvM7KwEtflt705yDe1sdvKWDm6/1Jq5I5+IxOhSP5FuwsxuAqqdc79tUm7E/tajKWlYE2b2NnCpcy7u7VkTsP1SYKxzriIZ2xfpDtTzF+mGzOzrZrbEzB4gdr+AXc3sITOb791f/IZGdd82s3FmFjCzCjO71bux0LtmNtCrc4uZXd6o/q1m9r7Xg/+WV97TzP7mrfuU973GtaPNT5jZH83sLTP7zMyO8coLzOwxi92v/cOGedu99t7p7ediM7u40eYuN7OFXvkenf6BinQzCn+R7ms08Ihzbn/n3Gpi9xGfAOwHHGVmo+Os0wd4w8VuLPQusZnG4jHn3EHAdKDhH4mfAOu8dW8ldiez5sxsNOx/a6PyYcBhwPHAQ2aWB/wUqHPO7QOcA/zFO6RxEbGZ3vZzzu1L7Ja0DdY75/YHHgaubKEdIlkpkOoGiEjSfO6c+6DR6zPN7AJif/eDif1z8EmTdWqccw23sl1A7Faq8cxuVGc37/m3gdsAnHMfmdnSFtp2ejPD/rO8wxPLzWwVMMrb7h3edpea2Rrg68TmfL/LORfxljW+wVbj9h3bQjtEspLCX6T72tLwxMxGAZcBBznnKszsCWL3qWiqrtHzCM1/RoTi1Il3K9L2anoSkmthuxanfoN47RMRj4b9RbJDb6AK2Ozd8W1iEr7H28BpAGa2D7GRhfY61bsz4x7EDgGsAN4EzvK2uzewK7EbnrwMXGRmfm9Zv07vgUiW0H/EItnhQ2JD/EuAlcA7Sfge9wKPm9li7/stASqbqTvTzGq85+udcw3/jJQQC/uBwDTnXJ2Z3Qs8aGYfE7vj27le+YPEDgssNrMw8EfggSTsl0i3o0v9RCQhzCwABJxztd5hhpeBUc65cBvXfwJ4xjk3N5ntFBH1/EUkcQqB17x/Agz4UVuDX0S6lnr+IiIiWUYn/ImIiGQZhb+IiEiWUfiLiIhkGYW/iIhIllH4i4iIZBmFv4iISJb5/zfGQyf+svUxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_llk(np.array(train_elbo), np.array(test_elbo), test_iter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "default_view": {},
   "name": "Untitled",
   "provenance": [],
   "version": "0.3.2",
   "views": {}
  },
  "kernel_info": {
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.5"
  },
  "notify_time": "30",
  "nteract": {
   "version": "0.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}