stochastic_segmentation_networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#  Stochastic Segmentation Networks: Modelling Spatially Correlated Aleatoric Uncertainty\n",
    "\n",
    "Monteiro et al. 2020, accepted at NeurIPS 2020.\n",
    "\n",
    "![caption](files/covmat.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "import tqdm.notebook as tqdm\n",
    "\n",
    "plt.rcParams.update({\n",
    "    'font.size': 8,\n",
    "    'text.usetex': True,\n",
    "    'text.latex.preamble': r'\\usepackage{amsfonts}'\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "MARGIN = 8\n",
    "DIM = 32\n",
    "\n",
    "rng = np.random.RandomState(7)\n",
    "\n",
    "def build_input(batch_size, shape):\n",
    "    h, w = shape\n",
    "    img = np.ones((batch_size, h, w), dtype=np.float32) + 0.5*rng.randn(h, w)\n",
    "    line_pos = rng.randint(0, h - MARGIN - 1, size=batch_size)\n",
    "    for i, lp in enumerate(line_pos):\n",
    "        img[i, lp:lp+MARGIN] = 0.9*rng.randn(MARGIN, w)\n",
    "        img[i, :lp] = -1.0 + 0.5*rng.randn(lp, w)\n",
    "    return img.astype(np.float32), line_pos\n",
    "\n",
    "def build_output(inp, line_pos):\n",
    "    batch_size = inp.shape[0]\n",
    "    offsets = line_pos\n",
    "    output = np.zeros_like(inp)\n",
    "    line_pos = rng.randint(0, MARGIN, size=batch_size) + offsets\n",
    "    for i, lp in enumerate(line_pos):\n",
    "        output[i, lp:] = 1.0\n",
    "    return output.astype(np.float32)\n",
    "\n",
    "def produce_data(batch_size, dim):\n",
    "    while True:\n",
    "        xs, line_pos = build_input(batch_size, (dim, dim))\n",
    "        ys = build_output(xs, line_pos)\n",
    "        yield np.expand_dims(xs, 1), np.expand_dims(ys, 1)\n",
    "        \n",
    "def plot_im(xs, cmap, title=None):\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
    "    xs = xs.squeeze().T[:16]\n",
    "    ax.imshow(xs, cmap, interpolation=\"none\")\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(title, fontsize=40, pad=20)\n",
    "\n",
    "def sample_bernoulli_im(im, n_samples=1):\n",
    "    im = np.concatenate([1 - im, im], axis=1)[0].transpose(1, 2, 0)\n",
    "    n_ims = np.tile(im[None, ...], [n_samples, 1, 1, 1])\n",
    "    sample = (n_ims.cumsum(-1) >= rng.uniform(size=n_ims.shape[:-1])[..., None]).argmax(-1)\n",
    "    return sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "inp, lp = build_input(1, (DIM, DIM))\n",
    "out = np.concatenate([build_output(inp, lp) for _ in range(100)], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Some data is messy \n",
    "\n",
    "For the same image $x$ we have multiple ambigous labels $y_{1,2,\\dots}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADoCAYAAADR9T5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAARU0lEQVR4nO3da3CU93XH8fPsStqVhG4ISYBkS0gCDDYXKwICxDZObDdOHV+Cce20aZq0NU1eJJPJOJ1Ob5OZtC5+0ek0TTKk07RjO2mLnTa0tZ0W29CaiQ3mWhNuRlzEXcJi0QVppd19+oLVjGoL9reO8ZHa7+eNsfTjPKvbT8+IPfoHYRgaAODDF/F+AADw/xUFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAMakEQfB4EAQ7gyDoCILgYhAE69/1+sogCJ7Lvr4jCIJNQRC0eT1e4FooYEwK2WLdaWYtZvaJMAxbzOwjZvZIEAQd2Uyzme00syfDMGzJZszMdmZfB0woFDAmi1fMbH0Yhr8bhmHCzCwMw6Nm9qSZNQdB8LiZPWdma8Iw3GVmlr07viv79+9670jAFwWMCS9bri+HYfj9cV59NPvf9Wa2Y7R8s9rH/Pnl6/X4gPerwPsBANcSBEGlma0zs1lXiSwZ8+f173rdGjN72K6U91EDJpiAU5ExkWXvfj8ShuHaq7x+k1358UIiDMOqD/XBAb8g7oAx0e0wsw3XeP3ojxn4EQMmHQoYE9q7fqb7v2R/PFGZ/d9NH8bjAT5I/CMcJrOxz2zgDhiTDgWMyezu7H8T/CMbJiMKGJPZ6B0wd7+YlChgTErZn/+Obre9eY0cCxiYsChgTFZji/Wq/1Bn731uMDBhUMCY0LK/fOfxcV71K2P+vONqf9fMnr8uDwz4ALCIgQkrCIJ1ZvaN7P+uCcPw+ezLK83sYvbl4y5gZDOv2JVf3JO47g8WeB+4A8ZE9vCYP499lsM6M1ubfVnlVf7uc3blt6IlrssjAz4AFDAmstFnN4z9DWfrzGxq9hfzrB3zMsv+uS37ayufG71jBiYqfgSBCS0Igm/YlaJN2JW73fVhGD415vVtduWOePQZEbvsyp3vtf5hDpgQKGAAcMKPIADACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOCvIJT6mYFk6tbcqZixXqxxz1X9ZypcWBPNNMu36sIC1PLAhHpNyFy3F9ZoH2NhUXZuSZgyPa99SCPL71Vqa7pVy6qFieORCWSrlIoH0sS4e65GtnYtq1LdEjz4yUl0u5IJ2SZ5p4XFg4on1umpmFpWVaLtA/QSLDQ1owjy/hQHzbM0Pitc0sKCoSk3p/hcPa+37Pqa4LYRjWvPvleRXw1Nome+Iv3syZa6rTi+1ne7RPyKULC+WZgfhFO6tS/wKryZyVcj/YdbM8s3ZaVMrdNFP8LmVmB8+USLmKUv2T7IGe70u5gcYF8szXh5dKubKY9gm+dP935GsPztGuHW58Vp4Z//g9Ui6a0L6ZmZlZRvs6Sp05o49cskqbWah/My0+c0gL5lHqQUr7uCcP7pdnFjU1aUHx/W5mNtx5UspVPvHtE+O9nB9BAIATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE7yeh5wMpmxo8cGcuamlotPdDezu5dpz/frS+rPWz3Xoz2/tivxnudFX9XcGVOk3KLZ8kjb+7b2fMOXNurP83zw4VYp90Dir+WZe1sflXLlBbk/N0YN9InPga48JeUOLf4N+done6ukXOeS2+WZa17Srl/afKM8M9LQJOWiVdrbY2YWprWvt6Gy6fLMospaKRdJ6p8f6f2HpdzwJz8rz4wd36MF1cUSM0v2XJKz4+EOGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADgJK9V5Hg8YrNbc68Zj+Rx7NV/7tZWUoeH9aEfXazlRlL6IVXf+cF5KfeHX9KPcpnZpp1RVVraIs9UHWz9jJzdd7JCys2boX0szczu2//HUu5M3W9KuRsHD8jXnlqhraA3VegrvuXpZVKuf9at8szSrqNSLlM/S54Z3b9Dyk2t0VbAzczOz10l5Sp7O+WZvXf9upQryOjn4akylxJyNv4ZcRX6yafHfTF3wADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOAkr024gqhZbWUmZy6d0TfMMuJZm4vna1tjZmZbd2iH6nWf7ZVnVkwrk3IZy/3+GfXiznIxqc88cXxQyh05oW23mZmtWdEl5d7s1A5nNDMLVj4h5aZZj5QrPvu2fO3TLXOkXOvuZ+SZI+e0TcnIDTfLM082fkzK1WxYJ88sXrRIymWmVMoza/e+KOXCmpnyzMKSainX/cTX5ZlVX35MyqV6++SZhT2n5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxW5MJKxurLcq65Hu3Mf3DmqtbFQyv3o2cPyzPtXa6umJSX6mz+vRctuOaivYXed11YeVyyZIs9sqteyN1ZrK8tmZjOSx6Tcls1peebCz2mHlyZS2sr0qQZtzdTM7Kk/PSLl1n75a/LMe+qel3LBYEKeecO5l7WZCxbKMy2lHW4bjAzLI4NYXMql9++VZ3bPXCXlmn/vK/LMMKl9zgeB/jVsPRf07Di4AwYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAneW3CDY5E7K2Tubfc+gf0QySnVWnfAx797Gx55o9/rG1urbyzSZ55SjuX0l7fclSe+fgXGqTcsfP6Zk5FqXbKaW38ojyz226Ucre2T5dn9iaTUq5nICblyoq1DS8zs2e+PiDlSs/9UJ6Z7jgk5aLl6kGsZgML7pByBS/+SJ7Zs/qrUm5Gx3/JM8NC7WN0+faH5Jmt2/9WyqVaF8gzR4q19328uVWfOa1ezo6HO2AAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgJO8VpGLCkJrrM298pkY0MeWxLS15f/Y0ivPnDVXW4kN9Y1pWzpHO9Cvv19b2zUzO3xaez8VaeeWmpnZgunnpVws0A/lLB/slnL3L9DWi83MBjPaoZxm2qGcN8Xelq89EorXjuifx8nln5JypR275Jmlp/ZLudSqe+WZda/+nZTrXvWr8syp3doadtnpffJMq50hxdKF2oGgZmbJDU9LuaL5+q89COsa5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJ3ltwpmZRYLchz4umdkpz9tyVNscW7akUp5588yElDt0Xt3GMut8p0TKHT+ibaKZmVW110m5c936gZObR7QtwDn1w/LMtoILUu7sUI08s73z76Xc9qLfknJl9TfI145FR6RcukH/8qh+81+lXChueJmZWeIdKTYyXT9E0pbcLcWmbf1HeWTyVu3w0MEXN8ozYw88ouUS+tdbtEH72uhr/6Q8s/hyj5wdD3fAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnea0iZzKB9Q7m/isH0vXyzIOHL0u51hZtFdjM7IXtU6Tc5ctD8swj+05LueZ5+qrpgf0XpdxXV+sHaG47OVPKLQu3yjML+vqlXNuBl+SZHbf/jpR7cPBn2kBtu9jMzN64vFTKzT3zT/LMVPN8KRc9qh20aWY2MrdNm/nqP8szj/zyH0m5m6f9XJ4Zu6StA8dvmSfPTEa1k2jDE3kcxnrbfVKufPfL8szMoNZfV8MdMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxMuHZr1DwY5cz/dpB0maGa2fEWtlLvQk5ZnplIZKdcyKy7PHBEPuxwe0R/nn3zqgJT7q23a5paZ2ZzG3IemmpntiK6UZzbWnpNyDScPyzMberVNq22R26Tcx/r+Rb52e/UeKTdUP1ueGbuovY9shnYIrZnZcLxcyvX9/Jg8c/Zt26RcWFwqz4z0agdTphrnyjOjKfHQ2Mh1uIcMtf4wM0uKh5yafWvcl3IHDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwktcqchiapYRN27qZ2gqlmdmBQwNSbvu/75Rn3vPoCikX0879MzOz5lnaamZVmT7z1f7l2syK3Ovfo4ZTWm7Fzj+TZ25e9PtSrmfhb8szCyPayvasqLbiu7n9D+Rrr/zWXVLurU+vk2e2n/22FqzVD22Nv6atV5d+/KPyTDuhrb8HMX1NPyytkHLRjn3yzMsL75ByJXk8zsKkdrhsZpZ+eGjBlo1ydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhohGzKcW5c3cu0Xt9b0eRlPven98kz9x+QnuzOjpH5JmJniEpt2LJFHlmEGgHaG78ob4FuPrz7dq1G5rkmQvKjki548P6gZPJtLaGWBXRDnzs/LeD8rUXN2sHgrZ1/VSeGVRUSbmh3frHMj5bOxQ0faFLntm9RTuUs2bpLfLMSIu2CZecv0yeGe/VNiDDmY3yzIKENjNZo8+MzdF7aTzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnOS1ihwEZvGiTM5c76B+2uXhA+9IuWlVtfLM3n5txbemWn+cp45flHKzqrVrm5ldHBL2us1swXJ93TE5rF1/e/Mj8sx4eljKhaF+eGhJYVLKlQ1dkHJ3tpbI1w4z2uN8ZmC1PPPz8aelXHTVvfLM4UB7nJHp+upszfw2bWaf9vluZpYu1tbvYxc65Zmp8hop11OtrWubmVVu/K6Uiy2OyjPDrrNydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhEpdGbOMLuQ8A/NJjpfLMzz2oZfcel0daRlxGe/21U/LMT9zdIOX2nZZHWs+l3FuFZmaNTfr7c2hIe+OXnt0gz7y06RUp1/nFv5Rnzgi0d9SUY7ul3NbpS+RrLy89I+Vmr10ozwy/uUbKFQ4NyjM3139Ryi0qPyDPDE3brisu1TfhSo7tlXL9e/9bnhmNaYf11s0/Ls9MlWibp5fqF8gzC7a+JmfHwx0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcBJXqvIdVMD+9qv5T7IctfJMnnmvgMDUq59kTzSKkvTUi4M6+WZp86mpFw8rn9PSyS0wy5/abl+0OexLm3dsqdeO5zRzCz6BW01c8/JSnlm//S4lKusPiLl5lVp68VmZu+Y9nFfLq4Xm5lZRlsr73tjmzxy1SrtYxn2xuSZXc/+g5Qru6VFnmnNWjZ5sU8eWb1qpZQb2KOvNxdVV0m5yoNb5ZnD94mfI9/8m3FfzB0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOMnvUM6BqP3kjfKcuVhM20QzM/v0HVEpt+OwtmlkZtYwd0TKPbRIP5TzqQ25324zs9tXTpVnts3Rvv+9dVw7oNDM7I653VIuUE8uNbOdg9rhlI+VvyDPLDh9XsplTh2Xcrsy98vXvj/4iZRLzdcP+uytuEHKVU9/XZ5p57XtvsEjR+WR0x+6V8qlKuvkmUPxKVJu8Cvatc3MBvuOS7mi0/rXsLqtaBHt4FIzs3SB/rU57qV+ob8NAHjfKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwEoShvpIaBEG3mZ24fg8HAP5PagzDsObdL8yrgAEAHxx+BAEATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABO/gehl8qxPJwTAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(inp[0], cmap=\"coolwarm\", title=r\"$$x$$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADoCAYAAADR9T5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAARU0lEQVR4nO3da3CU93XH8fPsStqVhG4ISYBkS0gCDDYXKwICxDZObDdOHV+Cce20aZq0NU1eJJPJOJ1Ob5OZtC5+0ek0TTKk07RjO2mLnTa0tZ0W29CaiQ3mWhNuRlzEXcJi0QVppd19+oLVjGoL9reO8ZHa7+eNsfTjPKvbT8+IPfoHYRgaAODDF/F+AADw/xUFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAMakEQfB4EAQ7gyDoCILgYhAE69/1+sogCJ7Lvr4jCIJNQRC0eT1e4FooYEwK2WLdaWYtZvaJMAxbzOwjZvZIEAQd2Uyzme00syfDMGzJZszMdmZfB0woFDAmi1fMbH0Yhr8bhmHCzCwMw6Nm9qSZNQdB8LiZPWdma8Iw3GVmlr07viv79+9670jAFwWMCS9bri+HYfj9cV59NPvf9Wa2Y7R8s9rH/Pnl6/X4gPerwPsBANcSBEGlma0zs1lXiSwZ8+f173rdGjN72K6U91EDJpiAU5ExkWXvfj8ShuHaq7x+k1358UIiDMOqD/XBAb8g7oAx0e0wsw3XeP3ojxn4EQMmHQoYE9q7fqb7v2R/PFGZ/d9NH8bjAT5I/CMcJrOxz2zgDhiTDgWMyezu7H8T/CMbJiMKGJPZ6B0wd7+YlChgTErZn/+Obre9eY0cCxiYsChgTFZji/Wq/1Bn731uMDBhUMCY0LK/fOfxcV71K2P+vONqf9fMnr8uDwz4ALCIgQkrCIJ1ZvaN7P+uCcPw+ezLK83sYvbl4y5gZDOv2JVf3JO47g8WeB+4A8ZE9vCYP499lsM6M1ubfVnlVf7uc3blt6IlrssjAz4AFDAmstFnN4z9DWfrzGxq9hfzrB3zMsv+uS37ayufG71jBiYqfgSBCS0Igm/YlaJN2JW73fVhGD415vVtduWOePQZEbvsyp3vtf5hDpgQKGAAcMKPIADACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOCvIJT6mYFk6tbcqZixXqxxz1X9ZypcWBPNNMu36sIC1PLAhHpNyFy3F9ZoH2NhUXZuSZgyPa99SCPL71Vqa7pVy6qFieORCWSrlIoH0sS4e65GtnYtq1LdEjz4yUl0u5IJ2SZ5p4XFg4on1umpmFpWVaLtA/QSLDQ1owjy/hQHzbM0Pitc0sKCoSk3p/hcPa+37Pqa4LYRjWvPvleRXw1Nome+Iv3syZa6rTi+1ne7RPyKULC+WZgfhFO6tS/wKryZyVcj/YdbM8s3ZaVMrdNFP8LmVmB8+USLmKUv2T7IGe70u5gcYF8szXh5dKubKY9gm+dP935GsPztGuHW58Vp4Z//g9Ui6a0L6ZmZlZRvs6Sp05o49cskqbWah/My0+c0gL5lHqQUr7uCcP7pdnFjU1aUHx/W5mNtx5UspVPvHtE+O9nB9BAIATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE7yeh5wMpmxo8cGcuamlotPdDezu5dpz/frS+rPWz3Xoz2/tivxnudFX9XcGVOk3KLZ8kjb+7b2fMOXNurP83zw4VYp90Dir+WZe1sflXLlBbk/N0YN9InPga48JeUOLf4N+done6ukXOeS2+WZa17Srl/afKM8M9LQJOWiVdrbY2YWprWvt6Gy6fLMospaKRdJ6p8f6f2HpdzwJz8rz4wd36MF1cUSM0v2XJKz4+EOGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADgJK9V5Hg8YrNbc68Zj+Rx7NV/7tZWUoeH9aEfXazlRlL6IVXf+cF5KfeHX9KPcpnZpp1RVVraIs9UHWz9jJzdd7JCys2boX0szczu2//HUu5M3W9KuRsHD8jXnlqhraA3VegrvuXpZVKuf9at8szSrqNSLlM/S54Z3b9Dyk2t0VbAzczOz10l5Sp7O+WZvXf9upQryOjn4akylxJyNv4ZcRX6yafHfTF3wADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOAkr024gqhZbWUmZy6d0TfMMuJZm4vna1tjZmZbd2iH6nWf7ZVnVkwrk3IZy/3+GfXiznIxqc88cXxQyh05oW23mZmtWdEl5d7s1A5nNDMLVj4h5aZZj5QrPvu2fO3TLXOkXOvuZ+SZI+e0TcnIDTfLM082fkzK1WxYJ88sXrRIymWmVMoza/e+KOXCmpnyzMKSainX/cTX5ZlVX35MyqV6++SZhT2n5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxW5MJKxurLcq65Hu3Mf3DmqtbFQyv3o2cPyzPtXa6umJSX6mz+vRctuOaivYXed11YeVyyZIs9sqteyN1ZrK8tmZjOSx6Tcls1peebCz2mHlyZS2sr0qQZtzdTM7Kk/PSLl1n75a/LMe+qel3LBYEKeecO5l7WZCxbKMy2lHW4bjAzLI4NYXMql9++VZ3bPXCXlmn/vK/LMMKl9zgeB/jVsPRf07Di4AwYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAneW3CDY5E7K2Tubfc+gf0QySnVWnfAx797Gx55o9/rG1urbyzSZ55SjuX0l7fclSe+fgXGqTcsfP6Zk5FqXbKaW38ojyz226Ucre2T5dn9iaTUq5nICblyoq1DS8zs2e+PiDlSs/9UJ6Z7jgk5aLl6kGsZgML7pByBS/+SJ7Zs/qrUm5Gx3/JM8NC7WN0+faH5Jmt2/9WyqVaF8gzR4q19328uVWfOa1ezo6HO2AAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgJO8VpGLCkJrrM298pkY0MeWxLS15f/Y0ivPnDVXW4kN9Y1pWzpHO9Cvv19b2zUzO3xaez8VaeeWmpnZgunnpVws0A/lLB/slnL3L9DWi83MBjPaoZxm2qGcN8Xelq89EorXjuifx8nln5JypR275Jmlp/ZLudSqe+WZda/+nZTrXvWr8syp3doadtnpffJMq50hxdKF2oGgZmbJDU9LuaL5+q89COsa5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJ3ltwpmZRYLchz4umdkpz9tyVNscW7akUp5588yElDt0Xt3GMut8p0TKHT+ibaKZmVW110m5c936gZObR7QtwDn1w/LMtoILUu7sUI08s73z76Xc9qLfknJl9TfI145FR6RcukH/8qh+81+lXChueJmZWeIdKTYyXT9E0pbcLcWmbf1HeWTyVu3w0MEXN8ozYw88ouUS+tdbtEH72uhr/6Q8s/hyj5wdD3fAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnea0iZzKB9Q7m/isH0vXyzIOHL0u51hZtFdjM7IXtU6Tc5ctD8swj+05LueZ5+qrpgf0XpdxXV+sHaG47OVPKLQu3yjML+vqlXNuBl+SZHbf/jpR7cPBn2kBtu9jMzN64vFTKzT3zT/LMVPN8KRc9qh20aWY2MrdNm/nqP8szj/zyH0m5m6f9XJ4Zu6StA8dvmSfPTEa1k2jDE3kcxnrbfVKufPfL8szMoNZfV8MdMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxMuHZr1DwY5cz/dpB0maGa2fEWtlLvQk5ZnplIZKdcyKy7PHBEPuxwe0R/nn3zqgJT7q23a5paZ2ZzG3IemmpntiK6UZzbWnpNyDScPyzMberVNq22R26Tcx/r+Rb52e/UeKTdUP1ueGbuovY9shnYIrZnZcLxcyvX9/Jg8c/Zt26RcWFwqz4z0agdTphrnyjOjKfHQ2Mh1uIcMtf4wM0uKh5yafWvcl3IHDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwktcqchiapYRN27qZ2gqlmdmBQwNSbvu/75Rn3vPoCikX0879MzOz5lnaamZVmT7z1f7l2syK3Ovfo4ZTWm7Fzj+TZ25e9PtSrmfhb8szCyPayvasqLbiu7n9D+Rrr/zWXVLurU+vk2e2n/22FqzVD22Nv6atV5d+/KPyTDuhrb8HMX1NPyytkHLRjn3yzMsL75ByJXk8zsKkdrhsZpZ+eGjBlo1ydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhohGzKcW5c3cu0Xt9b0eRlPven98kz9x+QnuzOjpH5JmJniEpt2LJFHlmEGgHaG78ob4FuPrz7dq1G5rkmQvKjki548P6gZPJtLaGWBXRDnzs/LeD8rUXN2sHgrZ1/VSeGVRUSbmh3frHMj5bOxQ0faFLntm9RTuUs2bpLfLMSIu2CZecv0yeGe/VNiDDmY3yzIKENjNZo8+MzdF7aTzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnOS1ihwEZvGiTM5c76B+2uXhA+9IuWlVtfLM3n5txbemWn+cp45flHKzqrVrm5ldHBL2us1swXJ93TE5rF1/e/Mj8sx4eljKhaF+eGhJYVLKlQ1dkHJ3tpbI1w4z2uN8ZmC1PPPz8aelXHTVvfLM4UB7nJHp+upszfw2bWaf9vluZpYu1tbvYxc65Zmp8hop11OtrWubmVVu/K6Uiy2OyjPDrrNydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhEpdGbOMLuQ8A/NJjpfLMzz2oZfcel0daRlxGe/21U/LMT9zdIOX2nZZHWs+l3FuFZmaNTfr7c2hIe+OXnt0gz7y06RUp1/nFv5Rnzgi0d9SUY7ul3NbpS+RrLy89I+Vmr10ozwy/uUbKFQ4NyjM3139Ryi0qPyDPDE3brisu1TfhSo7tlXL9e/9bnhmNaYf11s0/Ls9MlWibp5fqF8gzC7a+JmfHwx0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcBJXqvIdVMD+9qv5T7IctfJMnnmvgMDUq59kTzSKkvTUi4M6+WZp86mpFw8rn9PSyS0wy5/abl+0OexLm3dsqdeO5zRzCz6BW01c8/JSnlm//S4lKusPiLl5lVp68VmZu+Y9nFfLq4Xm5lZRlsr73tjmzxy1SrtYxn2xuSZXc/+g5Qru6VFnmnNWjZ5sU8eWb1qpZQb2KOvNxdVV0m5yoNb5ZnD94mfI9/8m3FfzB0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOMnvUM6BqP3kjfKcuVhM20QzM/v0HVEpt+OwtmlkZtYwd0TKPbRIP5TzqQ25324zs9tXTpVnts3Rvv+9dVw7oNDM7I653VIuUE8uNbOdg9rhlI+VvyDPLDh9XsplTh2Xcrsy98vXvj/4iZRLzdcP+uytuEHKVU9/XZ5p57XtvsEjR+WR0x+6V8qlKuvkmUPxKVJu8Cvatc3MBvuOS7mi0/rXsLqtaBHt4FIzs3SB/rU57qV+ob8NAHjfKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwEoShvpIaBEG3mZ24fg8HAP5PagzDsObdL8yrgAEAHxx+BAEATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABO/gehl8qxPJwTAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADoCAYAAADR9T5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAF5ElEQVR4nO3d0Y0TVxiG4f9EFICQch0JOlhIB9BBSAewHYBSQbR0AKmAbDqADhLTAZZyj0g6OLnwTFgc2+wGzHcGnkda2Tser+bq1dG/4+PWey8APr9v0hcA8LUSYIAQAQYIEWCAEAEGCBFggBABBggRYIAQAQYIEWCAEAEGCBFggBABBggRYIAQAQYIEWCAEAFmOK21R62119PPeWvt+gfOvzmd++IzXSJ8EtfSFwAXtdbOq2rde781/d6r6qSqbh142+Oqujn9wGJYATOM1tqjqnrbe3+89dLN1trJgbfenR5fHufK4DisgBnJ6bzyraraGj3c2PWG6Zx55WsEwaJYATOE1trDqvpt6/CPF57/seetdy88twJmUayAGcVpVd3fcayq6lXv/e897/t+ftJ7f3WE64KjsQJmFD/33tfzL9NoYZ77Pj3wvkvNf6c7K374qCuET0yAGULv/dD44dcDb50jvXP1O92idl5VZx9xeXAUAsyo5vHDy33jh9baxfnvi63XHk73BZ9W1fOjXCF8JDNghrM1fjg/cOq9+Unv/b0RRO/9WVU9m/6e+4MZkhUwI7rs+GFeAfvnG4skwIxoXtmuD9z9UPVulez2MxZJgBnRPDLYu7I9NP+FpRBgRnR9enx74Jy9819YCgFmROsPn1IPp0fzXxZLgBnRfOfDnV0vttZW9W6VbPXLYgkww5luIXtZVSfTHhFV9e+HKlb1/grZ/JfFEmCG1Hu/V5t9fh+31v6awntWVQ+q6vcLp+7bpAeG54MYDKv3/qSqnmwfb639Mj3d+yk5WAIrYBZl+lTbZT4lB8MTYIbRWrs7jRteH/geuHlTnfU0K4bFEmBGcl6buxt2fr/btPqdt5Tc3jsYFkeAGcn16fHVns3V532BT6+4+fr8d23Kw1AEmJGsq+p+7/329guttae12XznyWVGD9M440Vr7XVVrabDZ/PX11+8vQ1SWu89fQ1QVVXTN1b8VJsIr6djJ1X1S21Wrw92bNwOiyXADGXaZOesNmODG7VZFT+fbkmDL4oAA4SYAQOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQMiVvpKoteZzy3w1bt/+z6Zs8L+sVqs3vfdvt49faS8IAeZrYp8UPpXW2qr3fmf7uBEEQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMEDItSue/6aq/jzGhcBoWmvpS+DL8d2ug633/rkvBIAyggCIEWCAEAEGCBFggBABBggRYIAQAQYIEWCAEAEGCPkH9hvH0UDNbksAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(inp[0], cmap=\"coolwarm\", title=r\"$$x$$\")\n",
    "plot_im(out[0], cmap=\"gray\", title=r\"$$y_1$$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADoCAYAAADR9T5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAARU0lEQVR4nO3da3CU93XH8fPsStqVhG4ISYBkS0gCDDYXKwICxDZObDdOHV+Cce20aZq0NU1eJJPJOJ1Ob5OZtC5+0ek0TTKk07RjO2mLnTa0tZ0W29CaiQ3mWhNuRlzEXcJi0QVppd19+oLVjGoL9reO8ZHa7+eNsfTjPKvbT8+IPfoHYRgaAODDF/F+AADw/xUFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAMakEQfB4EAQ7gyDoCILgYhAE69/1+sogCJ7Lvr4jCIJNQRC0eT1e4FooYEwK2WLdaWYtZvaJMAxbzOwjZvZIEAQd2Uyzme00syfDMGzJZszMdmZfB0woFDAmi1fMbH0Yhr8bhmHCzCwMw6Nm9qSZNQdB8LiZPWdma8Iw3GVmlr07viv79+9670jAFwWMCS9bri+HYfj9cV59NPvf9Wa2Y7R8s9rH/Pnl6/X4gPerwPsBANcSBEGlma0zs1lXiSwZ8+f173rdGjN72K6U91EDJpiAU5ExkWXvfj8ShuHaq7x+k1358UIiDMOqD/XBAb8g7oAx0e0wsw3XeP3ojxn4EQMmHQoYE9q7fqb7v2R/PFGZ/d9NH8bjAT5I/CMcJrOxz2zgDhiTDgWMyezu7H8T/CMbJiMKGJPZ6B0wd7+YlChgTErZn/+Obre9eY0cCxiYsChgTFZji/Wq/1Bn731uMDBhUMCY0LK/fOfxcV71K2P+vONqf9fMnr8uDwz4ALCIgQkrCIJ1ZvaN7P+uCcPw+ezLK83sYvbl4y5gZDOv2JVf3JO47g8WeB+4A8ZE9vCYP499lsM6M1ubfVnlVf7uc3blt6IlrssjAz4AFDAmstFnN4z9DWfrzGxq9hfzrB3zMsv+uS37ayufG71jBiYqfgSBCS0Igm/YlaJN2JW73fVhGD415vVtduWOePQZEbvsyp3vtf5hDpgQKGAAcMKPIADACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOCvIJT6mYFk6tbcqZixXqxxz1X9ZypcWBPNNMu36sIC1PLAhHpNyFy3F9ZoH2NhUXZuSZgyPa99SCPL71Vqa7pVy6qFieORCWSrlIoH0sS4e65GtnYtq1LdEjz4yUl0u5IJ2SZ5p4XFg4on1umpmFpWVaLtA/QSLDQ1owjy/hQHzbM0Pitc0sKCoSk3p/hcPa+37Pqa4LYRjWvPvleRXw1Nome+Iv3syZa6rTi+1ne7RPyKULC+WZgfhFO6tS/wKryZyVcj/YdbM8s3ZaVMrdNFP8LmVmB8+USLmKUv2T7IGe70u5gcYF8szXh5dKubKY9gm+dP935GsPztGuHW58Vp4Z//g9Ui6a0L6ZmZlZRvs6Sp05o49cskqbWah/My0+c0gL5lHqQUr7uCcP7pdnFjU1aUHx/W5mNtx5UspVPvHtE+O9nB9BAIATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE7yeh5wMpmxo8cGcuamlotPdDezu5dpz/frS+rPWz3Xoz2/tivxnudFX9XcGVOk3KLZ8kjb+7b2fMOXNurP83zw4VYp90Dir+WZe1sflXLlBbk/N0YN9InPga48JeUOLf4N+done6ukXOeS2+WZa17Srl/afKM8M9LQJOWiVdrbY2YWprWvt6Gy6fLMospaKRdJ6p8f6f2HpdzwJz8rz4wd36MF1cUSM0v2XJKz4+EOGACcUMAA4IQCBgAnFDAAOKGAAcAJBQwATihgAHBCAQOAEwoYAJxQwADgJK9V5Hg8YrNbc68Zj+Rx7NV/7tZWUoeH9aEfXazlRlL6IVXf+cF5KfeHX9KPcpnZpp1RVVraIs9UHWz9jJzdd7JCys2boX0szczu2//HUu5M3W9KuRsHD8jXnlqhraA3VegrvuXpZVKuf9at8szSrqNSLlM/S54Z3b9Dyk2t0VbAzczOz10l5Sp7O+WZvXf9upQryOjn4akylxJyNv4ZcRX6yafHfTF3wADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOAkr024gqhZbWUmZy6d0TfMMuJZm4vna1tjZmZbd2iH6nWf7ZVnVkwrk3IZy/3+GfXiznIxqc88cXxQyh05oW23mZmtWdEl5d7s1A5nNDMLVj4h5aZZj5QrPvu2fO3TLXOkXOvuZ+SZI+e0TcnIDTfLM082fkzK1WxYJ88sXrRIymWmVMoza/e+KOXCmpnyzMKSainX/cTX5ZlVX35MyqV6++SZhT2n5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxW5MJKxurLcq65Hu3Mf3DmqtbFQyv3o2cPyzPtXa6umJSX6mz+vRctuOaivYXed11YeVyyZIs9sqteyN1ZrK8tmZjOSx6Tcls1peebCz2mHlyZS2sr0qQZtzdTM7Kk/PSLl1n75a/LMe+qel3LBYEKeecO5l7WZCxbKMy2lHW4bjAzLI4NYXMql9++VZ3bPXCXlmn/vK/LMMKl9zgeB/jVsPRf07Di4AwYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAneW3CDY5E7K2Tubfc+gf0QySnVWnfAx797Gx55o9/rG1urbyzSZ55SjuX0l7fclSe+fgXGqTcsfP6Zk5FqXbKaW38ojyz226Ucre2T5dn9iaTUq5nICblyoq1DS8zs2e+PiDlSs/9UJ6Z7jgk5aLl6kGsZgML7pByBS/+SJ7Zs/qrUm5Gx3/JM8NC7WN0+faH5Jmt2/9WyqVaF8gzR4q19328uVWfOa1ezo6HO2AAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgJO8VpGLCkJrrM298pkY0MeWxLS15f/Y0ivPnDVXW4kN9Y1pWzpHO9Cvv19b2zUzO3xaez8VaeeWmpnZgunnpVws0A/lLB/slnL3L9DWi83MBjPaoZxm2qGcN8Xelq89EorXjuifx8nln5JypR275Jmlp/ZLudSqe+WZda/+nZTrXvWr8syp3doadtnpffJMq50hxdKF2oGgZmbJDU9LuaL5+q89COsa5ex4uAMGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJ3ltwpmZRYLchz4umdkpz9tyVNscW7akUp5588yElDt0Xt3GMut8p0TKHT+ibaKZmVW110m5c936gZObR7QtwDn1w/LMtoILUu7sUI08s73z76Xc9qLfknJl9TfI145FR6RcukH/8qh+81+lXChueJmZWeIdKTYyXT9E0pbcLcWmbf1HeWTyVu3w0MEXN8ozYw88ouUS+tdbtEH72uhr/6Q8s/hyj5wdD3fAAOCEAgYAJxQwADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnea0iZzKB9Q7m/isH0vXyzIOHL0u51hZtFdjM7IXtU6Tc5ctD8swj+05LueZ5+qrpgf0XpdxXV+sHaG47OVPKLQu3yjML+vqlXNuBl+SZHbf/jpR7cPBn2kBtu9jMzN64vFTKzT3zT/LMVPN8KRc9qh20aWY2MrdNm/nqP8szj/zyH0m5m6f9XJ4Zu6StA8dvmSfPTEa1k2jDE3kcxnrbfVKufPfL8szMoNZfV8MdMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADjJaxMuHZr1DwY5cz/dpB0maGa2fEWtlLvQk5ZnplIZKdcyKy7PHBEPuxwe0R/nn3zqgJT7q23a5paZ2ZzG3IemmpntiK6UZzbWnpNyDScPyzMberVNq22R26Tcx/r+Rb52e/UeKTdUP1ueGbuovY9shnYIrZnZcLxcyvX9/Jg8c/Zt26RcWFwqz4z0agdTphrnyjOjKfHQ2Mh1uIcMtf4wM0uKh5yafWvcl3IHDABOKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwktcqchiapYRN27qZ2gqlmdmBQwNSbvu/75Rn3vPoCikX0879MzOz5lnaamZVmT7z1f7l2syK3Ovfo4ZTWm7Fzj+TZ25e9PtSrmfhb8szCyPayvasqLbiu7n9D+Rrr/zWXVLurU+vk2e2n/22FqzVD22Nv6atV5d+/KPyTDuhrb8HMX1NPyytkHLRjn3yzMsL75ByJXk8zsKkdrhsZpZ+eGjBlo1ydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhohGzKcW5c3cu0Xt9b0eRlPven98kz9x+QnuzOjpH5JmJniEpt2LJFHlmEGgHaG78ob4FuPrz7dq1G5rkmQvKjki548P6gZPJtLaGWBXRDnzs/LeD8rUXN2sHgrZ1/VSeGVRUSbmh3frHMj5bOxQ0faFLntm9RTuUs2bpLfLMSIu2CZecv0yeGe/VNiDDmY3yzIKENjNZo8+MzdF7aTzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4ATChgAnOS1ihwEZvGiTM5c76B+2uXhA+9IuWlVtfLM3n5txbemWn+cp45flHKzqrVrm5ldHBL2us1swXJ93TE5rF1/e/Mj8sx4eljKhaF+eGhJYVLKlQ1dkHJ3tpbI1w4z2uN8ZmC1PPPz8aelXHTVvfLM4UB7nJHp+upszfw2bWaf9vluZpYu1tbvYxc65Zmp8hop11OtrWubmVVu/K6Uiy2OyjPDrrNydjzcAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABOKGAAcEIBA4CTvDbhEpdGbOMLuQ8A/NJjpfLMzz2oZfcel0daRlxGe/21U/LMT9zdIOX2nZZHWs+l3FuFZmaNTfr7c2hIe+OXnt0gz7y06RUp1/nFv5Rnzgi0d9SUY7ul3NbpS+RrLy89I+Vmr10ozwy/uUbKFQ4NyjM3139Ryi0qPyDPDE3brisu1TfhSo7tlXL9e/9bnhmNaYf11s0/Ls9MlWibp5fqF8gzC7a+JmfHwx0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOKGAAcBJXqvIdVMD+9qv5T7IctfJMnnmvgMDUq59kTzSKkvTUi4M6+WZp86mpFw8rn9PSyS0wy5/abl+0OexLm3dsqdeO5zRzCz6BW01c8/JSnlm//S4lKusPiLl5lVp68VmZu+Y9nFfLq4Xm5lZRlsr73tjmzxy1SrtYxn2xuSZXc/+g5Qru6VFnmnNWjZ5sU8eWb1qpZQb2KOvNxdVV0m5yoNb5ZnD94mfI9/8m3FfzB0wADihgAHACQUMAE4oYABwQgEDgBMKGACcUMAA4IQCBgAnFDAAOMnvUM6BqP3kjfKcuVhM20QzM/v0HVEpt+OwtmlkZtYwd0TKPbRIP5TzqQ25324zs9tXTpVnts3Rvv+9dVw7oNDM7I653VIuUE8uNbOdg9rhlI+VvyDPLDh9XsplTh2Xcrsy98vXvj/4iZRLzdcP+uytuEHKVU9/XZ5p57XtvsEjR+WR0x+6V8qlKuvkmUPxKVJu8Cvatc3MBvuOS7mi0/rXsLqtaBHt4FIzs3SB/rU57qV+ob8NAHjfKGAAcEIBA4ATChgAnFDAAOCEAgYAJxQwADihgAHACQUMAE4oYABwEoShvpIaBEG3mZ24fg8HAP5PagzDsObdL8yrgAEAHxx+BAEATihgAHBCAQOAEwoYAJxQwADghAIGACcUMAA4oYABwAkFDABO/gehl8qxPJwTAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADoCAYAAADR9T5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAGhUlEQVR4nO3dUW4TRxzH8f9UHCCt1OdK4QYJvQG5AXADyA1APQEqN4CeoA03gJ6gDTcgVd8R5QbTB8+2xrWdhCb+beDzkSzb61lrn74aza7XrfdeAOzeV+kDAPhSCTBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAMzuttcettbfjcdJa2ztn/P4Y+2pHhwhX4lb6AGBZa+2kqs5677fH+15VB1V1e8tuT6pqfzzgxjADZjZaa4+r6n3v/cnKR/uttYMtu94dz6+v58jgepgBMyfH08y3qmpl6eGbdTuMMdPM1xIEN4oZMLPQWntUVS9XNj9Yev37hl3vLr02A+ZGMQNmLo6r6v6abVVVb3rvHzbs9/30ovf+5hqOC66NADMXT3vvZ9ObsbQwrfs+37LfxvXfsW58XIvli4Oq+lBVP/fen13B8cL/JsDMQu992/LDL1t2nSL90ey3tXavqvZ778dL2w6q6qS1dlxVh1tm1bAT1oCZqymcrzeFsrW2vP77amn7fi3i+9FMdyxRHNXipN2vV3q08AkEmNlZWX442TL0aHrRe19egjiuqhfrdhjLHC+r6uCcS9vg2gkwc3TR5YdpBrx68u1uVf0xZsLr/Dae73zCscGVEWDmaJrZnp2zTjvNYFdPwL2vqr3xWGf6zk2fw044CcccTTPXjZeVbVr/rarqvR+11va2xPvwvO+HXTADZo72xvP7LWM2rf9O2z5s2fdBVX1Ytx/skgAzR2fnD6lH4/lSs9hxedpeVT285DHBlRNg5mi68mHtSbLW2mn9O0u+8Cx2nJQ7qaoXa647hp0TYGan9/6iFmE9GPeIqKp/7vt7Wh/PkC9zA54pvsfnjoQdaL339DHAWuP2lNNPic/G42ktLjP7cQz7+iK/aGutPa+qEl/mRIC5ccYs+KAWv5I7usD4x1V1W3yZG0sQ3ChjHfciv5Kbxt+rDfE976+O4LoJMLPRWrvbWvtr/L/b3oZh09LD2Vgr3vp9VfX9lvg++s9OsEMCzJyc1OLqhrX/7zZmv/fG29V7B6+OPaiqozV/bzS5U36IQZg1YGZj/AFn1eIG7IdrPn9VixNwx9tmvyPUr2rzJWrfjO85XL4HMeyaADMbrbW3VfVk3TW64yqGR1X1bMusdho7naTbqvfePvVY4SoIMLMxTpj9UFX3p5npWEr4qRZLEg/9gILPiQAzK+PE2Y+1WAuerv/1N0J8lgQYIMRVEAAhAgwQIsAAIQIMECLAACECDBAiwAAhAgwQIsAAIQIMEHLrMoOXbhfIDXd4+J+7PQLX5PT09F3v/dvV7Ze6F4QAfz7cAwR2p7V22nu/s7rdEgRAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQIgAA4QIMECIAAOECDBAiAADhAgwQMitS45/V1V/XseBsFuttfQhwJfku3UbW+991wcCQFmCAIgRYIAQAQYIEWCAEAEGCBFggBABBggRYIAQAQYI+Rtwox83BvNKBQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(inp[0], cmap=\"coolwarm\", title=r\"$$x$$\")\n",
    "plot_im(out[4], cmap=\"gray\", title=r\"$$y_2$$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADrCAYAAABXYUzjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJ2klEQVR4nO3dQW4b1xnA8e/VDRIvAigqAiOrwlLWWcjpCaLeIHJvYN/ARk9gKDewfYJGuYGVE9TyDSwg6xY2b/CymDfWeEwydEzqG5G/HyBYJEfi0Is/n968GZZaawBw/f6SvQMAu0qAAZIIMEASAQZI8tfsHWC7lFL2xvfVWmfXvyc3m//H3WAEzLr9GhFvh1+llOPcXXpfKeWglPK2lHKQvS/zlFIexOj/MCKep+4UG2EEzCacR8Tp4PbLrB1Z4CAi9iLiKCIuc3dlrvOI+Ofg9umiDbnZBJhNuKy1nmfvxE1Va72MwRtDKeVN4u6wQaYgtlAp5biUUjf09SD79cG2MALeQrXW81LKSUTsR8RJRBxHxON4/8/t0/b444joR1j7EXEY3Z/oR+3fiIhnEfGi/cy9Te8/7AoB3lK11l/at89KKTUiZoP7opTyPCKe1FqfLfod7eDZaUQ8iIjXETHb3B7D7hHg3bE35/bSA1BtHvdem3Z42u6e2gE1uLHMAe+G2fDGYI3pbLzhPG2U/NNa9wgQYOYrpRyN7no6d0PgTxNgPtDmfi+G97WlUa9y9gi2kzlg5jmK6KYqRqe/nuTsDmwnAWaef8y7s42CgTUxBcF72gG6H6/heR6VUl63r7N5F58ZbX/Qtn2x6X37VKWUB6WUi7a/b0spT0eP77XX3L/+F3Pm3NkBAszY2aafoJRyFhF/q7Ue1loPowv+xR/82OPoTgyZ1IV9hlpYL6I7meWH9truRcT9Usrrts1BdK/1yeD1R0RcTPXiQGyOAO+2vYh3o8tHLRLHscElZ6WURxHxptb6ePTQwR+MAvvwTvkaE79GxNNa6+N+7rxN2zyJ7vU9iO4N7qTW+ioioo2O+9c22TcXNsMc8G47K6UMb1/G1anIjzb0nA8Ho77xdW/35/1A26YfHU5yCqLF9XzBmYX93PnTiHjWx7f5fvD9lN9c2AAB3m3n0Z2MsRfd6Osguj/1NxK5FqlfRnffH3y/6Cy74chwcpFqbxCnEXF3wSbDg5rj9dQn0U3BnDvIuXsEeLedDi8bOTjl+P7iH/kkD+PDpWwP27+vlnziw7uAjUaPU3E/In5esv/91MpsvP8tus4y3FECzDu11mellMvY3J/5T4ajvDZy7OO07Ey7qc//voyIn5c83k8zTHX/SeIgHO9pI+KNhGJ4NbZmONJeFrA+0lMc/UatdeHovb3J7LWbk5y/Jo8AM8/Gl6I1/fTD+ZKADed/b2LAJj1/TS4BZp7z6OYrZ/MeXMd61dH0w7Lgv/tstBv6MUf9/s8cZGNMgPlArfWy1vrVkk1el1I+9Wy5Vacf+hHkJKcfVjD1+WsSCTAfZRDe2Sf+qn5keLlk9UDE1Sj5xgVstH75v0u2cwLGjhJgVjZY7xrx6Z+M0Ydp4ch2y+Z/l43gXWt5RwnwbttbdcMW31+jC+fC+eE/8dzLPnL9Rsz/tovvzPu06H8Nvp/7hrXg5BR2hHXAu2FveKPWOhudgrxUG4n2B8ouYz0fznkZV6PgRfqoTXb+t5RyGu207VLKm36p3eiqcnPfsNo2DyPih+vYV6bHCHh3zObcXhrAUsqP7epeL6IL5r32c+v4YM4+6N/Pe7A97167OdnRb7x/6c7hKofT6OJ6GYv/0jiL7uSU2Ub2jMkzAt5CbWR1HN3FbfpTfw9Hfya/iYh/t5HwbHD/vejC/H1cheNVdPOUx9EdFPvkOct21t1JRByXUh70F7FpS9zOogtXfwBuyvO/59GN1IdXODuNiP3hmYWllNP+CnDtqm/Po7tymumHHVZqrdn7wJq1KYNNRutw0ZrWNnJ9WWt9OO/xOds/im6kuB9ddPvLNx7H1QG/r9Y5Shz8/5ysI4CD1zCL7k3raa31p8HjR9G9luGBxyerXteiXYR+Vmv1kVBbxgh4C7UDVqtP8iZqofrgYjSllOft24VnyU3FotcwePxVDA4oQs8cMJPTpiFWOUsObjQB5lqVUo7b56S9XvI5cP3Uw+WCC5zDVhBgrttZdPOkBzFnFUYb/fYrC8x5stUEmOu21/59teAgVL/C4uFEL74OayPAXLfL6FYf3Bs/MPiAyp9MPbALBJjr9ji69cfvph9KKUdt+dr96OI8/sRk2EqWoXGtaq2/lFJm0X0i815crf/9z3DtLOwCAWYT9ocrHMbreNs65Q+mILgyWiGyH+u5/gYTYwqCTfgxIt72XxO83m1/QaFJHuRrp4y/HXwdLf8JbiqnIrNW7bTb/dHdL6d+NtuUtPnx8RK9N1aFbB8BBkhiCgIgiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCDJR12Mp5TitLkt8cUXX6y87a1bt1ba7vPPP1/5d3722Wdr3e5jnnvVbW/fvr3y74RlLi4u/l9r/Xp8v6uh7ajDw8OVt/3yyy9X2u7bb79d+Xd+8803K213586dtT/33bt3V9ruu+++W/l3wjKllN/m3W8KAiCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkESAAZIIMEASAQZIIsAASQQYIIkAAyQRYIAkAgyQRIABkggwQBIBBkgiwABJBBggiQADJBFggCQCDJBEgAGSlFrr6huX8r+I+G1zuwOwlf5ea/16fOdHBRiA9TEFAZBEgAGSCDBAEgEGSCLAAEkEGCCJAAMkEWCAJAIMkOR3N7tFTOm+hgcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(out.mean(0), cmap=\"gray\", title=r\"$$\\mathbb{E}[y \\mid x]$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Probabilistic perspective\n",
    "\n",
    "#### Typical modelling approach:\n",
    "\n",
    " * Assume an underlying data distribution $p(y \\mid x)$,\n",
    " * construct a parametric model: $p_\\phi(y \\mid x) = \\int p(y \\mid \\phi) p(\\phi \\mid x)d\\phi$,\n",
    " * optimise $\\phi$ to minimise a divergence between $p$ and $p_\\phi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Simplifications:\n",
    "\n",
    " * Model only an \"output\" variable $\\eta$ (e.g. the logits before softmax);\n",
    "     - the model becomes $p_\\phi(y \\mid x) = \\int p(y \\mid \\eta) p_\\phi(\\eta \\mid x)d\\eta$\n",
    " * set $p_\\phi(\\eta \\mid x)$ to be a delta function, i.e. $\\eta = f_\\phi(x)$ is deterministic\n",
    "     - the integral reduces to the likelihood: $p_\\phi(y \\mid x) = p(y \\mid f_\\phi(x))$\n",
    " * set $p(y \\mid f_\\phi(x))$ to be a factorised categorical distribution\n",
    " * optimise KL-divergence between $p$ and $p_\\phi$ using MC approximation:\n",
    " \n",
    " Using these assumptions, we can derive an easy to implement model:\n",
    " \n",
    " \n",
    " $$\n",
    " \\begin{align}\n",
    " &\\text{argmin}_\\phi D_\\text{KL}\\left(p(y \\mid x) \\mid\\mid p(y \\mid f_\\phi(x))\\right) \\\\\n",
    " &= \\text{argmin}_\\phi \\int p(y \\mid x) \\log{\\frac{p(y \\mid x)}{p(y \\mid f_\\phi(x)))}} \\\\\n",
    " &= \\text{...steps left as an exercise to the reader...} \\\\\n",
    " &\\approx \\text{argmin}_\\phi -\\frac{1}{N}\\sum_{n,i,c}y^n_{i,c}\\log{\\text{softmax}(f_\\phi(x^n))_{i,c}},\\quad (x^n, y^n) \\sim p(x, y) \\quad \\text{for} \\, n=1, 2, \\dots, N\n",
    " \\end{align}\n",
    " $$\n",
    " \n",
    " **For a good solution $\\phi^*$ we will have $p(y_i \\mid f_{\\phi^*}(x)_i) = p(y_i \\mid x)$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![caption](files/gm_det.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "class MLP(torch.nn.Module):\n",
    "    def __init__(self, inp_size, n_hidden):\n",
    "        super(MLP, self).__init__()\n",
    "        self.layers = [torch.nn.Linear(inp_size, n_hidden[0])]\n",
    "        for cur, nxt in zip(n_hidden, n_hidden[1:]): \n",
    "            self.layers.append(torch.nn.Linear(cur, nxt))\n",
    "        self.layers.append(torch.nn.Linear(n_hidden[-1], inp_size))\n",
    "        self.model = torch.nn.Sequential(*self.layers)\n",
    "\n",
    "    def forward(self, xs, as_probs=False):\n",
    "        orig_shape = xs.shape\n",
    "        out = xs.view(xs.shape[0], -1)\n",
    "        for layer in self.layers[:-1]:\n",
    "            out = torch.nn.functional.relu(layer(out))\n",
    "        out = self.layers[-1](out)\n",
    "        if as_probs:\n",
    "            out = torch.sigmoid(out)\n",
    "        return out.view(*orig_shape)\n",
    "    \n",
    "    def log_likelihood(self, xs, ys):\n",
    "        params = self.forward(xs)\n",
    "        logprobs = torch.distributions.Bernoulli(logits=params).log_prob(ys).sum((1, 2, 3))\n",
    "        return logprobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def run_step(xs_t, ys_t, opt, net):\n",
    "    opt.zero_grad()\n",
    "    loss = -torch.mean(net.log_likelihood(xs_t, ys_t))\n",
    "    loss.backward()\n",
    "    opt.step()\n",
    "    return loss.detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fbeee2d2366c4576a4f00570314091fa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=1000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADrCAYAAABXYUzjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAR5ElEQVR4nO3dX27cRrrG4bdiW1acIOHYmABBAgwi3+RazrnMVeQdWJkVRNqBDa9goNmBMivIKDuQsoKJvANrMNfnjN0XdmLHf+pcsEqiqO7WVxTZH1v6PYAgqZtks7u63y4W+ZEhxigAwOJ94L0CAHBVEcAA4IQABgAnBDAAOCGAl1wIYcN7HXBxi2zHro/Fe61/BPASCyFsSdr0Xo9pQghrIYTnIYQ173UZO4d23EyPuaj5MMN17xVAN/lDG2O8770uM6xJqiStSzryXZXx8mjHGON2CGE/hKAY449Dz4fZCOAlFEJYl/QoxnjXe13GLIRQSXos6UG6aSJpdyzh4dyOm5IOQwi/xhifLGA+TEEAL6c9SdveKzFmKXwP07/3YoyTNIa5F0KoYox/91u7Y27tmF6PbUm/hBC+ijFOhpwP0zEGvGRCCDuSnsQYD7zXZeR+UT0Mst0I333VwyL/47li0jjaMT32r6q3EgafD2fRA14iqVf3UBJDD3OkTft16TgslP9P9he+Ug0ja8dtSU9DCLsxxpKx+q7zoYEe8HL5h6QD3vDnypv1x69TGnLYlHR3BGPAo2nHtA4HknYWMR9Oowe8JFKv6YGksR71MCbfp9+nNu9jjD87rMspI23HXZ2MjU8WMB8SesDLY0vShLHf+dJxx1X693DOpF5G147pi2miky+uQefDCQJ4eWyr1aPDVM2x3jG+XmNtxwN1OyKj63wQAbwUUq9uTc47j5bE8ab9GMZYm0bejvuS1jtULnadDyKAl0WuwR9jz2lsvkm/x/hajbkd8zqVnu+h63wQAdy7EMJGCOFp+jlMh0Q1718PIeyl+56mv9dnLS+5L5X36EIIDxvrspd2AM2bfi1NO8Ye2lQhhJ0QQsw/OhmC2GjePpLn1KkdsxDCVuN98zyEsNu6v0rtnNt83/DeUmudinYOdp0PSYyRn55+VG9ePpW0lv4/lBQlraf/d1Rvsq015tlN0zyYs9znkg4L12VP0k7j/yjp6Tnz5HWJPbwWG+c9r55e8yq97muqjy6I6WercfuapGoE74/idmw8x8P0/qka77XnuU0b7731xnz76bVYMz7O0/PeI33Ox0+kB9yzR6pr+3OvIG+ePc69lRjj/cb9ijFuq96TPK+HWqnghDYhhIeSnsUYH7XuWjunRzTmTeSpYoyTGONRPNurPMi3p5+Jx/q1VOp2YqJfVJ/D4lF+Hun5/k11m26p/sLdjOn8DOn9ltvTOjzwRHWQl+o635VHAPfr+3j6WNMq/X4gSVMCMTsO6vYdjZ0bRVVKKdjzMqrGfbenzZCmyY81hs31Lo5LjKcEsquO7ZjPlnYQpxeP5GXtSmqfHOebxt/WL9Sj1rpadZ3vyqMQoychhAc6+0bPb8hJMxCneJZ+T+up5GX817geW5LaBQfN4zR/nTFr87GXpgfcknv3YzxLV1E7SsdfijuSvpoxSfOcFrut+zZVf/GXVNzldVtT2RdF1/muPAK4P/cl/dS6LfdCzit9zb3SaT2IKv2eGNdjW2dP7p3D/8mcTfFm73GMAWaRX+9ZXzJmIYSN2G+xRJV+Twrm+V7SP+e0Wf7CmbTbLIVu6Rnf8uNUC5rvymMIoj97zeGHVkXWeZv0+YNUTbkvh/OzKfdN87dmjyf1ovLy272kpqUb/21Kz7NK//ZRAdf3FSpK21Gqv0hmDVtJwx1yN3WYasD5riwCuCdTeksbc+5rOx6q6GE95g0//HPOrGPefLdojnleqAecwnximG43DT0NIsY4c4ul9YXT15h9yZdDH/NdeQTwcPJxkXMDrbXjYojxszz8cDDnw9wc/13WHXDHR3f0MISSzx18RmuH5mTOfUO7DGP2Vx4BPBzrJn2zBzVv2uLNu9bww96cSZvlu8v6Yc5j2H18id2f8zo8TkUO+bjbo1R881y2k9L0tZme22zS4xEfXYZJLjLflcdOuAEUjv/+tfF3eyeedPKmrqbcdx7r8EP+sljW4QdpQUMoMcZHIYSflIoiVAfhM9WXPZoXhBdpx2mGGLOv0u/Jgua78ugBD8M0/tvqoT6Zsel8kd5N7iWdV4iQ12FZe7/SyTj6vy6ykEbPdqbUTvuqe36V6h2w57VTb8NLrWO2Zz7f1tBSCcaCF4QAHoZp/Fd1uWw2a293/uB2uXxN/pDOXI/LMP7bqu4r7gGncyhspWOodyXdTeeY2G2P64YQHoQQckDvqm63zTQsMa/K8CLt2NZss3nPd95RL9PkdSv9sug635XHEMQw8gfkvL3xeQfZj7N6yrG+oKTUrdSzSr/n9Uwuw/hvpyMgGoUOT2O6SnII4V4umklfTv8OIXzX2jq5H2M8SiXft2OM2+f1Nru0Y/pC0JQquOaw1dTnO6Mg5zxr6fEmC5rvyiOAe9Ya/535YUsf3jXVwwPnndD6iU6HjNXRvHVIci98mcd/76XfE2sIpN7qqfMntMUYD0IIP6o+F8Of0m3NULurNO5p/PIyt2O6avLD9Pez/LiNSxpJM55vmmZb0neWx2r4Rt3eB13nu/IYguhfsye0Nq0+PvWWdlQH5L32/VMcqNvOm3zkw9QPfQjhsLHcZe39SoUVcKlN8glunjRuX9fZIo59SdWMHu6u5u/cbCtpx+bRMc1N+x3V4Xo0Z1l7qgtyJgXrprS8Lu+DrvNdeQRw/453fKmuptpvhnDq+e5L+jnGeNf4IdlP8xbtVEmbrgeqr1hwPN6czvt7qNMf7KUc/03y2Kv1OeypPltcu1R3Q2eDJA/fnDnf7bxCiRlK2jGvR/MMZzuqhzx+VBq+Srcp/b2e2nUvFl6AtDF+XfQ+6DofagRw/46PKEgfnE1Ju/lE2qqPV70XYzSXujY2b4tPeh1jvK90msx0Eu98XtkfdHoP+oXPn+ChtYVx7mZwCr91Td9BdXfK0Qy5d20+ic4sJe2YhqUeSdppvHf+m983aVn3VH+5Pk33P5b0w4wzp53nr611HHo+iDHgXk07pWMK4T6uFvCzOl72JfX0zpyYJYTwj/TnzCq5JdCsgLOEQP7iswZGHiLqK2DM7Tir3Rr39/XeUlqn0p12F5kPogfcN+vhQV3squ7tVH0sLPUcLVVyY5cDyBoCuWrrVE839Yz3W7dVqndSzjpGu4te27EPjePRiw5b6zofThDA/RqiPFTSce9uoiknbW/LpbFp07SaMVkeOzzquMk6FnmIYFoV4TR52KVdEnxfZ3u5+TX6ocN6TVXSjgu0pfo9W9rL7zofEgK4X0Of0vGRThdvzLKn09dKOyX1fvNe9r5PuziIvONwysU111WHgLUHnKdrDwOsNYdh0k7LLdXH/Pa9NWNtx0V5rG5fMl3nQ0IA98RaHnoRqaf6rHlEwwxV+j1r0zlvMm4v0cnX91SH7UY+iqBxKkhzDz5tmWyq3rm1npZTKR3tkKricoXbvSF6dwXtOLi0Ds8KvsAuNB9OC7G+qikuqFHGKtUf3EGCLYXPXozxT3Omear64qBnPhwpXLYk/X3ONer6Ws991YdRXfhDms42VqnefP8qVZbtS/pm3msxZ3lrOjkWOxfPHKkemvhp6GCxtOMipNd1s/SLput8aIkjuDTzMv+o/vBuqD6BS74k+kPVvbVqoMfcVV1EMOv+B6oLCtYat+Uig+ca+FLx6fF6vSx9es6H+TXVySXo13ta9iBtdZF2XMDj70jaWdR8/Jz94TC0C5hz5qzjHVzq5+Qrp8T63AOHs65bFmP8OYQw0cml7m+ndfkpni0+WArpOe9IOkznVDjS9ON2u7gdHQ7DO68dh5R64BsxRksl5oXnw3QMQSypFKy/qN4MHN1ZqPoeghhSCGE3nn8+jqEeu9KC2zF1HPYkfVfyxdN1PszGTrgllT4A3+mktz02R6rHa0e9k2/G+R8Wxqkdd9QtRLvOhxnoAeNKS+fm+HmMWxG4/OgB46q7Q/jCCz1gAHBCDxgAnBDAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgJOis6GFEKjacJDO/mXywQe279QbN26Yl7m6utr7Mj/++GPTdCsrK6bprl+3v5Wtr+e1a9fMy7ROW9KWuDwODw//L8b45/btnI7SkTUsS8Ll1q1bpum++OIL8zK//vpr03SfffaZeZnffvutabovv/zSNN2dO3fMj33z5k3TdJ9++ql5mVVVmaaztrlkD2tCffxCCP+ZdjtDEADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATjgNeAiXHeb579840XUnRhHWZ1mNhJfvxsL///rtpurdv35of21pYUuL9+/em6UqOA7ZerYbjgJcXPWAAcEIAA4ATAhgAnBDAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBNKkXs2RFmotcxVKrt8kdWHH35oms5aNixJv/32m2m6zz//3LxMK2s5cEnZcMn146woMb786AEDgBMCGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEyrhema9kKJkr3QqWaZVyUUsrdOWLNN6UVBrNVrJRUatyyypQBwClXCXHz1gAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATSpEdWUuMSy60aS2zLbmIpLUkd4hy4CGW9+bNm96XaX2NSl536/uDkuXlRQ8YAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATAhgAnFAJtwTevXvX+zJLqqdWV1dN062srJiXaa3u++OPP0zTvX792vzY1vUsucjoEBdYpcLt8qMHDABOCGAAcEIAA4ATAhgAnBDAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQimyoyFKTa2lriVltlbWsmFJevXqlWk664UxSy6gaX2N+r5wqER5MU6jBwwATghgAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOqITrmXelk/XxS6q8hqhGs66nZ7VgyWOXXGyz72V6v+fQHT1gAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATSpEvmffv35umu37d3vQvXrwwTXfnzh3zMldXV03TraysmKa7du2a+bGtJb4ly+z7saVhLgqKcaGFAcAJAQwATghgAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAmVcD0b4uKMJRVR1uota8WcJL1798403evXr83LfPnypWk6axVeVVXmxy6pAvTERTkvP3rAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBMCGACcEMAA4IQABgAny1GTeUlZS0iHKBt+9eqVeZm3bt0yTffJJ5+Yl2ktmb5586ZpOuvzluyvZ8kyrUrKhikxvvzoAQOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATKuGWQMmFPq3VUysrK70/fknFnnWZ1guSlly41Dqtd9UaF+W8/OgBA4ATAhgAnBDAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBMCGACcUIq8BK5ftzfT27dvTdOVXHDSWupaUjJtLQd++fKlabqqqsyPbX0+JSW+1jLskmVaL1yK5UUPGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJxQCbcESqrWrBfbLKnIunHjRu/L/Oijj0zTWavBrOso9X9BUKAr3mEA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghAAGACeUIjsquYhl3968eWOe9sWLF6bprBcElewX27x9+7ZpupLnc/PmTdN01gttSsOULVvfHyUl4BgXesAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghEo4I2u10RDVbSWVTq9fvzZNV3KhT+tzWl1d7X2ZJdVofRuiwsyz+hHjQw8YAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATAhgAnBDAAOAklJRGhhD+V9J/hlsdALiU/hJj/HP7xqIABgD0hyEIAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcPL/ZqblINCBR5gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "producer = produce_data(10, DIM)\n",
    "net = MLP(inp_size=DIM*DIM, n_hidden=[64, 64, 64])\n",
    "opt = torch.optim.Adam(net.parameters(), lr=3e-4)\n",
    "\n",
    "N_ITERS = 1000\n",
    "pbar = tqdm.tqdm(range(N_ITERS))\n",
    "\n",
    "for i in pbar:\n",
    "    xs, ys = next(producer)\n",
    "    xs_t = torch.as_tensor(xs)\n",
    "    ys_t = torch.as_tensor(ys)\n",
    "    loss = run_step(xs_t, ys_t, opt, net)\n",
    "    pbar.set_description_str(f\"loss: {loss:4.2f} \")\n",
    "\n",
    "ys_pred = net(torch.as_tensor(inp[:, None, ...]), as_probs=True).detach().numpy()\n",
    "plot_im(ys_pred[0], cmap=\"gray\", title=r\"$$p(y \\mid f_{\\phi^*}(x))$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* What we see is a distribution, not a sample from it.  \n",
    "* We can obtain a useful sample by applying argmax/thresholding operation (mode extraction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAADrCAYAAADe1RnPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUSElEQVR4nO3d73HjRprH8aevNgDcbN3rq6Mz0MxGYE0Gkh3BShloaiLYojOgHYGXk4HGEaw5GUiue79nMoO+F3iabEJAAw3i3yN9P1UszZAA0QTQ+OFPo+G89wIAgCX/MXcBAADIRXgBAMwhvAAA5hBeAABzCC+gJ+fc9dxlEFlOOXCZKZdj32lNPV4K4QX04Jy7E5FbytHMObdyzu2dc6u5y7J0MyzHW53m0sdr9Jchvwx4C8KGxnv/kXIkrUSkEJErEXmetyjLNcdy9N7fO+cenXPivf95qeOlEF5ABufclYh88t5/RzmWzzlXiMhnEbnRtw4ishlqA3qpmZfjrYjsnHO/e++/LXi8WoQXkGcrIvdzF0KWU47F0uDa6X/fe+8Peu1l65wrvPc/zVe6o9mWo86PexH5zTn3P977wxLHa8I1L6Aj59xaRL55779SDhN+k/LU5X0UXI9Snsr825wFE1nGctRp/y7l0elix6vj6B4KaKd78XsR+c57P9v1m6WUo4soLG69918mnvaV6FGX997pew8istZB7uc8dbik5agNap5yyzL1eFUceQHd/CIiX+fe0CyoHEsXTsUd55OeJryVcqM59zWvxSxHLcNXOQX7Iser4poX0EL3km9EZO7WhYsohxE/6N+zU3JTHwHWWehy3MjpWuBhweMdceQFtLsTkcMCrjEtpRyLpqelCv3vLjHoXBa3HDXUD3IK/UWOFyO8gHb3UtmDn8lSyrF0V9G/lzi/lrocv0q/lo9TjycihBeQpHvxKykbHrz5chhxPB23hGtKsYUvx0cRuerRI8rU44kI4QW0CX2yzb2nvJRyWPBB/y5xXi15OYYy5fZDOPV4IpIRXs65B+fck762etExNfxKh13iHsbgnHM3zrmNvh71te4wnx50fu50fu2iz270s/B540KuDLvV6d9En99p2cJ33UTvx9P30TiF/oZ1NMxWmyHX/Y749+9y+zLrMw91Or7h9WI8LWfdsE0tnz6K5O3BO+euo7qyq84v59xVZZ7XztNLy6HTevX1VteR47KU02nD68oyXsJv6rUcA62vYb3ZO+c2lc8LXc5hmT92WLekUqashiRTjxd/QetLyrvA19H/vYg8tYyz0eF8l2lYfulvva68V+h88yJykxj3Tsffx/NVRB5E5E7/fZ2alzqdJxFZVaa/09c2LkN4X/99U5m+j6a5rZnWToe7it57jP8fld/XfceQ81BOfefdxb+hZfgnHWav5VwlyrUP86rj7wj3sKzq5peUzYMfK8tqkypzn3JE68Vs9TZabxt/1xAvXabhdNxNtA7cRe+vRKQYsxwdy5q9HKPfuNP1p4jWtb2cthlh3avWTZ9axyvTeWpbR5Ywnve+Pby0cm8q74WV46qlUF5EHudeYcZ86cqU2viFjdd1y/eEiv6kG+OH6LO7aJ4XlfE2TctCV2bfpbJEw3r996ZluH00/aJh2BdBN/I8vIp+Q+OwWqEfu2zMJCOAo/kR7yisw3foZ+uG8ULwNs3L3HLMXm9lovCqTDMOr04b7ClfucsxGm8nujNbs5xDUO/kPLg20bx4MW7DdLbSY8dl6vG87xZeT5X/F20biMowD30KVvmuUV4Droxne9KVz8MGNRkg0Tx7qttw6Ep60/Dd+8T3hg1/6wYkWmbb1PyR0wbueHTYMFyoWLVBOPQ8rEyzdm9OyvDttOcrp6CuDZyGcfaV/8cbkMb5IKcjzBfT6lmOWeutft8c4RV2FvxU08woW/Zy1PHumsaR87Cu7qzsos+6Hnmtc4afazzvffqal16zqN7UF7fL/71h1PjaTK8Lk3ru1ku5RzrKq3q+uEcZj+fW5dRr9Rl/6j35qu16g1pJzeMjvPc/+Zc3WP6of1PzOCyjHxPDVP3p0zcOhvJ99umeCsJwja2Jhp6HvuxF4auIrJxz25pBNtL9uUmh3P/XZWC9jlhdFuE7Dt77VLPgP/Vv3XXN3HLMVm8XIKxPF/daPoKs5ShyvKF5LSL/aBgk7qOxuj27FZFPUnZK3PUaWyhbbgvAqcdr7WHjXl5W9FABvyU2cMcZ6nt2fe/LjjTfi8i7PuN31FSJu3qW0wY6VdmfpVw476S8Ma9N3Ua3TqiofyaGedK/OStH24Xtg/5tm39huNQyHGMe3orIHyJy45y7CaGvOyubjIpc6N+26QUfReTXynuh5Vtbd0RhHtUtp9xyzFZvFyDM70vrtjjnrv2wNxIX+veQMc4PIvLPxDIL24BDdZnpep7bc36YTrHw8VrD6x9xRde9gDCzUkctgzQHXXoF0hUq5zk8Rcfh5r43pev0Ly7nGPNQd3xupQzhrXPuO9FQqDl6TQmBkto5iG3jjZ077+mhbYcg1Kui5rPccsxab+eiv7PQ/w7Rs8atDDsvcpejSBnC/0x8PtZtAX0PGiYbLxleNRU9PvWQmqFLPnQflZ4Gu5bzDfK7yt+kjCODbzqt1FFVKEfOyt11+k/tg+QbaB5+dc79JOU1sK2+937IctZNs/LWdeKzquPpxQHK8Vbr7Yfo3xcdeWkQHjoMt5HyGvUofSamduArYT3UbQA5wTrneNk3KYdTD1+bDmMr9yIt4b6K0UX3Q+2l7C36Wcqno97rdY7eC6hF2Iv+kBgmLI/q6axFGWMeeu8/SbkhvpJ5jibC/SvJMHDnPQyMcdT9Vurt8X6mAc7ahMe5vFC57npIfDa213CNsrecm5QLOa0cqWsycdcsr36G6lHCH1Lu4f/de//ee/+lpcHDIKJz2kV8Q3KlbCsR+XnJp2BHnodhD/yh682aNfqeCul6Gi5edqlhs8vxxuptuGY3xA7Ax8R8+Kw3AIf7qp71xvS9dOtodqjr+GGZHTLO1rTpc2pzjvGyHonS9dRDqLAXbyzH3ou5dOOoK284t55a2cf0Lynn9S/OuUMog26ot1IG12IfFz/mPNSjiXdSXrvYSvlk3f/M+IpQoYoe08653hW3BK07Qu5dDpmh3s5oktOe3vtPzrlfRW8YljJE/pT2Vn2XLMc6Y1yjLPTvYeHjZYVXSPnnlo1+WIEumqHu9KTR0TjnLt2whz3Zn3M3us651UB7S5+99+81rH50zsUP4btd8hGXGmUehibG4TqXc+6LlK0Pt977rk3lL1k+na53VY6MvjUsr0vKMWm9nVk4/fqvS74kOqJq5L3/pt1NhfrWpRXrYKeEdb1p/b0XtJhc/LWvnPAKMyp1AXGw8+ZLbyqvYdGlBZdI/W+4d879ekm4aCV7Fjme4196UJ0ZeR5uJWou7r2/dc49SRlgdy33pwVhY5PTGjLodL1LyhtQg08jlGPSejuXyinh7HqgYRCOUm+lPBUYjqw+xcGvp+jX+tpIuR2518YbqR3GS5ZjVbzMUr93kzm9MGxu0E49XlZ4Ffo3lZCDnjdf+FHD8SJ7RougoRVy/uwia0aZh8650DVStULcSnmKcuOca30Eu+5AnZUzQ9i4tO0ghT33xiPPC8tR6N/J6u1MerU0jG4CftIb3MU59z6ckdFg/8M5931lHf3ovX/Wde2d9/7eJTrOFum3HPWGc6nZ2YpPNdf+3oab1duEW0oOCx8vq7Vhl2QMe5FLDp3BtVyb+0FO53OrK+1BLqCVaaUVyLSh5qFuQP5W13RZ51e4abPrjeDfJN2a84XK9a5U7yIP+vlzh9PX2eVQb6XehtsgDl03hHq0tpNyR6f2Zl4N85+lvF4a3vsS7fjE9xB+7bDD3Xk56pHfRsqdrfgJEYWcGvnU/l4d5l6ae+Zo8kH6rQdTj5cVXqGy1854Vz7Ko9D/Wt17yxH/xtpHf+gK9FFOK1ARfTzUNa+fRWStj0cIj+AIr0d9PMJdZku7oU7VFpW/VYPOQw2uR0mcgtTm889SdjXVpXuwr5J/5Bzvga9czcP2tKxrLUuX+8/6lEPk7dTbrJ41dJn8JmVwfYveD4EWe5SyRW/dkdVG0g1hqnKWY9wKNd5WrKUMpufEd22lvFn9kFE20e/rsx5MPV57x7z+vBPF0L3+XfReaC0WOhZt7Pjztb2kPGW3l5rOR/WzrZw6Ag6d7hY6z9aV4eOe4zt3ZCrnj/hoez01LRs57+CzsVNWOT2GIXSmWyTKFXdKW9uT+RDzUIeLp5UqfyHnHZauJfG4DDl1Ltt5nY7qQnhCQPVxNVmPi+lbjmjcWeqtTNgxb5dlXxl+JzWdN0vNI3Lk1DF0Voe6ly7HaJ2uPqVgW/mudaWstT3Qd5jeVdeyzTnecfweP/BBK+M+WvmvogrpmzYEr/GlG8MHOX921rZaiaKN2F6i3p8rG9K6V+rRHjdyeibVi3ku58+6qt1IVd6vvnaVaaXKGW+cU8O92LhcMg+bppOYZ4+JsjX1ap+14ZLTzsQmKvej/ranUGd6rGu9N6Bz1FuZ7nleqy71paZcdevii57/5bRjeXFP+7nLMVpuYd2pqxOP+lnvdUu/a52qO0sZ7zj+gCtQ2Ai/6ud3TflKbUyiMMk5Igh7Or0e/ta1nENvBOd+6Qah62NUimhDOuhGO6ccGd85Wr2dMLyOO1Ydh089A68uvBqHX8pyHHBd6PusscnGC6/c7qFq6fnjLnfxI4NPn6/+RcrufjqfL/anBgurC3qbqPveQ+r/r8BGuj/SpmsT5rHL0eoV1dvQWrJry7pwTffsmnN0zTR+r5DyyKvpHrw+Bl2OQ4juN8x6TNTU48Vawyt0e6KNAYqGwdb699l3u38GF4has/W5JyeM06fZ9ZukOwgHEfncYfAxuuzJLscbq7ehsUbX/jvDTb3Vhkkf5WXjgTCP/t6jXLUy16ep3Em5zuY2nph6vKMuR17hgvlKajZ4uiENrWK69lyAy4R7di65afUwTFHejE/S0CKyYuzHinQtx6urt865VWhFW/noSsoNYdcjrzBctfXgyp/fjHwn5bz+OOBRV9B1OU7ls/QL6KnHO+lwXvLFxfvK5y9aMvEa/6Xzfd9jvF2f8XgdG2I0rudyfr1rkIv7fcqhw8xeb2Xga15y3rjpWt8L17uyGrLIqbHTVbTsNvG/dT4Pcp2r73Kc4iVliGZfB596vBff03EG1658crqQeXETUl7ZK0Chla/zhV85tSwcrUK+5pdujPeJz+PbHcbc6CXLocPMXm9HCK+9ft9etFFQ3504HXeldWKtf8M9gtuhynzpcpzipfOzzy0Yk4734ns6TOhG93jiptDhXoL9FAuZV+OyKbSi7aWhubwOd6cbs0dpaA7Oq/M830ilRZpuBK/l/H67B60ntctkjHJUPp+93o4QXhstfxH9xkF2FPS7R1lWlyzHCaa/7rMTM/V4da/Wvg2991+ccwcpH6deSHmR81lEfvUNXapgGr48P38btYjaat9pVd+kfFzDYbLCvVK+7MNuF3rrTvRAfmwMIcN0xJosR83nr67e6m9ei8hO1/NnEfnOD9Mw5t0c9aNtOY5JW1de+8wnjE89XuP3aRoC6EjD4DcpexAf48nHpsrRJGp6fuu7N6aYhXNu42d67t0cy1F3urYi8n1OaE89Xsog93kBb4lWvu/ldHT1psuR8Cxlq9ZFd/jb0J/hZGZajmvpFyRTj9eIIy8Ab5r27v9liUevaMaRF4C37q8Elz0ceQEAzOHICwBgDuEFADCH8AIAmEN4AQDMIbwAAOYQXgAAcwgvAIA5hBcAwJzWXuWXzjnHXdYYxPv3g3R2DZiw2+3+7b3/r7nL0Zf5HjYILwzFel0Acjjndt77D3OXoy9OGwIAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOYQXAMCcv8xdgAH8W0T+d+5CwD7n3NxFAKb033MX4BLOez93GQAAyMJpQwCAOYQXAMAcwgsAYA7hBQAwh/ACAJhDeAEAzCG8AADmEF4AAHMILwCAOf8Pk56fUeoJdgsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(ys_pred[0] > 0.5, cmap=\"gray\", title=r\"$$y = \\mathrm{argmax}(p(y \\mid f_{\\phi^*}(x)))$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Alternatively we can also sample from $p(y \\mid f_{\\phi^*}(x))$ at the cost of label inconsistencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADrCAYAAABXYUzjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQm0lEQVR4nO3dYVLjRhrG8ae39gDObO3nrZgbwOQEw9xgSE4QuMFQOcEWuYEnJ5gwN4A5wcbcALb2ezL4Br0f+m0jhCxLQvbbgv+vygXYli27xaNWq7sVYowCAOzf37xXAABeKwIYAJwQwADghAAGACcEMEYRQjj2Xgc8zz7LcOh7vbTtjADGs4UQTiWdeK9HkxDCPIRwH0KYe69LyRzK8MTec1/LFenv3iuAacv/uDHG997rssFc0kzSoaQ731Upk0cZxhjPQghXIQTFGD/terlSEcAYLIRwKOk8xnjgvS4lCyHMJP0i6YPdtZK0KCFAnMvwRNIyhPBHjPFmD8sVhwDGc1xKOvNeiZJZ+C7tz6MY48raMS9DCLMY469+ayfJsQztuziT9DWE8H2McbXL5UpEGzAGCSFcSLqJMV57r0vhvio1g5xVwvdKqVnkB88VK6EM7b3/UDpC2PlypaEGjN6sVvdREk0PLezw/lBaB4by3+Zq7ytlCivDM0m3IYRFjLFPO/3Q5YpBDRhD/Cbpeqob/R7lQ/v192RNDieSDpzbgIspQ1uHa0kX+1iuJNSA0YvVnD5IKrXXQ0l+tJ+PDvFjjF8c1mWt0DJc6KFdfLWH5YpADRh9nUpa0fbbzvodz+zPZctTPRRXhrZTWulhp7XT5UpBAKOvM9VqdGhUbest7fsqtQyvNaxHxtDl3BHA6MxqdXM5njyakPXhfQntrFnhZXgl6XDAqMWhy7kjgNFHHodfYu2pNG/tZ2nfVcllmNep73wPQ5dzN8kADiF8DCHc2u3STiq0PX9uzy1xr78TIYTjyne0tC5R1ccP7btbVr7Hw02vZ95L/Wt0r6G8QggXIYSYb3pogjiu3l/AZxpUhlkI4bSyzdyHEBa1x2dWxrm8rzpsV6qtU6+Tg0OXK0KMcVI3pZE7F5W/o6TbLcss7HnRe/339B3NJd1KmtvfS/v8h/b3hdJh27zhO/rQ8rr3kpZTKi+lWlHr5xrpO5/Z9z5X6mEQ7XZauX8uaea8bfQuw8rnW9q2M6tsZ/e5PCvb3WFluSv7HuYd3+d22/Yx5nLet0nVgEMIHyV9izGe1x6ab9nLlnzYtQvnSuP7c80gf+5fco0lxvi+8rhijGdKZ5Pbaqgz9ZjQ5jWVV4xxFWO8i09rltf5frutPNavYqZhkxJ9VZq/4jx/Bvus/1Yqz1Olne1JtPkZbFvLZdm1eeBGKcj7Grqcq0kFsNJwzvXZzlpQvGlawJ6TC2b0wz87lP9gh/yzrQt0e83nbkg/xsd9TWf284MkNQRitg7qlnXqNVKptPLak/UQ44ZAdjOwDPNsadexeeBIfq2FpPrkOG8rv3fdmd7V1rWrocu5mkwA20ZQ78Be7fv3x4ZFq3ve0WpUuY1V6bDsUiks7vu0eW16XaWNeejyH/T0c+aNclUNxAbf7GdTbSW/xl8d16Oo8tqzXP6lzdTVqwyl9Q7xQqmm26Q6n0V9uz1ROho76rEjyuvWN0iHLudqMgGs1M+vXsA5TG5aDu2qtZFR/iEsXK5UCbbKw8dKU+V9HPjy72uvN2T5z7X7ck1k29DXXCtt2ohn9nPVcT2KKS8H+fvetJPpbOQrQMzs56rHMj9K+r2lvPLOZlUvL2ty+bVnOeb3mfVY5jnLuZpSAP+7uhe1PXMu/LYa46jtiVa7ze93Lum7GON3McagNLFJnl7wYuAZ72M9b+Lwy2rzQ21E1rb1yd/nrOGxHM7fGh5rUkR57Zt9zpn9OcYIuDGvUtG3DKW0E9nUZCXtrrtdYxPVDpdzMZkAjk/Hz1cPZ39vWXTsw8E88ceB7d1X+QHb458rBfGdUhek267tUva8Q0n/Gbpy8enw0uOWx+qaavRD16OU8tq3arvns2rAFuarDs9bWNPT6GKMG49Wajubsdrr++wcxljO1WQCuEE+nL1u2UCqh2/P3kCs9nusdKZ3Yy3VgvhAqQ10rjRlXpd/kEulQ7kxJ2vJfSNbA622k9jFiaO9l5eTdfv/CE0oee7gJ2onNFctj+3SS2ivdzXJAK4dzl62PLU6HHSMDeRYqf2yU0DGGPNJCCl171o0/XPkzutKn6ntcG+Irof01R1E23N7H+I5lpeH3IY9xk7sfcv38Evl6OpW0p2dGL7X9olpxjpMz+W1GrG3x5Bmkucs52qq01F2PZzN4TPW4ewPenqCq1WM8dcQwo1S8JxKOg0hfFLlMjVKn2emFO6jzRHbs/33p8rvTZ8xb9izhse28SovD3tpQokxnocQPssGRiiF4Te19zh4Thk22UV7/cx+rva0nKupBnDe827r2J7/GcbaQObqGcBSqs2FEL5X6sx+qBTEddca/7Lgndp/azXUmw2Hzs+p4XiVl4fclDO4HV9a7zxv254TY7yxE725eWfblSFGa1qq9dfe+FlDCMcDj2ZeRVvwJJsg9FDwG2sZO2pP/KyB4WAjpY6UwuiL0rrnGf1PbGTaaqT1zDq1/+rxDmFTE0j+5x1yCRuv8tqrWv/v3jVga4o6tW6OC0kHNsfEk6YrG/yTA3qhVG4n1iyxqR/6c8qwrlpebZ+1b5/2vG59dxZDl3M11RrwzH627e1Gb0+MI1zB1tZlXzW8/E+y7Wx8rkF92vRdxXRBSWlYR/eZ/dxreTkY1AOiMtjhNm9jIYSjPGjGdk7/DSG8qx2dvI8x3lmf8zcxxrO2fsNDytB2BmpoGqs2WTV+1g2DcbaZ2/ut9rScq6kG8J22b0S5Vjfl9sTBau2/G78r++edKzUPbJvU+kaPQ6ar11JeR/Zz1TUIrLb6aA6FOmvC+qTUhPWd3VcNtgNZ22eHnVfnMgzpqskf7fdv+T3DwyWNpA2f1Z5zJuldl/eqeKth28DQ5VxNtQkin0lv3JBCCEs9hM9Ua1PPVa0JzZv6Iltt6UIpII/qjze41rATOK+lvHqNgLMyyZPc3FTuP9TTQRxXkmYbargLtZ/crOpThtWeMdVD+wulcL1rea1LpcE4q47vlc00bBsYupyrSQawHQ5dK82Cv26/DGke2aUebyyTbE8cwfrEl9LJvatqCFvN90rSlxjjQcd/lCtbttfw2FdUXrnttetnuFSaLa7etHWsp2GSm2+ezHnbNliiQZ8yzOtQneHsQqm545Os6cruk/1+aGV6GXv2Z6+0XffaBoYuV4JJBrCUplOUTbsY0sTQea7Sn/X4rOyzx+NP1LpHgf3znEhaBJtMW6lL3ZH1Ve6kcnjbe+Lrl15etSOMrYfCFoDVYe1VBw29GXLtuvNEOk36lKE1SZ0rDavP281feZux1zpS2rHe2uO/SPp5YHfKn2rruOvl3E21DVjS+qTYkxNjIYTf7NeNo65esqYpHS2Ex7hiwBcNvPTLCy+v6gi4LkGQd3xdQyM3EY0RMp3LcFOZVR4fa7uSrdOQUaBDl3M32RrwJpX5FKT2UVcvWdcuQkMslGo8szFe7AWVVw6hrkGQR249qulazfiqdt9M6STlpj7afY1ahmOo9EXv1W1t6HKlmFQA56GWdrgz2/C03B51N+aosonZxRBRSeva3UoNk7bXvbLyyk0EXQfq5GaX+rDg93pay83f0c8D1uuJPmW4R6dK22vfGv7Q5coQC7guUteb0vWn8rW2Dhsen7c9/lpuSiOootKJkF28/qmk+ymUl0a+Jpyt81LSVe3+2OU7afjsp7X7L2t/n9rzjj3KcF8321Z6l9HQ5Uq5TaoGrIcuL5sOxfJhyNmGx1+8rkNEnyOmmuq3ao+GDWb28yWVV5406Tj3JKjMdNe5Bh/TkcmJ0gmuQ3udmay3g42KyyPcjuLINbweZbhztg7rfsa7Xq4o3nuAnnu7W23Y2+nhSroX+1yn0m56qDHt9ChAqWZ5X3p5afwacK7V3+vh6sBX276Llteby64crYdLWy3s911fyXlrGe7jZt9l7xr+0OVKurmvQM8v/IPS4V/1cuq50/qkD0VG+G7m9g+Vmx+i0iimwxwUO3jPhdIggmLLawcBvLD1n1U+4yg7O3vtnZTV0DLcw/tfDNkJD12utNukuqHFGL+EEFZ6uHT6G6WzyJ/jCPM0TFXLzFnrE1waZwKWR2Kae2C5acarl1he9pkvlK77J9l3G8c52fkm7rkb3rYy3CVrwjmOaZKqnS9XomB7E2AQC9av2nKVEC+Vbl0nsfC2whDCIm6fj2MX7zvTnsvQKg2Xkt712ekMXa5UUzsJh8LYP8E7PdS2S3On1OWq6JN8G+Z/2AunMrzQsBAdulyRqAEDBbC5Ob6UeBSB3aEGDJThH4Tv60MNGACcUAMGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATAhgAnBDAAOCk12xoIQRGbWAUR0eTn8gK6Gy5XP4ZY/xn/f5eI+EIYIyFEZh4TUIIyxjj2/r9NEEAgBMCGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATiZ1WXqMp2f/7x2uyTj6rCN9kFEKasAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghAAGACcMRX6lpjC8uA+GF2OKqAEDgBMCGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAE0bCOeo6est71NpU1hOYGmrAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBMCGACcEMAA4IQABgAnDEV2NJWhu13XcyoXxnxpnwfTRQ0YAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATAhgAnDASDqPpM7LPc5QZI9xQCmrAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBMCGACcEMAA4IQABgAnDEXGaPoM8eXCmAA1YABwQwADgBMCGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHDCSLgXpuvIsV1cQHMqF+UESkENGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghKHIL0yf4cCerwmAGjAAuCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghAAGACcEMAA4YSTcBPS5gGXXUWveF8WcwnpykVHsGjVgAHBCAAOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4AThiJPwFQutDmV4bhTGAaN14EaMAA4IYABwAkBDABOCGAAcEIAA4ATAhgAnBDAAOCEAAYAJwQwADhhJBxceI4yY4QbSkENGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghKHIL0zXYba7uCgngH6oAQOAEwIYAJwQwADghAAGACcEMAA4IYABwAkBDABOCGAAcEIAA4ATRsK9MJ4j3Pq8NxfGBKgBA4AbAhgAnBDAAOCEAAYAJwQwADghgAHACQEMAE4IYABwQgADgBMCGACcMBQZRes6vJmhzZgiasAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJwQwADghJFwr5T3yDFGuAHUgAHADQEMAE4IYABwQgADgBMCGACcEMAA4IQABgAnBDAAOCGAAcAJAQwATvoORf5T0v92sSLYr65Dgb1NZT2BLf7VdGdgrD0A+KAJAgCcEMAA4IQABgAnBDAAOCGAAcAJAQwATghgAHBCAAOAEwIYAJz8H5tQUZlIDpYGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fd_sample = sample_bernoulli_im(ys_pred)\n",
    "plot_im(fd_sample, cmap=\"gray\", title=r\"$$y \\sim p(y \\mid f_{\\phi^*}(x))$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Modelling dependence in $y$\n",
    "\n",
    "Before:\n",
    "![caption](files/gm_det.png) \n",
    "\n",
    "Now:\n",
    "![caption](files/gm_sto.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Complications:\n",
    "\n",
    " Discard two simplifying assumptions of the older model:\n",
    " 1. Logits are independent\n",
    " 2. $p(\\eta \\mid x)$ is a delta function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Both can be tackled by letting $p(\\eta \\mid x) = \\mathcal{N}\\left(\\eta \\mid \\mu(x), \\Sigma(x)\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "However, $\\Sigma(x)$ will become intractable if one has to model all $\\frac{1}{2}(H\\times W\\times C)^2$ correlations!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Introduce low-rank parameterisation for $\\Sigma = PP^T + D$ where $P\\in \\mathbb{R}^{(H\\times W\\times C)\\times R}$ and $R \\ll H\\times W\\times C$.\n",
    "\n",
    "$R$ is a hyperparameter controlling the rank of the covariance factor matrix $P$ and $D$ is a diagonal matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This modifies the objective slightly:\n",
    "\n",
    "* the model becomes $p_\\phi(y \\mid x) = \\int p(y \\mid \\eta) \\mathcal{N}\\left(\\eta \\mid \\mu(x),\\Sigma(x)\\right)d\\eta \\approx \\frac{1}{M}\\sum_i p(y \\mid \\eta_i),\\quad \\eta_i \\sim \\mathcal{N}\\left(\\mu(x),\\Sigma(x)\\right)$\n",
    "$$\n",
    "\\begin{align}\n",
    "&\\text{argmin}_\\phi -\\frac{1}{N}\\sum_{n,i,c}y^n_{i,c}\\log{\\text{softmax}(\\underset{\\eta^n}{\\underbrace{f_\\phi(x^n)}})_{i,c}},\\quad (x^n, y^n) \\sim p(x, y)\\\\\n",
    "&\\rightarrow \\text{argmin}_\\phi -\\frac{1}{N}\\sum_n \\log{\\left(\\frac{1}{M}\\sum_j\\exp{\\left(\\sum_{i,c}y^n_{i,c}\\log{\\text{softmax}\\left(\\eta^n_j\\right)_{i,c}}\\right)}\\right)},\\quad \\eta^n\\sim \\mathcal{N}\\left(\\mu(x^n), \\Sigma(x^n)\\right);\\quad (x^n, y^n) \\sim p(x, y)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Lastly, we need to learn $\\mu(x)$ and $\\Sigma(x)$. Apply reparametrisation trick:\n",
    "\n",
    "$$\n",
    "\\eta \\sim \\mathcal{N}\\left(\\mu(x), \\Sigma(x)\\right) \\iff \\mu(x) + \\sqrt{\\Sigma(x)}\\epsilon,\\quad \\epsilon \\sim \\mathcal{N}\\left(0, I\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "class StochasticMLP(torch.nn.Module):\n",
    "    def __init__(self, inp_size, n_hidden, cov_factor=1):\n",
    "        super(StochasticMLP, self).__init__()\n",
    "        self.cov_factor = cov_factor\n",
    "        self.layers = [torch.nn.Linear(inp_size, n_hidden[0])]\n",
    "        for cur, nxt in zip(n_hidden, n_hidden[1:]): \n",
    "            self.layers.append(torch.nn.Linear(cur, nxt))\n",
    "        self.mu = torch.nn.Linear(n_hidden[-1], inp_size)\n",
    "        self.P = torch.nn.Linear(n_hidden[-1], inp_size*cov_factor)\n",
    "        self.D = torch.nn.Linear(n_hidden[-1], inp_size)\n",
    "        self.layers.extend([self.mu, self.P, self.D])\n",
    "        self.model = torch.nn.Sequential(*self.layers)\n",
    "\n",
    "    def forward(self, xs, as_probs=False, n_samples=1, use_mode=False):\n",
    "        out = xs.view(xs.shape[0], -1)\n",
    "        for layer in self.layers[:-3]:\n",
    "            out = torch.nn.functional.relu(layer(out))\n",
    "        if not use_mode:\n",
    "            cov_factor = self.P(out).view(xs.shape[0], -1, self.cov_factor)\n",
    "            cov_diag = torch.exp(self.D(out)) + 1e-3\n",
    "            out = torch.distributions.LowRankMultivariateNormal(\n",
    "                self.mu(out), cov_factor, cov_diag\n",
    "            ).rsample((n_samples,)).permute(1, 0, 2)\n",
    "        else:\n",
    "            out = self.mu(out)\n",
    "        if as_probs:\n",
    "            out = torch.sigmoid(out)\n",
    "        return out.view(xs.shape[0], n_samples, *xs.shape[1:])\n",
    "    \n",
    "    def log_likelihood(self, xs, ys, n_samples=8):\n",
    "        params = self.forward(xs, n_samples=n_samples).view(xs.shape[0], n_samples, -1)\n",
    "        logprobs = torch.distributions.Bernoulli(logits=params).log_prob(ys.view(xs.shape[0], 1, -1))\n",
    "        logprobs = torch.logsumexp(logprobs.sum(-1), dim=1) - np.log(n_samples)\n",
    "        return logprobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def compute_cov(self, xs):\n",
    "    out = xs.view(xs.shape[0], -1)\n",
    "    for layer in self.layers[:-3]:\n",
    "        out = torch.nn.functional.relu(layer(out))\n",
    "    cov_factor = self.P(out).view(xs.shape[0], -1, self.cov_factor)\n",
    "    cov_diag = torch.exp(self.D(out)) + 1e-3\n",
    "    return cov_factor, cov_diag\n",
    "\n",
    "StochasticMLP.compute_cov = compute_cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba21e3d17d4b4aeca4e4377d653fa6d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=20000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "net = StochasticMLP(inp_size=DIM*DIM, n_hidden=[64, 64, 64], cov_factor=2)\n",
    "opt = torch.optim.Adam(net.parameters(), lr=5e-4)\n",
    "\n",
    "N_ITERS = 20000\n",
    "pbar = tqdm.tqdm(range(N_ITERS))\n",
    "\n",
    "for i in pbar:\n",
    "    xs, ys = next(producer)\n",
    "    xs_t = torch.as_tensor(xs)\n",
    "    ys_t = torch.as_tensor(ys)\n",
    "    loss = run_step(xs_t, ys_t, opt, net)\n",
    "    pbar.set_description_str(f\"loss: {loss:4.2f} \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Sampling from $\\mathcal{N}\\left(\\eta \\mid \\mu(x),\\Sigma(x)\\right)$ results different probabilities for the output. Applying argmax on each gives samples $y_{1,2,\\dots}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADrCAYAAABn7V3CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAARH0lEQVR4nO3d73HbRhrH8WdvUgDOM/f6ZpgOZKeCyB1ISQVHdSCPK8jIHcipIMd0QKWCmO6Aurn3OVMd7L3As+JqiT+7IEBSD7+fGY5sEgssF8APILBcOu+9AADs+tuxKwAAmBZBDwDGEfQAYBxBDwDGEfQAYBxBDxyIc+7y2HUQOZ16nIOhbT32OiLogQNwzs1F5Jp6tHPOzZxzG+fc7Nh1GdG1tvmhyjX6bqwZAWgWwtV7/556dJqJSCUiFyLyeNyqjMN7f+OcWzrnxHv/eepybQh6YELOuQsR+eC9/556nK1rEVk55754778eoNwOLt0A01qIyM2xKyGnU4+z471/krrt/3DOVVOXa0LQAxNxzt2JyFfv/QP1OG/a9l9E5OMhyqUcY90A49MzsI2IfO+9P9r15lOpRw7tabKU+j7C78euz9j0JvNaCtfF0HIxzuiBafwqIg8nEK6nUo+zp+vgQUTuDlEuxs1YYGR6Fn0lIsfuZXMS9cAL9yKycM5Veg1+6nIiwhk9MIW5iDydwDXxU6kHlF6SehKRnw5RLiDogfHdSP1R+9hOpR546UGG9YAaWo6gB8akN85mUt9UPPt6oNFSRC4GfAN4aDmCHhhZGKPk2GfSp1IP7ArrpHQ8m6Hl8oPeOXfrnFvrY9HXgV/HrVg75zijaOCcu3LO3etjqY+7jHa91fZfafuuknkuotdbN4hk2oUu/yp6fa51C/O6ip6Pl++jMpW+h7tomoV+K7PpfcTvf1U6tseQNtTl+JbHTjmtZ9O0bT0g3os895TIfR+X0b61StvLOXeRtHljm+5bD12Wuf1ct9nwnpbpGXG03a6iNl5mtPEg0Topukk+tFwo3PuQ+lt1d9H/vYise8rc63Q+Zxnn9NC2uUyeq7SdvYhcdZSda/lNvB5E5FZE5vrvy6621+WsRWSWLH+lj0Vch/C8/vsqWb6PlrloWNZKp7uInlvG/4/q75vmMWYbynYslXn8HnqmX+s0G63nrKNem9BWme8j9JGeNbWX1F3qlsm6uu+q85B6RNvF0fbzaLttfV8D57kSkSpql030/7m2/1VSbjV2XZL5r/vadtRyGTO+FZH75Lmwc1z0VMiLyHKKhnqtD91xu4IibGCXPfMJO8Vag+s2em0eraMqKXfftu40dHxOQETTev33fc90m2j5Vcu0OweFidvwInoPrdNq0C7b6p1Mm32witojPqjehXnoa3ct5cJBqq0tS+tx9P18oqBfycuDZNg2bvXRWO+oLpux6pLMfyEDDo6Dy2XMeJ38v+rbOZJpbksrtWcDVlM9Rqqfl+QMLXk9hE9n2EZtvG7aWHUjTs9SwrxbN14pOJOJ1vGiq32iMHj+1NEyXTirbzxojN2GyTIbz5KkPlBlnRnL9qDWGM4tZTbJ/5/PkLvaQbafXHaWNbAeR9/PZeSgb1p30ba47lqvkpzIjFGfZP53Q+Y9tFznNXq9Zpp+FTnux/mlpWh8bfggN4P0OpuX7Uez0R/Oufs96xiu+V1KfQlkh9+OUneROZDRTBqGdPXef/K7XyP/Wf92rZOwTn/umCb1zXd/iSPU76PvHnI1TNfaq2DsNvTef5K6PWbOuUXDJPeSP357qPf/cibW+x7pugjzePLed3Wl+6Z/m+7DlNbj1eznhW5E5LfkudA2b0Tkx46y8TZYjVinIKyb0h40g8r1fTP2RnY38rDxfe3YuX8I//B7Dq+Zy3v/5Jx7K/UKnErbBp/rUbZh1rVjPEq9It9I/SWJPk0B1SSE5LeOadb6t2RD6rsR96R/+9ovTNe1Dqdow2sR+Y+IXDnnrsIBUg/s9z7/hmalf/uWF7yX3SB6p3/7xiAPbdS0nkrr8Wr280KXEoV5cnP1Q8/JyXO7TvTewrKrQ5TrC/pf4o1cz45CY3Wd3WZ17XLO3YrIY8OZ5yAnurE90w2rZDzwKnO6Y49jkrv8ves5RRvqScK11AeshXPue9EdvXDbDOHbdSCNLXz0rVXtDRLq23fwDPth1fBaaT0m3c+PKA3zd9G//91T9q3+fWqbYKT8GnpiWlSuM+gb3kD8ca6rocJG0hi8ukHfSf3R+yR/1uxQ9CzjUl6G15vkb6eCM86vuqyus/VQj5KdN3f56/5Jyo3Uhg/OuU9SX7Nf6HNvu0vtx+8OTXDZ8Vrq+RLPCPWYZD8/toY2DN0Suz6lBK0HsZHyK/cgPEq50i9MhY9zD20NlfTdXiavzbW/bdO1s7MR9dvdSD264KPUZx83el126EbQJ5ydveuYJqy/k14/U7Sh9/6D1KF1Icc5S30Ooq6Jkn7gU3ya22s/P2G5Vxoq2R5Il9Hzrza/skevTD7OdV0Tfu7Mnx5R9UbcZ52fpR8AzqZnn39I/ZH7oONue+8fw1lrfC06qdtMRD6f8mWwidvwi2h3VefcbwPbYejH8dxLIfFN6K5pi+sxxn5+inSbqfS/fQemxpvMI+dX6eW1vcqVDFOc+3EuNNLBgyKzl8pgGR/3OunGEb7J+v5IO8ifUq+bX51zT6EOuiMspA75k/3JuSnbUM9S30j9cXwh9cHk7wWzCDtfNWDZJdfn4x5RTWeWg+shr2A/HygO775OAeEg9lRwWbRUFZZxiHIlQR/e/GNP4IWzgYOGmNv+ks6Uy9g3BMMZ0ufSgHLOzUba6D56799qsP/snAvv51Hqs+NT33EnaUPdfu7CdXnn3O9S98JZeO9zr8Pus36yrs8nZ9xfW9bXPvU46f18D6Ncn5/AQa7VlwR9+KjSGgTHvG536t0rNVhzejKINL+Hmz0uJYQ6PPe51/mceqi/MHEbLiS6sea9v3bOraUO+3lP//8gBGxJr6Ag6/q81N96Dj5MUI+T3s/3EOrcuQ+3XZ+fQFg3pQflQeVKgr7Sv11HkqNetzvxs9Gsfrm6oVUT1aGSbVC+RpO0oXaTa+ovfy31ZaJ751zvz/HpycaLehbICiLZ3iht/USzZz0q/Xuy+3mppP/8oOvzEwjdd58OUa6k103OESScbZxy4B5dz72En2R7/S3dUZ9kDxqOMw22V22sNtSz0x+abuhqe33S/+Z+Ke2rdPdq2pFcn+/6VvCtvv6YcQmxuB7K4n5eEt7h/seU1+dF6nUzpP0GlSsJ+rChN248rh4ut9L/nvxR/gjiNmkcjlfD672I/KJPVdHLY12j/ywid865TTQsbngsdWjaeeEQrWNdLquSv6lR21BDfikdl4G0y+Wj1MMp5AyB8SDln8jiIJo19ejQut5pXXL69w+ph4jN/Tx8g7fvvoPI4a7PVwOXMaxc4YA6S6kH1JlHz4VeEGGQpd5RA6Nyo45Ud+oPqS+bbJret762kO0gamHgpUrb6i6ZPh6hMrsN5eWwu32Pddu6lLqLX++AVlr38J5bR4DUesUDejWOmDhGG+p08bK66l/JdqA3LzpyZsf7CANz9e4DUZlFVNcL2R1CumgI56H1iMqOtp8XLnf00St1vmF76RzkTV4OZNY6+F4ybXFdZTvoXlH7DS3nvS8L+mijW2vjhRV/EW2Mvm0nGKuhXvNDg+NWXo79vkjDJtrhNxKNYpiETtOja7jdK9mOqb6zjuTlWO2NO3TyfPpYJcvqqmccZF3T7YTwPm3YtpyONlt21K1t9MzeUEmmDwfe+6jeS31v67CPDdjWiuqRlB1lPy9c5lRBv5Tk4NkyXXzy1DftPkF/17XNjV3O+wFB31GJEEBZ41Lv01A8Wtu06ngtBG/JmWY4gyj+oYOSeo4dGMd+aCjmDm1cReEydsBl16NgnkX7eeG8Jwn6wvbykjEG/Z5Bv5KCT2f7lvO+Z5jiXHpNMefbdJiQ777++KvUX2nPvr7ntzcjZ2P+rFpaz556v0b3kj/MdHx9fuybmyX16HUG+/nk1+ej70EUDXk+tFzQG/Su/j3LcOOuapks/H7mo8/rb4wDinp1DOkXHMqc5ZAVQ+jB9ElEPmZMPtm3MEvqce77eeE3k/cxl3pdlx5MhpYTkbxeN+Hm1kwadnZtoDD2xlmPRHnCQp/ofb5A8zROVc7GB2npGZSY+iwytx7nvp8fqv/8RxH51wHLiUhe0Ff6t+3r1uGjxE3L633z5UxxYnpp5EFejmOS60b2OJM4V3rG+01/valR8i3MP49VD1Xp37H389ei9JNVpX+z80vXwTdfOAjf0HIvZNwA2PmF9Oi10EUt686+1EfNcAc87r2w1uc7uzTxGP6QesPcSMHNOdnenCru8cHjeXvfdLwe9/KYrI376qHTjLaf71nPo9yMlW0XzNabnfvmly5jSHfXQeVezCNjIVey+0vqF/rc5hgrhcfAlV2H/UI6uljqdPNo4x39h5HP6aEheZ88N9PQiAPjVverxnUyRT2S14++nx8r6OVl//lJli31/Y3iA+XQcumjd6wb7/3vzrknqX9irZL6W5CPIvKbr39YGa+Ery/hXOt6nEu9Tpsm/Soib7293jAH572/0W8dX/r6V6xm0vxLW883OmXYvZSiejS8fs77+Tep3+uTn+D3IfRbzZe+8BfLhpZrnJceNQBMRIPzD6mHgT7a7/ueSj3aRENSHPQHeaakB/aFiPxYcuI0tFybUfrRA2inO+qPsj1rP+t6dHiUuneXpZu9dzIsrIeWa8QZPQAYxxk9ABhH0AOAcQQ9ABhH0AOAcQQ9ABhH0AOAcQQ9ABhH0AOAcQQ9ABjXO6jZUM45vnJrxNu3e4+pBCDTarX6y3v/jzHnOdkQCAS9HQyTARyOc27lvX835jy5dAMAxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8AxhH0AGAcQQ8Axn034bz/EpH/Tjh/HIhz7thVAM7JP8eeofPejz1PAMAJ4dINABhH0AOAcQQ9ABhH0AOAcQQ9ABhH0AOAcQQ9ABhH0AOAcQQ9ABj3fx4+b8V0TZDmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADrCAYAAABn7V3CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAR3UlEQVR4nO3d33HbSLrG4bdPbQAoV53rrUNnIHsjGDkDaSaClTKQyxFsaTKQJ4IZOgN6IpiRM5Cmzv2upQx6L/C12ILwr0FApD79niqWbAlNNhvAC7DRaIYYowAAfv3PvisAAFgWQQ8AzhH0AOAcQQ8AzhH0AOAcQQ88kxDC8b7rIB1OPV6DqW099zoi6IFnEEI4k3RKPbqFEFYhhLsQwmrfdZnRqbX5c5Vr9be5nghAuxSuMcYP1KPXSlIl6UjS7X6rMo8Y43kIYRNCUIzx89LluhD0wIJCCEeSPsYY31KPV+tU0nUI4c8Y47dnKPcEXTfAstaSzvddCR1OPV6dGOO96rb/PYRQLV2uDUEPLCSEcCnpW4zxK/V43azt/5T06TnKNQXmugHmZ2dgd5Lexhj31t98KPUYw0aabFRfR/iy7/rMzS4y36hwXUwtl+OMHljGL5K+HkC4Hko9Xj1bB18lXT5HuRwXY4GZ2Vn0iaR9j7I5iHrgkStJ6xBCZX3wS5eTxBk9sIQzSfcH0Cd+KPWAsS6pe0k/Pke5hKAH5neu+qP2vh1KPfDYV00bATW1HEEPzMkunK1UX1R89fVAq42kowl3AE8tR9ADM0tzlOz7TPpQ6oGn0jopnc9marnxQR9CuAgh3NhjPTSA3+atuAkhcEbRIoRwEkK4ssfGHpcj2vXC2v/a2ve68Zzr7O+dG0Rj2bW9/kn29zOrW3quk+z3+evHrExl7+EyW2Ztd2W2vY/8/V+Xzu0xpQ3tdWLH40k5q2fbsl0jID5IDyMlxr6P42zfum62VwjhqNHmrW26az3stdzt57bNpve0aZ4RZ9vtddbGmxFtPEm2Tooukk8tlwoPPlTfVXeZ/T9Kuhkoc2XLxTGv8Zoe1jbHjd9V1s5R0klP2TMrf5evB0kXks7s38d9bW+vcyNp1Xj9a3us8zqk39u/TxqvH7PXXLe81rUtd5T9bpP/P6t/bHuOOdtQ27lUzvL3MLD8jS1zZ/Vc9dTrLrXVyPeRxkiv2tpL9ZC6TWNdXfXVeUo9su1ib/t5tt12vq+Jz3ktqcra5S77/5m1/0mj3PXcdWk8/81Q285absQTX0i6avwu7RxHAxWKkjZLNNRLfdiO2xcUaQM7HnietFPcWHBdZH87y9ZR1Sh31bXuLHTimIDIlo3276uB5e6y1686ln1yUFi4DY+y99C5rAXtpqvejWVHH6yy9sgPqpfpOexvlx3l0kGqqy1L67H3/XyhoL/W44Nk2jYu7NFa76wud3PVpfH8a004OE4uN+KJbxr/r4Z2jsYyF6WV2rEBq6UeM9UvqnGG1vh7Cp/esM3a+KZtY7WNuHmWkp67c+NVwZlMto7Xfe2ThcHDp46O5dJZfetBY+42bLxm61mS6gPVqDNjbQ9qreHcUeau8f+HM+S+dtD2k8uT15pYj73v55o56NvWXbYt3vStVzVOZOaoT+P5L6c899RyvX301mfavBU5H8f5Z0fRvG/4WS4GWT9b1Paj2eyPEMLVjnVMfX7HqrtAnojbWeqORk5ktFLLlK4xxp/j09vIf7KffeskrdOfepZp+h77b+JI9fsU+6dcTct1jiqYuw1jjD+rbo9VCGHdssiVxs/fnur9nzEL23WP5rpIz3EfY+wbSvfdfrZdhymtx4vZzwudS/q18bvUNm8k/dBTNt8GqxnrlKR1UzqCZlK5oTtjz/V0I08b37eenfsf6R9xx+k1x4ox3ocQ3qlegUvp2uDHutU2zPp2jFvVK/KN6pskhrQFVJsUkt97lrmxnyUb0tCFuHv7OdR+abm+dbhEG55K+kvSSQjhJB0g7cB+Fcdf0Kzs59DrJR/0NIje28+hOchTG7Wtp9J6vJj9vNCxsjBvXFz9OHBy8tCuC7239NrVc5QbCvp/5Ru5nR2lxuo7u+0c2mWNfa56Qz1SXfFf7cxqJwe6sT2wDatkPvBq5HL7nsdk7OvvXM8l2tBOEk5VH7DWIYS3sh295VNRnxS+fQfS3Dpmd63aaJBU36GDZ9oPq5a/ldZj9v38QDTD/H32798Gyr6zn3n5JfJr6olpUbneoG/ZyPOPc30NlTaSR8FrH1VX+UdSa7h1COFc0ruBo6w79v6P9Ti83jR+9io44/xmr9V3tp7qUbLzjn39m+FFys3Uhl9DCD+r7rNf2+/e9ZfaTXw6NcFxz9+aHrp4ZqjHrPv5oWhpwzQsse9TSvLkIDZzfo09CM9SrvSGqfQGv3a9ocbY7U32+5XqRnp05LOz8A+qN9zfC+vzImXjdu9Uzy54q/rs49w2oqkbwZB0dva+Z5m0/ppdCgdliTaMMX5UHVpH2s9Z6kMQ9S3UGAe+xKe5yfv5gRv1CcQ+0aQ23tjvXnR+ldwwVWl7BO/rE34YzN84op6ro9/Rzki/qL54tshNCofC3t9fqs8c/xljfBdj/PIcn2SsnX+WVOU3RzXqtpL0+ZC7wRZuw3Qd4WKHbXHqx/GxXSH5uutbtrgeM+znB8nWZWX/HTowtV1knju/SrvXdipXckY/9uNcaqRmUBxL+qt5V1rmD/vZd7bZy87yFntMrVdWv5Xs5g1JHwr7f+fyh+p180t+VpY+gqoO+YP9yrkl29Da4422FyZLz9DSzldNeO2S/vl8RFTbJ6/J9dDu+/mhysN7aFBAOojdZ92ic+dXlV5j5PI7lSuZjz69+duBM6d0RGse5b/b36qOcuk5u/7eK2y/SWcxIYRdQzCdIX0uPQsKIawK+uL7fIoxvrNg/8n6FqW6C+D0kM/kzSJtaNvPZeqXDyF8UT0KZx1jHDu8cpf1M6p/vnHG/a1jfe1Sj13380O1U/+8lsuvZ+mrLwn6dCTrDIK+frsY44fQP2l+uvA1KWgOfXilBeuYkQxS+3s4DyH8uksQ29nIrfTQt3joof7Iwm24VjbEMMZ4GkK4UR32ZwPj/5MUsCWjgpJR/fOq73pOPi5Qj5328wOW6ty7D7f1z0uL5FdaN6UH5UnlSoK+sp99R5LefruBI+mP2vFLEg78bHTUuFzb0KqF6lBpG5Qv0SJtGEJIt/83d55T1d1EVyGEwa/js5ONR/UsMCqItL1Q2vmJZsd6VPZz8n5+aBr95lP65yXNnl9p+G7fc85WrqSPfswRJJ1tFAWuXRisJP2zpNxLNdDf/6O2HwObO+q9dmDhuLJge9HmakM7O/1HW1+/tVcaZTH2prRvKrzO1Oif77sr+ML+fjuiC7G4Hmax/XyPSu7gTdc/7ocO7MnE/Hqvae03qVxJ0KcNvXXjCfV0uZX9d/RR3jbydBHQ3Te/Z/I2aZ2O18Lrg6R/2a+q7M9z9dF/lnQZQrgL22lx02MT6qlpzwpHnMzVXVY1fjbN2oYW8hv1dAPZkMtb1SMqxkyB8VXln8jyIFq1XfCzul5aXcaM759SD2mh/XzP0h28Q9cdpMKbwHbIr2rsa8xSbsyEOHE7oc5G9YQ6Z9nv0iiINMnS4KyBjee81sAkVl4eqrtN0qyDbROOrbWdRC1NvFRZG182ls9nqBw9CZQeT7s79LjpWpeqh/gNTmhldU/vuXMGSKtXPqFX64yJc7ShLZe/Vl/9K20neouymTN73keamKtkH1hndT3S0ymki6ZwnlqPrOzs+3lhnWedGjjbXnonedPjicw6J99rlCnOL20n3Stqv6nlYoxlQZ9tdDfWeGnFH2UbY+zaCVqe66q0kV76w4LjQo/nfl83wybb4e/yNmqETtujb7rdE23nVH+yjvR4rvbWHbrx++bjuvFaffXMg6xvuSchvEsbdr1OT5tteurWNXvmYKg0lk8H3qus3ht7bzdpH5uwrRXVo1F2tv284DWXCvqNGgfPjuXyk6fB2SE1Mb9kM1A+V7kYJwR9TyVSAI2al1ot81/z2HkdVD1/S8FbcqaZziCKv+igpJ5zB8a+HxaKY6c2rrJwmTvgRtej4DmL9vPC514k6AvbK2rEHPS75Fc6cD5XuRgHpikey/qpxtxNl5Y/kfQ2tlxQmuPGpNcq9vc//qL6lvbR/XtxezFyNecdy816DtT7JbrS+Gmm8/75uS9ultRjUOl+/gKNnSJhcn5l90EUTXk+tVwyGPSh/j7LdOGu6lgsfX/mbRwYb5yNcuhqpKLvDcWwbFTHlDHPqcyUoXqvkh1M7yV9GrF4212Yz16Puffzl2bsnckz5NeZpg0jn1pO0rhRN+ni1kotO7s1UJp7o/cOQjsr/BDrkQxtpg45Qr80JnqXG2ju56nKq/FR405alp7qd2w9ZtvPX6jBIZgz5dcnTRtGPrWcpHE3TFX2s+t26/RR4rzj75IeDUP62jFM7Y3qxl50atjXKNY30HxVPb68dAqHc+14I9trFGP8HEL42HdXbeMuzD/alnmOepjKfu60n79gvZ+s5sivUH+T1/dYOIx8arlHRlwAePIN6dnf0hC1wSv7Gh4t0jv6gcduD9U78p0KLs5pe3GqeMQHj4eLi3c9f89HeSzWxkP1sGVm2c9nqOdeLsZqOwSz9WLnHPllrzFluOukcvljTNfNR0mf8ps4QghHduPEj6onwur6KPMg1lPJhqHHiPpgglhf8Pw/SbfWF3vR1RdrN0yl8edvo88zuMXF+lPQb80zwBDCyvp68/3m2Par6rnq0TDLfv4SNfrnW7+HYdf8CiFcatpEfJPKNQ123cQYv4QQ7lV/i0ql+iPKrWb6+j88Hwv70+yi0drmRGn6plf4bV9LiDGe213Hx7H+FquV2r9p6+FCp6ZdSymqR8vfX/N+/l31e72PC9ydbwf141j4jWVTy7U+l300ALAQC87fVZ8V7+37fQ+lHl2yKSlOlwjcfcj69n8oOXGaWq7LLOPoAXSzHfUHbc/aX3U9etyqHt3lqavwUtPCemq5VpzRA4BznNEDgHMEPQA4R9ADgHMEPQA4R9ADgHMEPQA4R9ADgHMEPQA4R9ADgHNj5qOfJITALbd4Nd6942sUMI/r6+t/xxj/d87nXGwKBIIerwlTiWAuIYTrGOP7OZ+TrhsAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcI6gBwDnCHoAcO5vCz73vyX9/4LPDxyMEMK+qwA//j73E4YY49zPCQA4IHTdAIBzBD0AOEfQA4BzBD0AOEfQA4BzBD0AOEfQA4BzBD0AOEfQA4Bz/wXibHExP8GsjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ys_pred = net(torch.as_tensor(inp[:, None, ...]), as_probs=True, n_samples=200).detach().numpy()[0]\n",
    "plot_im(ys_pred[0] > 0.5, cmap=\"gray\", title=r\"$$y_1 = \\mathrm{argmax}(p(y \\mid \\eta_1))$$\")\n",
    "plot_im(ys_pred[1] > 0.5, cmap=\"gray\", title=r\"$$y_2 = \\mathrm{argmax}(p(y \\mid \\eta_2))$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verfiy that the average of all samples resembles the average of all ground truth labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADrCAYAAABXYUzjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAASOUlEQVR4nO3db3Lbxh3G8Wft/HuJutM3SWY6oXMCuT1B5BtIyQkq3UAan6BD34DOCRrlBlJOUNM3kDp93dbiARxtX2DXXMH4s6AI/kDo+5nhSJQAYgmSDxbL3YXz3gsAsHtPrAsAAI8VAQwARghgADBCAAOAEQK4hXOuqN6sy4Txc84dTmk7+2IfP68EcLvfJN2mN970aOOcO5F0vKPNHYftDc45N3PO3TrnZrvYXl9hP9xWbj+bFirDZ9YF2ANXkubJ/bdWBckRjvqHkmaS3nnvr2xL9HjE8PXev9zF9rz3p865S+ecvPdvBt7cTFIh6UDSzcDb2sSVpHS/z5sWHBMCuNvNvoSYc+5M0itJ5woHDufc3Hv/wrZk0+ecO5B07r1/vuNNH0taOufeeu/f7Xjbo+G9v1FyYHDOvTcsTrbRN0E45w6dc36g205O33bBOTdXedR/EWpDK5U14QOaTXbiQtLprjfqvV+F7f62D22euG/0NWDv/ZVz7ljSM5VH+0OVNbz0NGge/n8uKR75nkl6rvLU6SD8lKQ3ki7DOpOoGYba15mkN6EmIK33Q/V3OeeOvPe/7qp8UxcOfmbNPeEz8lbrsx/sidEHsCQlYfHGOeclrdIAcc79LOnvbe1goRY4l3Qi6VplDXEqXoWfi/iHUDNyzrki/J6aSyKAtyDUOs9UHuwtnUq6ds4tkoMwRm70TRANipr7rW867/1VaAs9VRlAB4OUzEZ8Lp/sg5rwlcqzA2zHz5KurEMvbL/6hTFGbh8DeJXeSdq9VtUF64Ra8uutlsjeTGoM23vC/iqGLc7jEPblkcYTegtJR7QF7499DODeQhtpalG74OMwmS8eR+BEZXPYKHrJhGa5laQfjYuCTJMP4ND2u0z/Fk7XHmuXnZ1/Uz9hpypP+8fkSrzGe2PyAazQPlpzWnas8Zw67kToJzzKkUz7JowIm6nsUTMmlyq7HvI674HHEMB/rfuj9/7G+ouTXQnDSBd6ZAecgcW+1WOsAUvr8mHE9qIb2qaSL0mG3MaBylO+VfLnQtK19/516CP6j6ZRSqGMJyqHUa6S9d+p7Fq3aljntubv1cubxAEZZxnLStJr7/15skwsVxyGOvPeu6QMsftbrA3ehDLfe66h5v08WfaZpEWf4bPOuSOth5rG2l3jPgrrLNXc2+UP1fVCOesOUvf2S/BS+ticlSXsz/g4N5JO0/WTfRrDswjLneeOcvPe3zjnYvmGHp68NWHfnKp8zs8k/eK9P03+X6jscZL2+MneL6Plvd+rm8rgOUvuF5K8pMOaZS/D/7ykYoNtLVUGRdP/F03/V/khWoZtHzQsc1R9PjX/O8koZ9ZzVBlcvnzZs57/UXiOt+l64bldNOyve883vAYHleXOwnKfPEbLfj6s/K1QOfrMSzpqWK9Q+YE9SZ9Dx/LXYZnbUM5Zy/tw2eO9FN8PRbL+bXI/9k8/qqy3bCtzw7auVVYAtvm5O+xbjszHLcJznCf7Yhb2zXVy/7rmfeVbXp/L3PeX5c28ABu8YE0BfJS8WGfJB2meE04N22oM4PCBaf0AhgCrDeCwfu2BI1nmID6Hju1khWqyrzqXraz3MbjD7037JC53G+4vmva7Og5OyXLzpg9Z5XEa92NlX3bt88twa32/qMcBJCnnrKbcZ+F22bBeDL7bHtu66PsaZzzmUAG8VE0lQ+uD9ElYJg3fRfJa1lZQ9iWAp9QGfBFOq6+1Po18ruG+JDmV9I+2Bfy6W9A9odlioXLocGMboi9Pr84lnYVTcBP+/mn2XA3DXf16QpQinM4vfXPf5Ljvur6xP5O0aPlS6W9JuRol+1Jq6IYYtvHMe/+ypdxKypLV/BCXr+zHIvw8lfSTb55BLa5T9Phi7aZSzlEKzQ5Xvr4pKj7vhaTqREN/SX4fWxt8L1MK4CuVw2vjCzJT+YErBtregaQ/ZixX9waJ85R2finmvY+DRsYyt+n7tnDS+oPzquGDVV2uMSSS/tuHamjLTz6YB10DEMK+vJI0c85d1CyyUN5cvrHM/8tYVqo/WMfHeCbph4xtSfnv5Viu0QZweK3mkv7esEj65Xn1gHms8rP9wu/5F+lTCuC59/441F6cyjf9iYYLrpXKmmnXkOZ/pnfC8gcqO/DnvnneqawBjWEQRdcZxSr87Jo3OS7XNiz6Jrm11XTifswZYn0ctn2UnlWEXiK58ygU4ecqY1mpPIB8PBhV3jPnHQe0jyHq879wio9XZC5v4UeVX7StGv4f99Gq+rx92YPpdY/9MVpTCuB7Qu3rpYZ7E8YP1NI5t2gK4po3yk/hZ58jdwyzXV1poU1uuR9cM/Her7z3z8Mt58NW5Dym1vvxInTROwz/y52gKAZ97pyz1ZBNT6F/6Vg3zti3aluowZjn/Hir9pnb4j7a6yaGLpPuhubLafquNECfSO/9eQjdQ5U17ZPQ/edKZS3xqiE0Yln6BNSqsq6l3HJfD7HxZJ+ns489q/xsFd4Xr1W2L1+Evw02NWlNO39s733XUfuVNutvPPrJyNsOqJX5SsY20GWrJlsDTtS19W1F+OLkVPdDKU57uQyXiykqq1Xv58hta5ykcIHFuXMuXucr9gE99WVf0d6B48t+ve9UnuruupaVFarhvRObICYdRBVpRWPSNeDHEMBXKtuRVnX/fOg3xd77N768DM1zlWH8RutAPlQZxEWyyphPC0cn1Hj/pbK2+jfv/Qvv/a8ZNcccsWknpy2/Tu/XMmynCHe7QnXTIOrbRDI28Qyhz/cke2nyARwa7P/Qssj1Jl28qpf5Cdt5E2plz1W+iW5U1mDSLwI/divqsbmc3hYP5pw7GVPXpVCWpcp99bJHG23OYx9qfZUVqbwCdq4YbMUGm07fN11fVG4aREVcr8c6YzLWYd5bN/kAbpME72qD1VubNkK73wuVgZsGfPzQ9ak9pUNvh2R9VYequI9b+0vXaTuQxC5QodfMryq7LxYNXdPqPKRWNnT7b2rvasCVZpd/tiw3hu9DHuzRBnDSD1Ha7FLznR3jwwdsHrYXl40f8j6nvHHZ1oEfU5J015O652+uO5idtjQrXCjpUeK9P1Y4UGZ29YsBvMkBKwZH63vuge2/sVz7ePqeBmtbhWMSc3pPKYCL3AXDm/s3lW/wxvbhDDlNFzfSehRUqMnFkUqdR/HKB/FBk6ukz7Nh0MJsRG1uWf1f+17hI4zQq+vvGwO5bdRdLM+qWsbMbacHhKHaf6UeV0ixEpq76g52PyW/1x6kwnqTuKbhPgZwkd7p+yYLofcvrWfvekjgvOpe5ON2UnH4bc4VbOM2Trf0gWobgVZs4fG3rmOE249aNyFVn9MqvRNe+7/WtSWHkI+jDnOaIt7pfn/eHH1CNQbRJl9E/UUjvuBAmCFwofJglw6GKbSu1NRWjMIyp2oeQbdX9jGApU/bbFfqqI04547C9ISXKkPoRVhvk+aHKM550OZclaANteBzSYdtp7whMOLl5jtrv5nXAovhc6/2HdbN+aBvqxdHUflZlQZU7T4KZX6p9Ycxfax7tfmwLy/VcuoauqbdqBzW3HWKe6X+B6w4vPZm4PbfYsP1diU9c6zOMxK7dRYN616oZQrSvWM9G1DXTeULcaTyQxinoIuXl4+3a62nD0z/vgjrpFMRLpP/N86m5NczNTXN/OWTbdRNJ1mEbbdNZxmnnPxk1jCtZ4P65LFrHsPnLJuUK+6vdIapCzVP7Ze1Da2nEfRqmVEslCGd0appus6D5PGq0zQehDIXWs/ydh1+nynMIBeWS7fVVv5C61nK4vtsVvc8tJ4drHUWtso68bl0zW6Xzj7XOR1pzX7pVa7Mx93abGjJ63GU/G2uMHtZsq158v8DNcyc1rCNvZgNzbwAPV74oW5dUx02BfBF8vtRCIPqrXWqxbBuoTJsl+F2qfX8qG1lu2h5Tl0HliJ8CK7jG7WurB3bWFaef9Y+7liu6UCW7p+4b88qyx1ofWBZtG2rZb9ctpTtk9dCGWFa8/jXba9rWO6kbbsd687bnuMWPodbmY5S6yljl+Fn3esZ99d17uepsq9HH8AuFBY1QpPFW5/MzL8PnHOFH8kpWrUsYyrbQ4VuazO/5WHM4XGPVLaDtvVhr1t3qbKJY6vzhiRNOLHr3qg55y5V7r8xzJ/SaF/bgNFiTAFXLcuYyrYFC2VMg7mBjdp/Qzlikwv2AAEMbMiXX6aulNcbJkvoAleEu337/56orPWN+Qs4JAhg4GHO1dBLY0MP6f/7SusrhGAPEMDAA/iye+D7LU6Wv9H8D2H77/ehfRZrBDDwcKfKuLxUpk37/8Y+tNgjBDDwQKHN9ZeMwRutKu2/2fN+hJFlvScsgr1JXxED2BXv/alzbumcO3xAEL5XOQpslduUELqHHW67Kxx2gwDu9iztZjSxblTYrh8k/eacO+7TfhuF91b2DGuhxjxX+1WVH42aCx+sbEqSjyaIbnG48K2k26nMQ4rtCwH6g7bXHtxlLumHHVUKblQG2ign+QlfQt4mt02ucLJzjIRrEaYPrE4+85ZaMDAu4WygOiHXez/yS9cTwABghCYIADBCAAOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMNJrMh7nHMPmJuKrr77KXvbLL7/MWu6LL77IfszPP/98q8v1eT655ezzfIA2y+Xyv977P1X/zmxoE/P06dOs5Waz6rD5Zt9//33Wct988032Y3799ddbfczcMkrSd999l7Xct99+m/2YQBvn3L/r/k4TBAAYIYABwAgBDABGCGAAMEIAA4ARAhgAjBDAAGCEfsATk3uJqSdPtn/svbu7y142t5y5yznnsrcNjAU1YAAwQgADgBECGACMEMAAYIQABgAjBDAAGCGAAcAIAQwARghgADBCAAOAEYYiT0zukNzcIb599BnebFnO3Ms2AUOjBgwARghgADBCAAOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGGAk3MUNcxHKI0WjbLmef0W2///579rLAkKgBA4ARAhgAjBDAAGCEAAYAIwQwABghgAHACAEMAEYIYAAwQgADgBECGACMMBR5Yoa42GXusnd3d9mPmSv3+fQZWg2MBTVgADBCAAOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMMJIuInJHbXWZyRc7kUsnzyxO573GYXX5wKewJCoAQOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARAhgAjDAUGZ0++8zubTLE8OYPHz5s/TGBTVADBgAjBDAAGCGAAcAIAQwARghgADBCAAOAEQIYAIwQwABghAAGACOMhJsY51zWcn0uYpn7mH0u9Jn7mEPoU05gSNSAAcAIAQwARghgADBCAAOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGGIr8SPW52GWfYcvbljtsuE8ZLS8yCqSoAQOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARhgRNzBAXnMy9gObTp0+3vu1cfUa3WY7sA1LUgAHACAEMAEYIYAAwQgADgBECGACMEMAAYIQABgAjBDAAGCGAAcAIAQwARhiKPDG5w4aHGI7b5zFzl932csCYUAMGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARAhgAjBDAAGCEAAYAI4yEQ6fcC332uSBo7og9YMqoAQOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARAhgAjDAUeWL6DAfe9mMOMbz46dOnW39MYCyoAQOAEQIYAIwQwABghAAGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARRsJNzJMnecfU3OX6LGs5Cq+PPs8dGBLvRAAwQgADgBECGACMEMAAYIQABgAjBDAAGCGAAcAIAQwARghgADBCAAOAEYYiT8zd3V3Wcn2G4w5xsc1cQwxF/vDhw9YfE9gENWAAMEIAA4ARAhgAjBDAAGCEAAYAIwQwABghgAHACAEMAEYIYAAwwki4ickdtZY7Yk7KH402xIi5IZ4PF+XEWPBOBAAjBDAAGCGAAcAIAQwARghgADBCAAOAEQIYAIwQwABghAAGACMEMAAYcX0ueuic+4+kfw9XHACYpD977/9U/WOvAAYAbA9NEABghAAGACMEMAAYIYABwAgBDABGCGAAMEIAA4ARAhgAjBDAAGDk/x/iP6dsu+NmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(ys_pred.mean(0), cmap=\"gray\", title=\"$$\\mathbb{E}[\\mathrm{softmax}(\\eta)\\mid x]$$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAADrCAYAAABU+oCOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAASjUlEQVR4nO3d33HjyLXH8dO3HADuVvnZZW4GmnUE1mQgzUZgKQNNTQQubQbajWCXm4FmI/BoMqBcfvcdMoO+DzgtNiGg0Q3iABD5/VSxOCMCRLOBxg9/G857LwAAjO1/5i4AAOA0ETAAABMEDADABAEDADBBwAAATBAwwEDOucu5yyCynHKcgynreui0ph4vhYABBnDO3YjINeVIc86tnHNb59xq7rIca4a6vtZpLn28Tn8a88uAcxBWNN7795Sj10pEKhG5EJHneYsy3Bx17b2/dc49OufEe//zUsdLIWCAAs65CxH56L3/nnKch5nr+lpEnpxzX7z3Xxc8XisOkQFl1iJyO3chZDnlOAez1bX3fqfT/sM5Vy11vC4EDJDJOXcvIl+9958px3lYQl3rtL+IyKclj9fG0RcZ0E+35rYi8r33frZzCUspRy69MulR6vMXv89dnhJLqmu9SGJTWpapx2tiDwbI84uIfJ57RbOgcpyDxdS1luGziNwvebwmTvIDPXRL9kpE5r5qbBHlOAcLresHEVk75yo9V7LU8V6wBwP0uxGR3QLOeSylHOdgcXWthxh3IvJhyePFCBig363UhwvmtpRynIOl1vVnGXZF29TjiQgBAyTpyc6V1Ceqz74c52Dhdf0oIhcDekaYejwRIWCAPqF/prm3ZpdSjnOw5LoOZSrtN2zq8USkIGCcc3fOuY2+1n034WgfRBvn3BK3AkbnnLtyzj3o61Ff9xn1dKf1+aT19dT4znX0eedMbgy71ulfRZ/faNnCd11Ff4+n76NxKv0N99Ewa72zue13xL//qbRfoyF1qNPxHa9X42k524btulrmvcjLVTW5v+MyaitPzfpyzl006ry1To8th07r5Nut/sZtvLwnhl27/r7RBtV1NI2baN5unXMPjc8rLUeYL48Z818aZSq6+GDq8eIv6H1JfSfrffR/LyKbnnEedDifM423/NLfetn4W6X15kXkKjHujY6/jetVRO5E5Eb/fZmqS53ORkRWjek/6WsdlyH8Xf991Zi+j6a5bpnWkw53Ef3tMf5/VH7f9h1j1qHs+7m6iX9Dz/AbHWar5VwlyrUNdZX5O8L9A6u2+pL6ss/Hxrx6SJV5SDmi5WLWdhstu52/bYRphHl+l1mWg+V3jLqOlq0nncdVtDxso3Ydlo9m+/Gp5bAxnU3ffFzCeN77/oDRBvjQMUNTMyk04kerBWsJL12YUiuosIK57PmesPBvdIV5F312E9V51RjvoWte6MLscxpLNKzXfz/0DLeNpl91DPsqjIzr8CL6DZ3DaoN+7Cp3Y9jskIzqIw7z+/Ad+tl9x3ghHLvqsrQci2i3YhwwchgafctHmBd+zHneWE5vOuaF13b8JIfh8hCV/9W4HdNZ9/2GJYznfV7AbBr/r/pmaGOY5FZF5g+srF4jlM1LY4u08XlY6SVX8lGdbdoaty6kVx3fvU18b1g59zbwaJ6tU3Uj+5XQy15Wx3ChYbWG1dh12Jhm6xaX1AGZtXUq+zBtDYWOcbaN/8crkM56kP2e2qtpDSzH7O1Wv9M6YLJCQ4cNbaFz/g+pax3vpmscqY8StC4DUZl81/Kf+M1Zw881nvc+fQ5Gj6E3u3eIr4n+0jFqfK5g8IkyPVbppd66M3k1j48Wli8cN72UeiF6xe97JL3oO/6tVtLSrbn3/if/uquNH/U9VcdhHv2YGKbpm0/fWBXK98mnu/UOw3Ue7x67Dr33P0ldHyvn3LplkAfJf6ZHKPf/5Qys5wCa8yJ8x857n7rc85u+t51nKy3HrO12YqHMOT3/hmUt9duK6lrk5abMexH5Z8cgf4v+3VzfXIvIRxF55/PP+YSylV7ZNfV4vXfy38rrxhgaydfESuilQv0RXT5773fOuXci8t3Q78jQ1dhyPMt+JZpaaJ+lnjnfSX3jUp+2FWOb0GC+JYbZ6HvJwtF3gnen7311F4ZLzT+LOrwWkX+LyJVz7ioEs25MPBQ05Erf+6YXvBeRXxt/+0Hf+56vEeqobT6VlmPWdjuxnNCQxkn01PJd6fuuoAwfROS3RL2Gae+a9arL4k8F04rLVi18vN6A+WfcGDWpQ2WltvxHu8xvyQu6LlAlz4ioMoebu++j3OkfXU6LOtQNk2upVyRr59z3oivulr3AlLDSTwV4bO2jO7/1SqVQ3r7QDu2qavmstByzt9spNK6q7KvfeNjUhlFpXYfv+y3xedjIGLteh254TzZeMmBaGmO8m52q0LAwLzYcrOiW0qUcrjS/a7wnFWxhf9VppfZOQjlKFu7c6W/6Byk3Uh1+ds79JPU5mbX+7d2Y5WybZuNPl4nPml4OpY1QjnNpty+XzmbUbxj2uefwb7HURrCGe6X/HevS75Lwm3O84hstw272566ZVLhVcRKi+0W2UvfA+iz1E/Bu9bj74BnUI2yN/pAYJsyP5qGbRbGoQ+/9R6lXlhcyz1Z5WKklV9iNezIs9l5Ptd2WnH+Za+/sLZ7XGk3JjZaV7LdwUucISrYq3jzd2v631FvK//Dev/Pe/z72VlKb6Pht1XaTmZZtJSI/L/lQo3EdhsMhd7k3s7UYekghd6UWz7vUsMXlOPF2O/b5l9hY531Dve4Kjkz0GXIYb47xirrrz93NLtmq6JV55dVgx6zEdMsz3Hn/fqaG+S+p6/oX59wulEEb1VrqcFnso3Ut61C3yr+T+oT3WkT+EJH/LfiK0KCqAdMuOf8SX+HXtqc5uBwyU7u1dsT5l77l65i6Tk17zHVDpe+7hY9XFDC5xzCztipyuP0T5cw4545ZAYctwp9LV4zOudVIWzSfvPfvNFB+dM6F3/Is9VMEl77CMKnDcOloOO/inPtd6qvK1t773MuUj5k/WedfGnsYXzvm1zHlmLzdTiTuuqTvasaXQ5VxHWhIfWnUy2iHKHXehsOf/0oMdzlww2rx52JKAiZUVOqE1qjHcZd8mbKu0HOuzBFpL/+tc+7XYwJAt5KfRV5ONC49TA4Y1+Faokt1vffXzrmN1CFz03P/ThBWNiVXuQVZ51+kvkEv+GhQjsnb7URCmXcZRyHCsM22/lHq81Px+MfUddd0RdLLwUPh9MKwpWE49XhFAVPpeyrFRj+Ou+At8JcTswVXkYytkv0K+i0yqUPnXOgmpdkgrqU+HPfgnOt9FK5u4ByUs0DXSq0p7HF27sEdWY5K3ydttxMIy31yHjbCcxP9vZL6zvSD8YfUtd7YKi0bLfGhz9bloOOm2D7hkvvdwscruoosJ73C1thSQ8FEz3miD7LfQmoutDs5gq6UV7pCfdPGqkNdofyt7X4Xra9wU1vuzaxfJX2V3iuN8y+pXgzu9PPnjMO0xeVQJ9duG6FRJYar5HDPOP59H6R7xZ5d167uhftB6o2WuPfySvYXb7TuZekwt9LdA0CXH2TYvJp6vKKACQ2yteJd3c18pf99K1tBx4h/Y2u39LoAvZf9AlRFH491DuZnEbnXbsFD9/Dh9ajdgt8UXkE11iHJqvHeNGod6ornURKH2/TS5Wepu53J6Sbos5TvgcYrwFVb1/Ba1nstS879OUPKIXKa7TY+/9Jav+oPqes4iNvbtXRful9S1/EVgPH330sdHs+J71pLfVPsLnNaQSXD5tXU4/V3dukPOz0L3UrfRH8LVwGFzvp6ezU9lZfUu+lbaenMTz9by75jzdCRZaV1dt8YPu4xObtjQDnsfr7vtemaN3LYIV9nR4ey7348dFBZJcoVd/TY2oPvGHWow8XTSpW/ksMOBu/1u7p+R+isMXuZjtpC6Bm7+SiFokcZDC1HNO5s7TYq92idXUbzLyw3T/H802k+hWnKvrPG8PiLG0l3PJpd19LyqAWd3rrxXfFjEy6ko+fljOld5JZtzvFexh/wA++0wWyjBfQiajS+q7Ge4ktXWHdy+OyVdXMlF61otvHC3VjZtb1S3c5fyf6ZJq/qXA6fldK6Imn8vfl6akwrVc54BZoa7tXK/5g67JpOos4eE2Xr6s35YAWRsUyEwH+Iyv2ov20T2syAZa2oHI1xZ2m3YhMwzY2D8Dykjf62h+a8jH7/RjKCvaSuo+8O87dtuX2Mpz9k/ut33aeW76WM9zL+iDM9rChP+vkvU75SDT5a4ZdsWYetkUEPD8otp8WKaub5sJb8Lv6raAU4ahf1JeUo+E7Tdjt2wETfN3ireq66HnF+DX1WzWTjhVdpVzGt9Bhozt3CKODTx2Z/kbrrj+xjo35/knt1xF3tbd+7S/3/BDxI/uMWci9NtS5Hrzfabqe64m3Uuh5DdM9U0SNGph4v1hswrn6+eDiBXHUMFk6kPfu8+wtwhOgqpSH3LIRxhlzyepZ0RbYTkU8Zg1t0DVJcjhNut5P0OFA4z6dyI/VyVRqsU4/3ImcPJpxkXUnLSklXduFKitw7pHGccE/DMTfe7cYpytn4KB1XujVYd6qYW45Tbbdhj2uKzltz63oqn0TkH29gvBc5AVPpe1c3FmH36bbjc4xMD0F9lsN+pnLdypFbJedIt/C/hZvq2uR2DWJdDlXp+8m028b9L+bLb0Fdm9MyfPNlzzOafLxXMk7ybKTjBJ3sL9EbdGULr+EvqVcgWyk4ESn7K8YGXcFy7i+p9062ic/jS83N6rivHDrMItqtjHiSX/aXG/ulzPMJy7GVYZeoTzpe85WzB/NRRD7FNzM55y70Bq0PUneo2NWHEoz4ei/mryLyrMfa77qOteuNluH+ke/9G9liXRpf7/X91rxB0zm30q3ruB1cajuppipHwym227AHc9xWdYHMujalvQUM6Qx20vHa9PZF5r3/3Tm3k/rRs5XUd3k/i8iv3vvSZ0ljRBoy1zpfbqSeR22DfhWRd/70ru6anPf+VntJuPT1UzNX0v5kz5cT6DJOx4nJcrR8fort9ovUhyBLu1Y5Sl9dW9INl0tf+DTWqcfr/D7dHQKQSVfYf0i9F2DxBMo3VY6UqPuea3/s8fwZzVHXuvGyFpG/l2wcTj1eyij3wQDnRBvf3+Wwn6uzLUePZ6mvWHzTh2Vnqut7Gbayn3q8TuzBAABMsAcDADBBwAAATBAwAAATBAwAwAQBAwAwQcAAAEwQMAAAEwQMAMAEAQMAMNHb2eUcnHN0L4BJvXs3St9+wKyenp7+673/89zlCBbZVQwBg6ktsR0ApZxzT977H+YuR8AhMgCACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAIAJAgYAYIKAAQCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJggYAAAJggYAICJP81dgA7/FZH/zF0InA/n3NxFAMbwl7kLEHPe+7nLAAA4QRwiAwCYIGAAACYIGACACQIGAGCCgAEAmCBgAAAmCBgAgAkCBgBggoABAJj4f5eKcmlUTdi9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mode_pred = net(torch.as_tensor(inp[:, None, ...]), as_probs=True, use_mode=True).detach().numpy()\n",
    "plot_im(mode_pred > 0.5, cmap=\"gray\", title=r\"$$y = \\mathrm{argmax}(p(y \\mid \\mu(x)))$$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Extra: explore the covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cf, cd = net.compute_cov(torch.as_tensor(inp[:, None, ...]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAC5CAYAAAAWJ+UfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAKtElEQVR4nO3dy2/WZbcG4KctpQcKUqDBolUUCUYTw9CBA6f+1Q6deRpIiDHmAwFrgdaeW2z7jfZk75re67PdK9n7uoZksd7D79ebN+S9+0ycnJwMAP73TXY/AYD/rwQwQBMBDNBEAAM0EcAATQQwQJNLpeFLl04uX7585tzMzEy8c2pqqvIUIru7u9HcX3/9de6PXflaXzo7Pz8f75ybmzvXuTHGWFhYiOYuXcpvp8nJ7N/+N2/eRHPpNR9jjMPDw2huZ2fn3HdW7rnj4+NorvLzll6jiYmJeOfe3l40V7nn0sevvJ/pbOW1p/fx1tbWy5OTk6X//uelAL58+fK4f//+mXP37t2Ld169ejWaq7wpP/zwQzS3uroa70wf/+joKN6Z3hAPHz6Md37yySfR3GeffRbv/OKLL6K5GzduxDvT6/78+fNo7vvvv48f+8mTJ9HcN998E+9Mn+fa2lq8c2trK5p78OBBvPPmzZvR3PT0dLwzfe8r91z6wWxjYyPemb73lX/Qkg+kY4zx9ddf/+u0P/dfEABNBDBAEwEM0EQAAzQRwABNBDBAEwEM0KT0PeDZ2dnx8ccfnzm3vr4e70y/D7q/vx/vTL8Yvrm5Ge+8cuVKNFf5/mT6Je50rqLyXce0ZFApYqSPn76fF/F7rSs707JK5bvnaRGj8h35tLBy/fr1eOf7778fzaXfvx5jjOXl5Wgu/R7uGHmhqVLuqHzv/zQ+AQM0EcAATQQwQBMBDNBEAAM0EcAATQQwQBMBDNBEAAM0EcAATUpV5Ddv3kRHr9y9ezfemdYDHz9+HO+8fft2NPfy5ct4Z1oLrdRX03p1ei7aGBdzxl5aha7UxdO650WcYZa+nkpd+9WrV9Fc5Xy/tDY8Ozsb70xnK/dcenbeRdwf6a8IGGOMg4ODc50bI8+Fv+MTMEATAQzQRAADNBHAAE0EMEATAQzQRAADNBHAAE0EMECTUhNucnIyOkQzbcaMkTeY7ty5E+988eJFPJu6detWNFc5kDQ9cLJy2GXaxKs0eNJmUOWAxPQww7QRVWkApi2vynuUXstKG+winmf6s1lpmKWvaXFxMd75+vXraK5yYG16H1dagJVDeE/jEzBAEwEM0EQAAzQRwABNBDBAEwEM0EQAAzQRwABNBDBAEwEM0KRURZ6amoqqyJVDAtMqX+Wwy/Swzcrz3N7ejuYqh0Om9dlKNTKtZlYqnGnFuLIzPfAy3ZlWmyuPXdmZ1lwPDw/jncnP2hi1gz5/+umnaO7hw4fxzp9//jma+/LLL+Od6T1XqemnB6dWft7Sqvzf8QkYoIkABmgigAGaCGCAJgIYoIkABmgigAGaCGCAJgIYoEmpCXd0dDQ2NzfPnNvY2Ih3zs3NRXMLCwvxzpWVlWju0aNH8c60NVc5IDFteVUOu7yIgxxTaRusMps2ICtNyfQQycr7nrbrLuKA1dXV1XjngwcPornr16/HO5eXl6O5Z8+exTs//fTTaC5tC44xxtOnT6O59FqO8c9/jnwCBmgigAGaCGCAJgIYoIkABmgigAGaCGCAJgIYoIkABmgigAGalKrIExMTUU1vb28v3pnWPXd3d+Odac31999/j3emFc708M4x8ipyenjnGOd/2OUYeSW2cjhk+prSQw8rldD0INjK4Yzp66ncH1euXInm3n333Xhnan19PZ5Na/pLS0vxzvRXFKQH8I6R18CTX7fwX9Ln+Xd8AgZoIoABmghggCYCGKCJAAZoIoABmghggCYCGKCJAAZoUmrCTU5ORi2itL00Rt6yunXrVrzzl19+ieY++uijeGf6mo6OjuKdaXuqsjNt+1RaXunOynVPX1PaWqscdpmamJiIZ9P3s3LYZdqaq1zLX3/9NZq7d+9evDNtIW5tbcU704ZbpSV6eHgYzVUO+qw0Sk/9+//obwPwHxPAAE0EMEATAQzQRAADNBHAAE0EMEATAQzQRAADNBHAAE3K/c2kelepB6bVzD/++CPemdYDK4eHpgcPpoeMjpFXSCvPM634Vg6xTKW18jHy9ymtj1Zq0On7ntagx8jvucrz3N/fj+Yq98d7770Xze3s7MQ733777XN97DHGePXqVTRXqYun16hSL67Uq099rH/0twH4jwlggCYCGKCJAAZoIoABmghggCYCGKCJAAZoIoABmpQP5bxy5cq5PoGZmZlobnd3N96ZHtBYadctLy9Hc5U2WNoYnJubi3deRMsrbc2lbcEx8teetpIqrydt4R0cHMQ70+teacK99dZb0Vzleabv+7Vr1+Kd6eN/99138c70nr9z5068c3FxMZq7iIbs3/EJGKCJAAZoIoABmghggCYCGKCJAAZoIoABmghggCYCGKCJAAZoUqoi7+/vj0ePHp05t7KyEu9MD99bWFiId66trUVzH374YbwzrSdWDvRLD5ys1LDTCmfleab15ouoi6c110rFNz3IMT3gtPL4lUNbU1evXo1n09eUHpY7xhjffvttNPf555/HO9P6+8bGRrwzPeR0c3Mz3jk/Px/PnsYnYIAmAhigiQAGaCKAAZoIYIAmAhigiQAGaCKAAZoIYIAmpSbczMzM+OCDD86cW1paineurq5Gc5WD8tIWzZMnT+KdN27ciGfPW+Ug1LS1Vml5pc2xtNlXmU2bfek1HyNvAVZaa+ls5XlexLVMn2flWqatxqdPn8Y709d++/bteGeaIZUW4Pb2djx7Gp+AAZoIYIAmAhigiQAGaCKAAZoIYIAmAhigiQAGaCKAAZoIYIAmpSry5OTkuHbt2plzL168iHe+88470VzlwMf08L30scfIDw+dmZmJd6azlZ2VCul5q1Q40/pqetjlyclJ/NhpHbiys1IHTqX3ceWap1XktH4+xhiLi4vRXOWw3mfPnkVzr1+/jnemlf4///wz3jk1NRXPnsYnYIAmAhigiQAGaCKAAZoIYIAmAhigiQAGaCKAAZoIYIAmpSbc8fHx2NnZOXNufn4+3pkeyrm2thbvnJ2djeZ+++23eGfaokmbW5XZSgswbc1Vmk6pSiMrbRClc9PT0/Fjp88zPRiy8viVa5leo4ODg3jn+vp6NFdpX6bvU+VnOP15e/nyZbwzlR4EO0btHjmNT8AATQQwQBMBDNBEAAM0EcAATQQwQBMBDNBEAAM0EcAATQQwQJNSFfnw8DCq7969e7e0M5Ee/DfGGJubm9FcpTZ8EYcuTk5m//5VasMXcShn+viVCmdaF68cjHnej12pNycV/TFqB5dub2+f+87nz59Hc5WfjbTSX7mW6a8zqBygmR4Em9agx/jn96dPwABNBDBAEwEM0EQAAzQRwABNBDBAEwEM0EQAAzQRwABNSk246enpsbS0dOZc2goaY4ybN29Gc+nhjGPkBx9WDtRLm3CVZkzaNqo04dL36fj4ON6ZSptbY+StufS1V15P+r5X7o/0EMtKUzFt4qWNyjHGuH//fjS3srIS7/zxxx+jucePH8c7v/rqq2hub28v3rm1tRXNVa5R5Z4/jU/AAE0EMEATAQzQRAADNBHAAE0EMEATAQzQRAADNBHAAE0EMECTiUp1dmJiYm2M8a+LezoA/ye9f3Jy8j9+j0MpgAE4P/4LAqCJAAZoIoABmghggCYCGKCJAAZoIoABmghggCYCGKDJvwHTjUL61HaW0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_im(cf.view(1, 32, 32, 2).detach().numpy().sum(-1), cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "cf = cf.detach().numpy()[0]\n",
    "cd = cd.detach().numpy()[0]\n",
    "cov_mat = cf @ cf.T + np.diag(cd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADrCAYAAACICmHVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABACElEQVR4nO192Y7jSLaka4mlli6gG7jofpuP6P//h/6FeWwM6uUiMysipJDmoWBKk8nO8YWkQpElAwRJpNPXsx8nuToej+WOO+64faw/ugN33HFHG+7MescdnwR3Zr3jjk+CO7PecccnwZ1Z77jjk+DOrHfc8Umw7Sn88PBwfH5+XqovFppaWq1WF8eWxGq1mv260TpHkM3Vj5K2czRxPB6vOs9z4eXlpex2O9vxLmZ9fn4u//73v0//eYIwMTpxbsLcxGr59/f3i3NYAJyLsF6vy+FwOH1PwXbrp2i9XpfNZhOeG7lutVqV9dobOxnhYYxKoMfjsez3+5ApMY+6ZofD4WL+UeZ4PJ7NKc8xznM9GE+2Dtvt1q6ZHtO5YZrjvr6/v9u54P+bzeaCntbrdZXBXb1TaYzxn//8JzzXxawZlKlKqWtBJ/2y8qvV6makJYjH9YeJncEE75CNLWqrRiz7/d4ej/qBurROLqtMEzEkCyX03605zqEefLtjo+A2uS43B0y/DrXzS6GbWVVyAk76tZhZmaQuxU9IpJlcmZayozgcDtUFZdQWNyPI7NpsniNGjq7R46zlXH0t68XHe+YLtBFpVlzTO9ctdUb18nzMrVVrmKxZebAsPYEWM1gJocd8VS3eKu1aytbKZJo+MmlxjbM6avVl/cy0dXTOzbGWd1qUr8s0cCtqa+3am6ppR1GbjyUxmxlciic2ZiImfreomamXMaNKbW2T6+Fy7n8P02fM9fDwEF4TEdp6vU6JMAu2ZX5pD7PiOPuz7hpdy9VqVQ6HwxDhsvXDY4ra5t9TmJZN7Mg6cDQb0dfSGGJWdDIzIfSYMkA2UMdgfCwLMME8cf0ZIaRRYhgxg9WsdnPm6mQ/k/t7OBxSrZsxcY0I0RcmdCeIWghar+sNMDm0WEXoV2bGO5r9qCh6N7MeDoeTJOQgC5hXB8TMA4DAnLbNJgv1bLfbbi2odblzar63EqyiFvHNosjRuYxIMMe4Fn3fbDbl/f09FFKZht/tdvac1gdtCoZCPzNzvlWIon88dvyuCVENKL2/v58de39/v9CsWdAQbWvd10Q3s/Jg8NuZoVzeLUYkgV10mEPsWp9eo1FHJz17sNlszsaLumCyRn7p4+OjPY60jiPazWZzwXAtiEzTUkp5fX09ldPzzkLB8ch6AdFzXTznnAZx6TN3nOfKmcGq3bRd10fuH4QIl+XUDSsWjkRHZi/36+YDTCrp+LhqVjWp9LxOmGoQx4hqMvF/9SOdueYIiOvib65LFyZbqEgzQcpH50aldsTYUT8yUzfTxkDESCxs8IHmV6Gnx7O6M9QCmu545CpF0eCI3m/eZ41QC3bUzrWcd2DmgsQcMZN729xsNqHVkPmDkQnHLkZvXxyOxzy10OrLtvrOeq52jV6XBSBrLoy7Zgo+yi/N0M2sOrk1RP6q1pcRjhKLI0AcY1M5MsEzDRltCijlUjvzd3QNQwlaNcKoUHEmG7cTQctFWsX91/ZLGU9jRIweWUDZdT1tsnuVlWn9vzSGNkWw5gCR4bhuP3Pah01M9hfAYJlWyoIzHCCITOaR8UbEAH+Wx1+K91mZoJHWcT5rtE3RlQey7YZZICba3ZRtUczSM5mvC2QxDN2eyGutVhuvsbpUrm6MNTrPEW18o77IZ54iYEcwrFkZPFnqC7mFUT9V6+iB+saZRo3aZs3iglMMFUbaDyBrv5ZLVU2XmZBs6rprI9/0eDyGzJoxpApAvU6FKcco0DctAwGlG2xcu2zZOGbRtUS/+Bz3K8KI1bg0ZvFZe82BGkPMhWvsLqlFLLUsE5mOXaPLLVLbBXBwnDWDi2yyJVBKn8CM/EtnjtfqjYSkWglz4hZ90hpmYdalgjgtPpRrO/I7eyKAI4vpGELPZ/VPJaBoHaJIZinxRooaWgJJLLQyGlGBlwWwsmM9NOiE5TVN2hHMpln1/5IMPIdUbOlfrZ2MGKP8a7alMMrz4X8El/t0fdLfzmxFNH0UkWbNouCjtDIHjUXWnXNFHHPffIAp+u/8JWeiqYZxGgj18XmYkAhaudxptCkC/dG+RmjJs0b7blerVZpnzXbsZMxaCzA5oB+O6HT/L69bBF47ZwZHWlI3e0Qxi8wM1qi92+4aaeeIxlqsNzf2qMySmC3ABEmt0VyN+qI816fnognQgEGUcpnDV62lbrJyKKuoSWKNNHLZLCilASa+ruXmc3c8m79sK6Wud4v2QVQ901ycQ8c5Z020aFqdyygI6n5HZa6BIc3qJjYySTKTuMaUer0GRkrxO5hKuSSoOSbWSfYeM6ymubg+7W+rtusZZ9T3LH2mGk77GG3Cz9AqWFm4uD62jp1pZpQulnLzMsy6KUIJxkUE8VvriXKs7vqadpszClzbMB4hS31Ei5ylTCLBkDFPKecEzvN+PPq9wbX+oy+AM9dxjPPQGZxgjUxsRzstTKPxhBYlEtV7bY0KTH5ShNNqtclzTFzz1fi424TN5SICqZmgrT5iC1oin67+Wg62F7qBpda/Us7vSNFrM4GS5Z91wwP/HtVQmXbM6sw0q/PHI2vn2hi6RS4jKGXEzKTi+qLAQvTfEV4LMWrQIro+q2MOjCz8qLUQEV5mBustZagHqR1XV4+5G/12/eb/vGat8xHFGnjLKIQQ03ZrnOVaGLpFjqFSqsfX4vrYFI4QSTgXjVRkPl1N40bIgmHR7qBS8udCjQqI3iBSra5s3SJhifuMe9EzZjVnHVqCgDWfuiVKXJunuTHLkyJYM+Ecm7cZg7CEzvwyF2ksJdbMLhDmBEnNJ80ISXOmHFzRx7pwPdljSkd85JrQqJmtLkWSCZssipz1MZrLLLrccn1Le+561awR/am1qALjWgw7aVNEZPvXwt6RGVWDCyyMBoBKqQdKWrSEWyi3C4j96ZHIbQZHSLXIs7s++j8HWtbXMVh2bSZs53RhPtpXBYZ81uz/EuB0TE0bursnSpk/qlvKnxsOooX8448/utvKNGQpYwQ4cg00fKSdIk3IT9XoQY0ZkGONBGML9PopdX0UZr35vIaapFQ4BquZikxkEbH1mlTZHTQ9x2v11c7VoLeZlVI34x1qAbdaWieK/EZQS8RZXDW/UvvZ4pP2CpaRa+bEpNRNjRFqxNzCOLpwU0xm/J9rwrMAQ8vGgdZr0FaEbDzsH+v2TK5b68juq9V+acyBYxm9G0e43toxl74DWtvMBHrvtUtjtrtuar8ZysyOwSOp3suovdf2INO40cPIwKxO42R3vESCIbqvU8tofzONl23m702r8bEaAyizufJZ0BC/nZmufXKMpm6IKogeRgbmDEDNfj9rJt10kE5TRpPK35lEU63t+joXooVwfh0i2VFOz+X4VOJHRJKZmZmGbInC1qwkRUvkfzRgEwldpT34uFm/oqhv1p5+jwYgRzH8dEMXtYNGcekMhqZaeLJB1HyXBr5HNOSIj9qKSBMej8fy9vZmrzke4yc0ZE/PzwgjY9bHx8c0mtpSHwslvPGNr+cUSItL45iN/UtnCUzVUCwoMzjTusWdi1Brs2dck/YG6zE0zv+zfZ+ZdOO9wsrYGeE6s2ip0HumWTONFi3uw8PDELNmZnAWoc3uq3X3uQJRRmBUIOJ6/XZlGJkg1tRV5Ea0uFvOPG9VHHPdxFDKALOyVuAFfH9/L9vtNpSMwPF4vLgflW+RwkQg6usCTPpQaD2f+RoaeMqwWq3Ky8tLOhfRZEOzakTycDikWjfaWJD5s5vN5sJaQf+fn59DguL3ogKY90hL1rbtRW5MRtTZBgwOVuncoG7W8MfjsTw9PV3cP8vjxDxvt9uy3+9P9IT1zOiX30bBtJv1P7Ngsoi1YpHnBtfOa8cis5rRKgHVNHHfUQTQtTmiLVSKM6PVAkxMUNy23mPK43DBD1wLt0T9LLZgtA+80aXXZ43WoJcwW6D01LNWPJ5oLaJrXBBqFD35/y5m7YnmZdAFixhxFHP6pyMmNIizlEvNyprLXVer1yGrTyU/fmf5ambuyL9ksNWg5fVYS86aBYRLNfWiNq9RPj/qq1owUwRQSx4a6GLWzDfAb7eYKpHUN6lJOdePSLNmWqYWpHJCI1to1oSqqSNTd7/fl7e3NyuRD4fvr1qE6YYxZJpVg3pcRu8pdUERt5GCTVstE0WxYY4z+KXWWTQ4Wg/df15DFPCr+bXqqkTXOetw6kaWVsyiWZUZssALg8uxbzEXoonG71azvaUNbcvdgXI8Hk9+UlSX+uB8LiIKvOjKSXzWoCrAHPPpOdbEuD7SrDqn6DNfl+0o4vZwrFezOgEbMTqPLTvvoNbSnCmaCJN91sw0ixYnM4NbJGgtl4dvF2Cq1dPSFoBXKgKshaBRtG0ENvgangsmeJ4L9ywl1O2YFddyqoUBrRuNKyL6bD6YMfHNmgjMp2kRt1ZTMBr9z1JH6G9moYxgMTPYIZpcNil08rIAU6s2azGD+Ty0BWtwJa6sHQfNi7KGUe3JzBppVkfMOJ6ZwdGCg5Gjc5kWZ7OV+59Bd1OxZo2Ed22OM0SbH6bUx98uEDki8DMsZgZnfp6G76F1Ms3qGApQAlSN6QiaJTWIns9rGRfUUNQ0q55H/Xo/K4+r5gdrXaXE75FZr7+/7xX95fl8enoK24tywY7BWxjAaSK3SWYqgQNuL3IpueXFfUMdXJ9qVnXVuI3aJpAWLJpnLcX7OMw8rL14MnnQUW5OJZke0+POD400gjLKVKKJtAI0mkO0OLW0BgezFNnN7MzIPdfpb/QvynXid21DQ+SOuDXlutVVyOZ+ChwD8hzMST+9GAowYeLYL9MyWUAB57A5IvIFauYuH2eTUM9HZnEkbVtNqShIdThcPmmBNSQ/AFz7oT4wm2Y1ZlXNwGPt8bXclkL8zgiUhTBfEwW5pkLpsYdRXVQX/7Ocq9L3yPt0FYtuiqhpVnyzlFJzIkKkWdUM5jJ6nsHla/5GFIhq6Sv/j+poEQJgPC2bPaw70p7Q8NG57PWZ/DvLiTMyzcpCs9cXzuDM4NbrSokfxs5wlmF07QhaTeEhn7WmWbEgUWQtW6BIs3L7zgTi8sp0el2NaZw2j/rqrArOl+q5Hs3KyDSr87NQt25waI1iRhHRyBrh9vg/fFbkjHW8bk1c3zJhjzGPMH6WRorKcFzk5jVrC6b6DoALErnAUKRRHfPy+azPvYvfIq2dT8T/NeXTAo7COnchignUxtF7ThG5R1GZUYxoVr22pZxaGDetWVsbbh1Ey8L1alanOdx1U4kW4IlmbRSZwW7XkwbBcKzFDIb2irRRtjc40nQapNMxZFAtyYIksnj0dy/m1qzufBZ3mYqralZ0XM1k/R059U4C60uIovJMXFgwJXglPNeettOy8JG5dDwe7b2p6Jve3cHXj9x1k/mlbtsg/0e/dOwR86zX57fPYd6zfpSSP31CtyKyYFHh6uIBeu8zAGZCHzVlyPSo6aUMNTO5B4vlWUfRQvhuESJNWDvuzrtySpRKJE6DMyItWANL62z7nTJQpFmdtsyQpcv4f0ZITOiOuXQcLoaB3xqY0p1gylROGOtcOPB1btwj6F37KVjkZcpR9KzXTMkmIjIzM9O6NrFRnRFUUisibaVBmCiP12ImZgKK+6ltcbnIusjQwtAtcBtmasjMxjkY51p7fXsxy3ZDRc3M7SUId9yZRey36PlM02bt1zRrFOxyggPldXcVX8N33fD57JEvGLtexxrNlR8N+qiP3ROk0d/KGLXoc+ZW9Qa/onFGjMrXfAQzz7KRv0dj1jRly+JHkV8wgTJL5KNG0eA5IpSoj0013gEUmZ6jmtWVrxFvTXC61Jua0LUtd5mvyXVEPrwLWPE5thZ6aHFEA+u4r82wV33IdwummFXsizos6V/UfOLIn4z2nro6s7bdtXwsysPW6ta+Aq1M4SLO2g8+pu249uZiklGGnUuY9+JmfNbeSWj1WXsw6rNGqPmsESNxgItNfLfHFpolYzzdSVZjUvVtuY8YC4A+txC+Bohu0S8s5Qf2WRmjWisyCXvbajH9roGob6V8vyHcaV0OWDHDlnKeItE6uW6NqvK+3Kx/tbFEGhJwfjqXzxjU7c/lerOyS6GVYa9Jb7M+kT9yvlsCFqqFWpFpBxdscsezUP6UABMD/2v3s04JMGngCn48m5f8m9dN519N0loswZnifDy71gkZtKkCQsvObUmh3ai8W9drMezNmMFLYOlJHEkhgHAjn81FdfG/ZgYrYTuzm39H5rAL2rQE/SLm1/FlmlT/16LBq9X3m++1bK3Pf0mftWaqfMTgQDy6IFk0WMvVFpM1a9YO/2emizaIR9HurM5sTDXNqkAfe7bwZcwPaLomSt+0alYnjNSfzrYD9tLlRzJqKTNrVhcw4eNZsMTVN0efGK0mSytDR5p1tVpdvA+U28XD1ByDZ28VjzQFCFKZAP3QlAP3yTGKBq20j6XkgkZNcVeej6n5rtrTpZBwLnOfWjVrFEyL1tZZL6PoyVPPuimixWdVcw7HenxWleI1n1ODOVFZvkbbcYh8VuzxdRbH8Xg8bcrXMY/6rDBTUYb7H+3XhdbUetW0jOAIOvNbp5iktWgz+uBeBhZpVrZUWvcG6zU4NoVhF90brI8CUemkktwRWCbpnXRVvLy8XDBhLVrJbUTCQY+BaCM8PT2F7UXPYHLCSs/3HC+l7h9n5x3YB65BUzrRPudszK3vgtVyWCMcx9pyW85SWa1Wp/XBtWBYHQ/TC9NNK5Nl4+6NbM++KWJOk/YWc12MjKBbCbC3PYeagJp7HiNCHW1n6v5i149aEKslNdNjjV0DN7eD6Q6PJUytUeaam/k/QijfuiJwmJ1ZNeLoEIXwPxtqecfo2DUl9Ij5rBHVUm5zV88UIRSNR+k382evPSeLataMaW9t4a+BlrD/KCOP+LoZambkj4osah3tZb4WFnmsS8uxHwFZxHSUSebO431kXnBpTKGrlmtvTWANP4o06njr9jSU+cxoSW043IoZXMshO9zSmi1pBkd0/KnMYLzWjzedsy+G4/qqP2C1+v4Ga51s3ta22+1OeS32cdfr9dneWs6xvb+/n71FXVM5WhegC8Pf/PoJLcdv2NYx/vTTTxdzh3l7enqy16zX64sNDtpHhygHezwey8vLi90GeDweT49E1Q0M+/2+vL6+XvQR5zSfy9dFmzoyYI5dKo7HwhtNMJdRum2/359F5JXpOD2F6DFom3OzmqfVNpmel0Y3s+qGAU1fKEG78DeYxu3g4Q3pqAffuI6JGPWAgJRJ3UYN3czACW5uR3fj6Bg3m83Za+uBzWZTHh8fLxYV556fn21OeLvdnvJ/PTtl8NKqaKNFxKxZnWBk7aPOLTOZvjyLy0eW1uFwKI+Pjxdtudym1tuTJnSKJZsDp1ld+meqm9GzqWLeBJdBa3L9WljKh2u9DS2aD8fYeo371Nq8pbkHbsmUZtziXDEWz7NGe0XnQq8wmJPAI7PaaWTWFE4zcTn3u9YHro/7oBZGS31LQ83RFkRm/pwMNnd9c2Px1I3bvsXQY1N2wUQmGPvSfE4xuvsq2mKnZizv4XX95XNcL9yJyCTlj2PWyKeKGED73gsVzGo1tG4NdXVyP5dmLHWdlkDPGBY3g//KUK3fu9fW1ddaj4sjtOBWNcuPnIJqxaypGw4o8W1a+HbE5oI3bNbhO4oWcgBJzVJE+jRoxXVrwIH71bJ31EWDNaLI87Hdbm00GOcQbOmBBtdQNyK+rClZG2fptxrTRprN7f7hIOFcmDMW4rIZpZw/F9plDOZAj4UwFA3ebrdn4fuHh4dTaP3x8fFEIKWU8vj4eBowzMDdbndxg3Epf4b+Udcff/xxOo+UzNvbW9lut+Xt7e00aTi22WxO16MuNgHVHHQmGgsApIIeHh7O7upAufV6XX766adT1JcZYrPZlH/9619n8wYGeHh4KH//+9/PhAHafX5+Lj///LP1713kEXXsdruLdEopfzLx77//fkp1KTN/+fLFzsGXL1/OBCVf59pC2ZeXl/L29nYxT5h7fekWrttut2W73Z7WFG3iGNYC9wHznKibBZpA6gZriTnAWm6327Lf7y9e08LvuoVAQ5nj8XhWL9Owxix4fJFCUiFaw2Sf9RZMrN7gifpNkYZoqTdLBbjb/zi/p1aH5g5bmJXPuWQ+tL87H5nVUTQ6QhRM4/98nPPw6he23FeqiIKXiAG0XNt6HnPxEcGoT+mzOrNqxKe5dT9orv4tPc5eoo0EvB5vUQRZ20sFhT6Kbj4ls2IRVXP1YgnJ6BayhWg0qqspF/ddS8vw8ZbIZkaEtQhxNpe1eZ77ftYl8ZECfvE8ay9DRCF9twPJ5Tb5P75d2J/PqY/RY+JowKuUy61w2ALpTF2uBzuiXOpDfSD1rwEX2OPUD37zO12jIJsbZ4RozVydztzt1aJLmaDRHKjL0uN6zdXXbmZlYnCdiCZU7fwewqgNuEWzqDPPwYyo7ihvC6hm1/Hp+N0nO4/2wXAarUYZjlbimwNDTjC4raA6Fz1Eiet4XhnqM6OszuMUzMUUcwivJTAUDebFUALjTf6lfA9+8Le7XtsAEGnja/gYonPv7++nvamIIEODoD+qiVm7IBihGjYjVqRoWGuV8n2Pr9s4/vT0dIr46iaN5+fn8uuvvzYxK44fDofy+vp6Gjsz6/v7e3l+fj7bs4vz7+/v5e3t7Wy+MV7uu2pqx1i8Ljyv+M3RYF4HtKc3ULRChZ6O3/VRBRj/1tQjBCGi0dpelOJq1aYcbW5JBw1pVg1388Z2MA6nKtwklHIekeTrlSnx2e/3ZbVaXaRuwLA4X8r3jei73e5E1PzGcZ4YCJD9fn86h/qiXU9gPN6Uj8jj09PTKT3DbaxWq/Lzzz+Xf/7znxeRytVqVX755Zfy22+/nQiCCcQ9pQ/j+vbt20lYcZQVT//nNBbmV9su5TvzcAqMBQbGofSAMvv9/pS6Qd/AoEiTqFUTWSrcpprwOIaUHY6B7nh8YAik4QD0a7PZnN19xGkZThm5VA6YzIEFUWSFctkWzPLKx+w/H89MiqzDqhVxzGmbDC4Xxv+5TE2rct8Bl6rR/2yCsv+DY3yOtYAjCr4W88EWjBKKaoJo7dwa9pjEpfS7JnNhinkajbsFc26UiDD7y5SjxahFGrP/mvzGsVLaJxMEr5oCv535U6tbx8eMrlYD2lETHUwAbQDrQM3gKOiB+vDh/rBb4AJyOhesDXqYSk1KN+6au1M73rrOU5jfXRvVF+WKR9FSx+eJmVeQmVBaLtKGU9Ar0dVv0uMZA2REXutHxAjRp2VsKpDdd4vgzqyaGiPNraF7MAcNtdTxwzCrSqaIcFRrtO5UqqGFuPR8RIBK2Fo2MzFb2s3a1E90jRMmkYCJ/uM3LCdOa7UIqBYB5hAxxqgJ3RsYi+qoYdE86xT/QdFjarRuQnA+q5bJUPPx2LxUoo5SMeqzluLdAFyjfm0UqXd95utrfXdzE5nnmWZls5/LfDZ8xEaOIWZVwuIAB4ehV6vvUUxN3agPCWjondM0iOixT8i+IRO3A0ehmaA5WMUbrGuEtF6vz57zw6mb5+fns7Egavn8/Fx++eWXk8/Kc/rrr7+Wv/3tb2f7d3FOI/AgfPhOiFzzpvPdbndK3Wgk8/39vXz79u1izg+HPx+zgnVTgYQ2OdKLcSDyzOVQJyLSWEeOBqsPj/rUn29ZE6YFpku+NhPmHLtQIcvrzHuvlY7Zn23FIqkbdE7/g6BAlOis3qGgpg6bEFy3BmzwzYGTUspZkIbTPAi5c36PiYDTIpgoDf27B5GxNtlsNqd0AK5Zr9fl8fGx/Pzzz5ZZf/311/Lbb7+dUkWob71el7/97W/lH//4x2l+uL/64Dkex9PT0+lhZRz02O125bfffjvdBcXMut/vT3c2YV3AXF+/fj17B4wSH+dTMber1cqmbjjtxgKG24sCbrzm+K8b8yPLQedK15IZTfPhzKRYO4yby3PbKmD1t0OvC3aV12fMaQ5HaJG62bU1kziDLhiOKRHVgjjOZ9MPm5Jsrah5rXVxP52miPqt56cgM51bEQUSr4FrtuWw+Iup3Hn1yVoQBXAyJnVlnK/E0rwlkANEO1gi1JjVleFjXI/+1vKRr5mdax1Dz/E75kM3s7YwhjvH5smIFnRBkFLKhSYCdLdLxAwR87g2FWzGsNnuxslmOvuQXJ5zsM4sdGYw1+nytuyDs38fmWuuHTaDM6FZgwqmGlydLuo/p6DIxpEJ8tGIcE/fu2u/NckaTW5P9HjJiGSk3TLCdeciLezaueYaZX2a2qa7dmkzuLW+pVyFDD+Mz1rKfEw3h8+q52ua3PmqHP2NzGX1V9UHdnXxN/fb7brKxuWQ+f9z+Ky9/ZkTS9DxogEmRyQ9DfP1mengzGY25xAVhOmHaChvseON+PgNgtFosKYHNBocEXGkJTlqiDZxDOktDRZxNJ19y9Xq+7OB9Bp9NhGnXJhxVAhEhKdmLebPRex5/ji6quciQYFv3MBRSjk9d4lNbxVIOIa54nkGXfFHaUojukwXmWXQq01b0jE9AmDS6zOATGLWjuG6bFBuryrnWTmfF+2HdT4k1695VkdkNTjfnM85f9CV4ToUzq+OxuYEnfto+9F8A0ykSAWtVqszQcOMwYyg51EH6nFpEdTDTOZiFY7ZSvEuUeZ+cJvO2lD65VTUKP5Se4NvBa0Cq1ZHFr3uiVhnbcyBj/IZl4RaFbeCO7POjEhitwR9nIbIjun5rM6o3anMMYXpRwTbNRioFoP4KHx4b0YXu0akI3Xh9xQCjggwGqfTkplJm12T1Rm1O4dmzsaXAfN8axpMNeuS2YIezBINrhF3dn40OpgR80hdka/Zg4zpnc8WnVe/jY9pe611ap8iYZcFn2pw/niroHL/PwqsWeGzTvVJ58Aws7qQfFa2Fn3MEBGqfkdE7cq5OvR3KxxRuuQ9Ryg1aMMBMswXR2J5wwN/a+QTUXI9zkE57a8KqlGf2EVSnRDh86MafknG1vmPAqC9e3unoptZlQDYZGDpo9FaMFH29IJSLjUbMyHfMIBjSF3w6xfwn49xisMxuEYvOcoZLQjSOxgvyr2/v5+9IgT1rVar8vr6Wr59+3bqE49zs9mUp6cne9eNvjoCxw+HQ/n69evFRv7j8c+7bvCMKggBnOc7ZLifx+P5jRIKXTN+UJuL+mbzrMygH715gGkG/eT14fG4D9fFqT+eS92VpjTuBNocvu0id90ossilk9AsYSOzSaGaL9qsrovCBKLEE905wW25RchSGjxud8cOGBQErswKBkKfOBesmpClP67jNnGc88WOyXRsTls4y0iZSeeSj/G5SNtmGsr1R5WFIupHK3gtozmZ4jKNYGhvsOsgFi8alBJcy+LUdgS5/8x4fD4q6/wT/K+9JwV9dWPLTCfcSseagW/TY2blsUSalZ/BxP2qPYPJmbyRGewET4TMV9Zy7vcU9DCPK5vR6Ef71MM7mBSsqbRs5C/WpGJmijkiw2++VpPabiFA9Eysrf5INDb1n9kcq+1g0rmEiezahumLvrIZzNYHn2cByv3k49xGq2aN5pjnVf3UFuHdA/a/W8tGx6Zq57kxa+pGB9Nj2jhMuYuBd7ro+RZJ33v7Wy8ignbn+BpXPqo/0hy3CmeRjPb3oxlrCcy6kT8jkJ7JG5GKI4sTMe0InObQfalMjKztuA4+pxqQywEadGn5cH+5Prd+I8ziLCr37X5D03NwrQdThdFUBbMkZvdZuZyaPPw9BZG5qVFH/qimVIJis9Sdd3DEjw/8Uh17FmDa7XanKLL6rG7eYb7vdruzABM+UYCJ/Vk3Ht7YruOsrUsUUHJr4ua2l1lG6Slzc3oFGurrRW/utptZ+ZUM8Af5fSavr69lt9ud5QjZbwRRIszORI26UX61Wp3djYEHeeEhYOzjIVXz9PR0dg0TPpvEmHANMAG4E0Sv4c/Ly8uZUEKdu92u/P7771Y4PD8/nyK+qnmfn59Pr88o5Tx1w4KP/cjD4VDe3t7Cd93897//PWNirvt///d/z+YAZb58+XJ6pQSvI9ZI5xDlsEY8n/wKCtaaOt88rxxv4HllKwJCSvsA2tG31fNrNjAOzBGw3+/Pxq1pHBzD86T0egA0dzwem1y5VuthKM/KE+QmWaFaFsf0POpXRDt6+LemcZxk5/q4HVzLRMTXMhPouJiY+IFreLcLzrHQiJgV+VkNkkRBE9asms9mLc6alXO4amrqeBi9mkP9f7TNQhvHWZvr/Gb9iGiR23emN8pG6+rgUkUj2pTRa+Jf5ebzUnzOdQqiRdDz/N8tzFyBJJX+eoyFUeQSKDFwORV2XMYl69UKcESemYE9hMxYcvP7VOZoASsOPnYLmJ1Zaz6Hm4iIEEuZfgNvLeqrbUVM7dBC8PwfgAbWcbEZh/8slCLNyiaatq++Kv+OmJ+FzBREcQLdyRQJ3owutM9zIrMSuUyrkJ9DGZRyxacbZhHA7LoW0zoLPHA59b9QP5uGbLLVEAVJItMb3/DfUBZ9hA+udavPqu1wYI/zqCij4+PrXf2RH58hmi8XpIEQ0jWpbULRNjQfPYdWdGPvEfgOPYyd4Wq3yM1tBo/CLXjEdKPQ8WW+VctxPuc+UduRf9XqH/bgGiZqS5tOoE3FpzeDM43mUgKsQaJzDpF0j9IBfF61k5pfri3WRDXtWjMldS8u7w1GOe5v9spHtiJ4LhHM0tdLIMDEUU0+17LfNdou2YMWM7g1wJRZWXNBx565OTVBENH8KIbyrLphHBMNYuRFR5qHTTF+xi3f9YDzADMNXtMAM9KdR8ge5zXvCyjhqxnlhIvbtsiRW27jcDiUL1++WEJ9e3s7/ddo8Ldv385eaeHm3h3Th8HhONIzzKx8/uvXr2dzjf68vr5WidWBH3yG+QMzugfBaZ96XCGmH/XV1TVQoV3K9zt09FrUy5FqjqZrfMCZ7pHL59ATkLvqkyLmMpWmbEOsHZvSxhS0BlIiM7m1DofWsTqB2tpmtvaq8fG7VvetmKfXwvBb5HrLsXlZi97W/Fv1v6K7UhhqQmowBWVciqVlfO5cZJa7wI4zv11ATNuOAlk8vkzSZ/nnaKw92wEjvzoKMmqeO0PUpxZk7tWtYkh9XFOi9ZgRPYs7cv4aWFqjK4G747V5GNGsGUY1a9SnJcrfAhbdFBFpyB/FfGkdh2p0Dao5YmVE/qPWqWayixi7tjjIE42LtWHvRnstU/uPY7cgPG8Jw8yqAQ0gI+DDof6S2RYG6DV5SunTHGijR6i4AI6eVxNcv12QCP2NzEl8XMTXMTBfC7TO5YifzIGmubGE0J8zejs3hqLBbqcMfEd97QRHKpmwdAO/hvA5vM/PXNIPJDBHgxExZh85Y6YoGqz+oBIHIop6fr1eX0SDgYeHh7O9wdze4+Nj+fr1a8jwylS41m3UR/kvX75caHTU/fr6elFXKaW8vLxcREvZF3bpl1LK2VvgNZqeMQHXo/46EO3N1thCpPH5pn5ci/JqbajFgnNMu61muhuz9rFVOAzddaNpB3R+t9udhdJ5gK7DOAcmPxwOFw/x0geF6ZMUmGn5nH5nxBLtYNJ0DcrgGtxhBKgmjNrit6UzUXPqQ5ERR5Z/1LtnWq7b7/fl9fX1YhwsYJyl8vDwcJG64Tqi+eenW+hNG6WcryFbVfrNZXkej8fjxV03zJQAaJv7j7KZO6HoYUIXYY/w4Q/5BnrM6d40DJ9nrRyBN0fMaRJd27xywaTa2PXaXukPZMSs5Zx5H2HpwNAtx1OGfFaWOHMQdCadsv/KUBoAqfm2Gg1FXa0aLNN2kX+nfWKiVosla9u148x8va4lksu+s/avFapZayZwrU+K1v5Eggnxk6yemkXAZa6BxTTrEoPoDVKMSOG5+t1Sj9M+mkNW5nKBI62rVav1jqHXR2vVynrNSIQ5QstmjM+C4U0RtQhfpqFc2Rrg37J/U8rlpohMOylc6qJmHrfWzUTH0hvHnS8HvxXg37W2nQZzmjy6FuXRFtep5TI4TT9nhHUkExDV0yNEsjLXYvqb2RQxupi6mR/Heq5H+7Xr5vSXbkWqO03Iv+fyWbPxKoN/pM96y+jWrJx/VL9QHzSNMpFfwvXUtEYp3yN8mrzXBeb3xfB3dD8k1xX9VyByHRExR4m1zUgoKLO0arWsr9E9ojXzMDuP5yxpXRrJdj60Q5Z71wi9btBHpFfb435FFpM+BAC0XVt7R+PXwNBDvjlHqj6WSmI3UZhgXghMNCYC5xlqArUEkRju6Yv4nSFqw70ig8/1IhNaGZNPIRZlhlJyvxEPrdNrai4EC1iXA0W7YBb8xnXK0OzysBukATH8x/VuY07E5Aql1RpaaLOHfmfZbqiE7wZSG+Ccfk0G59c51Igvu7ZFOjOYMHuioYxewaXX8nfNn9Mctv5Wf7xmIWhus2Z6Iy9/q5jLr1bM/n7WlhD9SJ60lNvzV6LgQmtEU6/BYzuj89l/aAunNXpR67vTxqrhdOcT96eWrlGt6NrvRYvwnQtL0emQGcySlRcl2upVqy+DW1w1hd01rn6O+CqhObQs8Ig1oHWyZo1MdTCzI+IpWsb59q3957y0iwf0RM7xjX7UNHtEa63rESkMHr+O4RqWX4bhvcE8EBxT4lZGcT5iLenOxyJprkSRmU8cMR6JHjMyc5d9Lz3OhMZlsveicmAP9QBzS/LMDHYbQwCXNhpJ/2i90fkoXeXgNLyrU2m6p0/c1s2bwW43UW8dHwHVLHPliCNfTMtEKRO+JmOglsAOl51KSDUzmNEyj04IzW0G92IKPd6MGQyN4aLAGinjDfoKlMf1h8Ofr1/gaDDqAHjPLq7DljEXWVYGcW8P53FFTBtNPh5U5vDy8mKPI+WDOeA+8tvI3XxFmioTBNGNAaXEKZPtdnu62cCd0wAT+sY3UkSaS/vKdMCCS8cDWmLr7Hj8/voKTgPy9U7olVIunqvE9Ksuib63iP87AbWUbzzp9RlMzDDT8A4aDJ7fH8oTBwaNTBgmThclxAft8kJn6YhIYjpzignJIRJEh8Ph9GA0N3f80DTGfr+/uI59wix4FzErv3tGkTGyu92thozZSvnOYGzSOyHP5fW/oxcnMDVYF43BtcH7BqK4hhMwDllchJVPizaeTQT0mI4oz98OWSrDTXKG3pTIaBQ1G1ckOFTTROf1XK2euRBFvFuOtSATNj3/GSPrp3CBwF7MaRLP9lgXSL1IAyjBqYR16DEnahPZOtGjPjeDmZLBFgBrFp2bqF/ZOW2PLY/adVFbWfBGLSoen/5WsDmM31qHBu9U+9aYcU5hxf3+yDqHmNXtAuFzjrBc4CarB2ghNCYsJ31bTbipYL9J24yijOobR5I4MiuV2VUIcH0tEVEcd8zjrokCS9l8zmkBRO6La7PXappDK7ak/lr7M8Ss2eI4SR6Zcy5gwqhFGMHsYA74yDhXynd/xgU+dAEjLZIhMncjEwo+vIMLZrEVkm0wiAjU7aXF72hOnGZ1/9mC6NEQqlm1XTcetcJcPzgGkrlJrQJC/U2+rtXPrOXqF9esWaOucyqp8a0BJoUGmIDMp+PAkgaY+Jvrcr/xv8cXxvUcVdR6mCGZAWqaVc1CHKtp1ugl0vw76ie3kRH/qH/P3zUzmM9pm1G2YSqUlmttOswZGR42gzNfKNJSIEjVLM4sZLj/mARN3QAcaUOd/HJj1BMJmda9ttF5F6E8Ho+n99lw+/jOUjduHnBMtSHPwXa7DfsY3SwRaTAw0BTz0PUTQor74zbdMANz2k9pTWlT23RmPveFj2tqp2X7ZO98LGoGo5HWc0yUUV4q26rIZhPgFtLtONEFYmLQhcfvHmnoynL+1x3Xu0w00MKaDn13/jj3QTUOM7AzV3HOmeRRWzWidNrQleExjWpo0Iwz2XvrcnAWV69v3oKrmMGRBnR7RBUjA3Q+DBZLTV+nbbIJ1z3Oo6hZBNw/ZUiMwzFK5pfiuCvDdbpzo4Tm/E0cj4hPTVtnUqrQVMZmLcrCmvsVCbapzIs2W+lDd3fNgVk1qzNHndnrNE4NrZpVNZQjDsUUzcrXKOG4BVO/lAkfmyx6NauWxX8QOKBCyfmsLcj8fw5aOYbmvmm/eCODc0F0P7dqVjcP3OZUhh21uObaKzz76zNaBhRpnLnhiE81aLQDqQcur6ompi4elwMhuYgvo8UM1vKqWaPfcyELULHmc0wcBXCcFafnW/v1EfhwzRqhRTMtMXE1govM4EiCj0xwzf9Sza9aIAveKFGzee80DK7JIppLEnBkersyU+B81UhwTR3viNVVyg1r1lL6I2RzDMS16YIYaG+ukHrm67pxqXleSrnQPC2IzHunsVyEtbZ7DHX1rCNraxUo2lfV8Gricr8jczvagLMUltqg34puZtW3jsPXgBR/fHy8kKq424UZB2F6vpZNxyi6ezgcbHqjRhRokyPAUaqGiRp3mTjGe3p6unh4GPD4+HjBJFqP5oLdA+e472xW95iBkckZ/cc6OD8Q/XCMCH88EpD4VjeAo+N8zj05g5n64eHh4k6cUr7fvcT94Ne+YE33+33ZbDansvxsKSeo0J8WAdGS+mtNDwLDO5h0QaIdO25TQqvErgWi1B9y1zHBqQDQ+pigWvqYnVf/E7+dGc5zFNW53W5tmoUDa45JdMyRrx6ZjrUyrj/8XZv33u14bo7cmJwg0d1ctXaXiqWMYtKmiAzsT/HGBD6fXcvQCKBDNPmZzzRSn/YzGge0jM5TpNFLyTUr2nP/IyGmfVytVmf3f7LgVMGbrU8UnHKMhP+8eUW1tprETlCq1o0QRbbdZoaWem4Js+wNdgwRaVb97RCdzxgrMn0ZKDOXxOzRrHw8KlfTrM6iKaXd4nCIzHNNweA3vqOoL/cBY1LN6kxblNX8Ka5VlyEbvwL1Rs8XbrEqbgGzaFZejMh8jHzKmgRnKFFFEx4xrpPObEJFvuUIWiSzEnL2HOKa76P/exLyUZmoLy54pGgVzlGkOkNmUbS081mxSDS4lEszOEJNgvVGJK+JTHP1El7tOt694861HIvaV8bTtiK3RM3WTLM6oahgq8FZENr+KPT6KJNwa5glwASNpdvl1uv1WaQPyHwbrle/M3NPER1rDTC1oGa2RoIqG7v2l32+qK3a6yeiaLC2kfXR9ctZSNE53XHG/YsCZzWs1+uzV2Ao7fAYMqWRjX/EyloqKDXLo0hL6ZfuyogOEcG0Sml3LPJbsoBLhprZ2rNwTlOPSPgel6NXSyFwNidc+oXbw7djev527k8UIXdoiQVw3XOgZy4nBZhckn2OgXCgYW5kxOuCGjWMaLssABYR0uiusIzZszmO2hv1/SJrSZk02kjSEpQcNZszpTCFDuem4UkBJtWsTrI5n0hNIQbqyF5MxSaKCxpFeTxnfmbM0eJP82+uryUFo/OSBZiy1MXcRJExZGS26nUsxPl/ZPLjmFvLnnXguvh89hgi9dd1DK69FqXUwuxX2RTB35EWijYCYMHZ/NFvDeGztu15XUTkl7Zi5Lm77jpe3IgZai5BL7OyUIyYJ8KIhtdXbaIethgyi8URbs0NiOIGzFAu0Ol8UsyxzpVew9diTA5TIvEOk6PB0a4RPhaZNvrdkgaYAib4aHNCT+qmxSeP+hFdo8+qZYEXEUfWVmbFZIi2RtbMe6c5NRjpUnLKAECLUKkFw7Q/PXVw+9kTTq6BbgekZWLmqm+0XpXmfNxpWtU+116EDK4vYFyOwPdcP9Kuc0m4L6MA8/KHjztkrs3cmDK2uXO7Q9FgRmY2lJL7IUtNeq2uuZhSGSWT4qppRtpy8wuhNGKKtbgTkbXET5F0kX31yyPLZSRw5t5sH0HTN65Oja9oKqi2E22KGdyDoVc+Ai6gpAzIBM2D54nWrWQ8qW5io6gp6sgCRzWz2/lcgJqnPDb1kbKcHvuzrcGKLM+aBaZ2u501Z9EXl8Ji5ndBOXeMr62ZlBrE0ffqMCNHBM9z7NwDXUdlQuRoVYgqjbo2df6uhcmatZR6xNVplSzPmk1IRgxK0L0aThkuqqOlLt7vqte0+n3an8ynd/MYuQP8wiWHmi/O57MIuEZZVTjgd9QPFUJ8vabZ8FvBbakS6NkkwWjJFCyBbmZFPrWU+t7RTPKoto3g8nJZYj7b65vdzwqohs+iwTVEwmJEAGQE0irde9Jn0Tl+U5zrp74MirWZrjn7qU47stDm8lq3uxbHIwaOAmEOkeC+NsN2e8DXVPtzIEof1cpfC7cwn3ONuTaWmpBXaO5+aYwEMq+JIbXR0km3MC0msML5VZmWUbNRAxE1s5LLZ34il43OOX/7cDh05Ylb2soeDK4EH/n6Pf3oIdLIZNY61d/k8lF7vN2QrysldkN4DG7DRZRFcHXcvM/aAue06/kpJkTkw2bR2Kx8rf4IWTtZ4j8KgmXz0rqjSsfmdpbVhGW2cUH7zf9d/tcJ6Gj+W5hkFE7Y1W6yaHHjroluZs2keCl9u4sUutiYEN16+PT0dFFuKUnn8rB63MG9TBkY6WPWVmbGtvrHS0CtIT2G3zWaGoG+NoTbjdZuCu1eA5OiwTwBGbM4DRL9ZwnI2oTvlqm9JnJuRIySEVl07ng8hkRRM/kiRCmkUv4k2ix67rBardLAWqZt+C333J9MoO73e6v5NBoeaV7N2fKYo3woB6AyTa+0ifparDi0P1e+dVKe1R2PIqAj9c+dVG6pr2fLYUbQ2VMRs7YjjGrWueH8PMaIdorM/1r01TEk+uh+MzJN3/K/labnpOHFnhQxB6LUyqj51qqRUc49CrMVNQ3qkC3sFDM4slxa0LO5fnRdeuICGnxijatmd88+78+AWZnVTUzrZLlyWR6wFB8pznDNyF3t6Q0ONQ05el2EWiDO1T/3zQt6bWQK9woZoIVhPwtDD6du3MTxnk3A+W7OnNDwu3tjN76zHVMZwYwSdbb3s9efrWncEU2eRYqz9FN2TXbPKs6pr6ha3Gk9x3x4IVe2yQFBxlq01tFmVCef5+szuGj4tTDErNniz+WzRneb1KKwGeaWoKPbznQvLJAFmEaJYuS61nFFNxa0/mZGZtcjqj+6pZGR0U/L2reM+yPSNqVcwWedc2BTGW0JUycaXxR8uoUdS8Do2sxtjo+uS891n8HMrWFob/AcjvuUekYJfiQPO7cUdfVlJpqWcbg2IbYG6lo2oJRSLu5CKuVyt1tUd7QRJEM037V1+Gjfdih1g8VS3/H9/d2+0CiaVJhb6quUcnmfJBN57VlF7hpFS8DieDyWh4cH60vpLXKKSLPy/DlkO5+itrK7Z7K5gl+tRMhb+NQvVbAA7Hn1B0NvmYzgdiBFzJrtwsK1qBP91px+zT/F+WsxcDez6mNHSvk+aXqXg/om2YYJx7BR7nbkflAt66KPKq1VUmv+LtPUtf26EaZocpfEbyEmd54JV4NBWfsjGN3BxLRYyuVjcCLwenI5ra9lO+Q1NW03s/JOomhvZba7g88rU9QCFrg+W9zsFrkRZBsfooBaKXEUdlSzjm5EUAJkZFHuSMDWmCAzWyO0mNVOU6ow78knw6JjtGh41bTXNI2HNKuCmYIH4wgl8iucZi3FEynefRpdE6F3N9GIRiqlvgF8t9ul7fa0tQQy02+UITPGYeHr3Bg2N1u3+bkNEgymY5dl4Gv0+o/yXWeJBjszUc9NqcOh9z7V3rIoP7qraG5M2XAwN3rztiP1Of9wtI0RoZttWfyoINNQNDhCTctlET4NZGjAI9poodJv7n2yWX1ZICIiKphaNXdBkdWZEU/mMmQaMjKtneXTGoWNkJnBkZbPNvm7zTkMrAG/JwfHM0QW082awTphfBwR0t4dHuz/uh1J6uNmknDu/aC8gLqw0U6q9XodPuR7tYrvaFE/SgkyAu/+UoEIk9ulRDQKiuNZ8CwLsmSbOkbTUny9zjf3u8XXdNHkSINrndG2x2vmzYdukYtC2iMmSsv2Lv3t9qqyZgXRaz8zH4TPuzJ6/nCIH3jtCLM159gSgVTiiYJBukYuWKdgJnDnswjoaHS011Uq5TJtg/Zb6+FAaYSWbMY1MfsOJpcacYxdOwZoKghlp5i8jglH/Fn+js674yO5VCCKgLZoNJcrduhJh3H5jOlatJ4KtJo1MYfLE0V01U35aEYtZYIZnJlJtUG2MiqX54Vjk5cJZUreK9Ki3HZth4vmnRXoZ6SRawSqqSh8q5vA9dWsg6gfqk0y/w/14JZClz7rTaNwH6McfBaYakndOF+brRUXV6hZK0ti2AyuHcvKunIZMrNMzeS5Ux9OSOC3e9MdzjkCxXzg4dLOLM1SPhoxr/mKmQ9ZO9cKl4JTjKbBnHB0KUJGa9CuZew18/fa2naWJ0U4vzIrX0qsXR3BRs69w1ybIbTtUi4DTPzNWK3ijQ+bzWZ4K6KLhoLpauN01k62mQL19ZjDrj5nobQis2BccG+UeXhde8zxm9esNWQSu5Q680bm4xRMYdqetjMzGPVkO4p6x8njyjRrZu5Gx0eIMJvnqWu4JDRNiO8oM/FRY5mFWSPNF5nGtYinHm8xp7n9yFx0fXZ1slTVckqQvYxQk8ZsPnN/3Fw6c1ivzXYctfisDqrVMnMev53vid89eVZ9KAF/63WZnx31O4uz6LmbN4O5gzzxzFgZk2kd2Tm3GLpgShTZhLt2ItPSCZ8osqo4HA5lv99bU6oWNR2JcGbCLmPWXiaunc+sI+c24Pfoc65a7qoBonVU5ouY0QXdro3hJ0W4qJjTAK2OfuaTOo3RYo60SMJIy7csSCa58ThMFQqZX5qldbI+ZcdH7v6p3RgQtTfKdLWobSl+vXvSNzVNW/PRo358Gp9VtWovRkwK1kxsXimTKVNr/S0LPZrLgzbT8TkTFXD+Xqt/lJmtWR9GzeAWos+OtdaXlXE73loDQ4prbRecimFmdVrDHXdauJT+tA4TUJY7LCXORyp6F6k1GozjjljYR1MNm22KqKWkaqZwKXnUXq+ZC3Mygq4jWyu9gZ+Wa1xk2J2/Fia9mMrlvLIkPY5lvlLtFjmYkVmAqZRywRC17WWRaa/9498ji7Va5S9abjEJGTWzvWbSRtB6WxghEza9zNSCyIUppZ5CZKssqiOa24/yXbuZtabVpiKrjyUpiL7HVG3ZSZP9V7TkKWv96bnrI6pz7juNHCKrRlHz+aagJQrfW9ecdS6Nm34ifw1ZauajUeuT26Orgkf/9zLlLc7LHeP41Mz6mdEi0ZU578z318adWW8Md4a8I8KdWRdCb060lGUYteXJHb3nUC/72qNprgx3wXWOO7MuhJHN30u8d3apfayjgRnH1HembMOdWRdC65ZKxdwR3yxVlSE7P0WDRrfQZXBR878ig38KZnWbHK6FuYhCc9HX6sfI0yBqjFzb29xrHremuXrK/4gM/SmYteUumltDlrP9qFusol1ncyNLN11L6H4GGunFp2DWHw0fxazXaLtHS/6IDLUkrmdP3nHHHZNwZ9Y77vgkuJvBC+GaQbC5cTdPbxN3Zl0Id4K/Y258Smb9DKmbkR1MU9pzcDdot2BKEGr0Zv076vhUzIqQf/bwrbmeb9SaWohurB99vlHrjfotaL3zx2F0U4T2WR/b01vfX30jBONTMWv2lAZd8KkbKVrLO6LOCN3dOL8kRph86qYIoPWxPb2bIv6q+FTMCozk8q692C2PZ7nlm89H5+vOVMvhUzLrZ0b2pI07g9yRYXZmVd/iR/Y1Rn0+nHfPkZobIyb3tW/ju6MNq57FXK1W/6+U8n+X684dd/zl8X+Ox+P/uBNdzHrHHXd8HD7vNps77viL4c6sd9zxSXBn1jvu+CS4M+sdd3wS3Jn1jjs+Ce7MescdnwR3Zr3jjk+CO7PecccnwZ1Z77jjk+D/A5ItWEwm7OB1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(cov_mat, cmap=\"gray\")\n",
    "plt.gca().set_xticks([])\n",
    "plt.gca().set_yticks([])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}