discussion/examples/cross_gpu_logprob.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "cross-gpu-logprob.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ICfV2dqVHpjh",
        "colab_type": "text"
      },
      "source": [
        "The objective of this notebook is to demonstrate splitting a log_prob and gradient computation across a number of GPU devices. For development purposes, this was prototyped in colab with a single GPU partitioned into multiple logical GPUs. \n",
        "\n",
        "*Note*: Since it runs on a single GPU, performance is not representative of what can be achieved with multiple GPUs. Usage of `tf.data` can likely benefit from some tuning when deployed to multiple GPUs.\n",
        "\n",
        "**Needs a GPU**: Edit > Notebook Settings: `Hardware Accelerator` => `GPU`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Dk5b7sR75zkU",
        "colab_type": "code",
        "outputId": "b704d781-ddc4-4cb7-992b-be29ecd0da21",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 122
        }
      },
      "source": [
        "%tensorflow_version 2.x\n",
        "import numpy as np\n",
        "import tensorflow as tf\n",
        "import tensorflow_probability as tfp\n",
        "tfb, tfd = tfp.bijectors, tfp.distributions\n",
        "\n",
        "physical_gpus = tf.config.experimental.list_physical_devices('GPU')\n",
        "print(physical_gpus)\n",
        "\n",
        "tf.config.experimental.set_virtual_device_configuration(\n",
        "    physical_gpus[0],\n",
        "    [tf.config.experimental.VirtualDeviceConfiguration(memory_limit=2000)] * 4)\n",
        "gpus = tf.config.list_logical_devices('GPU')\n",
        "print(gpus)\n",
        "\n",
        "st = tf.distribute.MirroredStrategy(devices=tf.config.list_logical_devices('GPU'))\n",
        "print(st.extended.worker_devices)"
      ],
      "execution_count": 42,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]\n",
            "[LogicalDevice(name='/device:GPU:0', device_type='GPU'), LogicalDevice(name='/device:GPU:1', device_type='GPU'), LogicalDevice(name='/device:GPU:2', device_type='GPU'), LogicalDevice(name='/device:GPU:3', device_type='GPU')]\n",
            "WARNING:tensorflow:NCCL is not supported when using virtual GPUs, fallingback to reduction to one device\n",
            "INFO:tensorflow:Using MirroredStrategy with devices ('/job:localhost/replica:0/task:0/device:GPU:0', '/job:localhost/replica:0/task:0/device:GPU:1', '/job:localhost/replica:0/task:0/device:GPU:2', '/job:localhost/replica:0/task:0/device:GPU:3')\n",
            "('/job:localhost/replica:0/task:0/device:GPU:0', '/job:localhost/replica:0/task:0/device:GPU:1', '/job:localhost/replica:0/task:0/device:GPU:2', '/job:localhost/replica:0/task:0/device:GPU:3')\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tpRKvEJU7cJP",
        "colab_type": "code",
        "outputId": "2df00801-9f2e-4aae-802a-2f3e9ceff0eb",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 153
        }
      },
      "source": [
        "# Draw samples from an MVN, then sort them. This way we can easily visually\n",
        "# verify the correct partition ends up on the correct GPUs.\n",
        "ndim = 3\n",
        "\n",
        "def model():\n",
        "  Root = tfd.JointDistributionCoroutine.Root\n",
        "  loc = yield Root(tfb.Shift(.5)(tfd.MultivariateNormalDiag(loc=tf.zeros([ndim]))))\n",
        "  scale_tril = yield Root(tfb.FillScaleTriL()(tfd.MultivariateNormalDiag(loc=tf.zeros([ndim * (ndim + 1) // 2]))))\n",
        "  yield tfd.MultivariateNormalTriL(loc=loc, scale_tril=scale_tril)\n",
        "\n",
        "dist = tfd.JointDistributionCoroutine(model)\n",
        "tf.random.set_seed(1)\n",
        "loc, scale_tril, _ = dist.sample(seed=2)\n",
        "\n",
        "samples = dist.sample(value=([loc] * 1024, scale_tril, None), seed=3)[2]\n",
        "samples = tf.round(samples * 1000) / 1000\n",
        "for dim in reversed(range(ndim)):\n",
        "  samples = tf.gather(samples, tf.argsort(samples[:,dim]))\n",
        "\n",
        "print(samples)"
      ],
      "execution_count": 43,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "tf.Tensor(\n",
            "[[-4.534 -5.856  3.606]\n",
            " [-4.527 -5.875  1.671]\n",
            " [-4.269 -4.346  5.697]\n",
            " ...\n",
            " [ 2.158  5.71  -2.926]\n",
            " [ 2.302  6.658 -3.491]\n",
            " [ 2.632  5.67  -4.854]], shape=(1024, 3), dtype=float32)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A1qBcZgSIN6q",
        "colab_type": "code",
        "outputId": "c8159097-5a70-44e5-ce2b-f1a6a72e78c6",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 119
        }
      },
      "source": [
        "print(loc)\n",
        "print(scale_tril)\n",
        "print(tf.reduce_mean(samples, 0))"
      ],
      "execution_count": 44,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "tf.Tensor([-1.0574996   0.24829748  1.0737331 ], shape=(3,), dtype=float32)\n",
            "tf.Tensor(\n",
            "[[ 1.1475685   0.          0.        ]\n",
            " [ 1.9094281   0.5724521   0.        ]\n",
            " [-1.1899896   0.49813363  1.5088601 ]], shape=(3, 3), dtype=float32)\n",
            "tf.Tensor([-0.9953702  0.3626416  1.0675195], shape=(3,), dtype=float32)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3uIDBMnndNP-",
        "colab_type": "text"
      },
      "source": [
        "### Single batch of data resident on GPU."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BY2R0E5ndPiz",
        "colab_type": "code",
        "outputId": "75eb688e-57b3-41ba-ec34-f4132f332d59",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 476
        }
      },
      "source": [
        "%%time\n",
        "\n",
        "def dataset_fn(ctx):\n",
        "  batch_size = ctx.get_per_replica_batch_size(len(samples))\n",
        "  d = tf.data.Dataset.from_tensor_slices(samples).batch(batch_size)\n",
        "  return d.shard(ctx.num_input_pipelines, ctx.input_pipeline_id)\n",
        "\n",
        "ds = st.experimental_distribute_datasets_from_function(dataset_fn)\n",
        "\n",
        "observations = next(iter(ds))\n",
        "# print(observations)\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "def log_prob_and_grad(loc, scale_tril, observations):\n",
        "  ctx = tf.distribute.get_replica_context()\n",
        "  with tf.GradientTape() as tape:\n",
        "    tape.watch((loc, scale_tril))\n",
        "    lp = tf.reduce_sum(dist.log_prob(loc, scale_tril, observations)) / len(samples)\n",
        "  grad = tape.gradient(lp, (loc, scale_tril))\n",
        "  return ctx.all_reduce('sum', lp), [ctx.all_reduce('sum', g) for g in grad]\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "@tf.custom_gradient\n",
        "def target_log_prob(loc, scale_tril):\n",
        "  lp, grads = st.run(log_prob_and_grad, (loc, scale_tril, observations))\n",
        "  return lp.values[0], lambda grad_lp: [grad_lp * g.values[0] for g in grads]\n",
        "\n",
        "singleton_vals = tfp.math.value_and_gradient(target_log_prob, (loc, scale_tril))\n",
        "\n",
        "kernel = tfp.mcmc.HamiltonianMonteCarlo(target_log_prob, step_size=.35, num_leapfrog_steps=2)\n",
        "kernel = tfp.mcmc.TransformedTransitionKernel(kernel, bijector=[tfb.Identity(), tfb.FillScaleTriL()])\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "def sample_chain():\n",
        "  return tfp.mcmc.sample_chain(\n",
        "      num_results=200, num_burnin_steps=100,\n",
        "      current_state=[tf.ones_like(loc), tf.linalg.eye(scale_tril.shape[-1])], \n",
        "      kernel=kernel, trace_fn=lambda _, kr: kr.inner_results.is_accepted)\n",
        "samps, is_accepted = sample_chain()\n",
        "\n",
        "print(f'accept rate: {np.mean(is_accepted)}')\n",
        "print(f'ess: {tfp.mcmc.effective_sample_size(samps)}')\n",
        "\n",
        "print(tf.reduce_mean(samps[0], axis=0))\n",
        "# print(tf.reduce_mean(samps[1], axis=0))\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "for dim in range(ndim):\n",
        "  plt.figure(figsize=(10,1))\n",
        "  plt.hist(samps[0][:,dim], bins=50)\n",
        "  plt.title(f'loc[{dim}]: prior mean = 0.5, observation = {loc[dim]}')\n",
        "  plt.show()"
      ],
      "execution_count": 45,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "accept rate: 0.8\n",
            "ess: [<tf.Tensor: shape=(3,), dtype=float32, numpy=array([42.702564, 56.667793, 43.04328 ], dtype=float32)>, <tf.Tensor: shape=(3, 3), dtype=float32, numpy=\n",
            "array([[57.469185,       nan,       nan],\n",
            "       [34.60636 , 43.471287,       nan],\n",
            "       [13.848625, 30.807411, 72.34708 ]], dtype=float32)>]\n",
            "tf.Tensor([-0.18598816  0.7396643   0.4074543 ], shape=(3,), dtype=float32)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAV+0lEQVR4nO3debgcVZnH8e8PgolssiSywxXZwyPo\nMIAIGgTZd1CBCEFhEOdhxnFDEBeWOIIbwig6DEKCLKLIvigBhACCGhSIkTUYIRBIEAgEIus7f5zT\npFJ039t1l+57b/8+z1PP7a46VfX2qdPVb51TfVsRgZmZmZk1b4l2B2BmZmY21DiBMjMzM6vICZSZ\nmZlZRU6gzMzMzCpyAmVmZmZWkRMoMzMzs4qcQHUISbMk7dgP21go6WdNlj9c0gJJIWm9vuy7m31s\nJ+mBgdi29Q9JXbkNjGh3LFVIGi/p+nbHYWaDkxMoq2rPiDik9iR/OP5W0kuS7i8maRHx04hYdiCD\niYhbI2LDgdzHUCZppKRzJD0v6UlJn++m7GGSXs9Jb20a18Jw26ZekhcRF0TETu2MqyeSTpY0XdJr\nkk7ooawknSrpH3k6VZIKy0PSi4Vjf3Zh2XWldvGKpOl19vGhvJ2JhXkjJZ0m6QlJz0o6U9JSheUb\nS7pJ0nxJD0vat7TNpfM6T+cyU+vs922S7pM0uzBvu1LMtYu5/fPyTSX9Jm+34T9ElLS+pH9KOr80\n/2BJf891drmklQrLbs7r1Pb7QIV1z5c0J79nH5R0RKPYrL2cQFlfXQT8GVgZOB64RNKYVuy4rz0a\nkpbsr1gGsROA9YF1gO2BYyTt0k35OyJi2cJ0cwti7LOcHHTi+exh4BjgmibKHgnsA2wGvAfYE/h0\nqcxmhWP/5gd3ROxabBfA74BfFlfMSdHpwO9L2zwW2ALYFNgAeB/w1bzOCOAK4GpgpRzj+ZI2KKx/\nVl62cf77uTqv7UvAvOKMfHFVjHkPYAHw61zkVeAXwOF1tlf0I+CPpdc6Fvhf4BBgFeAl4MzSekcX\n9r9hhXW/BXRFxPLAXsBESf/SQ4zWBp14wul4+YrwB/mK8In8eGRh+d6S7s5XQDMbfeDmk9z7gG9E\nxMKI+BUwHdi/D7HNknScpL/mq9VzJY3Ky8ZJmi3py5KeBM6tzSusv3G++ntO0gxJexWWTZL0Y0nX\nSnqRlFCU93+zpImSfpevHK+StLKkC3J9/FFSV6H8RpKmSHpG0gOSPlZYtrukP+f1HlOhh0CLejwm\nSHo0XwUf39t668YE4OSIeDYi7gP+Dzisv3ciaQlJX81X1XMlnSfpHaVin8rtbY6kLxbW3VLStFxP\nT0n6fmHZ1vlYPCfpHhV6xPKx+qak20kfQl+SNK0U1+ckXZkfNzweQK1X47l83N+v1CN3W2Fb2+Tj\nPz//3aYUy8mSbpf0gqTrJY3udYU2KSImR8R1wAtNFJ8AfC8iZkfE48D36EVbyO1/O+C80qIvANcD\n95fm7wmcERHPRMQ84AzgU3nZRsDqwGkR8XpE3ATcTkoukLQRKYk4MiLm5TJ3leJ5F/AJUuLRnQnA\nJRHxIkBEPBARPwVmdPNaDwSeA24sLRoPXBURUyNiAfA1YD9Jy/UQQ4/rRsSMiHg5l408vbuJ7VqL\nOYHqTMcDWwObk65Gt2TRFeGWpBPjl4AVgA8CsxpsZyzwSEQUT9735Pl1STpW0tU9xDce2Jl00tig\nFlu2KukqdB3S1Wpx20sBV5FO4u8E/gO4QFJxiO9g4JvAcsBt1Hcg6QS+Ro7hDuDcvN/7gG/k/S0D\nTAEuzPs7EDhT0iZ5Oy8Ch5LqcXfgM5L2Ke1rW2BDYAfg65I2rhdQrrfnGk0N1lkRWI10TGq6PT7A\ne3My96Ckr6n5Xr7D8rQ9sC6wLPDDUpntSb1hOwFf1qLh3tOB0/MV97tJvQJIWoPUszKRVPdfBH6l\nxXs4DyG1g+WAnwAbSlq/sPxg0vGB7o/HB/PfFXKPwR3FwJWGWK4hffivDHwfuEbSyqV9fZLUFt6W\n430LSWt3dywlHVxvvX4wlp7bwlSlod5LixcKJYcCt0bErNoMSeuQkqKTGqyj0uM16yTYxeWb5sdb\nAn8HTsztcrryEFzB/wBfARY22F7tvXoAMLlRmTrrLE96PfWGvRery4iYCbxCOl/VfCvHfLsWHwrv\ncV2lIcuXSMnoHODaZuO21nEC1ZnGAydFxNx8RXgi+YqP1J19TkRMiYg3IuLxiChfUdYsC8wvzZtP\n+jCrKyJOiYg9eojvhxHxWEQ8Q0p2Diose4PU4/VyRJRPmFvnmE6JiFfy1ezVpfWviIjb82v7Z4P9\nnxsRMyNiPnAdMDMiboiI10jDFu/N5fYAZkXEuRHxWkT8GfgV8NH8Wm+OiOl5X/eShjs/VNrXibn3\n7h7SSXWzegHleluh0dTgddTuPyseo+6Oz1TSB9c7Sb2IB5ES6WaMB74fEY/kq+rjgANLCdiJEfFi\nREwnJaS14/IqsJ6k0RGxICLuzPM/AVwbEdfmOpwCTAN2K2xzUr5ify0frytq282J1EbAldD08Whk\nd+ChiPhZ3tdFpA+3PQtlzo2IB3O7/AXpAuUtIuLR7o5lRFxYb71+UH6/zgeWld68D+pDQBepzp4A\nrm6QQB8KTCrNOwP4Wj72Zb8GPitpjKRVgf/M85cGHgDmknoPl5K0U45j6VxmTVKbnE/qqToamFy7\n0FC6X2rJiLish9e+H/A0cEsP5YpOBn4aEbPrLOvp3Pdl0oXEGqQhyKskvbvJdYmIf8/PtwMuBV7G\nBh0nUJ1pddJVXc3f8zyAtYCZTW5nAbB8ad7yNDec0J3HCo+LsQHM6ybxWR14LCLeKK2/RoNtN/JU\n4fHCOs9rick6wFalnqDxpF4yJG2ldIP9PEnzgaOA8rDOk4XHLxW23R9qH2bFY9Tw+OTk5285wZhO\nuvo+oMl91WtTI0j3eNQ0Oq6Hk66+789DY7UEex3go6X63ZbUq1Zvm5B6m2qJ2cHA5RHxEjR9PJp9\nfbXXUGxbA3ksURqSrt2UvF0vNlF+vy4PLIhIvyifh5ReiYjngM8C7yLdd1SMYVtS+76kMG9PYLmI\nuLjBfr9Juk/ybtK9U5eTkuanIuJV0n1Zu5Pq7wuk5LOWtCzMZSfm2G4BfgvslHuVvs2ihKw7E4Dz\naq+1J5I2B3YETmtQpNtzX0T8PiJeyBd6k0nDkrs1s25NHq68jZREfqaZuK21nEB1pidIH041a+d5\nkD6Qmh1vnwGsWxr334xu7ilo0loNYoN0P0AjTwBrafGbidcGHm9y/aoeA24p9R4sGxG1k92FpN6P\ntSLiHaQhJjXaWHckfUVv/UbRm1O9dSLiWVL3f7FXq8rxiQrx1mtTr7F48ln3uEbEQxFxEKnn61TS\nFxGWIdXvz0r1u0xEnFKKsWgKMCZ/AB7EouE76P549NQuyq+v9hoer1O2W3kIr+GxlDS+3noRMTYW\n3ZR8a9X9ko57lbZQ7/hPAC4t9TTtAGyRh/6eBD4O/JekK3LcCyPi6IhYIyLWBf4B3FW70ImIeyPi\nQxGxckTsTOq5+UPe9r0N4oI0HNwF3Jr3eymwWo6jq1ZY0lrAON56z1Z3xuVtP5q3/UVgf0l/yssX\nq0tJ6wIjgQcbbK9Yl1XXHYHvgRqcIsJTB0yk+5h2zI8nkq4Ex5CuwG8jXeFBuufgOdJJcQnSFfZG\n5W0Utnsn8F1gFLBvXndMqUwA61WIczrpqmulHNt/52XjgNml8m/OI9138gjpWz9L5WUvFOKfVHud\n3ez/ZuCIwvOJpGGi2vMdgYfz4+VIvRCH5P0tBfwrsHFePheYUKjXucD5+XlXrpcRjfbdT8f9FNKw\nxYqkoZk5wC4Nyu4KrJIfbwT8hTRcWozvhAbrHgE8ROq1WJbUQ1F+rReQhmbG5rrYKS//RK3N5Pr9\nJ/B2UsL1JOl+uCVzGxsHrNldfQE/JiVSc0v1293xWBp4HdigUP4w4Lb8eGVS2z6Y9IH28fx8dIN2\n8+a6A/y+XirXy4W5rY4iDWnVK3sU6R6+NUg9ajOAo/KysaQhxyXz8fsBaXhtqcL6bycNNX24tN3l\nSL1SteliUs/NSnl5bX8iDbM/Vjv2efl7ctxLkxKVvwEjC6/vYdKN1iOAD5Df0/l5cb/7kRLdVYt1\nQLo/amqd+lDe7yak9jmqsN+lS9v+LqlNjynU1/OkIbZlgPOBn+dlK5Da7Kgc43jS/XcbNLFu7V7K\nZfOx2Dmvu9dAtyVPvXj/tTsATy060IsnUKNI9yzMydMZwKhC2X1JV34v5JPXzuVtFMp2kT48FuYT\n7o519v1mApVPZtf1EOdxwF9JH1CTgaXzsnF0k0Dl52NJCcP8vI19C8sm0Y8JVH6+Ienm4nmkK+ub\ngM3zsgNICdYLpHuxfkjrE6iRwDn5hP0U8PnCsrVJwwlr5+ffzWVeJCWiJ7H4B+hM4CMN9rME8HXS\nh+M80ofCiqXXeiTpA+5J4JjCuueTkpkFpA/1fQrLtsrH85m83WsK8datL9IHUwA/Ks1veDzy8pPy\nPp4jfdAfRiEJIg0f3pXb1l3Att20m8XWHcD39SQWfVOrNh1WqIcFhbIiDXk9k6dvA8rLPkx6/76Y\nj8XlwPqlfR2U609NxDSx8Lz2RZSX8j7Gl8p/B3g2H//rKF1skd7Td+TYFntPd3cuKMy/Hzi8zvxa\nuyxOsxps+4RiW8nzDgYezXFdwaKEcQzp3x68kNvSnZTeNz2se0te73nSxeS/DXQ78tS7qfbmMeuR\n0j+DWw24LCImNFH+k6Qr0VHAJhHxSBPrzCJ9EN3Qx3CtH0laE/hFRGzTY2Ezsw7gBMoGFSdQZmY2\nFPgmcjMzM7OK3ANlZmZmVpF7oMzMzMwq6tOPsVY1evTo6OrqauUuzczMzHrlrrvuejoixtRb1tIE\nqquri2nTpvVc0MzMzKzNJJV/geBNHsIzMzMzq6ilPVBmZu3Udew1TZWbdcruAxyJmQ117oEyMzMz\nq8gJlJmZmVlFTqDMzMzMKnICZWZmZlaREygzMzOzipxAmZmZmVXkBMrMzMysIidQZmZmZhU5gTIz\nMzOryAmUmZmZWUVOoMzMzMwqcgJlZmZmVpETKDMzM7OKnECZmZmZVeQEyszMzKwiJ1BmZmZmFTmB\nMjMzM6toRE8FJJ0D7AHMjYhN87yVgIuBLmAW8LGIeHbgwjQza52uY69pqtysU3Yf4EjMbLBqpgdq\nErBLad6xwI0RsT5wY35uZmZm1hF6TKAiYirwTGn23sDk/HgysE8/x2VmZmY2aPU4hNfAKhExJz9+\nElilUUFJRwJHAqy99tq93J2Z2eAznIb6htNrMWuFPt9EHhEBRDfLz4qILSJiizFjxvR1d2ZmZmZt\n19sE6ilJqwHkv3P7LyQzMzOzwa23Q3hXAhOAU/LfK/otIjMb9jxc1HeuQ7P26rEHStJFwB3AhpJm\nSzqclDh9RNJDwI75uZmZmVlH6LEHKiIOarBoh36OxczMzGxI6O0Qnpl1CA8V9V2zdQiuR7Ohwj/l\nYmZmZlaREygzMzOzijyEZ4aHqcw6kd/31hfugTIzMzOryAmUmZmZWUVOoMzMzMwq8j1QZtYvqnxV\nv7+36XtUGhuI49Kf/C8ebKhyD5SZmZlZRU6gzMzMzCryEJ7ZENHfQzHDaThksA9TVTGcXovZcOYe\nKDMzM7OKnECZmZmZVeQhPLMO5aEi6412thu32eFvKH3z1j1QZmZmZhU5gTIzMzOryEN41hKd+A0y\nDzeYWT1DaZjKGnMPlJmZmVlFTqDMzMzMKvIQXg+G0+809Xe38VD4Nk5/HxMPy5l1Hg+5WT3ugTIz\nMzOryAmUmZmZWUV9GsKTtAtwOrAkcHZEnNIvUfXBUOhqHezfSBtOw1TD6bWY2eDW3+ebgTh/9fct\nGp08bNnrHihJSwI/AnYFNgEOkrRJfwVmZmZmNlj1ZQhvS+DhiHgkIl4Bfg7s3T9hmZmZmQ1eioje\nrSgdAOwSEUfk54cAW0XE0aVyRwJH5qcbAg/0PtyWGg083e4gBiHXS32ul/pcL/W5XupzvTTmuqlv\noOtlnYgYU2/BgP8bg4g4CzhroPfT3yRNi4gt2h3HYON6qc/1Up/rpT7XS32ul8ZcN/W1s176MoT3\nOLBW4fmaeZ6ZmZnZsNaXBOqPwPqS3iXpbcCBwJX9E5aZmZnZ4NXrIbyIeE3S0cBvSP/G4JyImNFv\nkbXfkBt2bBHXS32ul/pcL/W5XupzvTTmuqmvbfXS65vIzczMzDqV/xO5mZmZWUVOoMzMzMwqcgKV\nSfqOpPsl3SvpMkkrNCi3i6QHJD0s6dhWx9lqkj4qaYakNyQ1/KqopFmSpku6W9K0VsbYDhXqpdPa\ny0qSpkh6KP9dsUG513NbuVvSsP3ySU/HX9JISRfn5b+X1NX6KFuviXo5TNK8Qhs5oh1xtpqkcyTN\nlfSXBssl6Yxcb/dKel+rY2yHJuplnKT5hfby9VbE5QRqkSnAphHxHuBB4LhygQ79+Zq/APsBU5so\nu31EbN4h/6ukx3rp0PZyLHBjRKwP3Jif17Mwt5XNI2Kv1oXXOk0e/8OBZyNiPeA04NTWRtl6Fd4X\nFxfayNktDbJ9JgG7dLN8V2D9PB0J/LgFMQ0Gk+i+XgBuLbSXk1oQkxOomoi4PiJey0/vJP1fq7KO\n+/maiLgvIobKf49vmSbrpePaC+n1Tc6PJwP7tDGWdmvm+Bfr6xJgB0lqYYzt0Invi6ZExFTgmW6K\n7A2cF8mdwAqSVmtNdO3TRL20hROo+j4FXFdn/hrAY4Xns/M8gwCul3RX/vke68z2skpEzMmPnwRW\naVBulKRpku6UNFyTrGaO/5tl8gXcfGDllkTXPs2+L/bPw1SXSFqrzvJO1InnlGa9X9I9kq6TNLYV\nOxzwn3IZTCTdAKxaZ9HxEXFFLnM88BpwQStja6dm6qUJ20bE45LeCUyRdH++ahiy+qlehp3u6qX4\nJCJCUqP/k7JObi/rAjdJmh4RM/s7VhuyrgIuioiXJX2a1Ev34TbHZIPXn0jnlAWSdgMuJw1zDqiO\nSqAiYsfulks6DNgD2CHq/4OsYfnzNT3VS5PbeDz/nSvpMlI3/ZBOoPqhXjquvUh6StJqETEnDy3M\nbbCNWnt5RNLNwHuB4ZZANXP8a2VmSxoBvAP4R2vCa5se6yUiinVwNvDtFsQ1FAzLc0pfRcTzhcfX\nSjpT0uiIGNAfX/YQXiZpF+AYYK+IeKlBMf98TR2SlpG0XO0xsBPpJutO14nt5UpgQn48AXhLT52k\nFSWNzI9HAx8A/tqyCFunmeNfrK8DgJsaXLwNJz3WS+m+nr2A+1oY32B2JXBo/jbe1sD8wpB5x5K0\nau3eQUlbknKbgb8QiQhP6Xz1MGls+e48/STPXx24tlBuN9K39GaShnLaHvsA18u+pHH2l4GngN+U\n6wVYF7gnTzNcLx3dXlYmffvuIeAGYKU8fwvg7Px4G2B6bi/TgcPbHfcA1sdbjj9wEulCDWAU8Mt8\n/vkDsG67Yx4k9fKtfC65B/gtsFG7Y25RvVwEzAFezeeXw4GjgKPycpG+wTgzv3e2aHfMg6Reji60\nlzuBbVoRl3/KxczMzKwiD+GZmZmZVeQEyszMzKwiJ1BmZmZmFTmBMjMzM6vICZSZmZlZRU6gzMzM\nzCpyAmVmZmZW0f8Ds9XTtDkl2RkAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAW6klEQVR4nO3debgcVZ3G8e8LYVPCEhIigSzs6wgi\nsgnIEB6FqCwCkUUII4jI4IaPiOggIrLOKDKAwAACAhFEkVURkMiiLIFh3wNhCSHsS4QBA7/545yG\nupXue7uTvt19+76f56nnVtU5VXVOnVpOnTrVVxGBmZmZmdVvgXYnwMzMzGygcQXKzMzMrEGuQJmZ\nmZk1yBUoMzMzswa5AmVmZmbWIFegzMzMzBrkClSXkTRd0tZNWMdbkn5dZ/ytJc2W9N78bruXbYzJ\n21iwP9ZvzSEpJK3S7nQ0QtLmkh5udzrMbGBxBcpq+XxE7FmZkPQTSfdKmiPp8GLEiLg2IhYHnuqv\nxETEUxGxeES821/bGMiUHCvppTwcK0k14m6ZK7uzC8OkVqe5XcqVvIi4MSJWb2eayiSNl/SQpDcl\nXS9pbI14y0qaLOlZSa9JulnSRjXinlXOu6Rxkq6S9Iqk5ySdJGlIDltN0qWSXpD0sqSrJa1eWHYR\nST/P235F0imSFiqEzy4N70r67yrpOiyna+vCvGGSLszH8ouSzpe0RD15lnRoabtv5eN9eA7/T0mP\nSnoj7+O9SunZStKdkl6X9Lik/QphkvQDSU/l8N9U0pXDz5b0Tmn7CxbCJ0p6MG/7AUk7FMLWyfv4\nRUlz/UCjpPMkzczbfUTSvtXK2VrHFSir12PAwcCVrd5w5YLeruUHiP2AHYB1gY8Cnwe+2kv8Z3OF\ntDKc04pENkO3t0LmG/3vgf8AhgFTgQtrRF8cuB34eI57DnClpMVL69wMWLnK8qcAzwPLAesBnwIO\nyGFLAZcBqwMjgduASwvLHgJsAKwDrAasD/ywElg8voCPAG8Bvy2la2VgF2BmKV1HAksDK+Z0jwQO\nryfPEXFUadvHAlMi4sW8/D9I58eSwCTgF5I2zelZCLgEOC2HfxH4maR187J7AXsCnwRGAYsB5Urh\ncaVz69287uWB84CDgCWA7wIXSFo2L/dP4CJgH6o7GhgXEUsA2wFHSvp4jbjWChHhoYsGYDqwdR5f\nBDgBeDYPJwCLFOJuD9wFvA5MA7Ypr6PK+s8DDu9r23WkcwrpgnBb3v6lwLAcNg4I0oXkKeCGwrwh\nOc4o0sX9ZVLl7iuFdR8OXJzT+jqwb5Xtn026efwRmA3cTLrInwC8AjwEfKwQfxTwO+AF4AngG4Ww\nDYG/A6+SbgQnAQsXwgPYH3g0xzkZUJPL/W/AfoXpfYBbasTdEnhmPrb1lbzPX85lMKqU128AjwMv\nAscDC+SwVYC/Aq/lsAsLy60BXJPX+TAwsVRWvwSuIt38vgc8ByxYiLMjcE9f5ZGPpcjrmU26QfbY\nH8Ca+fh8Fbgf2K6UlpNJDxJvALcCKze5LPcD/laY/jCp8rFGncu/Dny8MD0E+F9SxTqAVQphDwIT\nCtPHA6fVWO+wvPwyeXoqsEshfHfg6RrLTsrHhErz/wRMoHTtIJ2XBxSm/x24ut48F+Yrb3dSL8te\nBnwnj4/MefxQIfx2YLc8fjHw3ULYpsD/VeLn4+PIGtvZCHi+NO8FYJPSvFWA6KOMV8/H9sTe4nno\n38EtUN3tB8DGpCfLdUk3lh8CSNoQOJf0FLQUsAXpItYUknaXdE8f0fYCvkx6+p0DnFgK/xTpZvaZ\nKsv+BniGVLHZGThK0laF8O1JF7ulgPNrbH8iaX8MB94m3XTvzNMXAz/LeVkAuBy4G1geGA98S1Il\nXe8C387LbZLDD6CnzwGfIN3EJtbIU2W/vdrLMKZGXtbO6au4O8+rZVlJsyQ9kV/DfLiXuMX0bUWq\n+E4klduTpLIo2pHUMrE+qRy+nOf/BPgzqWVhBfKTe972NcAFwLLArsApktYqrHN34KfAUOAXpArQ\nVqXwC/J4zfKIiC1ynHUjtQ70aNnJLRCX53QuC3wdOL/46iqn78c5H4/ldFXVR1keUmOxHmUZEf8g\nPeD0Vp6V7a0HLJzTVfFt4IaIqHY+ngDsKulDuYVkW1KlppotgOci4qXiJkvjK0hassqyk4BzI9/9\nc1p3Ad6OiKuqxD8Z+JykpSUtDexEqlTNpUaeKzYnlePvaiy7GOm8vB8gImYBk4F/k7SgpE2AscBN\npXwWxxcBVi3MO0DplecdknYqzJ8KPChpu7zuHUjXnb6uk8X0niLpTdID3kzSQ4W1S7trcB6aO9Cz\nBWoaPZ8uPwNMz+OnAT/vax1VwprZAnVMYXot4B1gQT5obVqpEF6ZNwQYTbpJDi2EHw2cnccPJ90w\netv+2cD/FKa/DjxYmP4X4NU8vhHwVGn57wO/qrHubwGXFKYD2KwwfRFwSJPL/V0KLRSkC3pQpaWL\n1NK2FukV/oqkVpmqrQ5Vlj2T9IqiMr046dXDuEJetymEHwBcl8fPBU4HViit84vAjaV5pwE/KpTV\nuaXwI4Gz8vhQUoVqbAPlUWyF2ZLcAkW64T5HbjXL8yZXjvmcljMKYROAh5pclmcWz40872Zg7z6W\nWwK4F/h+Yd5oUsViyRp5XxO4g/QAEzl/1Y6ZFYAZ5JaYQhncDIzIx9SteR3LlZYdm4/PFQvzhpJa\nZCvHzXR6tkCNAq4F3svDNRRadXvLc5V9eXYv++wcUoVRhXmfB2blfTKHnq3b+wKPkK5HS5Jar4Lc\nikR6aFiGdJ2aQGql/GRh+X1ILZ9zgDeBz1ZJU68tUKRr5Gakh7+FmnnseWhscAtUdxtFaiGoeDLP\ng3RhndbyFPX0dGH8SWAhUqtBtfCiUcDLEfFGafnl61i2aFZh/K0q05V+JGOBUcXWA+BQUnN/pbPt\nFUqdcF8HjirlA9JNueLNwrqbZTbpZlKxBDA78hW3KCKei4gHIuK9iHiC1Ldtp3K8GnocUxExG3iJ\n2vu+eMwdTHpiv03S/ZIqLVNjgY1K+3cP0k252johtTZ9QdIiwBeAOyPiSai7PHrL39MR8V4pD8X8\ntbosydNvVIkLvN+Scjnpte3RhaATgCMi4rUqyyxAqjz8nvSacDipVe3YUrwRpBa5UyJiciHop6RX\ng3eRXiH/gVSZLp5HkPoM3ZSPtYrDgV9HxPQaWbqIVFEZSsr7NNLDWz15roR/iNS/qmr/PknHk/pv\nTaycJ5LWILWo7kVq1VobOFjSZ/NiZ5Eq1FNIrVbX5/nPAETEnRHxUkTMidSydj7p+ESpk/xxpAr7\nwqQW9jNyC1rdIuLdiLiJVKn9WiPLWnO5AtXdniXdnCrG5HmQbkjVOpW20ujC+BjSxffFwry5bv7Z\ns8AwSUNLy8+oY9l58TTwREQsVRiGRsSEHP5LUpP6qpE6eB5Kz2b+uknaQ3N/vVQcar3Cu5/0mrZi\n3TyvHkH914Iex1R+/bYMPfd9uVyfhfcrbl+JiFGkDu6nKH0R9jTw19L+XTwiijeHHuUZEQ+QKjbb\n0vP1HcxfeTwLjM6Vi2IeZtSI36s+yvLQGov1KMu8j1emRnnmSuQfSDfx8ocD44Hjc2WyUvH7u6Td\nSX2axgAnRcTbkV7N/YrUclJZ99KkytNlEdHjVWVEvBURB0bE8hGxEqkifUep8gmpMlKuxIwHvlFI\n12jgIknfy+HrkVpF/5Er6aeW0tVbnit2JPWpm1IOkPRj0rHz6Yh4vRC0DvBIRFydHzAeJvV32zbn\n+b2I+FFEjIuIFUhlMoPax0fwwbG3HqllfGpez+2kVrt5/emXIbT/Gj6ouQLV3SYDP5Q0In/Zcxgf\nPMWdSXrPP17SApKWz09fVUlaSNKipGNmiKRFNf9fQ31J0lr5SfEI4OKo42cKIuJp0hPv0TkdHyU1\njZ/X+5Lz7DbgDUnfk7RY7r+wjqRP5PChpE6ss/M+nOenwog4P3p+wVMeav1UxLnAQbkcRwHfIb2O\nmYukf5U0Vslo4BgKX1cpfYpddVk+6B+yXr6JHQXcWmpJ+G7uuzIa+Cb5CzJJu0haIcd5hXRzeQ+4\nAlhN0p75OFtI0ickrdnH7rogr38Len7d1Vd5zAJWqrHOW0mtSgfndGxJeqVT7udVlz7K8qgai10C\nrCNpp3zOHUbqIP9QOaJSn62LSS2mk6pUXlYjVcbWywM5P5dE+irtCeBrkoZIWorUV+mevO4lgKuB\nmyNirv5alWMtH0cbk74a/FEpzqak1rvflhYfT6qsVNL1LKkidHIOvx3YN59vi5E61lfS1VeeK+bq\nd5WX/z6p0r119OzPBalFbVWlnzKQ0leCnytse5iklXPYWqR+kkdU0iBpZ0mL52vqp4EvkV7zVfK0\neaXFSdLHSK+MK+tWLu+F8/Si+Ryr/HTDrnndCyr1v9wNuK5G3q0V2v0O0UNzB3r2gVqU1DF7Zh5O\nBBYtxN2RdPK+Qeon8ZnyOgpxzybd8IrD3r1sew/g/l7SOYWeX+FdDgzPYeMofHFXbR6p+foK0hPm\nNGD/QtzDgfP62E9nU/hahtS3YUphehVgTmF6FKny8Bzp5n9LIa9bkFo8ZgM3kiqDNxWWLfc76bHt\nJpW7SK8HXs7DcfTs1zEb2DyPH0R6Yn6T1PpzIj37k11Hod9HlW3tn/f5y7kMViiEBR98hfcS8F/k\nr+VymmbktEyj51eDq5Oe9F/Iy/0FWK+3/UVqPXkPuLI0v6/y2J90PrxK6gy/JT2/wlubD74WfADY\nsZfjpseyTSzPrXMe3iKdK+MKYacCp+bxT+V9/mbOb2XYvMZ6y8fienn9r5Bafy8CRuawSfT8YrEy\njCns5+l52w8De1TZ3mmkV3V1X7fy9Iqka8JL+Tj7E6lFsa48kyptc4p5Le2Dt0vLHloInwjcR7ou\nPkN6pVn5knS1nNc3SS2gB5XWfWM+bl4nfQiwayn8QNK19g3SOfKdQtg45r7GTs9hI/Ix+Wpe9730\nco56aM2gXDhm71P6VeblSE+pk+qIP570lcsipE7r1/exCJKmkCo5Z8xncq2JJC1MuvB/NCL+2e70\nmJl1qsHwA4PWoGjwV5kj4jrSzwXYABcR75C+zDIzs164D5SZmZlZg/wKz8zMzKxBboEyMzMza1BL\n+0ANHz48xo0b18pNmpmZmc2TO+6448WIGFEtrKUVqHHjxjF16tRWbtLMzMxsnkh6slaYX+GZmZmZ\nNcg/Y2BmVjLukCvrijf9mM/2HcnMupJboMzMzMwa5AqUmZmZWYNcgTIzMzNrkCtQZmZmZg1yBcrM\nzMysQa5AmZmZmTXIFSgzMzOzBvl3oMzMrG7+jSyzpM8WKElnSXpe0n2FecMkXSPp0fx36f5NppmZ\nmVnnqOcV3tnANqV5hwDXRcSqwHV52szMzGxQ6LMCFRE3AC+XZm8PnJPHzwF2aHK6zMzMzDrWvPaB\nGhkRM/P4c8DIWhEl7QfsBzBmzJh53JyZ2eDgPkZmA8N8f4UXEQFEL+GnR8QGEbHBiBEj5ndzZmZm\nZm03rxWoWZKWA8h/n29ekszMzMw627xWoC4DJuXxScClzUmOmZmZWefrsw+UpMnAlsBwSc8APwKO\nAS6StA/wJDCxPxNpZu3jPjm1DYR9U28azawxfVagImK3GkHjm5wWMzMzswHB/8rFzMzMrEGuQJmZ\nmZk1yP8Lz2yA6PT+No30tRkIaex03ZSXenX6OWCDi1ugzMzMzBrkCpSZmZlZg1yBMjMzM2uQ+0CZ\nWcdynxcz61RugTIzMzNrkCtQZmZmZg1yBcrMzMysQe4DZdaAgfBbR9Z5BuNvNtn88/Wms7kFyszM\nzKxBrkCZmZmZNcgVKDMzM7MGuQ+UGe6jYmZmjXELlJmZmVmDXIEyMzMza5ArUGZmZmYNch8os37S\n6f/Hzf2+rD91+vEPAyONzTYY89xf3AJlZmZm1iBXoMzMzMwa5AqUmZmZWYPcB2oQ8bvvwWEg9G1q\ndhoHQp6tuoFQdr52WjVugTIzMzNrkCtQZmZmZg1yBcrMzMysQe4D1QXa1YegP7Zbbx+CgdBvol7d\nlBezTtBN10TrXG6BMjMzM2uQK1BmZmZmDXIFyszMzKxBXdcHyr/XYWZmg027flut2ffSgXQPn68W\nKEnbSHpY0mOSDmlWoszMzMw62TxXoCQtCJwMbAusBewmaa1mJczMzMysU81PC9SGwGMR8XhEvAP8\nBti+OckyMzMz61yKiHlbUNoZ2CYi9s3TewIbRcSBpXj7AfvlydWBh+c9uQ0ZDrzYom11Gud98Bms\n+YbBm/fBmm9w3gdj3tuV77ERMaJaQL93Io+I04HT+3s7ZZKmRsQGrd5uJ3DeB1/eB2u+YfDmfbDm\nG5z3wZj3Tsz3/LzCmwGMLkyvkOeZmZmZdbX5qUDdDqwqaUVJCwO7Apc1J1lmZmZmnWueX+FFxBxJ\nBwJXAwsCZ0XE/U1L2fxr+WvDDuK8Dz6DNd8wePM+WPMNzvtg1HH5nudO5GZmZmaDlf+Vi5mZmVmD\nXIEyMzMza1DXVKAk7SLpfknvSar5qaOk6ZLulXSXpKmtTGN/aSDvXfevdyQNk3SNpEfz36VrxHs3\nl/ldkgbsxw59laGkRSRdmMNvlTSu9ansH3XkfW9JLxTKed92pLPZJJ0l6XlJ99UIl6QT8365R9L6\nrU5jf6gj31tKeq1Q3oe1Oo39RdJoSddLeiBf279ZJU7XlXud+e6cco+IrhiANUk/1DkF2KCXeNOB\n4e1Ob6vzTuroPw1YCVgYuBtYq91pb0LejwMOyeOHAMfWiDe73WltQl77LEPgAODUPL4rcGG7093C\nvO8NnNTutPZD3rcA1gfuqxE+AfgjIGBj4NZ2p7lF+d4SuKLd6eynvC8HrJ/HhwKPVDneu67c68x3\nx5R717RARcSDEdGqXznvKHXmvVv/9c72wDl5/Bxghzampb/VU4bF/XExMF6SWpjG/tKtx2+fIuIG\n4OVeomwPnBvJLcBSkpZrTer6Tx357loRMTMi7szjbwAPAsuXonVdudeZ747RNRWoBgTwZ0l35H8z\nM1gsDzxdmH6GDj4wGzAyImbm8eeAkTXiLSppqqRbJA3USlY9Zfh+nIiYA7wGLNOS1PWveo/fnfLr\njIslja4S3o269dyuxyaS7pb0R0lrtzsx/SG/hv8YcGspqKvLvZd8Q4eUe7//K5dmknQt8JEqQT+I\niEvrXM1mETFD0rLANZIeyk86Ha1JeR+Qest7cSIiQlKt3+UYm8t9JeAvku6NiGnNTqu11eXA5Ih4\nW9JXSS1xW7U5TdZ/7iSd17MlTQD+AKza5jQ1laTFgd8B34qI19udnlbpI98dU+4DqgIVEVs3YR0z\n8t/nJV1CejXQ8RWoJuR9wP7rnd7yLmmWpOUiYmZuvn6+xjoq5f64pCmkJ5uBVoGqpwwrcZ6RNARY\nEnipNcnrV33mPSKK+TyD1D9uMBiw5/b8KN5YI+IqSadIGh4RXfGPdiUtRKpEnB8Rv68SpSvLva98\nd1K5D6pXeJI+LGloZRz4NFD1C48u1K3/eucyYFIenwTM1RonaWlJi+Tx4cAngQdalsLmqacMi/tj\nZ+AvkXteDnB95r3U/2M7Uv+JweAyYK/8VdbGwGuF19pdS9JHKv37JG1Iup91w8MCOV9nAg9GxM9q\nROu6cq8n3x1V7u3uxd6sAdiR9A74bWAWcHWePwq4Ko+vRPp6527gftLrr7anvRV5z9MTSF81TOui\nvC8DXAc8ClwLDMvzNwDOyOObAvfmcr8X2Kfd6Z6P/M5VhsARwHZ5fFHgt8BjwG3ASu1OcwvzfnQ+\nr+8GrgfWaHeam5TvycBM4J/5PN8H2B/YP4cLODnvl3vp5SvkgTTUke8DC+V9C7Bpu9PcxLxvRuqv\new9wVx4mdHu515nvjil3/ysXMzMzswYNqld4ZmZmZs3gCpSZmZlZg1yBMjMzM2uQK1BmZmZmDXIF\nyszMzKxBrkCZmZmZNcgVKDMzM7MG/T+tPUoNSzp5WAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAASqElEQVR4nO3deZQeVZ3G8e+TBcKQAIMJYEKSFmUR\nGQWN4AzLRBaJiCI4ohJiZFhkXDkuDKIgMmziEQEHOERgAsOiCDKgjOMZkcgiW0BBWcKwhYSwCDGQ\nBWTJb/64t6G68m7V3W/eTvfzOadOv1W3llt1q+r91b33rVZEYGZmZmatG9bpDJiZmZmtaRxAmZmZ\nmVXkAMrMzMysIgdQZmZmZhU5gDIzMzOryAGUmZmZWUUOoAY5SY9J2r0f1vGipP9scf7dJS2TtLKv\n226wjUl5G8PbsX7rH5JC0ts6nY8qJO0saV6n82FmA5sDKGvVhyNiBoCkjSRdJmmRpOcl3Sxph+4Z\nI+LXETEaeLxdmYmIxyNidES81q5trMmUfFfSc3n4riTVmXdqDnaXFYaZqzvPnVIO8iLixojYspN5\nKpK0lqQr8oNMSJraZP4NJV0labmk+ZIOKKQdXSrnF3PZj83pp0paIOmFvOzRhWV3Li27LOfnYzn9\nk5Lm5XvCM5IulLReYfmLJT2Z1/2gpENK+d5N0gOSVki6XtLkQtoESVdLWixpoaTDS8vOytteKekz\npbRm+Zoj6aXCPtUMniVdUD5XGu2TpOmlY7UiL/+enN7wGm20TzYwOICy3hgN3AG8B9gQuBC4VtLo\n1bFxSSM6ufwa4jDgo8C7gHcCHwY+22D+RTkg7R4uXB2Z7A8aGrWQNwEHAk+1MO9ZwMvAxsB04BxJ\n7wCIiJOK5Qx8F5gTEc/mZc8HtoqI9YB/AKZL2i8ve2Np2b2BZcD/5GVvBnaMiPWBzYARwAmFfJ0M\ndOV1fwQ4oRBMjAV+BhxDuqfMBX5SWPZi4NG8Tx8CTpL0/kL63cDngLtqHI9m+QL4QmHfVgmeJe0E\nvLXGuuvuU0RcUjpenwMeKeSx2TXaaJ9sAHAANYRIWlvS6Uo1R4vy57UL6ftI+kN+mnpY0rRa64mI\nRyLitIh4MiJei4hZwFpAr5/a81PgyZJuz9u/WtKGOa0rP7kdLOlx4DeFaSPyPOMlXZOfUB+SdGhh\n3ccpPcFfLOkF4DM1tj9b0tmSfpmfFm+WtEk+Rn/JT8bbFeYfL+lKSX+W9KikLxXStpd0i6Ql+en0\n3yWtVUgPSYdL+r88z1nFJ89+MhP4fkQsjIgngO/X2u/+IOnQfMwX5zIYX5plL0mPSHpW0vckDcvL\nvU3Sb3PNwLOSflJY51aS/jevc56k/QtpsyWdI+m/JS0HvibpqWIgJWlfSffkz3XLQ9INeZG7c7l/\nQqlGbmFhXW/P5+cSSfdK+kgpL2dJulbSUkm3Sar1RdtrEfFyRJweETcBDWtcJa0LfAw4JiKW5WWu\nAWbUmFfAp0kPQN3bmhcRywuzrQTqNcHOBK7onj8iFhQCMXJeizV790bEX7tH89B9rPYD7o2In0bE\nS8BxwLvyeTAamAqcGBGvRMTdwBXAPxfWfVZEXAe8VM5ks3w1k+8xPwS+WGPdjfapbCZwUbzx7z8a\nXqON9skGiIjwMIgH4DFg9/z5eOBWYCNgHPA74N9y2vbA88AepMB6AulJtMc66mxjW9JFvn6DbR8A\n3NNgHXOAJ4BtgHWBK4GLc1oX6cZ0UU5bpzBtRJ7nBuBsYFTOz5+BXXPaccArpKe9YcA6NbY/G3iW\nVKs2CvgN6Yn308Bw0hPr9XneYcCdwLGkwHEz0pPlnjn9PcD7SE+6XcD9wBGFbQXwC2ADYFLO67Q6\nx+UAYEmDYVKd5Z4HdiiMTwGW1pl3KqnG4um8zz8A1m3x/No1H7d3A2uTvmhuKO3r9aRahUnAg8Ah\nOe0y4Jv5eI4CdsrT1wUWAAflY7hd3sbWhbJ6HtixsOzDwB6F7f4UOKpCebytdDwW5s8jgYeAo3NZ\n7wosBbYs5OU50vUzArgE+HGD49WoLI9q4XgvBKY2SN8OWFGa9jXg5zXm3YVUgzS6NP2oPD1I5/Wm\nNZZdNx+HqaXpO+WyCWA58IFS+tnAipx+V/e2gTOAc0rz/okUDI7J829USPsR8Psa+boJ+EyN6XXz\nRbr3/DmfYzfX2KevA2fUOlca7VNpnsmkwO0tVa/RevvkofNDxzPgoc0F3DOIeRjYq5C2J/BY/nwu\n8INm66iRth7wR+AbVZarMe8c4JTC+NakL/XhvBEsbVZI7542ApiYb05jCuknA7Pz5+MofKnX2f5s\n4EeF8S8C9xfG/w5Ykj/vADxeWv4bwH/UWfcRwFWF8SAHC3n8clr48qxY7q+RA+A8vnnermrMu0k+\n3sOAt5CC0XNb3M75wKmF8dGkYLWrsK/TCumfA67Lny8CZlH6ggY+AdxYmnYu8O1CWV1USj8BuCB/\nHkP6kpxcoTzqBVA7k5rNhhXSLwOOK+TlvELaXsAD/VmWpbw3C6B2Bp4qTTuU1ExXq+xm11mPSMHY\nd4rXVSF9BinYXuV8yukT8nW3RY204aSA5lvAyEJeTinNdzM5cCAFET8kBcvvBhYD82qsu2GwUStf\npOt5DOkBYCYpMHxrTptICqDXr3WuNNqnUvox5TKgxWu02T556NzgJryhZTwwvzA+P0+DdKN4uMrK\nJK0D/By4NSJO7of8LSjlbSQwtk560XhgcUQsLS0/oYVli54ufH6xxnh3H6/JwPjcpLNE0hJSDcXG\nAJK2kPSL3Kz0AnBSaT+gZ1+WFYV195dlpOC223rAssh35KKIeCoi7ouIlRHxKHAk6cm/FT3OqYhY\nRqqRqXfsi+fckaQv6ttz01h3k8xkYIfS8Z1OCvRqrRPgUmC/3CS9H3BXRMyHlsuj0f4tiIiVpX0o\n7l+7y7KKcrmTx4vXBpL+Bvg4hea7okh+Tzrvv1Njlpn0bI4qL/8EqW/Uj2ukvRapaXFT4F9azPd0\nUnC/ADiH1CdqIRXVyldE3BYRSyPir5H6/t1MCoQBTgeOj4jnm6y31j4V9WgqzVq+Rm1gcgA1tCwi\nfTl1m5SnQboxtdx3I39R/RfpJtaoc3IVE0t5e4VUrd6t3o1lEbChpDGl5Z9oYdneWAA8GhEbFIYx\nEdF90z0HeADYPFLn0qNJgUJlWvWXPOVhUp1F7yV1Tu32rjytFUHr94Ye51Tug/Mmeh77crkugtcD\nt0MjYjzpHDpb6RdOC4Dflo7v6IgofjH1KM+IuI8U2HyQ1Ox5aSG5L+WxCJjY3W+rsA9P1Jm/oSZl\neXTzNTT1IDBC0uaFabXKfl9SLc6cJusbQem+IGkiqZbuoqrLNkjvcb7m8+iteToRMT8i9o6IcRGx\nAykAvr3J9nubr+CN82M34Hs5+O4OlG9R4ZeNzdYtaUdSIH5Fad6+XKM2ADiAGlouA74laZzSr16O\nJT3JQapCP0jpp8TDlH42vFWtlUgaSboZvAjMLD2d98WBkrbOT8fHkzqoNn1NQUQsIPXnOlnSKEnv\nBA7mjX3rb7cDSyX9q6R1JA2XtI2k9+b0McALwLJ8DGs9kbYkSr/kqTHUe1XERcBXcjmOB75Kam5a\nhaT3S5qsZCJwCnB1IX22pJrLks6pgyRtm4Pqk4DbIuKxwjxfl/S3ed1fJv+6StLHJW2a5/kL6Ytr\nJal/2BaSZkgamYf3Snp7k8N1aV7/LqQ+UN2alcfTpH5stdxGqlU6MudjKunXUqvUrLSiSVmeVG85\npR+AjMqja+XzfJUgMFKH7p8Bx0taN3957wOU3+G2Sg1Svu4/m8tKkrYHPg9cV1p2BvC7iOhRY52D\n/Un582TgxO5llV598klJo/P1sifwqcK6rwK2kfSxvJ/HkvpMPpCXf7ukMUqvdDgQ+ABwWmHba+Xl\nBIzMx2dYC/naQNKeef4RkqaTzp/uXxZuQQpsts0DpPK/qoV9Kh7rK0s15NDkGm20TzZAdLoN0UN7\nB3r2gRoFnAk8mYczgVGFefcF7iFVmz/EG52iX19HHv9H0pfdClI1dPewc4NtTyf9yqZePueQ+i3d\nTvqy+zkwNqd1UegwXmsaqer8F6Sn6oeBwwvzHkfukN5g+7OBEwrjh1Dos0D61c6rhfHxpODhKdKX\n/62Ffd2FVOOxDLiRFAzeVFi23Oemx7b7qdwFnJqPx+L8WYX018sL+AqpRmUFqfbnTHr2J7sOOLTB\ntg7Px3xxLoNNC2kBfInUGfk50i+Nhue0U/N2l+XlDysstyVwLalz73OkTv3bNjpepJqhlcC1penN\nyuNw0vWwBNifQh+onP4O4LekTr/3Afs2OG96LNvP13GUhq6cdjTwy8K8G5Jqh5eT3sV2QGldE4BX\nWbUz9DBS4LA4H6sH87rLfXIeAA6ukccTSTXSy/PfWcCbctq4fAyXkK7vP5bPKWD3vO4XSfeDrkLa\nEflcWE7qEzSlxv2jfHymtpivO0j3vCWk63iPWmVQvnZb3KdROX23XlyjdffJw8AYlAvKrC6lF8u9\nmdTxdmYL8+9G+hXd2qRO69e3sMwcUpBzXh+za/1I6ef+dwPvjIhXOp0fM7OBYii8UND6KCq+lTnS\nu0s2aFN2bDWKiJeBZk1nZmZDjttTzczMzCpyE56ZmZlZRa6BMjMzM6totfaBGjt2bHR1da3OTZqZ\nmZn1yp133vlsRIyrlbZaA6iuri7mzp27OjdpZmZm1iuS5tdLcxOemZmZWUV+jYGZmbWs66hrW5rv\nsVM+1OacmHWWa6DMzMzMKnIAZWZmZlaRAygzMzOzihxAmZmZmVXkAMrMzMysIgdQZmZmZhU5gDIz\nMzOryAGUmZmZWUUOoMzMzMwq8pvIzWyN57djm9nq5hooMzMzs4ocQJmZmZlV5CY8M7NBzM2bZu3h\nGigzMzOzihxAmZmZmVXUtAlP0gXA3sAzEbFNnrYh8BOgC3gM2D8i/tK+bJq1l5s5zMysilZqoGYD\n00rTjgKui4jNgevyuJmZmdmQ0DSAiogbgMWlyfsAF+bPFwIf7ed8mZmZmQ1Yvf0V3sYR8WT+/BSw\ncb0ZJR0GHAYwadKkXm7OzMyKWm12NrP26HMn8ogIIBqkz4qIKRExZdy4cX3dnJmZmVnH9TaAelrS\nmwHy32f6L0tmZmZmA1tvA6hrgJn580zg6v7JjpmZmdnA18prDC4DpgJjJS0Evg2cAlwu6WBgPrB/\nOzNpNlBU6Xcy1F554GNTn4+N2eDTNICKiE/VSdqtn/NiZmZmtkbwm8jNzMzMKvI/EzYbZPxWdTOz\n9nMNlJmZmVlFDqDMzMzMKnITnpkNWJ1627bf8m1mzbgGyszMzKwiB1BmZmZmFbkJz6xNBsuv4QZT\nc9Zg2pf+5mNjVo1roMzMzMwqcgBlZmZmVpEDKDMzM7OK3AfKbA3hPiq2JhksfQDN6nENlJmZmVlF\nDqDMzMzMKnITnpmZdUyVpulWm/vcfGirg2ugzMzMzCpyAGVmZmZWkZvwrNfWhKp3/3JtYHK51Odj\nM/C4SdBqcQ2UmZmZWUUOoMzMzMwqchPeAOZq46GhU002bioy61++Zw8troEyMzMzq8gBlJmZmVlF\nDqDMzMzMKlJE9H5haRpwBjAcOC8iTmk0/5QpU2Lu3Lm93l4rOtmvY6D/tH6g58/MbCjo5L14oL/N\nfaD1I5N0Z0RMqZXW6xooScOBs4APAlsDn5K0dW/XZ2ZmZram6EsT3vbAQxHxSES8DPwY2Kd/smVm\nZmY2cPW6CU/SPwHTIuKQPD4D2CEivlCa7zDgsDy6JTCv99m1DhsLPNvpTFjbuZyHDpf10OBy7r3J\nETGuVkLb3wMVEbOAWe3ejrWfpLn12oJt8HA5Dx0u66HB5dwefWnCewKYWBjfNE8zMzMzG9T6EkDd\nAWwu6S2S1gI+CVzTP9kyMzMzG7h63YQXEa9K+gLwK9JrDC6IiHv7LWc2ELkpdmhwOQ8dLuuhweXc\nBn16D5SZmZnZUOQ3kZuZmZlV5ADKzMzMrCIHUNYySR+XdK+klZL8k9hBSNI0SfMkPSTpqE7nx9pD\n0gWSnpH0p07nxdpH0kRJ10u6L9+7v9zpPA0mDqCsij8B+wE3dDoj1v/875mGlNnAtE5nwtruVeCr\nEbE18D7g876m+48DKGtZRNwfEX6T/ODlf880RETEDcDiTufD2isinoyIu/LnpcD9wITO5mrwcABl\nZt0mAAsK4wvxzdZsUJDUBWwH3NbZnAwebf9XLrZmkfRrYJMaSd+MiKtXd37MzKxvJI0GrgSOiIgX\nOp2fwcIBlPUQEbt3Og/WMf73TGaDjKSRpODpkoj4WafzM5i4Cc/MuvnfM5kNIpIEnA/cHxGndTo/\ng40DKGuZpH0lLQT+HrhW0q86nSfrPxHxKtD975nuBy73v2canCRdBtwCbClpoaSDO50na4sdgRnA\nrpL+kIe9Op2pwcL/ysXMzMysItdAmZmZmVXkAMrMzMysIgdQZmZmZhU5gDIzMzOryAGUmZmZWUUO\noMzMzMwqcgBlZmZmVtH/A9iu/KOtB473AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "CPU times: user 14.3 s, sys: 1.57 s, total: 15.8 s\n",
            "Wall time: 12.9 s\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b4RG2YJQdHGA",
        "colab_type": "text"
      },
      "source": [
        "### Two batches of data per log-prob eval (2x slower)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eopN1CSgHrZY",
        "colab_type": "code",
        "outputId": "0872100d-9684-4cf5-aec1-7dbf87740996",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 476
        }
      },
      "source": [
        "%%time\n",
        "batches_per_eval = 2\n",
        "\n",
        "def dataset_fn(ctx):\n",
        "  batch_size = ctx.get_per_replica_batch_size(len(samples))\n",
        "  d = tf.data.Dataset.from_tensor_slices(samples).batch(batch_size // batches_per_eval)\n",
        "  return d.shard(ctx.num_input_pipelines, ctx.input_pipeline_id).prefetch(2)\n",
        "\n",
        "ds = st.experimental_distribute_datasets_from_function(dataset_fn)\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "def log_prob_and_grad(loc, scale_tril, observations, prev_sum_lp, prev_sum_grads):\n",
        "  with tf.GradientTape() as tape:\n",
        "    tape.watch((loc, scale_tril))\n",
        "    lp = tf.reduce_sum(dist.log_prob(loc, scale_tril, observations)) / len(samples)\n",
        "  grad = tape.gradient(lp, (loc, scale_tril))\n",
        "  return lp + prev_sum_lp, [g + pg for (g, pg) in zip(grad, prev_sum_grads)]\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "@tf.custom_gradient\n",
        "def target_log_prob(loc, scale_tril):\n",
        "  sum_lp = tf.zeros([])\n",
        "  sum_grads = [tf.zeros_like(x) for x in (loc, scale_tril)]\n",
        "  sum_lp, sum_grads = st.run(\n",
        "      lambda *x: tf.nest.map_structure(tf.identity, x), (sum_lp, sum_grads))\n",
        "  def reduce_fn(state, observations):\n",
        "    sum_lp, sum_grads = state\n",
        "    return st.run(\n",
        "        log_prob_and_grad, (loc, scale_tril, observations, sum_lp, sum_grads))\n",
        "  sum_lp, sum_grads = ds.reduce((sum_lp, sum_grads), reduce_fn)\n",
        "  sum_lp = st.reduce('sum', sum_lp, None)\n",
        "  sum_grads = [st.reduce('sum', sg, None) for sg in sum_grads]\n",
        "  return sum_lp, lambda grad_lp: [grad_lp * sg for sg in sum_grads]\n",
        "\n",
        "multibatch_vals = tfp.math.value_and_gradient(target_log_prob, (loc, scale_tril))\n",
        "\n",
        "kernel = tfp.mcmc.HamiltonianMonteCarlo(target_log_prob, step_size=.35, num_leapfrog_steps=2)\n",
        "kernel = tfp.mcmc.TransformedTransitionKernel(kernel, bijector=[tfb.Identity(), tfb.FillScaleTriL()])\n",
        "\n",
        "@tf.function(autograph=False)\n",
        "def sample_chain():\n",
        "  return tfp.mcmc.sample_chain(\n",
        "      num_results=200, num_burnin_steps=100,\n",
        "      current_state=[tf.ones_like(loc), tf.linalg.eye(scale_tril.shape[-1])], \n",
        "      kernel=kernel, trace_fn=lambda _, kr: kr.inner_results.is_accepted)\n",
        "samps, is_accepted = sample_chain()\n",
        "\n",
        "print(f'accept rate: {np.mean(is_accepted)}')\n",
        "print(f'ess: {tfp.mcmc.effective_sample_size(samps)}')\n",
        "\n",
        "print(tf.reduce_mean(samps[0], axis=0))\n",
        "# print(tf.reduce_mean(samps[1], axis=0))\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "for dim in range(ndim):\n",
        "  plt.figure(figsize=(10,1))\n",
        "  plt.hist(samps[0][:,dim], bins=50)\n",
        "  plt.title(f'loc[{dim}]: prior mean = 0.5, observation = {loc[dim]}')\n",
        "  plt.show()"
      ],
      "execution_count": 46,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "accept rate: 0.8\n",
            "ess: [<tf.Tensor: shape=(3,), dtype=float32, numpy=array([42.702564, 56.667793, 43.043278], dtype=float32)>, <tf.Tensor: shape=(3, 3), dtype=float32, numpy=\n",
            "array([[57.469193,       nan,       nan],\n",
            "       [34.606365, 43.471287,       nan],\n",
            "       [13.848625, 30.807411, 72.34709 ]], dtype=float32)>]\n",
            "tf.Tensor([-0.18598816  0.7396643   0.4074543 ], shape=(3,), dtype=float32)\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAV+0lEQVR4nO3debgcVZnH8e8PgolssiSywxXZwyPo\nMIAIGgTZd1CBCEFhEOdhxnFDEBeWOIIbwig6DEKCLKLIvigBhACCGhSIkTUYIRBIEAgEIus7f5zT\npFJ039t1l+57b/8+z1PP7a46VfX2qdPVb51TfVsRgZmZmZk1b4l2B2BmZmY21DiBMjMzM6vICZSZ\nmZlZRU6gzMzMzCpyAmVmZmZWkRMoMzMzs4qcQHUISbMk7dgP21go6WdNlj9c0gJJIWm9vuy7m31s\nJ+mBgdi29Q9JXbkNjGh3LFVIGi/p+nbHYWaDkxMoq2rPiDik9iR/OP5W0kuS7i8maRHx04hYdiCD\niYhbI2LDgdzHUCZppKRzJD0v6UlJn++m7GGSXs9Jb20a18Jw26ZekhcRF0TETu2MqyeSTpY0XdJr\nkk7ooawknSrpH3k6VZIKy0PSi4Vjf3Zh2XWldvGKpOl19vGhvJ2JhXkjJZ0m6QlJz0o6U9JSheUb\nS7pJ0nxJD0vat7TNpfM6T+cyU+vs922S7pM0uzBvu1LMtYu5/fPyTSX9Jm+34T9ElLS+pH9KOr80\n/2BJf891drmklQrLbs7r1Pb7QIV1z5c0J79nH5R0RKPYrL2cQFlfXQT8GVgZOB64RNKYVuy4rz0a\nkpbsr1gGsROA9YF1gO2BYyTt0k35OyJi2cJ0cwti7LOcHHTi+exh4BjgmibKHgnsA2wGvAfYE/h0\nqcxmhWP/5gd3ROxabBfA74BfFlfMSdHpwO9L2zwW2ALYFNgAeB/w1bzOCOAK4GpgpRzj+ZI2KKx/\nVl62cf77uTqv7UvAvOKMfHFVjHkPYAHw61zkVeAXwOF1tlf0I+CPpdc6Fvhf4BBgFeAl4MzSekcX\n9r9hhXW/BXRFxPLAXsBESf/SQ4zWBp14wul4+YrwB/mK8In8eGRh+d6S7s5XQDMbfeDmk9z7gG9E\nxMKI+BUwHdi/D7HNknScpL/mq9VzJY3Ky8ZJmi3py5KeBM6tzSusv3G++ntO0gxJexWWTZL0Y0nX\nSnqRlFCU93+zpImSfpevHK+StLKkC3J9/FFSV6H8RpKmSHpG0gOSPlZYtrukP+f1HlOhh0CLejwm\nSHo0XwUf39t668YE4OSIeDYi7gP+Dzisv3ciaQlJX81X1XMlnSfpHaVin8rtbY6kLxbW3VLStFxP\nT0n6fmHZ1vlYPCfpHhV6xPKx+qak20kfQl+SNK0U1+ckXZkfNzweQK1X47l83N+v1CN3W2Fb2+Tj\nPz//3aYUy8mSbpf0gqTrJY3udYU2KSImR8R1wAtNFJ8AfC8iZkfE48D36EVbyO1/O+C80qIvANcD\n95fm7wmcERHPRMQ84AzgU3nZRsDqwGkR8XpE3ATcTkoukLQRKYk4MiLm5TJ3leJ5F/AJUuLRnQnA\nJRHxIkBEPBARPwVmdPNaDwSeA24sLRoPXBURUyNiAfA1YD9Jy/UQQ4/rRsSMiHg5l408vbuJ7VqL\nOYHqTMcDWwObk65Gt2TRFeGWpBPjl4AVgA8CsxpsZyzwSEQUT9735Pl1STpW0tU9xDce2Jl00tig\nFlu2KukqdB3S1Wpx20sBV5FO4u8E/gO4QFJxiO9g4JvAcsBt1Hcg6QS+Ro7hDuDcvN/7gG/k/S0D\nTAEuzPs7EDhT0iZ5Oy8Ch5LqcXfgM5L2Ke1rW2BDYAfg65I2rhdQrrfnGk0N1lkRWI10TGq6PT7A\ne3My96Ckr6n5Xr7D8rQ9sC6wLPDDUpntSb1hOwFf1qLh3tOB0/MV97tJvQJIWoPUszKRVPdfBH6l\nxXs4DyG1g+WAnwAbSlq/sPxg0vGB7o/HB/PfFXKPwR3FwJWGWK4hffivDHwfuEbSyqV9fZLUFt6W\n430LSWt3dywlHVxvvX4wlp7bwlSlod5LixcKJYcCt0bErNoMSeuQkqKTGqyj0uM16yTYxeWb5sdb\nAn8HTsztcrryEFzB/wBfARY22F7tvXoAMLlRmTrrLE96PfWGvRery4iYCbxCOl/VfCvHfLsWHwrv\ncV2lIcuXSMnoHODaZuO21nEC1ZnGAydFxNx8RXgi+YqP1J19TkRMiYg3IuLxiChfUdYsC8wvzZtP\n+jCrKyJOiYg9eojvhxHxWEQ8Q0p2Diose4PU4/VyRJRPmFvnmE6JiFfy1ezVpfWviIjb82v7Z4P9\nnxsRMyNiPnAdMDMiboiI10jDFu/N5fYAZkXEuRHxWkT8GfgV8NH8Wm+OiOl5X/eShjs/VNrXibn3\n7h7SSXWzegHleluh0dTgddTuPyseo+6Oz1TSB9c7Sb2IB5ES6WaMB74fEY/kq+rjgANLCdiJEfFi\nREwnJaS14/IqsJ6k0RGxICLuzPM/AVwbEdfmOpwCTAN2K2xzUr5ify0frytq282J1EbAldD08Whk\nd+ChiPhZ3tdFpA+3PQtlzo2IB3O7/AXpAuUtIuLR7o5lRFxYb71+UH6/zgeWld68D+pDQBepzp4A\nrm6QQB8KTCrNOwP4Wj72Zb8GPitpjKRVgf/M85cGHgDmknoPl5K0U45j6VxmTVKbnE/qqToamFy7\n0FC6X2rJiLish9e+H/A0cEsP5YpOBn4aEbPrLOvp3Pdl0oXEGqQhyKskvbvJdYmIf8/PtwMuBV7G\nBh0nUJ1pddJVXc3f8zyAtYCZTW5nAbB8ad7yNDec0J3HCo+LsQHM6ybxWR14LCLeKK2/RoNtN/JU\n4fHCOs9rick6wFalnqDxpF4yJG2ldIP9PEnzgaOA8rDOk4XHLxW23R9qH2bFY9Tw+OTk5285wZhO\nuvo+oMl91WtTI0j3eNQ0Oq6Hk66+789DY7UEex3go6X63ZbUq1Zvm5B6m2qJ2cHA5RHxEjR9PJp9\nfbXXUGxbA3ksURqSrt2UvF0vNlF+vy4PLIhIvyifh5ReiYjngM8C7yLdd1SMYVtS+76kMG9PYLmI\nuLjBfr9Juk/ybtK9U5eTkuanIuJV0n1Zu5Pq7wuk5LOWtCzMZSfm2G4BfgvslHuVvs2ihKw7E4Dz\naq+1J5I2B3YETmtQpNtzX0T8PiJeyBd6k0nDkrs1s25NHq68jZREfqaZuK21nEB1pidIH041a+d5\nkD6Qmh1vnwGsWxr334xu7ilo0loNYoN0P0AjTwBrafGbidcGHm9y/aoeA24p9R4sGxG1k92FpN6P\ntSLiHaQhJjXaWHckfUVv/UbRm1O9dSLiWVL3f7FXq8rxiQrx1mtTr7F48ln3uEbEQxFxEKnn61TS\nFxGWIdXvz0r1u0xEnFKKsWgKMCZ/AB7EouE76P549NQuyq+v9hoer1O2W3kIr+GxlDS+3noRMTYW\n3ZR8a9X9ko57lbZQ7/hPAC4t9TTtAGyRh/6eBD4O/JekK3LcCyPi6IhYIyLWBf4B3FW70ImIeyPi\nQxGxckTsTOq5+UPe9r0N4oI0HNwF3Jr3eymwWo6jq1ZY0lrAON56z1Z3xuVtP5q3/UVgf0l/yssX\nq0tJ6wIjgQcbbK9Yl1XXHYHvgRqcIsJTB0yk+5h2zI8nkq4Ex5CuwG8jXeFBuufgOdJJcQnSFfZG\n5W0Utnsn8F1gFLBvXndMqUwA61WIczrpqmulHNt/52XjgNml8m/OI9138gjpWz9L5WUvFOKfVHud\n3ez/ZuCIwvOJpGGi2vMdgYfz4+VIvRCH5P0tBfwrsHFePheYUKjXucD5+XlXrpcRjfbdT8f9FNKw\nxYqkoZk5wC4Nyu4KrJIfbwT8hTRcWozvhAbrHgE8ROq1WJbUQ1F+rReQhmbG5rrYKS//RK3N5Pr9\nJ/B2UsL1JOl+uCVzGxsHrNldfQE/JiVSc0v1293xWBp4HdigUP4w4Lb8eGVS2z6Y9IH28fx8dIN2\n8+a6A/y+XirXy4W5rY4iDWnVK3sU6R6+NUg9ajOAo/KysaQhxyXz8fsBaXhtqcL6bycNNX24tN3l\nSL1SteliUs/NSnl5bX8iDbM/Vjv2efl7ctxLkxKVvwEjC6/vYdKN1iOAD5Df0/l5cb/7kRLdVYt1\nQLo/amqd+lDe7yak9jmqsN+lS9v+LqlNjynU1/OkIbZlgPOBn+dlK5Da7Kgc43jS/XcbNLFu7V7K\nZfOx2Dmvu9dAtyVPvXj/tTsATy060IsnUKNI9yzMydMZwKhC2X1JV34v5JPXzuVtFMp2kT48FuYT\n7o519v1mApVPZtf1EOdxwF9JH1CTgaXzsnF0k0Dl52NJCcP8vI19C8sm0Y8JVH6+Ienm4nmkK+ub\ngM3zsgNICdYLpHuxfkjrE6iRwDn5hP0U8PnCsrVJwwlr5+ffzWVeJCWiJ7H4B+hM4CMN9rME8HXS\nh+M80ofCiqXXeiTpA+5J4JjCuueTkpkFpA/1fQrLtsrH85m83WsK8datL9IHUwA/Ks1veDzy8pPy\nPp4jfdAfRiEJIg0f3pXb1l3Att20m8XWHcD39SQWfVOrNh1WqIcFhbIiDXk9k6dvA8rLPkx6/76Y\nj8XlwPqlfR2U609NxDSx8Lz2RZSX8j7Gl8p/B3g2H//rKF1skd7Td+TYFntPd3cuKMy/Hzi8zvxa\nuyxOsxps+4RiW8nzDgYezXFdwaKEcQzp3x68kNvSnZTeNz2se0te73nSxeS/DXQ78tS7qfbmMeuR\n0j+DWw24LCImNFH+k6Qr0VHAJhHxSBPrzCJ9EN3Qx3CtH0laE/hFRGzTY2Ezsw7gBMoGFSdQZmY2\nFPgmcjMzM7OK3ANlZmZmVpF7oMzMzMwq6tOPsVY1evTo6OrqauUuzczMzHrlrrvuejoixtRb1tIE\nqquri2nTpvVc0MzMzKzNJJV/geBNHsIzMzMzq6ilPVBmZu3Udew1TZWbdcruAxyJmQ117oEyMzMz\nq8gJlJmZmVlFTqDMzMzMKnICZWZmZlaREygzMzOzipxAmZmZmVXkBMrMzMysIidQZmZmZhU5gTIz\nMzOryAmUmZmZWUVOoMzMzMwqcgJlZmZmVpETKDMzM7OKnECZmZmZVeQEyszMzKwiJ1BmZmZmFTmB\nMjMzM6toRE8FJJ0D7AHMjYhN87yVgIuBLmAW8LGIeHbgwjQza52uY69pqtysU3Yf4EjMbLBqpgdq\nErBLad6xwI0RsT5wY35uZmZm1hF6TKAiYirwTGn23sDk/HgysE8/x2VmZmY2aPU4hNfAKhExJz9+\nElilUUFJRwJHAqy99tq93J2Z2eAznIb6htNrMWuFPt9EHhEBRDfLz4qILSJiizFjxvR1d2ZmZmZt\n19sE6ilJqwHkv3P7LyQzMzOzwa23Q3hXAhOAU/LfK/otIjMb9jxc1HeuQ7P26rEHStJFwB3AhpJm\nSzqclDh9RNJDwI75uZmZmVlH6LEHKiIOarBoh36OxczMzGxI6O0Qnpl1CA8V9V2zdQiuR7Ohwj/l\nYmZmZlaREygzMzOzijyEZ4aHqcw6kd/31hfugTIzMzOryAmUmZmZWUVOoMzMzMwq8j1QZtYvqnxV\nv7+36XtUGhuI49Kf/C8ebKhyD5SZmZlZRU6gzMzMzCryEJ7ZENHfQzHDaThksA9TVTGcXovZcOYe\nKDMzM7OKnECZmZmZVeQhPLMO5aEi6412thu32eFvKH3z1j1QZmZmZhU5gTIzMzOryEN41hKd+A0y\nDzeYWT1DaZjKGnMPlJmZmVlFTqDMzMzMKvIQXg+G0+809Xe38VD4Nk5/HxMPy5l1Hg+5WT3ugTIz\nMzOryAmUmZmZWUV9GsKTtAtwOrAkcHZEnNIvUfXBUOhqHezfSBtOw1TD6bWY2eDW3+ebgTh/9fct\nGp08bNnrHihJSwI/AnYFNgEOkrRJfwVmZmZmNlj1ZQhvS+DhiHgkIl4Bfg7s3T9hmZmZmQ1eioje\nrSgdAOwSEUfk54cAW0XE0aVyRwJH5qcbAg/0PtyWGg083e4gBiHXS32ul/pcL/W5XupzvTTmuqlv\noOtlnYgYU2/BgP8bg4g4CzhroPfT3yRNi4gt2h3HYON6qc/1Up/rpT7XS32ul8ZcN/W1s176MoT3\nOLBW4fmaeZ6ZmZnZsNaXBOqPwPqS3iXpbcCBwJX9E5aZmZnZ4NXrIbyIeE3S0cBvSP/G4JyImNFv\nkbXfkBt2bBHXS32ul/pcL/W5XupzvTTmuqmvbfXS65vIzczMzDqV/xO5mZmZWUVOoMzMzMwqcgKV\nSfqOpPsl3SvpMkkrNCi3i6QHJD0s6dhWx9lqkj4qaYakNyQ1/KqopFmSpku6W9K0VsbYDhXqpdPa\ny0qSpkh6KP9dsUG513NbuVvSsP3ySU/HX9JISRfn5b+X1NX6KFuviXo5TNK8Qhs5oh1xtpqkcyTN\nlfSXBssl6Yxcb/dKel+rY2yHJuplnKT5hfby9VbE5QRqkSnAphHxHuBB4LhygQ79+Zq/APsBU5so\nu31EbN4h/6ukx3rp0PZyLHBjRKwP3Jif17Mwt5XNI2Kv1oXXOk0e/8OBZyNiPeA04NTWRtl6Fd4X\nFxfayNktDbJ9JgG7dLN8V2D9PB0J/LgFMQ0Gk+i+XgBuLbSXk1oQkxOomoi4PiJey0/vJP1fq7KO\n+/maiLgvIobKf49vmSbrpePaC+n1Tc6PJwP7tDGWdmvm+Bfr6xJgB0lqYYzt0Invi6ZExFTgmW6K\n7A2cF8mdwAqSVmtNdO3TRL20hROo+j4FXFdn/hrAY4Xns/M8gwCul3RX/vke68z2skpEzMmPnwRW\naVBulKRpku6UNFyTrGaO/5tl8gXcfGDllkTXPs2+L/bPw1SXSFqrzvJO1InnlGa9X9I9kq6TNLYV\nOxzwn3IZTCTdAKxaZ9HxEXFFLnM88BpwQStja6dm6qUJ20bE45LeCUyRdH++ahiy+qlehp3u6qX4\nJCJCUqP/k7JObi/rAjdJmh4RM/s7VhuyrgIuioiXJX2a1Ev34TbHZIPXn0jnlAWSdgMuJw1zDqiO\nSqAiYsfulks6DNgD2CHq/4OsYfnzNT3VS5PbeDz/nSvpMlI3/ZBOoPqhXjquvUh6StJqETEnDy3M\nbbCNWnt5RNLNwHuB4ZZANXP8a2VmSxoBvAP4R2vCa5se6yUiinVwNvDtFsQ1FAzLc0pfRcTzhcfX\nSjpT0uiIGNAfX/YQXiZpF+AYYK+IeKlBMf98TR2SlpG0XO0xsBPpJutO14nt5UpgQn48AXhLT52k\nFSWNzI9HAx8A/tqyCFunmeNfrK8DgJsaXLwNJz3WS+m+nr2A+1oY32B2JXBo/jbe1sD8wpB5x5K0\nau3eQUlbknKbgb8QiQhP6Xz1MGls+e48/STPXx24tlBuN9K39GaShnLaHvsA18u+pHH2l4GngN+U\n6wVYF7gnTzNcLx3dXlYmffvuIeAGYKU8fwvg7Px4G2B6bi/TgcPbHfcA1sdbjj9wEulCDWAU8Mt8\n/vkDsG67Yx4k9fKtfC65B/gtsFG7Y25RvVwEzAFezeeXw4GjgKPycpG+wTgzv3e2aHfMg6Reji60\nlzuBbVoRl3/KxczMzKwiD+GZmZmZVeQEyszMzKwiJ1BmZmZmFTmBMjMzM6vICZSZmZlZRU6gzMzM\nzCpyAmVmZmZW0f8Ds9XTtDkl2RkAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAW6klEQVR4nO3debgcVZ3G8e8LYVPCEhIigSzs6wgi\nsgnIEB6FqCwCkUUII4jI4IaPiOggIrLOKDKAwAACAhFEkVURkMiiLIFh3wNhCSHsS4QBA7/545yG\nupXue7uTvt19+76f56nnVtU5VXVOnVpOnTrVVxGBmZmZmdVvgXYnwMzMzGygcQXKzMzMrEGuQJmZ\nmZk1yBUoMzMzswa5AmVmZmbWIFegzMzMzBrkClSXkTRd0tZNWMdbkn5dZ/ytJc2W9N78bruXbYzJ\n21iwP9ZvzSEpJK3S7nQ0QtLmkh5udzrMbGBxBcpq+XxE7FmZkPQTSfdKmiPp8GLEiLg2IhYHnuqv\nxETEUxGxeES821/bGMiUHCvppTwcK0k14m6ZK7uzC8OkVqe5XcqVvIi4MSJWb2eayiSNl/SQpDcl\nXS9pbI14y0qaLOlZSa9JulnSRjXinlXOu6Rxkq6S9Iqk5ySdJGlIDltN0qWSXpD0sqSrJa1eWHYR\nST/P235F0imSFiqEzy4N70r67yrpOiyna+vCvGGSLszH8ouSzpe0RD15lnRoabtv5eN9eA7/T0mP\nSnoj7+O9SunZStKdkl6X9Lik/QphkvQDSU/l8N9U0pXDz5b0Tmn7CxbCJ0p6MG/7AUk7FMLWyfv4\nRUlz/UCjpPMkzczbfUTSvtXK2VrHFSir12PAwcCVrd5w5YLeruUHiP2AHYB1gY8Cnwe+2kv8Z3OF\ntDKc04pENkO3t0LmG/3vgf8AhgFTgQtrRF8cuB34eI57DnClpMVL69wMWLnK8qcAzwPLAesBnwIO\nyGFLAZcBqwMjgduASwvLHgJsAKwDrAasD/ywElg8voCPAG8Bvy2la2VgF2BmKV1HAksDK+Z0jwQO\nryfPEXFUadvHAlMi4sW8/D9I58eSwCTgF5I2zelZCLgEOC2HfxH4maR187J7AXsCnwRGAYsB5Urh\ncaVz69287uWB84CDgCWA7wIXSFo2L/dP4CJgH6o7GhgXEUsA2wFHSvp4jbjWChHhoYsGYDqwdR5f\nBDgBeDYPJwCLFOJuD9wFvA5MA7Ypr6PK+s8DDu9r23WkcwrpgnBb3v6lwLAcNg4I0oXkKeCGwrwh\nOc4o0sX9ZVLl7iuFdR8OXJzT+jqwb5Xtn026efwRmA3cTLrInwC8AjwEfKwQfxTwO+AF4AngG4Ww\nDYG/A6+SbgQnAQsXwgPYH3g0xzkZUJPL/W/AfoXpfYBbasTdEnhmPrb1lbzPX85lMKqU128AjwMv\nAscDC+SwVYC/Aq/lsAsLy60BXJPX+TAwsVRWvwSuIt38vgc8ByxYiLMjcE9f5ZGPpcjrmU26QfbY\nH8Ca+fh8Fbgf2K6UlpNJDxJvALcCKze5LPcD/laY/jCp8rFGncu/Dny8MD0E+F9SxTqAVQphDwIT\nCtPHA6fVWO+wvPwyeXoqsEshfHfg6RrLTsrHhErz/wRMoHTtIJ2XBxSm/x24ut48F+Yrb3dSL8te\nBnwnj4/MefxQIfx2YLc8fjHw3ULYpsD/VeLn4+PIGtvZCHi+NO8FYJPSvFWA6KOMV8/H9sTe4nno\n38EtUN3tB8DGpCfLdUk3lh8CSNoQOJf0FLQUsAXpItYUknaXdE8f0fYCvkx6+p0DnFgK/xTpZvaZ\nKsv+BniGVLHZGThK0laF8O1JF7ulgPNrbH8iaX8MB94m3XTvzNMXAz/LeVkAuBy4G1geGA98S1Il\nXe8C387LbZLDD6CnzwGfIN3EJtbIU2W/vdrLMKZGXtbO6au4O8+rZVlJsyQ9kV/DfLiXuMX0bUWq\n+E4klduTpLIo2pHUMrE+qRy+nOf/BPgzqWVhBfKTe972NcAFwLLArsApktYqrHN34KfAUOAXpArQ\nVqXwC/J4zfKIiC1ynHUjtQ70aNnJLRCX53QuC3wdOL/46iqn78c5H4/ldFXVR1keUmOxHmUZEf8g\nPeD0Vp6V7a0HLJzTVfFt4IaIqHY+ngDsKulDuYVkW1KlppotgOci4qXiJkvjK0hassqyk4BzI9/9\nc1p3Ad6OiKuqxD8Z+JykpSUtDexEqlTNpUaeKzYnlePvaiy7GOm8vB8gImYBk4F/k7SgpE2AscBN\npXwWxxcBVi3MO0DplecdknYqzJ8KPChpu7zuHUjXnb6uk8X0niLpTdID3kzSQ4W1S7trcB6aO9Cz\nBWoaPZ8uPwNMz+OnAT/vax1VwprZAnVMYXot4B1gQT5obVqpEF6ZNwQYTbpJDi2EHw2cnccPJ90w\netv+2cD/FKa/DjxYmP4X4NU8vhHwVGn57wO/qrHubwGXFKYD2KwwfRFwSJPL/V0KLRSkC3pQpaWL\n1NK2FukV/oqkVpmqrQ5Vlj2T9IqiMr046dXDuEJetymEHwBcl8fPBU4HViit84vAjaV5pwE/KpTV\nuaXwI4Gz8vhQUoVqbAPlUWyF2ZLcAkW64T5HbjXL8yZXjvmcljMKYROAh5pclmcWz40872Zg7z6W\nWwK4F/h+Yd5oUsViyRp5XxO4g/QAEzl/1Y6ZFYAZ5JaYQhncDIzIx9SteR3LlZYdm4/PFQvzhpJa\nZCvHzXR6tkCNAq4F3svDNRRadXvLc5V9eXYv++wcUoVRhXmfB2blfTKHnq3b+wKPkK5HS5Jar4Lc\nikR6aFiGdJ2aQGql/GRh+X1ILZ9zgDeBz1ZJU68tUKRr5Gakh7+FmnnseWhscAtUdxtFaiGoeDLP\ng3RhndbyFPX0dGH8SWAhUqtBtfCiUcDLEfFGafnl61i2aFZh/K0q05V+JGOBUcXWA+BQUnN/pbPt\nFUqdcF8HjirlA9JNueLNwrqbZTbpZlKxBDA78hW3KCKei4gHIuK9iHiC1Ldtp3K8GnocUxExG3iJ\n2vu+eMwdTHpiv03S/ZIqLVNjgY1K+3cP0k252johtTZ9QdIiwBeAOyPiSai7PHrL39MR8V4pD8X8\ntbosydNvVIkLvN+Scjnpte3RhaATgCMi4rUqyyxAqjz8nvSacDipVe3YUrwRpBa5UyJiciHop6RX\ng3eRXiH/gVSZLp5HkPoM3ZSPtYrDgV9HxPQaWbqIVFEZSsr7NNLDWz15roR/iNS/qmr/PknHk/pv\nTaycJ5LWILWo7kVq1VobOFjSZ/NiZ5Eq1FNIrVbX5/nPAETEnRHxUkTMidSydj7p+ESpk/xxpAr7\nwqQW9jNyC1rdIuLdiLiJVKn9WiPLWnO5AtXdniXdnCrG5HmQbkjVOpW20ujC+BjSxffFwry5bv7Z\ns8AwSUNLy8+oY9l58TTwREQsVRiGRsSEHP5LUpP6qpE6eB5Kz2b+uknaQ3N/vVQcar3Cu5/0mrZi\n3TyvHkH914Iex1R+/bYMPfd9uVyfhfcrbl+JiFGkDu6nKH0R9jTw19L+XTwiijeHHuUZEQ+QKjbb\n0vP1HcxfeTwLjM6Vi2IeZtSI36s+yvLQGov1KMu8j1emRnnmSuQfSDfx8ocD44Hjc2WyUvH7u6Td\nSX2axgAnRcTbkV7N/YrUclJZ99KkytNlEdHjVWVEvBURB0bE8hGxEqkifUep8gmpMlKuxIwHvlFI\n12jgIknfy+HrkVpF/5Er6aeW0tVbnit2JPWpm1IOkPRj0rHz6Yh4vRC0DvBIRFydHzAeJvV32zbn\n+b2I+FFEjIuIFUhlMoPax0fwwbG3HqllfGpez+2kVrt5/emXIbT/Gj6ouQLV3SYDP5Q0In/Zcxgf\nPMWdSXrPP17SApKWz09fVUlaSNKipGNmiKRFNf9fQ31J0lr5SfEI4OKo42cKIuJp0hPv0TkdHyU1\njZ/X+5Lz7DbgDUnfk7RY7r+wjqRP5PChpE6ss/M+nOenwog4P3p+wVMeav1UxLnAQbkcRwHfIb2O\nmYukf5U0Vslo4BgKX1cpfYpddVk+6B+yXr6JHQXcWmpJ+G7uuzIa+Cb5CzJJu0haIcd5hXRzeQ+4\nAlhN0p75OFtI0ickrdnH7rogr38Len7d1Vd5zAJWqrHOW0mtSgfndGxJeqVT7udVlz7K8qgai10C\nrCNpp3zOHUbqIP9QOaJSn62LSS2mk6pUXlYjVcbWywM5P5dE+irtCeBrkoZIWorUV+mevO4lgKuB\nmyNirv5alWMtH0cbk74a/FEpzqak1rvflhYfT6qsVNL1LKkidHIOvx3YN59vi5E61lfS1VeeK+bq\nd5WX/z6p0r119OzPBalFbVWlnzKQ0leCnytse5iklXPYWqR+kkdU0iBpZ0mL52vqp4EvkV7zVfK0\neaXFSdLHSK+MK+tWLu+F8/Si+Ryr/HTDrnndCyr1v9wNuK5G3q0V2v0O0UNzB3r2gVqU1DF7Zh5O\nBBYtxN2RdPK+Qeon8ZnyOgpxzybd8IrD3r1sew/g/l7SOYWeX+FdDgzPYeMofHFXbR6p+foK0hPm\nNGD/QtzDgfP62E9nU/hahtS3YUphehVgTmF6FKny8Bzp5n9LIa9bkFo8ZgM3kiqDNxWWLfc76bHt\nJpW7SK8HXs7DcfTs1zEb2DyPH0R6Yn6T1PpzIj37k11Hod9HlW3tn/f5y7kMViiEBR98hfcS8F/k\nr+VymmbktEyj51eDq5Oe9F/Iy/0FWK+3/UVqPXkPuLI0v6/y2J90PrxK6gy/JT2/wlubD74WfADY\nsZfjpseyTSzPrXMe3iKdK+MKYacCp+bxT+V9/mbOb2XYvMZ6y8fienn9r5Bafy8CRuawSfT8YrEy\njCns5+l52w8De1TZ3mmkV3V1X7fy9Iqka8JL+Tj7E6lFsa48kyptc4p5Le2Dt0vLHloInwjcR7ou\nPkN6pVn5knS1nNc3SS2gB5XWfWM+bl4nfQiwayn8QNK19g3SOfKdQtg45r7GTs9hI/Ix+Wpe9730\nco56aM2gXDhm71P6VeblSE+pk+qIP570lcsipE7r1/exCJKmkCo5Z8xncq2JJC1MuvB/NCL+2e70\nmJl1qsHwA4PWoGjwV5kj4jrSzwXYABcR75C+zDIzs164D5SZmZlZg/wKz8zMzKxBboEyMzMza1BL\n+0ANHz48xo0b18pNmpmZmc2TO+6448WIGFEtrKUVqHHjxjF16tRWbtLMzMxsnkh6slaYX+GZmZmZ\nNcg/Y2BmVjLukCvrijf9mM/2HcnMupJboMzMzMwa5AqUmZmZWYNcgTIzMzNrkCtQZmZmZg1yBcrM\nzMysQa5AmZmZmTXIFSgzMzOzBvl3oMzMrG7+jSyzpM8WKElnSXpe0n2FecMkXSPp0fx36f5NppmZ\nmVnnqOcV3tnANqV5hwDXRcSqwHV52szMzGxQ6LMCFRE3AC+XZm8PnJPHzwF2aHK6zMzMzDrWvPaB\nGhkRM/P4c8DIWhEl7QfsBzBmzJh53JyZ2eDgPkZmA8N8f4UXEQFEL+GnR8QGEbHBiBEj5ndzZmZm\nZm03rxWoWZKWA8h/n29ekszMzMw627xWoC4DJuXxScClzUmOmZmZWefrsw+UpMnAlsBwSc8APwKO\nAS6StA/wJDCxPxNpZu3jPjm1DYR9U28azawxfVagImK3GkHjm5wWMzMzswHB/8rFzMzMrEGuQJmZ\nmZk1yP8Lz2yA6PT+No30tRkIaex03ZSXenX6OWCDi1ugzMzMzBrkCpSZmZlZg1yBMjMzM2uQ+0CZ\nWcdynxcz61RugTIzMzNrkCtQZmZmZg1yBcrMzMysQe4DZdaAgfBbR9Z5BuNvNtn88/Wms7kFyszM\nzKxBrkCZmZmZNcgVKDMzM7MGuQ+UGe6jYmZmjXELlJmZmVmDXIEyMzMza5ArUGZmZmYNch8os37S\n6f/Hzf2+rD91+vEPAyONzTYY89xf3AJlZmZm1iBXoMzMzMwa5AqUmZmZWYPcB2oQ8bvvwWEg9G1q\ndhoHQp6tuoFQdr52WjVugTIzMzNrkCtQZmZmZg1yBcrMzMysQe4D1QXa1YegP7Zbbx+CgdBvol7d\nlBezTtBN10TrXG6BMjMzM2uQK1BmZmZmDXIFyszMzKxBXdcHyr/XYWZmg027flut2ffSgXQPn68W\nKEnbSHpY0mOSDmlWoszMzMw62TxXoCQtCJwMbAusBewmaa1mJczMzMysU81PC9SGwGMR8XhEvAP8\nBti+OckyMzMz61yKiHlbUNoZ2CYi9s3TewIbRcSBpXj7AfvlydWBh+c9uQ0ZDrzYom11Gud98Bms\n+YbBm/fBmm9w3gdj3tuV77ERMaJaQL93Io+I04HT+3s7ZZKmRsQGrd5uJ3DeB1/eB2u+YfDmfbDm\nG5z3wZj3Tsz3/LzCmwGMLkyvkOeZmZmZdbX5qUDdDqwqaUVJCwO7Apc1J1lmZmZmnWueX+FFxBxJ\nBwJXAwsCZ0XE/U1L2fxr+WvDDuK8Dz6DNd8wePM+WPMNzvtg1HH5nudO5GZmZmaDlf+Vi5mZmVmD\nXIEyMzMza1DXVKAk7SLpfknvSar5qaOk6ZLulXSXpKmtTGN/aSDvXfevdyQNk3SNpEfz36VrxHs3\nl/ldkgbsxw59laGkRSRdmMNvlTSu9ansH3XkfW9JLxTKed92pLPZJJ0l6XlJ99UIl6QT8365R9L6\nrU5jf6gj31tKeq1Q3oe1Oo39RdJoSddLeiBf279ZJU7XlXud+e6cco+IrhiANUk/1DkF2KCXeNOB\n4e1Ob6vzTuroPw1YCVgYuBtYq91pb0LejwMOyeOHAMfWiDe73WltQl77LEPgAODUPL4rcGG7093C\nvO8NnNTutPZD3rcA1gfuqxE+AfgjIGBj4NZ2p7lF+d4SuKLd6eynvC8HrJ/HhwKPVDneu67c68x3\nx5R717RARcSDEdGqXznvKHXmvVv/9c72wDl5/Bxghzampb/VU4bF/XExMF6SWpjG/tKtx2+fIuIG\n4OVeomwPnBvJLcBSkpZrTer6Tx357loRMTMi7szjbwAPAsuXonVdudeZ747RNRWoBgTwZ0l35H8z\nM1gsDzxdmH6GDj4wGzAyImbm8eeAkTXiLSppqqRbJA3USlY9Zfh+nIiYA7wGLNOS1PWveo/fnfLr\njIslja4S3o269dyuxyaS7pb0R0lrtzsx/SG/hv8YcGspqKvLvZd8Q4eUe7//K5dmknQt8JEqQT+I\niEvrXM1mETFD0rLANZIeyk86Ha1JeR+Qest7cSIiQlKt3+UYm8t9JeAvku6NiGnNTqu11eXA5Ih4\nW9JXSS1xW7U5TdZ/7iSd17MlTQD+AKza5jQ1laTFgd8B34qI19udnlbpI98dU+4DqgIVEVs3YR0z\n8t/nJV1CejXQ8RWoJuR9wP7rnd7yLmmWpOUiYmZuvn6+xjoq5f64pCmkJ5uBVoGqpwwrcZ6RNARY\nEnipNcnrV33mPSKK+TyD1D9uMBiw5/b8KN5YI+IqSadIGh4RXfGPdiUtRKpEnB8Rv68SpSvLva98\nd1K5D6pXeJI+LGloZRz4NFD1C48u1K3/eucyYFIenwTM1RonaWlJi+Tx4cAngQdalsLmqacMi/tj\nZ+AvkXteDnB95r3U/2M7Uv+JweAyYK/8VdbGwGuF19pdS9JHKv37JG1Iup91w8MCOV9nAg9GxM9q\nROu6cq8n3x1V7u3uxd6sAdiR9A74bWAWcHWePwq4Ko+vRPp6527gftLrr7anvRV5z9MTSF81TOui\nvC8DXAc8ClwLDMvzNwDOyOObAvfmcr8X2Kfd6Z6P/M5VhsARwHZ5fFHgt8BjwG3ASu1OcwvzfnQ+\nr+8GrgfWaHeam5TvycBM4J/5PN8H2B/YP4cLODnvl3vp5SvkgTTUke8DC+V9C7Bpu9PcxLxvRuqv\new9wVx4mdHu515nvjil3/ysXMzMzswYNqld4ZmZmZs3gCpSZmZlZg1yBMjMzM2uQK1BmZmZmDXIF\nyszMzKxBrkCZmZmZNcgVKDMzM7MG/T+tPUoNSzp5WAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABlCAYAAACCyt1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAASqElEQVR4nO3deZQeVZ3G8e+TBcKQAIMJYEKSFmUR\nGQWN4AzLRBaJiCI4ohJiZFhkXDkuDKIgMmziEQEHOERgAsOiCDKgjOMZkcgiW0BBWcKwhYSwCDGQ\nBWTJb/64t6G68m7V3W/eTvfzOadOv1W3llt1q+r91b33rVZEYGZmZmatG9bpDJiZmZmtaRxAmZmZ\nmVXkAMrMzMysIgdQZmZmZhU5gDIzMzOryAGUmZmZWUUOoAY5SY9J2r0f1vGipP9scf7dJS2TtLKv\n226wjUl5G8PbsX7rH5JC0ts6nY8qJO0saV6n82FmA5sDKGvVhyNiBoCkjSRdJmmRpOcl3Sxph+4Z\nI+LXETEaeLxdmYmIxyNidES81q5trMmUfFfSc3n4riTVmXdqDnaXFYaZqzvPnVIO8iLixojYspN5\nKpK0lqQr8oNMSJraZP4NJV0labmk+ZIOKKQdXSrnF3PZj83pp0paIOmFvOzRhWV3Li27LOfnYzn9\nk5Lm5XvCM5IulLReYfmLJT2Z1/2gpENK+d5N0gOSVki6XtLkQtoESVdLWixpoaTDS8vOytteKekz\npbRm+Zoj6aXCPtUMniVdUD5XGu2TpOmlY7UiL/+enN7wGm20TzYwOICy3hgN3AG8B9gQuBC4VtLo\n1bFxSSM6ufwa4jDgo8C7gHcCHwY+22D+RTkg7R4uXB2Z7A8aGrWQNwEHAk+1MO9ZwMvAxsB04BxJ\n7wCIiJOK5Qx8F5gTEc/mZc8HtoqI9YB/AKZL2i8ve2Np2b2BZcD/5GVvBnaMiPWBzYARwAmFfJ0M\ndOV1fwQ4oRBMjAV+BhxDuqfMBX5SWPZi4NG8Tx8CTpL0/kL63cDngLtqHI9m+QL4QmHfVgmeJe0E\nvLXGuuvuU0RcUjpenwMeKeSx2TXaaJ9sAHAANYRIWlvS6Uo1R4vy57UL6ftI+kN+mnpY0rRa64mI\nRyLitIh4MiJei4hZwFpAr5/a81PgyZJuz9u/WtKGOa0rP7kdLOlx4DeFaSPyPOMlXZOfUB+SdGhh\n3ccpPcFfLOkF4DM1tj9b0tmSfpmfFm+WtEk+Rn/JT8bbFeYfL+lKSX+W9KikLxXStpd0i6Ql+en0\n3yWtVUgPSYdL+r88z1nFJ89+MhP4fkQsjIgngO/X2u/+IOnQfMwX5zIYX5plL0mPSHpW0vckDcvL\nvU3Sb3PNwLOSflJY51aS/jevc56k/QtpsyWdI+m/JS0HvibpqWIgJWlfSffkz3XLQ9INeZG7c7l/\nQqlGbmFhXW/P5+cSSfdK+kgpL2dJulbSUkm3Sar1RdtrEfFyRJweETcBDWtcJa0LfAw4JiKW5WWu\nAWbUmFfAp0kPQN3bmhcRywuzrQTqNcHOBK7onj8iFhQCMXJeizV790bEX7tH89B9rPYD7o2In0bE\nS8BxwLvyeTAamAqcGBGvRMTdwBXAPxfWfVZEXAe8VM5ks3w1k+8xPwS+WGPdjfapbCZwUbzx7z8a\nXqON9skGiIjwMIgH4DFg9/z5eOBWYCNgHPA74N9y2vbA88AepMB6AulJtMc66mxjW9JFvn6DbR8A\n3NNgHXOAJ4BtgHWBK4GLc1oX6cZ0UU5bpzBtRJ7nBuBsYFTOz5+BXXPaccArpKe9YcA6NbY/G3iW\nVKs2CvgN6Yn308Bw0hPr9XneYcCdwLGkwHEz0pPlnjn9PcD7SE+6XcD9wBGFbQXwC2ADYFLO67Q6\nx+UAYEmDYVKd5Z4HdiiMTwGW1pl3KqnG4um8zz8A1m3x/No1H7d3A2uTvmhuKO3r9aRahUnAg8Ah\nOe0y4Jv5eI4CdsrT1wUWAAflY7hd3sbWhbJ6HtixsOzDwB6F7f4UOKpCebytdDwW5s8jgYeAo3NZ\n7wosBbYs5OU50vUzArgE+HGD49WoLI9q4XgvBKY2SN8OWFGa9jXg5zXm3YVUgzS6NP2oPD1I5/Wm\nNZZdNx+HqaXpO+WyCWA58IFS+tnAipx+V/e2gTOAc0rz/okUDI7J829USPsR8Psa+boJ+EyN6XXz\nRbr3/DmfYzfX2KevA2fUOlca7VNpnsmkwO0tVa/RevvkofNDxzPgoc0F3DOIeRjYq5C2J/BY/nwu\n8INm66iRth7wR+AbVZarMe8c4JTC+NakL/XhvBEsbVZI7542ApiYb05jCuknA7Pz5+MofKnX2f5s\n4EeF8S8C9xfG/w5Ykj/vADxeWv4bwH/UWfcRwFWF8SAHC3n8clr48qxY7q+RA+A8vnnermrMu0k+\n3sOAt5CC0XNb3M75wKmF8dGkYLWrsK/TCumfA67Lny8CZlH6ggY+AdxYmnYu8O1CWV1USj8BuCB/\nHkP6kpxcoTzqBVA7k5rNhhXSLwOOK+TlvELaXsAD/VmWpbw3C6B2Bp4qTTuU1ExXq+xm11mPSMHY\nd4rXVSF9BinYXuV8yukT8nW3RY204aSA5lvAyEJeTinNdzM5cCAFET8kBcvvBhYD82qsu2GwUStf\npOt5DOkBYCYpMHxrTptICqDXr3WuNNqnUvox5TKgxWu02T556NzgJryhZTwwvzA+P0+DdKN4uMrK\nJK0D/By4NSJO7of8LSjlbSQwtk560XhgcUQsLS0/oYVli54ufH6xxnh3H6/JwPjcpLNE0hJSDcXG\nAJK2kPSL3Kz0AnBSaT+gZ1+WFYV195dlpOC223rAssh35KKIeCoi7ouIlRHxKHAk6cm/FT3OqYhY\nRqqRqXfsi+fckaQv6ttz01h3k8xkYIfS8Z1OCvRqrRPgUmC/3CS9H3BXRMyHlsuj0f4tiIiVpX0o\n7l+7y7KKcrmTx4vXBpL+Bvg4hea7okh+Tzrvv1Njlpn0bI4qL/8EqW/Uj2ukvRapaXFT4F9azPd0\nUnC/ADiH1CdqIRXVyldE3BYRSyPir5H6/t1MCoQBTgeOj4jnm6y31j4V9WgqzVq+Rm1gcgA1tCwi\nfTl1m5SnQboxtdx3I39R/RfpJtaoc3IVE0t5e4VUrd6t3o1lEbChpDGl5Z9oYdneWAA8GhEbFIYx\nEdF90z0HeADYPFLn0qNJgUJlWvWXPOVhUp1F7yV1Tu32rjytFUHr94Ye51Tug/Mmeh77crkugtcD\nt0MjYjzpHDpb6RdOC4Dflo7v6IgofjH1KM+IuI8U2HyQ1Ox5aSG5L+WxCJjY3W+rsA9P1Jm/oSZl\neXTzNTT1IDBC0uaFabXKfl9SLc6cJusbQem+IGkiqZbuoqrLNkjvcb7m8+iteToRMT8i9o6IcRGx\nAykAvr3J9nubr+CN82M34Hs5+O4OlG9R4ZeNzdYtaUdSIH5Fad6+XKM2ADiAGlouA74laZzSr16O\nJT3JQapCP0jpp8TDlH42vFWtlUgaSboZvAjMLD2d98WBkrbOT8fHkzqoNn1NQUQsIPXnOlnSKEnv\nBA7mjX3rb7cDSyX9q6R1JA2XtI2k9+b0McALwLJ8DGs9kbYkSr/kqTHUe1XERcBXcjmOB75Kam5a\nhaT3S5qsZCJwCnB1IX22pJrLks6pgyRtm4Pqk4DbIuKxwjxfl/S3ed1fJv+6StLHJW2a5/kL6Ytr\nJal/2BaSZkgamYf3Snp7k8N1aV7/LqQ+UN2alcfTpH5stdxGqlU6MudjKunXUqvUrLSiSVmeVG85\npR+AjMqja+XzfJUgMFKH7p8Bx0taN3957wOU3+G2Sg1Svu4/m8tKkrYHPg9cV1p2BvC7iOhRY52D\n/Un582TgxO5llV598klJo/P1sifwqcK6rwK2kfSxvJ/HkvpMPpCXf7ukMUqvdDgQ+ABwWmHba+Xl\nBIzMx2dYC/naQNKeef4RkqaTzp/uXxZuQQpsts0DpPK/qoV9Kh7rK0s15NDkGm20TzZAdLoN0UN7\nB3r2gRoFnAk8mYczgVGFefcF7iFVmz/EG52iX19HHv9H0pfdClI1dPewc4NtTyf9yqZePueQ+i3d\nTvqy+zkwNqd1UegwXmsaqer8F6Sn6oeBwwvzHkfukN5g+7OBEwrjh1Dos0D61c6rhfHxpODhKdKX\n/62Ffd2FVOOxDLiRFAzeVFi23Oemx7b7qdwFnJqPx+L8WYX018sL+AqpRmUFqfbnTHr2J7sOOLTB\ntg7Px3xxLoNNC2kBfInUGfk50i+Nhue0U/N2l+XlDysstyVwLalz73OkTv3bNjpepJqhlcC1penN\nyuNw0vWwBNifQh+onP4O4LekTr/3Afs2OG96LNvP13GUhq6cdjTwy8K8G5Jqh5eT3sV2QGldE4BX\nWbUz9DBS4LA4H6sH87rLfXIeAA6ukccTSTXSy/PfWcCbctq4fAyXkK7vP5bPKWD3vO4XSfeDrkLa\nEflcWE7qEzSlxv2jfHymtpivO0j3vCWk63iPWmVQvnZb3KdROX23XlyjdffJw8AYlAvKrC6lF8u9\nmdTxdmYL8+9G+hXd2qRO69e3sMwcUpBzXh+za/1I6ef+dwPvjIhXOp0fM7OBYii8UND6KCq+lTnS\nu0s2aFN2bDWKiJeBZk1nZmZDjttTzczMzCpyE56ZmZlZRa6BMjMzM6totfaBGjt2bHR1da3OTZqZ\nmZn1yp133vlsRIyrlbZaA6iuri7mzp27OjdpZmZm1iuS5tdLcxOemZmZWUV+jYGZmbWs66hrW5rv\nsVM+1OacmHWWa6DMzMzMKnIAZWZmZlaRAygzMzOzihxAmZmZmVXkAMrMzMysIgdQZmZmZhU5gDIz\nMzOryAGUmZmZWUUOoMzMzMwq8pvIzWyN57djm9nq5hooMzMzs4ocQJmZmZlV5CY8M7NBzM2bZu3h\nGigzMzOzihxAmZmZmVXUtAlP0gXA3sAzEbFNnrYh8BOgC3gM2D8i/tK+bJq1l5s5zMysilZqoGYD\n00rTjgKui4jNgevyuJmZmdmQ0DSAiogbgMWlyfsAF+bPFwIf7ed8mZmZmQ1Yvf0V3sYR8WT+/BSw\ncb0ZJR0GHAYwadKkXm7OzMyKWm12NrP26HMn8ogIIBqkz4qIKRExZdy4cX3dnJmZmVnH9TaAelrS\nmwHy32f6L0tmZmZmA1tvA6hrgJn580zg6v7JjpmZmdnA18prDC4DpgJjJS0Evg2cAlwu6WBgPrB/\nOzNpNlBU6Xcy1F554GNTn4+N2eDTNICKiE/VSdqtn/NiZmZmtkbwm8jNzMzMKvI/EzYbZPxWdTOz\n9nMNlJmZmVlFDqDMzMzMKnITnpkNWJ1627bf8m1mzbgGyszMzKwiB1BmZmZmFbkJz6xNBsuv4QZT\nc9Zg2pf+5mNjVo1roMzMzMwqcgBlZmZmVpEDKDMzM7OK3AfKbA3hPiq2JhksfQDN6nENlJmZmVlF\nDqDMzMzMKnITnpmZdUyVpulWm/vcfGirg2ugzMzMzCpyAGVmZmZWkZvwrNfWhKp3/3JtYHK51Odj\nM/C4SdBqcQ2UmZmZWUUOoMzMzMwqchPeAOZq46GhU002bioy61++Zw8troEyMzMzq8gBlJmZmVlF\nDqDMzMzMKlJE9H5haRpwBjAcOC8iTmk0/5QpU2Lu3Lm93l4rOtmvY6D/tH6g58/MbCjo5L14oL/N\nfaD1I5N0Z0RMqZXW6xooScOBs4APAlsDn5K0dW/XZ2ZmZram6EsT3vbAQxHxSES8DPwY2Kd/smVm\nZmY2cPW6CU/SPwHTIuKQPD4D2CEivlCa7zDgsDy6JTCv99m1DhsLPNvpTFjbuZyHDpf10OBy7r3J\nETGuVkLb3wMVEbOAWe3ejrWfpLn12oJt8HA5Dx0u66HB5dwefWnCewKYWBjfNE8zMzMzG9T6EkDd\nAWwu6S2S1gI+CVzTP9kyMzMzG7h63YQXEa9K+gLwK9JrDC6IiHv7LWc2ELkpdmhwOQ8dLuuhweXc\nBn16D5SZmZnZUOQ3kZuZmZlV5ADKzMzMrCIHUNYySR+XdK+klZL8k9hBSNI0SfMkPSTpqE7nx9pD\n0gWSnpH0p07nxdpH0kRJ10u6L9+7v9zpPA0mDqCsij8B+wE3dDoj1v/875mGlNnAtE5nwtruVeCr\nEbE18D7g876m+48DKGtZRNwfEX6T/ODlf880RETEDcDiTufD2isinoyIu/LnpcD9wITO5mrwcABl\nZt0mAAsK4wvxzdZsUJDUBWwH3NbZnAwebf9XLrZmkfRrYJMaSd+MiKtXd37MzKxvJI0GrgSOiIgX\nOp2fwcIBlPUQEbt3Og/WMf73TGaDjKSRpODpkoj4WafzM5i4Cc/MuvnfM5kNIpIEnA/cHxGndTo/\ng40DKGuZpH0lLQT+HrhW0q86nSfrPxHxKtD975nuBy73v2canCRdBtwCbClpoaSDO50na4sdgRnA\nrpL+kIe9Op2pwcL/ysXMzMysItdAmZmZmVXkAMrMzMysIgdQZmZmZhU5gDIzMzOryAGUmZmZWUUO\noMzMzMwqcgBlZmZmVtH/A9iu/KOtB473AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 720x72 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "CPU times: user 1min 8s, sys: 8.71 s, total: 1min 17s\n",
            "Wall time: 51.5 s\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R6II-MMbjl8k",
        "colab_type": "text"
      },
      "source": [
        "Sanity check logprob and gradients."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fSKHFj2f8w3V",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for i, (sv, mv) in enumerate(zip(tf.nest.flatten(singleton_vals), \n",
        "                                 tf.nest.flatten(multibatch_vals))):\n",
        "  np.testing.assert_allclose(sv, mv, err_msg=i, rtol=1e-5)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3cibdyv4kZj2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}