not-MIWAE-demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# not-MIWAE: Deep Generative Modelling with Missing not at Random Data\n",
    "This notebook illustrates how to fit a *deep latent variable model* to data affected by a missing process which depends on the missing data itself, i.e. *missing not at random*.\n",
    "\n",
    "We fit a linear PPCA-like model to a relatively small UCI dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preamble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:516: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint8 = np.dtype([(\"qint8\", np.int8, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:517: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint8 = np.dtype([(\"quint8\", np.uint8, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:518: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint16 = np.dtype([(\"qint16\", np.int16, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:519: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint16 = np.dtype([(\"quint16\", np.uint16, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:520: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint32 = np.dtype([(\"qint32\", np.int32, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:525: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  np_resource = np.dtype([(\"resource\", np.ubyte, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:541: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint8 = np.dtype([(\"qint8\", np.int8, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:542: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint8 = np.dtype([(\"quint8\", np.uint8, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:543: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint16 = np.dtype([(\"qint16\", np.int16, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:544: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_quint16 = np.dtype([(\"quint16\", np.uint16, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:545: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  _np_qint32 = np.dtype([(\"qint32\", np.int32, 1)])\n",
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorboard/compat/tensorflow_stub/dtypes.py:550: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.\n",
      "  np_resource = np.dtype([(\"resource\", np.ubyte, 1)])\n",
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "import keras\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "\n",
    "# import sys\n",
    "# sys.path.append('./')\n",
    "\n",
    "import os\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"6\"\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"serif\"\n",
    "plt.rcParams['font.size'] = 15.0\n",
    "plt.rcParams['axes.spines.right'] = False\n",
    "plt.rcParams['axes.spines.top'] = False\n",
    "plt.rcParams['savefig.format'] = 'pdf'\n",
    "plt.rcParams['lines.linewidth'] = 2.5\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data\n",
    "Here we use the white-wine dataset from the UCI database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-white.csv\"\n",
    "data = np.array(pd.read_csv(url, low_memory=False, sep=';'))\n",
    "# ---- drop the classification attribute\n",
    "data = data[:, :-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N, D = data.shape\n",
    "n_latent = D - 1\n",
    "n_hidden = 128\n",
    "n_samples = 20\n",
    "max_iter = 30000\n",
    "batch_size = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ---- standardize data\n",
    "data = data - np.mean(data, axis=0)\n",
    "data = data / np.std(data, axis=0)\n",
    "\n",
    "# ---- random permutation \n",
    "p = np.random.permutation(N)\n",
    "data = data[p, :]\n",
    "\n",
    "# ---- we use the full dataset for training here, but you can make a train-val split\n",
    "Xtrain = data.copy()\n",
    "Xval = Xtrain.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduce missing \n",
    "Here we denote\n",
    "- Xnan: data matrix with np.nan as the missing entries\n",
    "- Xz: data matrix with 0 as the missing entries\n",
    "- S: missing mask \n",
    "\n",
    "The missing process depends on the missing data itself:\n",
    "- in half the features, set the feature value to missing when it is higher than the feature mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ---- introduce missing process\n",
    "Xnan = Xtrain.copy()\n",
    "Xz = Xtrain.copy()\n",
    "\n",
    "mean = np.mean(Xnan[:, :int(D / 2)], axis=0)\n",
    "ix_larger_than_mean = Xnan[:, :int(D / 2)] > mean\n",
    "\n",
    "Xnan[:, :int(D / 2)][ix_larger_than_mean] = np.nan\n",
    "Xz[:, :int(D / 2)][ix_larger_than_mean] = 0\n",
    "\n",
    "S = np.array(~np.isnan(Xnan), dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build the model\n",
    "The model we are building has a Gaussian prior and a Gaussian observation model,\n",
    "\n",
    "$$ p(\\mathbf{z}) = \\mathcal{N}(\\mathbf{z} | \\mathbf{0}, \\mathbf{I})$$\n",
    "\n",
    "$$ p(\\mathbf{x} | \\mathbf{z}) = \\mathcal{N}(\\mathbf{x} | \\mathbf{\\mu}_{\\theta}(\\mathbf{z}), \\sigma^2\\mathbf{I})$$\n",
    "\n",
    "$$ p(\\mathbf{x}) = \\int p(\\mathbf{x} | \\mathbf{z})p(\\mathbf{z}) d\\mathbf{z}$$\n",
    "\n",
    "where $\\mathbf{\\mu}_{\\theta}(\\mathbf{z}): \\mathbb{R}^d \\rightarrow \\mathbb{R}^p $ in general is a deep neural net, but in this case is a linear mapping, $\\mathbf{\\mu} = \\mathbf{Wz + b}$.\n",
    "\n",
    "The variational posterior is also Gaussian\n",
    "\n",
    "$$q_{\\gamma}(\\mathbf{z} | \\mathbf{x}) = \\mathcal{N}(\\mathbf{z} | \\mu_{\\gamma}(\\mathbf{x}), \\sigma_{\\gamma}(\\mathbf{x})^2 \\mathbf{I})$$\n",
    "\n",
    "If the missing process is *missing at random*, it is ignorable and the ELBO becomes, as described in [the MIWAE paper](https://arxiv.org/abs/1812.02633)\n",
    "\n",
    "$$ E_{\\mathbf{z}_1...\\mathbf{z}_K} \\left[ \\log \\frac{1}{K}\\sum_{k=1}^K \\frac{p_{\\theta}(\\mathbf{x^o} | \\mathbf{z}_k)p(\\mathbf{z}_k)}{q_{\\gamma}(\\mathbf{z}_k | \\mathbf{x^o})} \\right] $$\n",
    "\n",
    "When the missing process is MNAR it is non-ignorable and we need to include the missing model. In this example we include the missing model as a logistic regression in each feature dimension\n",
    "\n",
    "$$ p_{\\phi}(\\mathbf{s} | \\mathbf{x^o, x^m}) = \\text{Bern}(\\mathbf{s} | \\pi_{\\phi}(\\mathbf{x^o, x^m}))$$\n",
    "\n",
    "$$ \\pi_{\\phi, j}(x_j) = \\frac{1}{1 + e^{-\\text{logits}_j}} $$\n",
    "\n",
    "$$ \\text{logits}_j = W_j (x_j - b_j) $$\n",
    "\n",
    "The ELBO in the MNAR case becomes\n",
    "\n",
    "$$ E_{(\\mathbf{z}_1, \\mathbf{x}_1^m)...(\\mathbf{z}_K, \\mathbf{x}_K^m)} \\left[ \\log \\frac{1}{K} \\sum_{k=1}^K \\frac{p_{\\phi}(\\mathbf{s} | \\mathbf{x}^o, \\mathbf{x}_k^m) p_{\\theta}(\\mathbf{x}^o | \\mathbf{z}_k) p(\\mathbf{z}_k)}{q_{\\gamma}(\\mathbf{z} | \\mathbf{x}^o)} \\right]$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inputs\n",
    "Let's first define the inputs of the model\n",
    "- x_pl: data input\n",
    "- s_pl: mask input\n",
    "- n_pl: number of importance samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating graph...\n"
     ]
    }
   ],
   "source": [
    "print(\"Creating graph...\")\n",
    "tf.reset_default_graph()\n",
    "# ---- input\n",
    "with tf.variable_scope('input'):\n",
    "    x_pl = tf.placeholder(tf.float32, [None, D], 'x_pl')\n",
    "    s_pl = tf.placeholder(tf.float32, [None, D], 's_pl')\n",
    "    n_pl = tf.placeholder(tf.int32, shape=(), name='n_pl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the noise variance is learned as a shared parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/ops/init_ops.py:1251: calling VarianceScaling.__init__ (from tensorflow.python.ops.init_ops) with dtype is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Call initializer instance with the dtype argument instead of passing it to the constructor\n"
     ]
    }
   ],
   "source": [
    "# ---- parameters\n",
    "with tf.variable_scope('data_process'):\n",
    "    logstd = tf.get_variable('logstd', shape=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encoder\n",
    "The encoder / inference network consists of two hidden layers with 128 units and tanh activation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py:3636: The name tf.random_uniform is deprecated. Please use tf.random.uniform instead.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x = keras.layers.Dense(units=n_hidden, activation=tf.nn.tanh, name='l_enc1')(x_pl)\n",
    "x = keras.layers.Dense(units=n_hidden, activation=tf.nn.tanh, name='l_enc2')(x)\n",
    "\n",
    "q_mu = keras.layers.Dense(units=n_latent, activation=None, name='q_mu')(x)\n",
    "\n",
    "q_logstd = keras.layers.Dense(units=n_latent, activation=lambda x: tf.clip_by_value(x, -10, 10),\n",
    "                           name='q_logstd')(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variational distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q_z = tfp.distributions.Normal(loc=q_mu, scale=tf.exp(q_logstd))\n",
    "\n",
    "# ---- sample the latent value\n",
    "l_z = q_z.sample(n_pl)                    # shape [n_samples, batch_size, dl]\n",
    "l_z = tf.transpose(l_z, perm=[1, 0, 2])   # shape [batch_size, n_samples, dl]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = keras.layers.Dense(units=D, activation=None, name='mu')(l_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observation model / likelihood function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_x_given_z = tfp.distributions.Normal(loc=mu, scale=tf.exp(logstd))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing model\n",
    "- first mix observed data and samples of missing data\n",
    "- feed through missing model\n",
    "- find likelihood of missing model parameters\n",
    "\n",
    "We have to expand the dimensions of x_pl and s_pl, since mu has size [batch, n_samples, D]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l_out_mixed = mu * tf.expand_dims(1 - s_pl, axis=1) + tf.expand_dims(x_pl * s_pl, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = tf.get_variable('W', shape=[1, 1, D])\n",
    "W = -tf.nn.softplus(W)\n",
    "b = tf.get_variable('b', shape=[1, 1, D])\n",
    "\n",
    "logits = W * (l_out_mixed - b)\n",
    "\n",
    "p_s_given_x = tfp.distributions.Bernoulli(logits=logits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluating likelihoods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/tensorflow/python/ops/nn_impl.py:180: add_dispatch_support.<locals>.wrapper (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use tf.where in 2.0, which has the same broadcast rule as np.where\n"
     ]
    }
   ],
   "source": [
    "# ---- evaluate the observed data in p(x|z)\n",
    "log_p_x_given_z = tf.reduce_sum(tf.expand_dims(s_pl, axis=1) * \n",
    "                                p_x_given_z.log_prob(tf.expand_dims(x_pl, axis=1)), axis=-1)  # sum over d-dimension\n",
    "\n",
    "# --- evaluate the z-samples in q(z|x)\n",
    "q_z2 = tfp.distributions.Normal(loc=tf.expand_dims(q_z.loc, axis=1), scale=tf.expand_dims(q_z.scale, axis=1))\n",
    "log_q_z_given_x = tf.reduce_sum(q_z2.log_prob(l_z), axis=-1)\n",
    "\n",
    "# ---- evaluate the z-samples in the prior p(z)\n",
    "prior = tfp.distributions.Normal(loc=0.0, scale=1.0)\n",
    "log_p_z = tf.reduce_sum(prior.log_prob(l_z), axis=-1)\n",
    "\n",
    "# ---- evaluate the mask in p(s|x)\n",
    "log_p_s_given_x = tf.reduce_sum(p_s_given_x.log_prob(tf.expand_dims(s_pl, axis=1)), axis=-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Losses for the MIWAE and not-MIWAE respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lpxz = log_p_x_given_z\n",
    "lpz = log_p_z\n",
    "lqzx = log_q_z_given_x\n",
    "lpsx = log_p_s_given_x\n",
    "\n",
    "# ---- MIWAE\n",
    "# ---- importance weights\n",
    "l_w = lpxz + lpz - lqzx\n",
    "\n",
    "# ---- sum over samples\n",
    "log_sum_w = tf.reduce_logsumexp(l_w, axis=1)\n",
    "\n",
    "# ---- average over samples\n",
    "log_avg_weight = log_sum_w - tf.log(tf.cast(n_pl, tf.float32))\n",
    "\n",
    "# ---- average over minibatch to get the average llh\n",
    "MIWAE = tf.reduce_mean(log_avg_weight, axis=-1)\n",
    "\n",
    "\n",
    "# ---- not-MIWAE\n",
    "# ---- importance weights\n",
    "l_w = lpxz + lpsx + lpz - lqzx\n",
    "\n",
    "# ---- sum over samples\n",
    "log_sum_w = tf.reduce_logsumexp(l_w, axis=1)\n",
    "\n",
    "# ---- average over samples\n",
    "log_avg_weight = log_sum_w - tf.log(tf.cast(n_pl, tf.float32))\n",
    "\n",
    "# ---- average over minibatch to get the average llh\n",
    "notMIWAE = tf.reduce_mean(log_avg_weight, axis=-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ---- training stuff\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.allow_growth = True\n",
    "sess = tf.Session(config=config)\n",
    "global_step = tf.Variable(initial_value=0, trainable=False)\n",
    "optimizer = tf.train.AdamOptimizer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choose wether you want to train the MIWAE or the notMIWAE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "loss = -notMIWAE\n",
    "# loss = -MIWAE\n",
    "\n",
    "tvars = tf.trainable_variables()\n",
    "train_op = optimizer.minimize(loss, global_step=global_step, var_list=tvars)\n",
    "sess.run(tf.global_variables_initializer())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do the training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0/30000 updates, 0.48 s, 21.12 train_loss, 21.61 val_loss\n",
      "100/30000 updates, 0.39 s, 20.42 train_loss, 20.17 val_loss\n",
      "200/30000 updates, 0.30 s, 19.16 train_loss, 19.21 val_loss\n",
      "300/30000 updates, 0.31 s, 18.26 train_loss, 18.44 val_loss\n",
      "400/30000 updates, 0.33 s, 17.41 train_loss, 17.78 val_loss\n",
      "500/30000 updates, 0.34 s, 16.48 train_loss, 17.30 val_loss\n",
      "600/30000 updates, 0.37 s, 16.70 train_loss, 16.89 val_loss\n",
      "700/30000 updates, 0.35 s, 16.92 train_loss, 16.59 val_loss\n",
      "800/30000 updates, 0.38 s, 16.20 train_loss, 16.22 val_loss\n",
      "900/30000 updates, 0.38 s, 15.79 train_loss, 15.92 val_loss\n",
      "1000/30000 updates, 0.35 s, 14.78 train_loss, 15.69 val_loss\n",
      "1100/30000 updates, 0.39 s, 14.25 train_loss, 15.48 val_loss\n",
      "1200/30000 updates, 0.37 s, 15.40 train_loss, 15.24 val_loss\n",
      "1300/30000 updates, 0.37 s, 14.81 train_loss, 15.07 val_loss\n",
      "1400/30000 updates, 0.37 s, 15.50 train_loss, 14.83 val_loss\n",
      "1500/30000 updates, 0.37 s, 14.10 train_loss, 14.63 val_loss\n",
      "1600/30000 updates, 0.33 s, 14.95 train_loss, 14.48 val_loss\n",
      "1700/30000 updates, 0.33 s, 13.51 train_loss, 14.29 val_loss\n",
      "1800/30000 updates, 0.32 s, 14.09 train_loss, 14.07 val_loss\n",
      "1900/30000 updates, 0.33 s, 13.76 train_loss, 13.94 val_loss\n",
      "2000/30000 updates, 0.42 s, 13.85 train_loss, 13.76 val_loss\n",
      "2100/30000 updates, 0.36 s, 13.07 train_loss, 13.58 val_loss\n",
      "2200/30000 updates, 0.35 s, 13.96 train_loss, 13.50 val_loss\n",
      "2300/30000 updates, 0.36 s, 12.50 train_loss, 13.31 val_loss\n",
      "2400/30000 updates, 0.37 s, 12.24 train_loss, 13.35 val_loss\n",
      "2500/30000 updates, 0.33 s, 11.59 train_loss, 13.07 val_loss\n",
      "2600/30000 updates, 0.37 s, 11.78 train_loss, 12.90 val_loss\n",
      "2700/30000 updates, 0.36 s, 11.92 train_loss, 12.84 val_loss\n",
      "2800/30000 updates, 0.36 s, 11.85 train_loss, 12.72 val_loss\n",
      "2900/30000 updates, 0.35 s, 12.07 train_loss, 12.64 val_loss\n",
      "3000/30000 updates, 0.34 s, 12.45 train_loss, 12.56 val_loss\n",
      "3100/30000 updates, 0.35 s, 12.11 train_loss, 12.58 val_loss\n",
      "3200/30000 updates, 0.35 s, 12.46 train_loss, 12.43 val_loss\n",
      "3300/30000 updates, 0.35 s, 12.40 train_loss, 12.37 val_loss\n",
      "3400/30000 updates, 0.35 s, 11.16 train_loss, 12.35 val_loss\n",
      "3500/30000 updates, 0.35 s, 11.91 train_loss, 12.33 val_loss\n",
      "3600/30000 updates, 0.35 s, 12.16 train_loss, 12.26 val_loss\n",
      "3700/30000 updates, 0.30 s, 11.99 train_loss, 12.25 val_loss\n",
      "3800/30000 updates, 0.31 s, 11.46 train_loss, 12.20 val_loss\n",
      "3900/30000 updates, 0.33 s, 12.72 train_loss, 12.15 val_loss\n",
      "4000/30000 updates, 0.33 s, 11.84 train_loss, 12.15 val_loss\n",
      "4100/30000 updates, 0.32 s, 11.51 train_loss, 12.13 val_loss\n",
      "4200/30000 updates, 0.30 s, 12.25 train_loss, 12.08 val_loss\n",
      "4300/30000 updates, 0.30 s, 11.74 train_loss, 12.04 val_loss\n",
      "4400/30000 updates, 0.29 s, 11.79 train_loss, 12.05 val_loss\n",
      "4500/30000 updates, 0.29 s, 12.74 train_loss, 12.04 val_loss\n",
      "4600/30000 updates, 0.34 s, 10.75 train_loss, 11.97 val_loss\n",
      "4700/30000 updates, 0.35 s, 14.03 train_loss, 12.00 val_loss\n",
      "4800/30000 updates, 0.35 s, 11.48 train_loss, 11.99 val_loss\n",
      "4900/30000 updates, 0.36 s, 11.93 train_loss, 11.92 val_loss\n",
      "5000/30000 updates, 0.35 s, 11.72 train_loss, 11.92 val_loss\n",
      "5100/30000 updates, 0.34 s, 11.51 train_loss, 11.96 val_loss\n",
      "5200/30000 updates, 0.35 s, 11.25 train_loss, 11.93 val_loss\n",
      "5300/30000 updates, 0.35 s, 11.69 train_loss, 11.94 val_loss\n",
      "5400/30000 updates, 0.34 s, 11.19 train_loss, 11.85 val_loss\n",
      "5500/30000 updates, 0.36 s, 11.04 train_loss, 11.83 val_loss\n",
      "5600/30000 updates, 0.36 s, 13.52 train_loss, 11.87 val_loss\n",
      "5700/30000 updates, 0.33 s, 12.29 train_loss, 11.85 val_loss\n",
      "5800/30000 updates, 0.34 s, 11.87 train_loss, 11.98 val_loss\n",
      "5900/30000 updates, 0.33 s, 11.53 train_loss, 11.85 val_loss\n",
      "6000/30000 updates, 0.31 s, 11.93 train_loss, 11.84 val_loss\n",
      "6100/30000 updates, 0.30 s, 11.99 train_loss, 11.90 val_loss\n",
      "6200/30000 updates, 0.32 s, 11.76 train_loss, 11.82 val_loss\n",
      "6300/30000 updates, 0.30 s, 11.58 train_loss, 11.79 val_loss\n",
      "6400/30000 updates, 0.33 s, 11.36 train_loss, 11.77 val_loss\n",
      "6500/30000 updates, 0.32 s, 10.63 train_loss, 11.88 val_loss\n",
      "6600/30000 updates, 0.34 s, 10.65 train_loss, 11.77 val_loss\n",
      "6700/30000 updates, 0.33 s, 11.25 train_loss, 11.77 val_loss\n",
      "6800/30000 updates, 0.31 s, 12.00 train_loss, 11.78 val_loss\n",
      "6900/30000 updates, 0.33 s, 11.67 train_loss, 11.74 val_loss\n",
      "7000/30000 updates, 0.36 s, 10.09 train_loss, 11.76 val_loss\n",
      "7100/30000 updates, 0.35 s, 11.45 train_loss, 11.79 val_loss\n",
      "7200/30000 updates, 0.37 s, 12.06 train_loss, 11.77 val_loss\n",
      "7300/30000 updates, 0.36 s, 10.94 train_loss, 11.76 val_loss\n",
      "7400/30000 updates, 0.38 s, 11.65 train_loss, 11.75 val_loss\n",
      "7500/30000 updates, 0.35 s, 11.33 train_loss, 11.70 val_loss\n",
      "7600/30000 updates, 0.35 s, 11.86 train_loss, 11.74 val_loss\n",
      "7700/30000 updates, 0.35 s, 12.85 train_loss, 11.73 val_loss\n",
      "7800/30000 updates, 0.35 s, 10.48 train_loss, 11.71 val_loss\n",
      "7900/30000 updates, 0.35 s, 11.18 train_loss, 11.71 val_loss\n",
      "8000/30000 updates, 0.34 s, 11.15 train_loss, 11.73 val_loss\n",
      "8100/30000 updates, 0.34 s, 13.02 train_loss, 11.68 val_loss\n",
      "8200/30000 updates, 0.35 s, 10.79 train_loss, 11.67 val_loss\n",
      "8300/30000 updates, 0.36 s, 13.34 train_loss, 11.68 val_loss\n",
      "8400/30000 updates, 0.34 s, 13.99 train_loss, 11.70 val_loss\n",
      "8500/30000 updates, 0.31 s, 11.74 train_loss, 11.70 val_loss\n",
      "8600/30000 updates, 0.28 s, 10.22 train_loss, 11.64 val_loss\n",
      "8700/30000 updates, 0.29 s, 11.77 train_loss, 11.68 val_loss\n",
      "8800/30000 updates, 0.30 s, 10.11 train_loss, 11.68 val_loss\n",
      "8900/30000 updates, 0.34 s, 12.05 train_loss, 11.68 val_loss\n",
      "9000/30000 updates, 0.29 s, 12.50 train_loss, 11.66 val_loss\n",
      "9100/30000 updates, 0.29 s, 12.39 train_loss, 11.66 val_loss\n",
      "9200/30000 updates, 0.30 s, 11.85 train_loss, 11.62 val_loss\n",
      "9300/30000 updates, 0.29 s, 12.91 train_loss, 11.63 val_loss\n",
      "9400/30000 updates, 0.29 s, 11.84 train_loss, 11.67 val_loss\n",
      "9500/30000 updates, 0.35 s, 10.99 train_loss, 11.59 val_loss\n",
      "9600/30000 updates, 0.34 s, 12.21 train_loss, 11.59 val_loss\n",
      "9700/30000 updates, 0.35 s, 11.39 train_loss, 11.62 val_loss\n",
      "9800/30000 updates, 0.35 s, 11.04 train_loss, 11.65 val_loss\n",
      "9900/30000 updates, 0.35 s, 11.22 train_loss, 11.63 val_loss\n",
      "10000/30000 updates, 0.35 s, 11.95 train_loss, 11.60 val_loss\n",
      "10100/30000 updates, 0.34 s, 11.24 train_loss, 11.64 val_loss\n",
      "10200/30000 updates, 0.35 s, 11.70 train_loss, 11.60 val_loss\n",
      "10300/30000 updates, 0.35 s, 11.33 train_loss, 11.62 val_loss\n",
      "10400/30000 updates, 0.37 s, 12.13 train_loss, 11.60 val_loss\n",
      "10500/30000 updates, 0.37 s, 10.70 train_loss, 11.60 val_loss\n",
      "10600/30000 updates, 0.31 s, 12.48 train_loss, 11.58 val_loss\n",
      "10700/30000 updates, 0.33 s, 12.51 train_loss, 11.58 val_loss\n",
      "10800/30000 updates, 0.35 s, 11.89 train_loss, 11.63 val_loss\n",
      "10900/30000 updates, 0.34 s, 10.81 train_loss, 11.59 val_loss\n",
      "11000/30000 updates, 0.33 s, 9.97 train_loss, 11.57 val_loss\n",
      "11100/30000 updates, 0.32 s, 11.09 train_loss, 11.56 val_loss\n",
      "11200/30000 updates, 0.30 s, 13.57 train_loss, 11.57 val_loss\n",
      "11300/30000 updates, 0.37 s, 11.21 train_loss, 11.56 val_loss\n",
      "11400/30000 updates, 0.37 s, 10.18 train_loss, 11.55 val_loss\n",
      "11500/30000 updates, 0.36 s, 12.49 train_loss, 11.56 val_loss\n",
      "11600/30000 updates, 0.37 s, 10.99 train_loss, 11.53 val_loss\n",
      "11700/30000 updates, 0.39 s, 12.05 train_loss, 11.52 val_loss\n",
      "11800/30000 updates, 0.38 s, 11.04 train_loss, 11.61 val_loss\n",
      "11900/30000 updates, 0.39 s, 12.42 train_loss, 11.52 val_loss\n",
      "12000/30000 updates, 0.40 s, 12.25 train_loss, 11.59 val_loss\n",
      "12100/30000 updates, 0.35 s, 10.90 train_loss, 11.52 val_loss\n",
      "12200/30000 updates, 0.36 s, 10.75 train_loss, 11.56 val_loss\n",
      "12300/30000 updates, 0.35 s, 10.97 train_loss, 11.54 val_loss\n",
      "12400/30000 updates, 0.30 s, 10.65 train_loss, 11.53 val_loss\n",
      "12500/30000 updates, 0.28 s, 10.38 train_loss, 11.59 val_loss\n",
      "12600/30000 updates, 0.30 s, 11.06 train_loss, 11.54 val_loss\n",
      "12700/30000 updates, 0.31 s, 11.24 train_loss, 11.52 val_loss\n",
      "12800/30000 updates, 0.28 s, 12.97 train_loss, 11.55 val_loss\n",
      "12900/30000 updates, 0.29 s, 10.63 train_loss, 11.53 val_loss\n",
      "13000/30000 updates, 0.31 s, 13.13 train_loss, 11.50 val_loss\n",
      "13100/30000 updates, 0.32 s, 11.46 train_loss, 11.54 val_loss\n",
      "13200/30000 updates, 0.31 s, 13.50 train_loss, 11.60 val_loss\n",
      "13300/30000 updates, 0.31 s, 11.45 train_loss, 11.48 val_loss\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13400/30000 updates, 0.31 s, 11.72 train_loss, 11.53 val_loss\n",
      "13500/30000 updates, 0.33 s, 11.52 train_loss, 11.53 val_loss\n",
      "13600/30000 updates, 0.32 s, 11.57 train_loss, 11.52 val_loss\n",
      "13700/30000 updates, 0.33 s, 11.00 train_loss, 11.51 val_loss\n",
      "13800/30000 updates, 0.28 s, 11.71 train_loss, 11.50 val_loss\n",
      "13900/30000 updates, 0.30 s, 10.65 train_loss, 11.49 val_loss\n",
      "14000/30000 updates, 0.29 s, 10.99 train_loss, 11.51 val_loss\n",
      "14100/30000 updates, 0.31 s, 12.16 train_loss, 11.50 val_loss\n",
      "14200/30000 updates, 0.36 s, 13.33 train_loss, 11.48 val_loss\n",
      "14300/30000 updates, 0.35 s, 12.58 train_loss, 11.49 val_loss\n",
      "14400/30000 updates, 0.35 s, 11.10 train_loss, 11.48 val_loss\n",
      "14500/30000 updates, 0.34 s, 11.47 train_loss, 11.53 val_loss\n",
      "14600/30000 updates, 0.34 s, 10.57 train_loss, 11.52 val_loss\n",
      "14700/30000 updates, 0.34 s, 11.30 train_loss, 11.49 val_loss\n",
      "14800/30000 updates, 0.34 s, 11.46 train_loss, 11.59 val_loss\n",
      "14900/30000 updates, 0.35 s, 13.10 train_loss, 11.48 val_loss\n",
      "15000/30000 updates, 0.35 s, 10.55 train_loss, 11.46 val_loss\n",
      "15100/30000 updates, 0.33 s, 11.07 train_loss, 11.53 val_loss\n",
      "15200/30000 updates, 0.35 s, 10.54 train_loss, 11.54 val_loss\n",
      "15300/30000 updates, 0.34 s, 14.11 train_loss, 11.53 val_loss\n",
      "15400/30000 updates, 0.31 s, 12.01 train_loss, 11.48 val_loss\n",
      "15500/30000 updates, 0.31 s, 10.67 train_loss, 11.47 val_loss\n",
      "15600/30000 updates, 0.32 s, 11.49 train_loss, 11.49 val_loss\n",
      "15700/30000 updates, 0.34 s, 10.82 train_loss, 11.48 val_loss\n",
      "15800/30000 updates, 0.33 s, 11.69 train_loss, 11.51 val_loss\n",
      "15900/30000 updates, 0.32 s, 11.62 train_loss, 11.52 val_loss\n",
      "16000/30000 updates, 0.30 s, 10.77 train_loss, 11.46 val_loss\n",
      "16100/30000 updates, 0.36 s, 11.54 train_loss, 11.54 val_loss\n",
      "16200/30000 updates, 0.36 s, 12.74 train_loss, 11.50 val_loss\n",
      "16300/30000 updates, 0.37 s, 12.50 train_loss, 11.47 val_loss\n",
      "16400/30000 updates, 0.36 s, 11.43 train_loss, 11.44 val_loss\n",
      "16500/30000 updates, 0.40 s, 11.18 train_loss, 11.47 val_loss\n",
      "16600/30000 updates, 0.38 s, 12.55 train_loss, 11.48 val_loss\n",
      "16700/30000 updates, 0.38 s, 11.25 train_loss, 11.48 val_loss\n",
      "16800/30000 updates, 0.37 s, 11.14 train_loss, 11.46 val_loss\n",
      "16900/30000 updates, 0.36 s, 11.56 train_loss, 11.48 val_loss\n",
      "17000/30000 updates, 0.33 s, 10.68 train_loss, 11.49 val_loss\n",
      "17100/30000 updates, 0.35 s, 12.05 train_loss, 11.48 val_loss\n",
      "17200/30000 updates, 0.29 s, 11.99 train_loss, 11.51 val_loss\n",
      "17300/30000 updates, 0.32 s, 11.03 train_loss, 11.48 val_loss\n",
      "17400/30000 updates, 0.32 s, 10.50 train_loss, 11.47 val_loss\n",
      "17500/30000 updates, 0.29 s, 10.54 train_loss, 11.46 val_loss\n",
      "17600/30000 updates, 0.28 s, 11.98 train_loss, 11.44 val_loss\n",
      "17700/30000 updates, 0.30 s, 10.63 train_loss, 11.46 val_loss\n",
      "17800/30000 updates, 0.30 s, 10.81 train_loss, 11.42 val_loss\n",
      "17900/30000 updates, 0.29 s, 11.29 train_loss, 11.43 val_loss\n",
      "18000/30000 updates, 0.28 s, 13.02 train_loss, 11.47 val_loss\n",
      "18100/30000 updates, 0.30 s, 11.59 train_loss, 11.45 val_loss\n",
      "18200/30000 updates, 0.34 s, 12.23 train_loss, 11.63 val_loss\n",
      "18300/30000 updates, 0.35 s, 11.36 train_loss, 11.47 val_loss\n",
      "18400/30000 updates, 0.33 s, 11.09 train_loss, 11.47 val_loss\n",
      "18500/30000 updates, 0.35 s, 10.48 train_loss, 11.49 val_loss\n",
      "18600/30000 updates, 0.34 s, 10.94 train_loss, 11.53 val_loss\n",
      "18700/30000 updates, 0.34 s, 11.69 train_loss, 11.44 val_loss\n",
      "18800/30000 updates, 0.34 s, 10.81 train_loss, 11.48 val_loss\n",
      "18900/30000 updates, 0.34 s, 10.67 train_loss, 11.47 val_loss\n",
      "19000/30000 updates, 0.37 s, 11.37 train_loss, 11.51 val_loss\n",
      "19100/30000 updates, 0.35 s, 11.24 train_loss, 11.45 val_loss\n",
      "19200/30000 updates, 0.34 s, 11.76 train_loss, 11.44 val_loss\n",
      "19300/30000 updates, 0.32 s, 11.25 train_loss, 11.40 val_loss\n",
      "19400/30000 updates, 0.27 s, 11.19 train_loss, 11.43 val_loss\n",
      "19500/30000 updates, 0.29 s, 12.46 train_loss, 11.47 val_loss\n",
      "19600/30000 updates, 0.29 s, 10.96 train_loss, 11.43 val_loss\n",
      "19700/30000 updates, 0.30 s, 12.31 train_loss, 11.43 val_loss\n",
      "19800/30000 updates, 0.30 s, 10.63 train_loss, 11.47 val_loss\n",
      "19900/30000 updates, 0.31 s, 10.90 train_loss, 11.39 val_loss\n",
      "20000/30000 updates, 0.32 s, 10.80 train_loss, 11.42 val_loss\n",
      "20100/30000 updates, 0.30 s, 11.30 train_loss, 11.47 val_loss\n",
      "20200/30000 updates, 0.29 s, 11.04 train_loss, 11.42 val_loss\n",
      "20300/30000 updates, 0.30 s, 11.88 train_loss, 11.44 val_loss\n",
      "20400/30000 updates, 0.32 s, 10.56 train_loss, 11.44 val_loss\n",
      "20500/30000 updates, 0.33 s, 11.09 train_loss, 11.44 val_loss\n",
      "20600/30000 updates, 0.34 s, 10.56 train_loss, 11.48 val_loss\n",
      "20700/30000 updates, 0.36 s, 12.55 train_loss, 11.43 val_loss\n",
      "20800/30000 updates, 0.30 s, 12.01 train_loss, 11.44 val_loss\n",
      "20900/30000 updates, 0.32 s, 12.08 train_loss, 11.42 val_loss\n",
      "21000/30000 updates, 0.33 s, 13.42 train_loss, 11.42 val_loss\n",
      "21100/30000 updates, 0.34 s, 12.96 train_loss, 11.45 val_loss\n",
      "21200/30000 updates, 0.36 s, 11.21 train_loss, 11.44 val_loss\n",
      "21300/30000 updates, 0.34 s, 11.43 train_loss, 11.50 val_loss\n",
      "21400/30000 updates, 0.36 s, 12.21 train_loss, 11.43 val_loss\n",
      "21500/30000 updates, 0.37 s, 12.73 train_loss, 11.47 val_loss\n",
      "21600/30000 updates, 0.36 s, 11.61 train_loss, 11.49 val_loss\n",
      "21700/30000 updates, 0.34 s, 13.35 train_loss, 11.41 val_loss\n",
      "21800/30000 updates, 0.35 s, 10.78 train_loss, 11.43 val_loss\n",
      "21900/30000 updates, 0.37 s, 11.35 train_loss, 11.41 val_loss\n",
      "22000/30000 updates, 0.37 s, 11.90 train_loss, 11.51 val_loss\n",
      "22100/30000 updates, 0.37 s, 11.11 train_loss, 11.40 val_loss\n",
      "22200/30000 updates, 0.35 s, 11.83 train_loss, 11.40 val_loss\n",
      "22300/30000 updates, 0.32 s, 10.86 train_loss, 11.48 val_loss\n",
      "22400/30000 updates, 0.30 s, 11.11 train_loss, 11.40 val_loss\n",
      "22500/30000 updates, 0.29 s, 11.13 train_loss, 11.41 val_loss\n",
      "22600/30000 updates, 0.28 s, 11.61 train_loss, 11.41 val_loss\n",
      "22700/30000 updates, 0.29 s, 11.32 train_loss, 11.40 val_loss\n",
      "22800/30000 updates, 0.27 s, 10.27 train_loss, 11.45 val_loss\n",
      "22900/30000 updates, 0.31 s, 12.08 train_loss, 11.43 val_loss\n",
      "23000/30000 updates, 0.31 s, 10.58 train_loss, 11.40 val_loss\n",
      "23100/30000 updates, 0.36 s, 12.18 train_loss, 11.41 val_loss\n",
      "23200/30000 updates, 0.35 s, 11.56 train_loss, 11.40 val_loss\n",
      "23300/30000 updates, 0.35 s, 11.79 train_loss, 11.41 val_loss\n",
      "23400/30000 updates, 0.35 s, 11.01 train_loss, 11.44 val_loss\n",
      "23500/30000 updates, 0.36 s, 11.54 train_loss, 11.45 val_loss\n",
      "23600/30000 updates, 0.34 s, 12.40 train_loss, 11.42 val_loss\n",
      "23700/30000 updates, 0.35 s, 12.33 train_loss, 11.53 val_loss\n",
      "23800/30000 updates, 0.34 s, 11.29 train_loss, 11.43 val_loss\n",
      "23900/30000 updates, 0.34 s, 11.70 train_loss, 11.41 val_loss\n",
      "24000/30000 updates, 0.34 s, 12.01 train_loss, 11.44 val_loss\n",
      "24100/30000 updates, 0.35 s, 11.45 train_loss, 11.58 val_loss\n",
      "24200/30000 updates, 0.30 s, 11.36 train_loss, 11.47 val_loss\n",
      "24300/30000 updates, 0.28 s, 11.84 train_loss, 11.45 val_loss\n",
      "24400/30000 updates, 0.29 s, 10.89 train_loss, 11.48 val_loss\n",
      "24500/30000 updates, 0.30 s, 11.49 train_loss, 11.44 val_loss\n",
      "24600/30000 updates, 0.28 s, 11.06 train_loss, 11.49 val_loss\n",
      "24700/30000 updates, 0.29 s, 11.02 train_loss, 11.42 val_loss\n",
      "24800/30000 updates, 0.30 s, 13.28 train_loss, 11.39 val_loss\n",
      "24900/30000 updates, 0.29 s, 10.11 train_loss, 11.42 val_loss\n",
      "25000/30000 updates, 0.30 s, 10.70 train_loss, 11.40 val_loss\n",
      "25100/30000 updates, 0.35 s, 10.35 train_loss, 11.40 val_loss\n",
      "25200/30000 updates, 0.34 s, 11.19 train_loss, 11.42 val_loss\n",
      "25300/30000 updates, 0.35 s, 11.28 train_loss, 11.42 val_loss\n",
      "25400/30000 updates, 0.38 s, 11.57 train_loss, 11.40 val_loss\n",
      "25500/30000 updates, 0.35 s, 12.69 train_loss, 11.39 val_loss\n",
      "25600/30000 updates, 0.37 s, 11.78 train_loss, 11.61 val_loss\n",
      "25700/30000 updates, 0.36 s, 10.62 train_loss, 11.41 val_loss\n",
      "25800/30000 updates, 0.37 s, 11.21 train_loss, 11.43 val_loss\n",
      "25900/30000 updates, 0.36 s, 10.49 train_loss, 11.41 val_loss\n",
      "26000/30000 updates, 0.37 s, 10.90 train_loss, 11.38 val_loss\n",
      "26100/30000 updates, 0.35 s, 9.99 train_loss, 11.42 val_loss\n",
      "26200/30000 updates, 0.31 s, 13.20 train_loss, 11.43 val_loss\n",
      "26300/30000 updates, 0.31 s, 11.01 train_loss, 11.41 val_loss\n",
      "26400/30000 updates, 0.30 s, 10.29 train_loss, 11.45 val_loss\n",
      "26500/30000 updates, 0.33 s, 10.59 train_loss, 11.38 val_loss\n",
      "26600/30000 updates, 0.39 s, 10.52 train_loss, 11.42 val_loss\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26700/30000 updates, 0.37 s, 10.83 train_loss, 11.49 val_loss\n",
      "26800/30000 updates, 0.34 s, 11.80 train_loss, 11.45 val_loss\n",
      "26900/30000 updates, 0.34 s, 10.30 train_loss, 11.42 val_loss\n",
      "27000/30000 updates, 0.36 s, 12.08 train_loss, 11.40 val_loss\n",
      "27100/30000 updates, 0.31 s, 11.63 train_loss, 11.43 val_loss\n",
      "27200/30000 updates, 0.34 s, 11.83 train_loss, 11.78 val_loss\n",
      "27300/30000 updates, 0.34 s, 10.94 train_loss, 11.45 val_loss\n",
      "27400/30000 updates, 0.31 s, 11.50 train_loss, 11.41 val_loss\n",
      "27500/30000 updates, 0.33 s, 10.76 train_loss, 11.38 val_loss\n",
      "27600/30000 updates, 0.30 s, 12.25 train_loss, 11.43 val_loss\n",
      "27700/30000 updates, 0.30 s, 11.20 train_loss, 11.45 val_loss\n",
      "27800/30000 updates, 0.28 s, 11.67 train_loss, 11.38 val_loss\n",
      "27900/30000 updates, 0.34 s, 12.55 train_loss, 11.46 val_loss\n",
      "28000/30000 updates, 0.37 s, 11.23 train_loss, 11.38 val_loss\n",
      "28100/30000 updates, 0.34 s, 10.60 train_loss, 11.39 val_loss\n",
      "28200/30000 updates, 0.34 s, 11.60 train_loss, 11.38 val_loss\n",
      "28300/30000 updates, 0.35 s, 11.23 train_loss, 11.40 val_loss\n",
      "28400/30000 updates, 0.35 s, 11.60 train_loss, 11.39 val_loss\n",
      "28500/30000 updates, 0.34 s, 11.14 train_loss, 11.39 val_loss\n",
      "28600/30000 updates, 0.35 s, 11.08 train_loss, 11.40 val_loss\n",
      "28700/30000 updates, 0.34 s, 11.94 train_loss, 11.37 val_loss\n",
      "28800/30000 updates, 0.35 s, 11.42 train_loss, 11.47 val_loss\n",
      "28900/30000 updates, 0.33 s, 12.49 train_loss, 11.38 val_loss\n",
      "29000/30000 updates, 0.32 s, 11.15 train_loss, 11.37 val_loss\n",
      "29100/30000 updates, 0.29 s, 12.22 train_loss, 11.45 val_loss\n",
      "29200/30000 updates, 0.29 s, 11.11 train_loss, 11.38 val_loss\n",
      "29300/30000 updates, 0.31 s, 11.12 train_loss, 11.40 val_loss\n",
      "29400/30000 updates, 0.29 s, 11.07 train_loss, 11.35 val_loss\n",
      "29500/30000 updates, 0.29 s, 12.29 train_loss, 11.44 val_loss\n",
      "29600/30000 updates, 0.29 s, 10.81 train_loss, 11.43 val_loss\n",
      "29700/30000 updates, 0.28 s, 12.47 train_loss, 11.39 val_loss\n",
      "29800/30000 updates, 0.29 s, 10.84 train_loss, 11.41 val_loss\n",
      "29900/30000 updates, 0.28 s, 10.28 train_loss, 11.40 val_loss\n"
     ]
    }
   ],
   "source": [
    "batch_pointer = 0\n",
    "\n",
    "start = time.time()\n",
    "best = float(\"inf\")\n",
    "\n",
    "\n",
    "for i in range(max_iter):\n",
    "    x_batch = Xz[batch_pointer: batch_pointer + batch_size, :]\n",
    "    s_batch = S[batch_pointer: batch_pointer + batch_size, :]\n",
    "\n",
    "    _, _loss, _step = sess.run([train_op, loss, global_step], {x_pl: x_batch, s_pl: s_batch, n_pl: n_samples})\n",
    "\n",
    "    batch_pointer += batch_size\n",
    "    \n",
    "    if batch_pointer > N - batch_size:\n",
    "        batch_pointer = 0\n",
    "\n",
    "        p = np.random.permutation(N)\n",
    "        Xz = Xz[p, :]\n",
    "        S = S[p, :]\n",
    "    \n",
    "    if i % 100 == 0:\n",
    "        took = time.time() - start\n",
    "        start = time.time()\n",
    "        \n",
    "        # --- change the following batch if you want a true validation set\n",
    "        x_batch = Xz  \n",
    "        s_batch = S\n",
    "        \n",
    "        val_loss, _step = sess.run([loss, global_step], {x_pl: x_batch, s_pl: s_batch, n_pl: n_samples})\n",
    "\n",
    "        print(\"{0}/{1} updates, {2:.2f} s, {3:.2f} train_loss, {4:.2f} val_loss\".format(i, max_iter, took, _loss, val_loss))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single imputation RMSE\n",
    "The *self-normalized importance sampling* approach for the MIWAE is described in this [paper](https://arxiv.org/pdf/1812.02633.pdf). This needs to be modified slightly in the MNAR case to account for the missing model, as described in the not-MIWAE paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def imputationRMSE(sess, Xorg, Xnan, L):\n",
    "\n",
    "    N = len(Xorg)\n",
    "    \n",
    "    Xz = Xnan.copy()\n",
    "    Xz[np.isnan(Xnan)] = 0\n",
    "    S = np.array(~np.isnan(Xnan), dtype=np.float32)\n",
    "\n",
    "    def softmax(x):\n",
    "        e_x = np.exp(x - np.max(x, axis=1)[:, None])\n",
    "        return e_x / e_x.sum(axis=1)[:, None]\n",
    "\n",
    "    def imp(xz, s, L):\n",
    "        _mu, _log_p_x_given_z, _log_p_z, _log_q_z_given_x = sess.run(\n",
    "            [mu, log_p_x_given_z, log_p_z, log_q_z_given_x],\n",
    "            {x_pl: xz, s_pl: s, n_pl: L})\n",
    "\n",
    "        wl = softmax(_log_p_x_given_z + _log_p_z - _log_q_z_given_x)\n",
    "\n",
    "        xm = np.sum((_mu.T * wl.T).T, axis=1)\n",
    "        xmix = xz + xm * (1 - s)\n",
    "\n",
    "        return _mu, wl, xm, xmix\n",
    "\n",
    "    XM = np.zeros_like(Xorg)\n",
    "\n",
    "    for i in range(N):\n",
    "\n",
    "        xz = Xz[i, :][None, :]\n",
    "        s = S[i, :][None, :]\n",
    "\n",
    "        _mu, wl, xm, xmix = imp(xz, s, L)\n",
    "\n",
    "        XM[i, :] = xm\n",
    "\n",
    "        if i % 100 == 0:\n",
    "            print('{0} / {1}'.format(i, N))\n",
    "\n",
    "    return np.sqrt(np.sum((Xorg - XM) ** 2 * (1 - S)) / np.sum(1 - S)), XM\n",
    "\n",
    "\n",
    "def not_imputationRMSE(sess, Xorg, Xnan, L):\n",
    "\n",
    "    N = len(Xorg)\n",
    "    \n",
    "    Xz = Xnan.copy()\n",
    "    Xz[np.isnan(Xnan)] = 0\n",
    "    S = np.array(~np.isnan(Xnan), dtype=np.float32)\n",
    "\n",
    "    def softmax(x):\n",
    "        e_x = np.exp(x - np.max(x, axis=1)[:, None])\n",
    "        return e_x / e_x.sum(axis=1)[:, None]\n",
    "\n",
    "    def imp(xz, s, L):\n",
    "        _mu, _log_p_x_given_z, _log_p_z, _log_q_z_given_x, _log_p_s_given_x  = sess.run(\n",
    "            [mu, log_p_x_given_z, log_p_z, log_q_z_given_x, log_p_s_given_x],\n",
    "            {x_pl: xz, s_pl: s, n_pl: L})\n",
    "\n",
    "        wl = softmax(_log_p_x_given_z + _log_p_s_given_x + _log_p_z - _log_q_z_given_x)\n",
    "\n",
    "        xm = np.sum((_mu.T * wl.T).T, axis=1)\n",
    "        xmix = xz + xm * (1 - s)\n",
    "\n",
    "        return _mu, wl, xm, xmix\n",
    "\n",
    "    XM = np.zeros_like(Xorg)\n",
    "\n",
    "    for i in range(N):\n",
    "\n",
    "        xz = Xz[i, :][None, :]\n",
    "        s = S[i, :][None, :]\n",
    "\n",
    "        _mu, wl, xm, xmix = imp(xz, s, L)\n",
    "\n",
    "        XM[i, :] = xm\n",
    "\n",
    "        if i % 100 == 0:\n",
    "            print('{0} / {1}'.format(i, N))\n",
    "\n",
    "    return np.sqrt(np.sum((Xorg - XM) ** 2 * (1 - S)) / np.sum(1 - S)), XM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the single imputation RMSE using 10k importance samples\n",
    "If you used the MIWAE loss use the imputationRMSE \n",
    "\n",
    "If you used the notMIWAE loss use the not_imputationRMSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 / 4898\n",
      "100 / 4898\n",
      "200 / 4898\n",
      "300 / 4898\n",
      "400 / 4898\n",
      "500 / 4898\n",
      "600 / 4898\n",
      "700 / 4898\n",
      "800 / 4898\n",
      "900 / 4898\n",
      "1000 / 4898\n",
      "1100 / 4898\n",
      "1200 / 4898\n",
      "1300 / 4898\n",
      "1400 / 4898\n",
      "1500 / 4898\n",
      "1600 / 4898\n",
      "1700 / 4898\n",
      "1800 / 4898\n",
      "1900 / 4898\n",
      "2000 / 4898\n",
      "2100 / 4898\n",
      "2200 / 4898\n",
      "2300 / 4898\n",
      "2400 / 4898\n",
      "2500 / 4898\n",
      "2600 / 4898\n",
      "2700 / 4898\n",
      "2800 / 4898\n",
      "2900 / 4898\n",
      "3000 / 4898\n",
      "3100 / 4898\n",
      "3200 / 4898\n",
      "3300 / 4898\n",
      "3400 / 4898\n",
      "3500 / 4898\n",
      "3600 / 4898\n",
      "3700 / 4898\n",
      "3800 / 4898\n",
      "3900 / 4898\n",
      "4000 / 4898\n",
      "4100 / 4898\n",
      "4200 / 4898\n",
      "4300 / 4898\n",
      "4400 / 4898\n",
      "4500 / 4898\n",
      "4600 / 4898\n",
      "4700 / 4898\n",
      "4800 / 4898\n",
      "imputation RMSE:  1.035147088721874\n"
     ]
    }
   ],
   "source": [
    "# ---- S has been permuted during training, so just reinstantiate it\n",
    "S = np.array(~np.isnan(Xnan), dtype=np.float32)\n",
    "\n",
    "rmse, imputations = not_imputationRMSE(sess, Xtrain, Xnan, 10000)\n",
    "# rmse, imputations = imputationRMSE(sess, Xtrain, Xnan, 10000)\n",
    "\n",
    "print(\"imputation RMSE: \", rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare to missForest and MICE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "from sklearn.ensemble import RandomForestRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/nbip/miniconda3/envs/python36/lib/python3.6/site-packages/sklearn/impute/_iterative.py:638: ConvergenceWarning: [IterativeImputer] Early stopping criterion not reached.\n",
      "  \" reached.\", ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "estimator = RandomForestRegressor(n_estimators=100)\n",
    "imp = IterativeImputer(estimator=estimator)\n",
    "imp.fit(Xnan)\n",
    "Xrec = imp.transform(Xnan)\n",
    "rmse_mf = np.sqrt(np.sum((Xtrain - Xrec) ** 2 * (1 - S)) / np.sum(1 - S))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "missForst imputation RMSE:  1.6231575357178838\n"
     ]
    }
   ],
   "source": [
    "print(\"missForst imputation RMSE: \", rmse_mf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "imp = IterativeImputer(max_iter=100)\n",
    "imp.fit(Xnan)\n",
    "Xrec = imp.transform(Xnan)\n",
    "RMSE_iter = np.sqrt(np.sum((Xtrain - Xrec) ** 2 * (1 - S)) / np.sum(1 - S))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MICE, imputation RMSE 1.4102763723316918\n"
     ]
    }
   ],
   "source": [
    "print(\"MICE, imputation RMSE\", RMSE_iter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the learned missing model\n",
    "There is a separate missing process in each feature dimesion, inspect each of them, plot as function of feature value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl83PV54PHP85tD933blnxjY3zbBAiEcCXhSnZzNCFHSY8tIdk2m8CW0JQmdEsTyA3dkgSSLknaHG2A5iAhpaQc4TYEzGVsbMuWbcm6z9Formf/+P00GgnJlmzNSJp53q+X4Hd8Z/TMWJpH31tUFWOMMWYyzlwHYIwxZv6yJGGMMWZKliSMMcZMyZKEMcaYKVmSMMYYMyVLEsYYY6ZkScIYY8yULEkYY4yZkiUJY4wxU/LPdQAn6+KLL9b7779/rsMwxpiFRqZTaMHXJDo7O+c6BGOMyVoLPkkYY4xJH0sSxhhjpmRJwhhjzJQsSRhjjJmSJQljjDFTytgQWBGpB24CNqnq6ZPcd4AvAIPAUuC7qvpkpuIzxhjzRpmcJ3EO8DNg8xT33w+Uqur1IlIJPCkip6pqPB3BRMLhdDytZ1rDj6ddbKz85A+QKa7P9HmOxZHjVzpTn3Y6Mc04bmNMxmUsSajqT0XkvGMUuQz4D69st4iEgdOAnemI5/tX381IfkM6nnp+0kTyUPC2rNXR/6ibrzT1OOGW0wSiiqCIxkEVR+OIxpFEHNEYTiKOozEkEcUXj+AkojiJME48jD8exomH8MWHcOJDOLF+fLEelDBxH0T8QtQHIwEhHBSGg8JQnjBQ4NBf4NBb5Ke9NEhHaZBQXj6O5uGQh48CAlJCUEoodCqozGugpqCBxcX1NFYUs6SikKVVhSypKLBkZMxJmE8zrmuBgZTzfu/aG4jIVcBVAE1NTemPLBuk1ASSu5qnfHZmeqdzX3yE/OFOCkNHKRxup3zwMKVHm8kPd09ZweorhH31wr562L1YeGmZMOgXuoFDcWAQtD9AfO8S4sPLiIeWUco6tjVVsW1pJe/c1MCSisIMvkpjFr75lCTagZKU81Lv2huo6h3AHQDbt28/oc83f8nLOOFXT+Shs2Smf91O9TJn8jwy7kjHPVbGrql7pohXUBAcr7zgjndwAB+KA/hR/ECAhARRgiQkj4TkgwQmjSTuy2OoeDFDxYvHXfdHByjrf53a9p1Udb9MMDqUvFcWgi37lC373PcjHIDnVglPrBWePkVQRxAnir9oP/6i/QBEYsU80rmNB18/nS/9pprzTqnhw2cs5YK1tTiO1TCMOR5RzdzfkF5z01dUdbt3XgQUqmqHiFwBnKuqnxjtkwCO2yexfft23bFjR7pDN9ORiMPIAIR7IdQNoS7i/R2EO48S7uom3N3LUNcAg4PKULyK3lgDvfFFDMRrUXxveDoRqF8U4JTFIeriLUT27Cb88itE9u3zmsbGxBvrOfz+c3hqXYAd7b+nuX8vOiGxRnpOZ6T9Mkjkc8bySr5xxWYaygrS+pYYM49N66+kjCUJEXkrcCVwMfBN4KvAnwAbVPVqb3TTF4EQ0ATcOZ3RTZYkFqBICDp3w+EdcGgHseZn6OxwOBo9hbboKbSMbGFEi8c9pLgyj00XNLLhrUtgZJihxx5j4D/+g8GHHiIxNFbbyFu9mkW33Exk1RIeanmIu3ffzXPtzyXva6yM4SPvIT60hrKCALe8dyMXr6/P2Es3Zh6ZX0kiXSxJZAFV6HgNXrsPXvkZ8SMv0RZZy76RM9k1fD4RLUoWLast4NwrTqFpXRUAieFhen78E7ruvJN4dzcAkp9Pw003UXb5ZQDs7d3LzU/fzJOtY39zjHRcSKTzbQB86qLVfOqiUzL1ao2ZLyxJmAXq0LPw9LfhpXuIxHzsDr+VF0Lvoje2KFlk5dYazv/IWvIK3T6PxNAQXf/vLjpvvx0S7kiuyj/9E2qvuQbx+VBVfrrnp3x1x1cZGu3n6PrvDLSfCcA3P7yVSzbk0Gg3YyxJmAWv9yD85q/h1Z8TVx8vhi7j6dCHicaDAFTUF3LpJzZSXjs2Ymnwd49x+NprSfT1uWU+9CHqP/c3yfsH+g9w5a+vpDvc7XbOt3+Yga71FAV9/OzPz2ZVbQnG5Ijc2E/CZLHyJvjAD+DKn+GrXs7mop/z4cqrWVbyMgA9bSF+essODr/Wk3xI8Tlns/zf/pXgqpVumR/+kN577k3eX1q6lG9e9E0K/YUoiq/uJ/gK9zIUifOxHzzL4Egss6/RmHnOkoSZ/1acB3/2X7DmUop8PVxa+DdsLf0FACNDMX5+6/MceKkrWTzY1ETTnXfiq6wEoO3GGxl+8aXk/XVV67j1glvxO37iGqVmxb0gUfZ2DHHDvS9m8pUZM+9ZkjALQ34pfOBf4K3XI6KcVfhPXFj9HRwfJBLKb+58ic5DY3MxAw0NLP7618HnQyMRDn3yk8S6xhLJmQ1ncv3p1wMwFO9kzSnPA/CzF46wt2Mws6/NmHnMkoRZOBwHzv8ruPwbAKz138c7au8EgehInF/+350M9owkixed8SbqrvtLAGKtrRz9+y+Me7r3nfI+VpWvAqA3eD/iG0IVvvPovgy9IGPmP0sSZuHZ/sdwwQ0ArNBfcXbdLwEY6h3hvttfIBIe61eouPJKSt52EQD9v/oV4V27kvd8jo9rtl0DQCg2xMpTHgfg7ucO0zEwlmyMyWWWJMzC9Jb/DWdcDcAm/S7rl+wGoLNlkMfv2ZssJiLUXnst+NwZ3R233jbuac5ZfA5nNJzhPlYeQgKdRGIJvvd4cwZehDHznyUJszCJwDu+CKvfgQi8JfpXLF7s1iBefuQwrXv7kkWDy5ZR9u7/DsDgf/0Xw88/n/I0wrXbrkUQ4hqnbulvAfjBkwcYspFOxliSMAuY48C7/gEKKnAkwXnO5/D53aHfD/3LLuKxseXRaz7xCSTgTrxrv/XWcU9zatWpXLriUgCGAy+Ab4i+4Sg/eaYlQy/EmPnLkoRZ2Erq4LKvAVAefZXtjc8A0H1kiN8/cDBZLLBoEeUf+AAAoSeeZOjJp8Y9zXtXvxeABHHq692mq7seb2ahTzY15mRZkjAL3/r3wGnvAWBL6BYqK92Fg3fc10xveyhZrPpjVyH5+QB03XnnuKfYVreN2kJ3+5LyWney3sHukA2HNTnPkoTJDpd9FYpq8EmM80pvByAeS/DUz8eGs/praih/z7sBGHriCWKdncl7jjhcvOxiAI6EX0b8bp/Gw7vHyhiTiyxJmOxQWOmOeAIawr9l9Sq3BrH32Xb6OsZqE6WXv9M9SCTov/83455itF9CUapqXwHgkd0d6Y7cmHnNkoTJHtv+CErdne62xdyhrqrw+/8Y65so2LKZwCJ3Ndn+X/5y3MPXVa5jaelSAPIr3K3Vn9rfRTh6zH2vjMlqliRM9gjkw7lubaIq/BTLGt0lwV99opWhPndynIhQepm7z8Tw888TOXQo+XAR4ZLllwDQl9iPBDoJRxM8vb87k6/CmHnFkoTJLps/AuVubWCrfhOAREzZ+dux4ayll1+ePO6/71fjHj6aJACC5S8A1uRkcpslCZNd/EF462cAaIg9RkPdMAAvPXyYkWF3clz+mlPIW70aeGOT04qyFaytXAtAcaW7cuwjeyxJmNxlScJkn40fgMoVAGzN/xEAkXCclx85nCwy2uQ0smcP4dd2j3v4RU3uWk8RpxXxDbD76CCtfcOZiNyYeceShMk+Pj9s+2MAloZ/QWWNOwv71cdbk5PjSr39r+GNtYmtdVvHnqqwGbAmJ5O7LEmY7LT5Q+AEEIG1Fb8HoPdoiPYD7p4TwSVLKNi0CYCBBx8c99D11evxix+AwlK3L+MRmy9hcpQlCZOdiqrhVLeD+pSBO5K7+e5+qm2syFveAkBk375xE+sK/AWcWnUqAMVl7uin373eSSJhS3SY3GNJwmSvrR8FoEhbWdLgTqjbs+Moibi78F/h6acni4Z27Bj30M21mwEY0gMgEfqGo7T0hDAm11iSMNlr+VuhYhkAaxx3qOvwQJSWXT0AFGzamFwZNvT0M+MeuqV2C+Au+OcrcGsTr7b2ZyJqY+YVSxImezkObL0SgBWRe/C53QzJJicnPz/ZLxF6ZvIkAeArOADAq60DGJNrLEmY7Lb5I+D4CTrDLK9xawT7nu9IbnFa+Ca3yWlkzx5iPT3Jh1UXVLOkeAkAxWVu5/WuNqtJmNxjScJkt5I6WHkhAGsS9wIQiyTY/4LbUX2sfolkk1NeM5CwmoTJSZYkTPbzRjk18jD5Be6lPTuOAlCweTOM9ktMbHKqc5NEnBBOsIOD3SEGbUtTk2MsSZjst+ZSEAefxFleewSAw7t6iEXjOAUFFKxfD0DomQk1iZqUfglvUt1rbVabMLnFkoTJfkXV0HgmAE1Rdw+JWDRB6x53Y6HRJqeRXbuI9/UlH7aifAUlwRIgtfPa+iVMbrEkYXLDaJNT7D8Rb2LdgVe6gJR+CVVCzz6XfIgjDptr3PkS/qJmwDqvTe6xJGFyw1p3raY8J0R9ldtkdPAlN0kUbNkCPh/wxn6JTTXuEFkJdIMTYpd1Xpsck9EkISIXicjtInKjiHx+kvvLReRuEbleRH4sIu/KZHwmi1Usg7oNADT5nwKgpy1Ef9cwvuIi8k87DXhjklhdsTp57OR1sKttwJbnMDklY0lCRAqBbwGfVtUbgY0icuGEYtcBv1PVm4FbgK9mKj6TA7zaRFP018lLLa+4u84VbnGblUZeew2NRpP3V5WvSh778o4yOBLjcK8tG25yRyZrEmcBB1R1xDt/DLhsQpmjQI13XAM8m6HYTC7w+iVq/PspyHeHsh7wmpzy1roL+mk0ysj+/cmHLC5eTNAJAuAE2wHrvDa5JZNJohZIbdDt966l+hpwhoh8Dfgc8P8meyIRuUpEdojIjo4OW+ffTFPdeihvQkRpKnY3Gjr0Wg/xWIL8tWuSxUZeey157HN8LC9bDoCTN5okrF/C5I5MJol2oCTlvNS7luou4Duqeg3wbuAnIlI58YlU9Q5V3a6q22tqaibeNmZyIrDqbQA0xd09JKLhOG37+giuXAl+d3Gn8Ku7xj1sRbm7y12gwP1xtRFOJpdkMkk8ASwVkTzv/GzgPhGpFJFS71oj0Ood9wCJDMdost3ycwFoDIx1UB98uQsnGCRvhZsMRnaNTxIry1YCoL5ecMLssgl1Jodk7ANYVUPAx4HbROQmYKeqPghcD3zCK/Zp4M9F5LPAN4HPqqptCWZmzzJ3o6ECZ4DaSvfDvuVVd2G/PK/JKZzS3ASwsnxl8tgJttPcNUQoYstzmNzgz+Q3U9UHgAcmXLsu5fh3wO8yGZPJMUVV7lDYoy+yOPgy7ZxJZ8sAkXCM/DVr6ecXxLu6iHV04PeaMkebm8Dtl4iFmzjYHWJtfelU38WYrGFNOSb3eE1ODZGHAVCFo/v7kzUJgPCusdpEY0kjfsf9e8rndV63dNswWJMbLEmY3OMliXr/y8lLrXv7yF+7Nnk+8tpYv0TACbCsdBkATp67euzBbtvK1OQGSxIm9yx9M4iPAmeAipIhANr29uKvqsJXUw2Mr0nAWL/EWE3CkoTJDZYkTO7JL4VF7jLg9UE3GbTt6ycRT5C/xq1NhHe9Ou4hoyOcJNADMmJJwuQMSxImN432S0QfByA6Eqfr8FByUl1kfzOJkZFk8fGd1x209FiSMLnBkoTJTaNJIvhK8lLr3r7k8hzE44zseT15b7QmAW6/REv3MKq20J/JfpYkTG5qPAN8Qcp8rRQE3RpD697eCctzjHVeLy1dik/c5cSdYDvD0Tidg5HMxmzMHLAkYXJTsBCWvAkRqM/fC0Db3j6Cy5YhQXdBv9TO64AvQFNpEzDWeW0jnEwusCRhclfTGQA0qLtEx2DPCIP9MfJWu3tITLU8x+gw2EPWL2FygCUJk7uWvAmAhsDYSKbWvb1jy3Ps3j2u32G089od4RSxEU4mJ1iSMLlryXYAagJ78fkSALS+3kfeSnejoURfH/GenmTx0Ql1IooT6LHmJpMTLEmY3FVUDRXL8UmM2qI2wB3hFFy2LFkk0tycPF5SsiR5LMFuW5rD5ARLEia3NbpNTvXyAgDdR4ZwFjclb0f2NyePlxSPJQkn0G01CZMTLEmY3LbkdABqcNdx0oQyIOXgc4e7ptYkqguqyfO526E4gW5a+4aJxhOZjdeYDLMkYXLbaJII7E1e6jgyTGDJYmB8khARFhUvco+DPSQUjvRak5PJbpYkTG6rOw38BZT52gj6owB0HBxI9kukJgmAxcVu8nAC3YAtGW6ynyUJk9t8AVi0BRGoyT8EQPvBAfKWLQcgcuAAmhhrUhrtl3CThFq/hMl6liSMafSanHQnAD1HhvA1LQNAIxGiR1qTRUdHOIkvgvhCttCfyXqWJIzx+iVq/e6CfomEMlTWmLyd2uQ02twEIDbCyeQASxLGTNJ53asVyeOp5ko4gR4OWZIwWc6ShDEl9VDWRJmvjYDPXdm1q1eQggJg6pqEE7SahMl+liSMAViyHRGlJtAMeCOcli4FxieJkmAJpcFSwG1u6glFGRyJZTpaYzLGkoQxkNzOtMZxF/vrPjKEb5m7oN/Uw2DddZ1sNViTzSxJGAPQsAmAWq9fIhFXhuvd1WCjhw+TiIxtMDTaL+EE3bkSrX3hTEZqTEZZkjAGkkkitfN6oMjrpFYlevBg8vroXAkJ9AIJ2ixJmCxmScIYgIJyqFhGua812XndmyhL3h7Zvz95PNrcJBJH/P2WJExWsyRhzKiGzW7ndfAAAN39Y78e40Y4lYwf4WRJwmQzSxLGjBptchKv8/poGKmoAibMlUhZMlwC3bT2W5Iw2cuShDGjFm0GoNobBpuIKdEVbuKINB8YK1a8CEEAdw2ntj5b5M9kL0sSxoyqdxNClb85eSnkjXBKrUkEfUFqCmsAcII91txkspolCWNGFVVBWSMV/kOIuCu/Dha5/Q/xri7iAwPJomMjnLrpD8cYsgl1JktlNEmIyEUicruI3Cgin5/kvojIJ72vr4nIP2UyPmNo2IRfolQEjgLQz9gIp+ihQ8nj5FwJb1+JNuuXMFkqY0lCRAqBbwGfVtUbgY0icuGEYh8BelX1NlW9BvhGpuIzBoAGt1+iytkDQO+QP3kr0tKSPB6bdT0AErUmJ5O1MlmTOAs4oKoj3vljwGUTynwYqPRqEl8ABjMYnzHJzusqv9tRHRpKEAkUARBtGatJjFvoL2D9EiZ7ZTJJ1AIDKef93rVUS4FSVb0NuAu4X0R8E59IRK4SkR0isqOjoyNd8Zpc5A2DHR3hBBCqXgVA5NBYTWJ0r2twZ15bc5PJVplMEu1AScp5qXctVT/wFICq7vbKNE4og6reoarbVXV7TU1NmsI1Oam4FkoWjR/htOhUYHxNor6oPnnsBHpptWGwJktlMkk8ASwVkTzv/GzgPhGpFJFS79qDwAoA75oPaMtgjMZAwyaKnG7yfEMAhMqaAIim9EnUFdYl50qIv8+am0zWyliSUNUQ8HHgNhG5Cdipqg8C1wOf8IrdAmwWkc8CXwc+qqr222cyq2ETIlDt2wfAQKAagMiRI2g8DrhzJaoL3OuONTeZLOY/fpHZo6oPAA9MuHZdynEf8LFMxmTMG9RvANxJdYcjG+iLFJIQBycaJXb0KIFFbn9EQ1EDHcMdbp9EtyUJk51sMp0xE9WvB8ZGOCVUGC5w+74ik/RLOIFeOgcjjMTiGQ7UmPSzJGHMROVLIa+U6sDY8uCjM6+jKSOcRpOE+PuABO39IxiTbWacJESkaLJhqcZkDRGo30Cl/xCCtzxHiTvILnVCXUNRg1vciSG+IeuXMFnpuElCRBwR+ZCI3Cci7cAuoFVEXhaRL4vI6vSHaUyG1a3HLxHK/a0AhKqWA+OHwY4mCQAJ9Nk2piYrTacm8V/ASuCvgHpVbVTVWuAtwJPAzSLykTTGaEzmJTuvvRFOhW5ndeqEuvri8XMlbMlwk42mM7rpIlWNTryoqt3A3cDdIhKY9ciMmUspndev8xbCUkTUX4BvqpqEv9dqEiYrHbcmMZogROQbIiLHKmNM1qg5FcSXHOEEMFS0iHh3N/FBd5JdRV4FeT53bqgT6OWo9UmYLDSTjutB4OciUgQgIm8XkcfSE5YxcyyQDzVrqPIfTF4aLHKbnEZHOInI2AingNUkTHaadpJQ1RuAHwEPicjvgGtxZ0sbk53q1lPi6yAg7of/kJckUkc4jc2VsKU5THaadpLw9n74M2AIqAE+qaqPpiswY+Zc/QZElEqvySlZk5ikX0ICvbQPjBBPaObjNCaNZtLc9NfA51T1POB9wE9E5IK0RGXMfJDsvHabnIaKF6OMn1A3miQc/wBxjdI5aBPqTHaZSXPTBaM1B1V9EbgEuCldgRkz5+rcYbCjNYmYv4CRvPJxS3OMH+HUx5FeGwZrsstxh8CKSNMxbv9pyv1eVe2fnbCMmQeKa6CkgaqR8SOcSifpkwAb4WSy03TmSXwPUGCy4a+j1xV3J7nvz1pkxswHdeup6nsqeTpYtIjo4YfReBzx+cYlCQn0Wue1yTrHTRKqen4mAjFmXqrfQMHrD1Dg9DCcqGCwaDGasmT4uJqEv5c2W+TPZJmZjG56eHQHORG5WkQ+JSLB9IVmzDwwYdnwIa8PInLI7Zco8BdQkVcBuOs3WXOTyTYzGd1Urqr9IrINdyhsBXBnesIyZp6o3wiMjXAKFdWTEGfS/a5tr2uTjWaSJKIi4geuBG5R1c8Dp6UnLGPmicoVECik0ksSCSfAcEHNuIX+UudKHLXmJpNlZrJ96W3AC0A+YzOti2c9ImPmE8cHteuoCo0f4TRuQl2xN1fCWwlWVZlimTNjFpyZzJP4PnAGsF5Vh0VkFfBE2iIzZr6o3+DVJLwNiIoWEZ1086EIw/FB+sOxuYjSmLSYzqZDyT+JVHVQVYe949dV9Y8nljEm69SvJyARynxtgJskRjuuAeqK6pLHNlfCZJtpbTokIn8xcVKdiARF5AIR+R7w0fSEZ8w84HVeV6YszxHv6iIx5C4ZPn6HOlsN1mSX6SSJi4E48CMRaRWRV0RkP7AH+CDwdVW9K40xGjO3atcBkhzhNJxfRdwJEjl0GIBF3sJ/AE6gh6OWJEwWmc5kujBwO3C7twNdNTCsqr3pDs6YeSGvGCpXUBX2Oq/FYaiogeihFvLXnEJVQRUBJ0A0EXU7r625yWSRmUymuwR4FHgIuENEzkxXUMbMO/UbqPI3J08Hixcn95VwxBk3DNaShMkmM5kncTvuRkNnAncAXxGRD6YlKmPmm/r1lPna8OHOgxgsWjzpvhKOrd9kssxMksRRVX1MVXtU9T+Bd+DuMWFM9qvfiCOJZL+EO8IpZRhs8WhNoseShMkqM0kSzSJyU8p6TVFgIA0xGTP/1Lt7S1QFvF3qiheN21ditPPa8Q/S1m+/FiZ7zCRJKPAeoMXb4/p13P2uV6clMmPmk5IGKKhMLvQXCxQz2D6AJtwJdqM1CYCeSAcjsfichGnMbJv2shyq+kEAEckH1gObvK/viMgKVW1MT4jGzAMibud1b3Py0mCwhlhHB4G6ujcMg23vH6GxsnAOAjVmds1k7SYgOSR2h/dlTO6o30DV3u8lT93O6xYCdXXjJtSNDoO1JGGywUyam4zJbfUbKXAGKJQuYHQYrNsvUV9Uj3ibN1rntckmGU0SInKRiNwuIjeKyOePUe7DIqIiYqvMmvmjYRMwtgHRYNEiot4aTgFfgKr8asDWbzLZJWNJQkQKgW8Bn1bVG4GNInLhJOVOBdZlKi5jpq16NfgLqA42AxAqrCecMsJpcYnbL2HrN5lsksmaxFnAAVUd3ZXlMeCy1AJeIrkO+NsMxmXM9Dg+b+a1W5NQx0fPkaHk7eQwWJt1bbJIJpNELePnVfR711L9PfB3qho51hOJyFUiskNEdnR0dMxymMYcQ8MmqlOW5+jp1bFbxSlLc/SFMh2ZMWmRySTRDpSknJd61wAQkUbcfbPfLyKjO99dIyLbJz6Rqt6hqttVdXtNTU06YzZmvIZNlPuP4Ki7sVB/opTEsLuvdXL9JknQOtg+5VMYs5BkMkk8ASwVkTzv/GzgPhGpFJFSVW1R1T9S1ZtV9WavzNdU1Ybamvlj0WZ8EqNU3WXC3ZnX7vIci4rH5kp0httIJHTSpzBmIclYklDVEPBx4DYRuQnYqaoP4u6X/YnRciJSIyI3eKfXicjiTMVozHHVrAVfkCrf6AinxUQOuMepcyUSvh66Q8dsNTVmQZjxZLqToaoPAA9MuHbdhPMO4Cbvy5j5xReAutOo6d7P3uFzieSVMfB6M6VvG1+TGF0Ntro47xhPZsz8Z5PpjJmphk3U5O9Pnrbvd/ffKgoUUeR3u91sQp3JFpYkjJmphs3UBMaSRFdHLHlcX1QPuDWJI33DGQ/NmNlmScKYmWrYRIHTT0HMXZ6jN5yfvNVY4nahib+XQz2WJMzCZ0nCmJmqXQeOn8pEMwB9wTrig4PAWL+EE+ihpXtoqmcwZsGwJGHMTAXyoeZUqn17AQgX1DCwy21+Gk0S4ovQ0ts1ZyEaM1ssSRhzIho2UZu/N3l69CV3DafUYbCHB49kPCxjZpslCWNOxKLNNJS8njxtb3ZHOKUOgx2ItxOKxN7wUGMWEksSxpyIRVsoCXaTF3GTQ3enu11pY8nYBo0S6OZIr3Vem4XNkoQxJ6J+A/iClI24s617RgoAKMsrS86VcIKdtNgIJ7PAWZIw5kT486BhE5XaDMCQU8bIcBSAxcVubcIJdtkwWLPgWZIw5kQtOZ1q3z73WBzaX3E7qleWLwPcJHHYkoRZ4CxJGHOilmynLn+s87rtJXdl2KVlTQCIv4+D3b1zEpoxs8WShDEnasnplBcfJRBx99LqaO4DoKnESxKiNPe1zFl4xswGSxLGnKiyRvIaqikZdBNBV1cCgKbSpmSRtuFDkz7UmIXCkoQxJ0oEaTqdcm+EU38kn2gkztKSpckiA/E2wtH4XEVozEmzJGHMyViynUodXRFW6GwZpCyvjHynGPA6r22uhFnALEkYczKWnE6tsydNFkw/AAATRklEQVR52ra3FxGhrtBdDdZGOJmFzpKEMSdj0RbKizvIC/cAcOSVowAsK3ObnJxgp82VMAuaJQljTkawiGBjI2X97nyJo82DqCqnVC4D3GGwB2wYrFnALEkYc5LyNm6lrM9NEsNh6O8Ms7xsGeAOg93bc3AOozPm5FiSMOYk+deeQ0U4ZVLdvr5xC/21DNpcCbNwWZIw5mQ1nUlNYD9OPAJA296+cXMlOsI2V8IsXJYkjDlZlSsoqA1QOtAMQOveXiryKghKIQChxFFGYjZXwixMliSMOVki5K9ZQ1mfO1+i+8gQ0XCcqjxvK9NAF6294bmM0JgTZknCmFmQt/3c5AgnVTi6v58lxW6Tkw2DNQuZJQljZkHeWe+idGBf8rx1Xx8rK9y5EhLoY19n31yFZsxJsSRhzCyQ2lWUlIQoDLUB7gin02pWuvdE2Xl071yGZ8wJsyRhzCzJX1qfnC9xdF8fS0vGRjjt7t4/1cOMmdcsSRgzS/LWbUgmiUg4TlmoLnnv0KAlCbMwWZIwZpbknXFRsvMaYGB/nCJfFQAhaaHP2wPbmIXEkoQxsyR/y5spDB1NLvbX8mo3TcWrAfDlHeH19oG5DM+YE5LRJCEiF4nI7SJyo4h8fpL7nxGRr3v//1cRWZvJ+Iw5Gf7qavzFQSp7XgXg8O4eNlSuA0CCXbzU2jGX4RlzQjKWJESkEPgW8GlVvRHYKCIXTihWDFyjqrcAdwNfzlR8xsyG/BWNVHa7SSIWSbAhsR5wRzg9e+SVuQzNmBOSyZrEWcABVR3xzh8DLkstoKp/o6qaEttgBuMz5qTlbTqdyp7XQN39rgvaqpP39vS+NldhGXPCMpkkaoHURtl+79obiEgQ+ChwwxT3rxKRHSKyo6PDqvBm/shfv5lAbIjSAXd58J7XI/goAKA1bHMlzMKTySTRDpSknJd618bxEsQ3gb9W1Ul/q1T1DlXdrqrba2pq0hKsMScib63bjVbZ7TYtdR0aohG3XyLiHLIRTmbByWSSeAJYKiJ53vnZwH0iUikipQAiUgB8G/iaqj4rIu/NYHzGnLS8Vatwiouo7NmVvLZ+5HQAnLw2drX1zFVoxpwQf6a+kaqGROTjwG0i0gHsVNUHReRLQDdwM/AvwHpguYgAFOF2YBuzIIjPR+G27cQeeRR/bJiYv4Cm/uVQA+LEeLJlF2csn7SV1Zh5KWNJAkBVHwAemHDtupTj92QyHmPSoWD7NgYffpjyntforNmM70g+VAuIsrP9FeDcuQ7RmGmzyXTGzLLC7dsBqPLmS0RDUDXkbme6f2D3nMVlzImwJGHMLCs47TQkPz85XwJgTc+bAeiK2hpOZmGxJGHMLJNgkILNmykId1EaOgDAqq71oBD3H6Y3FJnjCI2ZPksSxqTBaJNTbeuz7vlwCdVDSxDfME+32HwJs3BYkjAmDZJJov3Z5OzrVZ3bAHi8ZeecxWXMTFmSMCYNCjZthECA/JFeqiNuZ/Wqrq2gwnNHn5/j6IyZPksSxqSBU1BAwXp3cb+6jh0AFEfKaehfwYFhSxJm4bAkYUyaFG53m5cqDuzEIQbAqq5tJAKHeLH10FyGZsy0WZIwJk1G+yWC0SEa1O2HWNG1CSfh8G8v/3YuQzNm2ixJGJMmBdu2QSAAQH3Pc+61WDFL+k7lybYn5jI0Y6bNkoQxaeIrLqb4nHMAKNm7Gz/DAGxofQttkZ0kEom5DM+YabEkYUwalV56CQDO0BCr9CEAGvtOpSJczO8OvjiHkRkzPZYkjEmj4vMvQPLc1fEbO38PuLWHTUfO555d1i9h5j9LEsakka+4iOLzzgNA93WyMvg44I5y2nXopTmMzJjpsSRhTJqVXnopAInhEdYE3MTgUx+NhxcRiobmMjRjjsuShDFpVvzWc3EKCwEI9AUoyHe3Nl3Xfha/etlGOZn5zZKEMWnm5OdTfOGFAAzueIWtNS8AEIzn8+wDrx7rocbMOUsSxmTA6CgnDYdprFpJOL8FgPr9qzh8pH0uQzPmmCxJGJMBxWefja+yEoDenz/K4lp3F19/Isg9dzw4l6EZc0yWJIzJAAkGqfof/wOAyL59bK85h8NVTwLgtNWx92nbsc7MT5YkjMmQig9ega+6GgD56b8zXLOXEZ87uuk/f7iT6Eh8LsMzZlKWJIzJEKeggOqr/gyAyIEDXM4ZPNf4SwBi4RIe/65NrjPzjyUJYzKo/P3vx19bC0DTL3/KwbwR2ordpqaXdvp46Wc2JNbML5YkjMkgJz+fqo9dBUDsUAtXtTXywOq7CAX6AXjk14McetpmYpv5w5KEMRlW/gd/QGDRIgC23XcvNQeruH/Nd4hJFMXH/Xc107u/ZY6jNMZlScKYDHOCQRZ95Svg90Msxt880E5I23l45Y8AGEkUcs9XnuPocy/McaTGWJIwZk4Ubt1C3XXXAVDS38On/72U16t28MySXwEwHC/h3+88wr777p/LMI2xJGHMXKn4w48kZ2JvOHyYq+8t54WG+3lk+b8CCWKax69/4eeZW79DfKh/boM1OUtUda5jOCnbt2/XHTt2zHUYxpyQxNAQ+z/wASKv7wXgtUV+vvI+pSxyGpfs/iiSyAegItjKeZcVs+jt7wKRuQzZZI9p/SBZkjBmjsW6ujj0F59k+Dl3H+yOEofbL4e2mkbetecPyRtuSJZdUfYKmy9ooP6Cy5FA3lyFbLKDJQljFopEJELb52+k7957k9deWCb85Fw/9c4lbDhwPqrB5L26vH2sWxdh2blbKVzzJnCs5djMmCUJYxYSVaXn+9+n/evfQMPh5PXXFsMza2rwlb6Tur5NoP7kPSFOff4+GheHqV9VRd2W9QQb14PPP9m3MCbV/EsSInIR8B6gHVBV/dsJ9/OBrwCHgdXAzaq6+1jPaUnCZJtYRwcd3/o23T/+CU48Nu5eS3Up+5adS6z4bITSSR6doMzXRkVhHxUVMUqr8imuLqa4vprCulryq2pxSmrAH5zksSbHzK8kISKFwE7gNFUdEZG7gdtV9cGUMtcDCVX9kohs8O6/5VjPa0nCZKvIocMc/N4/0/2rX1PSdXTcPUXoK11OR80muqrWEyqsn/bzOhrCYQiHMI6MIEQRiSIS974SIHGQhPsp4iiQcPvLBRRFxPvcSH7MqHdPQPSNnz7T+ThK7ZB/Q/mF3eKRTu+94XoKiopO5KHzLklcCHxWVS/0zq8BlqjqNSllHvXKPOqd93tlphz/Z0nCZDtVpeP3L/LCv/6C0EtPU3u4mfLhyLgykUAR/aXL6StdRqiwgaHCOoYLalHHN0dRm0y54sY1VNUvPpGHTitJZLLhshYYSDnv965Np8y4JCEiVwFXATQ1Nc16oMbMJyJC7daNvG3rRgDi8QR7Xm3m0Mu7OLLneYYP7cbp7yZvqIPC0AEKOyNUReLkRQCnGPVXEPNXEAuWEPUXEQ0UEfUXEvfnE/PlE/fnkXCCxJ0ACV+QhPhIOH5U/JZkTEaTRDtQknJe6l2baRlU9Q7gDnBrErMbpjHzm8/nsHb9CtauXwFcOmU5VSWWUKLxBOFIjJHICJHwMCMjISKRMPFYzP2KR9FEgkQiTiIeA1UURRMJNB4HhUQ8gSbiSCJOPB5HEoom4pBQ0ASqQCKBoKiq1zrk/WomRn9FlXENF+N+cxNT3Rh/Omu/7dkz16Sk/M1pff5MJokngKUikqeqI8DZwO0iUgnEvCal+4CzgEe9PokXjtXUZIyZmogQ8AkBn0Nh0A/kA2VzHZZZYDKWJFQ1JCIfB24TkQ5gp6o+KCJfArqBm4Fbga+IyA3AKuBPMxWfMcaYN7J5EsYYk5um1eZm0zSNMcZMyZKEMcaYKVmSMMYYMyVLEsYYY6ZkScIYY8yUFvzoJm847YETfHg10DmL4cyW+RoXzN/YLK6ZsbhmJhvj6lTVi49XaMEniZMhIjtUdftcxzHRfI0L5m9sFtfMWFwzk8txWXOTMcaYKVmSMMYYM6VcTxJ3zHUAU5ivccH8jc3imhmLa2ZyNq6c7pMwxhhzbLlekzDGGHMMObVbuoh8HQgBg8Am4FOq2jZJuY8AW4A4sFdVv53muBzgz4C/Ay5Q1ZemKPckEPZO46O7/M2DuI65d3maYqvEXTl4H+5+6J9V1aOTlGsGmr3Tw6r64TTEMut7t2corj8CrmbsZ+q7qvqDDMRVD9wEbFLV0ye57wBfwP09XerF9eQ8iOs84BtAr3fpPlX9cgbiWunF9RywBOhS1f8zoUz6fsZUNWe+gJtSjj8D/MMkZZYAzzPWFPcMsDrNcW0BNuN+mK0/RrkbM/x+HTcuoBB4Hcjzzu8GLsxAbN8C3u8dvxP4wVy8Z9N5/cD1wHXe8Qbg0Qy8P9OJ64+AZZn8mfK+7/u8f7MdU9y/And/e4BKYDfgmwdxnQecNwfv1+nAf0s5fwXYlqmfsZxqblLVG1JOHdy/VCZ6B/Cseu827mZJl6Q5rt+r6vPTKLpBRD4jIjeKyGXpjGkGcZ0FHFB3IymAx4C0x+Z9jyem8T3PFZHrROTvRCQdW3hN5/UnY1XVF4FNIlKahlhmGhfAn4vI/xaRz3m1s7RT1Z8yfpviiVLfr27cms5p8yAugD/03q//IyKN6Y7Ji+sZVf1ZyiUHGJpQLG0/Y1nX3CQivwHqJrn1OVX9uVemHHg78N5Jyk1nL+60xDUNt6jq0yLiAx4RkQFVfWSO40rL+3W82CZ8336gQkT8qhqbUPZ67z0rBJ4TkctV9fXZiM8za3u3z7LpxPUwbpNJh4hcCvwbkNYmzGlK28/USXoF+DtVbRaR04AHRGSdqiaO98DZIiLvBn6jqrsm3Erbz1jWJQlVfcex7otIGXA78CfeXykTtePuijeqFLfanta4pvkcT3v/j4vIo8D5wEkliVmIa1r7kp+IY8UmIqPft9f7nj2TJIjU9ywkIs/jbps7m0li1vZun2XH/Z6quj/l9LfAz0XEp6rxNMd2PHPxfh2XqranHL/s/bHZyIkvCzQjInI+7u/8pya5nbb3LKeam0SkGvhH4C9Vdb+IvNe77ohIk1fsN8A2ERndteks4NeZj9YlIsu9/68VkdTtXFczux92MzIaFyl7l3vnZ+PuVZ5uo/uhj/ueqf+WInKhiKSuTbMK2DvLcUz6+kWkMqW6n4w1g3u3HzcuEfmiiIz+obga2D9XCUJEikSkxjtNfb8qcTfnfnmu4xKR60eb5Lz/B4E3DJZIUxyX4TaF/y+gXkTOytTPWE7NkxCR53BrT6M1iAFVfaeIbMbt+NzglfsIsB13dNNuTf/opgrgfwLXAj8AfqiqT3o/nM8DK3E78P4Rd4RDKRAArklnVXc6calqWETehtvp1wFENXOjm27B/StuJW6z0tHUf0vvl+VG4FlgEe7opi+mIZY3vP7RvdtV9WYRKcAdedKKm6i+oJkZ3XS8uP4XsB7Yj9vZeatmZhTRW4ErgYuBbwJfBf4E2KCqV3ujm76IOxKxCbhznsR1BW7b/yvAOuDHqpr2P4hEZBtu0+DoPs1FuJ8F68jAz1hOJQljjDEzk1PNTcYYY2bGkoQxxpgpWZIwxhgzJUsSxhhjpmRJwhhjzJQsSZisJSIfEpF/EJFvexORZvLYcm8BvHlBRD7hLVZoTEbZEFiTtUTkAeBjuAsUykwmionIMuAuVT0vHbGdCBFpVtVlcx2HyS1WkzBZSUTehTuL+JO4k6OKROSfROSvROQfReQdXrkzReTH3iKAP0yZeX8VsMxbTPFMEblZRB7yHvN+Een1js8SkedF5Lsi8lURafdqIR8XkW+IyA0i8uWUGfyj8X1SRPpE5F0iUigiv/C+xwoR+Xcvnh94kwMnvrYrUr7/m7zvf553fpqIfF9E/tKLaUU63l+TQzK55K192Vcmv4CH8JbCxp3B+1nvuAB3prYfd6bxcu/6e4Ave8fLgIdSnmvieXPK8Y3Al7zjjd5zvspYTf0uUpZ6Tnnco3hLPuOuJ+bgLlW/1bu2Ffi3Kb5n6vFdeEtY4y7H8Wbv+Dzg3rn+d7Cvhf2VdQv8GTOFjUCXiFzvnb+Iu9TJMO5y2Z3Actz1eE7EqwCqulNE/gBIAJ/xKhBR3KVUJvq/wF+IyM3Aa6qaEJEocIWIXOI9pmaSxx3LRuDtInIubjKcbDl8Y6bNkoTJFS8Abap6G4CI/CHQBXwb96/t74vI24EPeeXjgHhlNwMteB/04u4CNvHDO7Vz70VgWFVv9spvxU0UE92DuwNbDHd9LHA3j+lR1b8XkVOAM6Z4PQMiUqruIm5NKddfAO7xklUe8O6p3hBjpsOShMlK3qqZS3FrCbfjNjd9SURuAMqAfeouuf7PuH/NL8dd9nmjiGzHXcAwLCJfA15R1ee9tv8v4G6Z2iciV+MusX0u7oZQL6nqDlXd5Y2o+hruGv9VuB/+46hqVES+B9Srap93+W7gi94HfBB3JdcLgTVAmYh8UFV/hFsL+Qdxt7RN4G6G8wzwp8C1IrLfez3/PKtvrMk5NrrJGGPMlGx0kzHGmClZkjDGGDMlSxLGGGOmZEnCGGPMlCxJGGOMmZIlCWOMMVOyJGGMMWZKliSMMcZM6f8DzCzvo8UQN4YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6a919975f8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_W, _b = sess.run([W, b])\n",
    "x = np.linspace(-2, 2, 100)\n",
    "\n",
    "def sigmoid(logits):\n",
    "    return 1 / (1 + np.exp(-logits))\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "for d in range(D // 2):\n",
    "    logits = _W[0][0][d] * (x - _b[0][0][d])\n",
    "    psx = sigmoid(logits)\n",
    "    plt.plot(x, psx)\n",
    "    \n",
    "plt.xlabel('feature value')\n",
    "plt.ylabel('$p(s|x)$')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should illustrate that the probability of the feature value being observed when it is below the feature mean should be close to 1, while the probability of being observed above the feature mean should be close to 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Close the session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# sess.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:python36] *",
   "language": "python",
   "name": "conda-env-python36-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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}