04_continuous_normalizing_flows.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous normalizing flows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/homebrew/Caskroom/miniforge/base/envs/torch-mps/lib/python3.9/site-packages/equinox/_ad.py:753: UserWarning: As of Equinox 0.10.7, `equinox.filter_custom_vjp.defvjp` is deprecated in favour of `.def_fwd` and `.def_bwd`. This new API supports symbolic zeros, which allow for more efficient autodifferentiation rules. In particular:\n",
      "- the fwd and bwd functions take an extra `perturbed` argument, which     indicates which primals actually need a gradient. You can use this     to skip computing the gradient for any unperturbed value. (You can     also safely just ignore this if you wish.)\n",
      "- `None` was previously passed to indicate a symbolic zero gradient for     all objects that weren't inexact arrays, but all inexact arrays     always had an array-valued gradient. Now, `None` may also be passed     to indicate that an inexact array has a symbolic zero gradient.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "from functools import partial\n",
    "\n",
    "import jax\n",
    "\n",
    "from jax import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "\n",
    "import jax.numpy as np\n",
    "import flax.linen as nn\n",
    "import optax\n",
    "import diffrax as dfx\n",
    "import math\n",
    "from tensorflow_probability.substrates import jax as tfp\n",
    "\n",
    "from sklearn import datasets, preprocessing\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import trange"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.0, 2.0)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGiCAYAAADulWxzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPhElEQVR4nO3df3RV1Zk//ncguRci+UUMScCQBrSEFgVEgaR+JXyMEnRZow7jr1XBsWgd6FeE1oJTYaHjJ6NicepY0WkVnak/ylTsYBUHosEvGkCDqFCIBSLhRxJsyG80AXK+fzje9uznSbJzc29yc/J+rZW1vNtzzt0399ybzX728+wox3EcEBEREXnQoL7uABEREVG4cKBDREREnsWBDhEREXkWBzpERETkWRzoEBERkWdxoENERESexYEOEREReRYHOkRERORZHOgQERGRZ3GgQ0RERJ4V1oFOUVERLr74YsTFxWHEiBEoLCxEeXl5l+etW7cO2dnZGDJkCM4//3y88cYb4ewmEREReVRYBzpbtmzBggULsG3bNmzatAmnTp3CFVdcgZaWlg7Pef/993HTTTfh9ttvx0cffYTCwkIUFhZi9+7d4ewqEREReVBUb27q+cUXX2DEiBHYsmULLr30UvWYG264AS0tLXj99dcDbdOnT8ekSZOwZs2a3uoqEREReUB0bz5ZQ0MDAGD48OEdHlNaWorFixe72mbNmoXXXntNPb61tRWtra2Bx+3t7Thx4gSSk5MRFRXV804TERFR2DmOg6amJowcORKDBoUu4NRrA5329nYsWrQI3/ve9zBhwoQOj6uurkZqaqqrLTU1FdXV1erxRUVFWLlyZUj7SkRERH3j8OHDOOecc0J2vV4b6CxYsAC7d+/G1q1bQ3rdZcuWuWaAGhoaMHr0aBw+fBjx8fEhfS4iIiIKj8bGRmRkZCAuLi6k1+2Vgc7ChQvx+uuv49133+1ylJaWloaamhpXW01NDdLS0tTj/X4//H6/aI+Pj+dAh4iIqJ8J9bKTsGZdOY6DhQsXYv369Xj77beRlZXV5Tk5OTkoLi52tW3atAk5OTnh6iYRERF5VFhndBYsWIAXX3wRf/jDHxAXFxdYZ5OQkIChQ4cCAG699VaMGjUKRUVFAIC7774bM2bMwGOPPYarrroKL7/8Mj788EM888wz4ewqEREReVBYZ3SeeuopNDQ0IC8vD+np6YGfV155JXBMZWUlqqqqAo9zc3Px4osv4plnnsHEiRPxX//1X3jttdc6XcBMREREpOnVOjq9obGxEQkJCWhoaOAaHSIion4iXH+/udcVEREReRYHOkRERORZHOgQERGRZ3GgQ0RERJ7FgQ4RERF5Fgc6RERE5Fkc6BAREZFncaBDREREnsWBDhEREXkWBzpERETkWWHd1LO/uHzQnL7uAnUiOuVs0Xb6i7/0+bUi5Tmjx53rvnb5/l7vAxFRVza1r+uT5+WMDhEREXkWBzpERETkWQxdUcQLZYilJ9cywz/hvpZtuEkLVdmcZ9OvnlyLiCgScEaHiIiIPIsDHSIiIvIsDnSIiIjIs7hGh/qlvkiXNq+vrmcZnijbTtQH94TKtdQPrHGczZodwG49ji2u7SGiSMUZHSIiIvIsDnSIiIjIsxi6oogX0srIRhXhjthUFz49dpQ8cdvHoqmt4GLRFltW0WW/6i6UrzvuJaVfZohr+kTZh0SfPM/oQ0/YpscTEfU2zugQERGRZ3GgQ0RERJ7FgQ4RERF5FtfoUMQL5fYFJ7OSRFtTZoxoS7VICY+uaxFtp5XjzPU4AHBySlaX14+r+FI2KutvYPSjKWuo3bVCiLujE1Gk4owOEREReRYHOkRERORZDF1Rv6SGSizSvX318pjkjR/I8yzS0GvyRoi2OCU0VjVBhsYyNhx3PVZDakoI6vg02Y+Mje7UcS1MpV0r6YC8lvk71MJz4f7SYBiMiEKJMzpERETkWRzoEBERkWcxdEW9xjYkISoQK8doYSotzALjWjUTzxKHpB6Q/Tp4swxLjXnRHW5qi5NPp2VwnVKOM69vXhsAWpNkvxL3yn+b+OrdoarK2fI1jn5T/m60zC+z/8my0DOitYw0bTNT4z1Ss8Fe2mZ1LfOLiqEsIrLFGR0iIiLyLA50iIiIyLM40CEiIiLP4hod6jW2O1zbrL9QqxInybUpbUbadvLHynnKep+zjsnnNHcTz3j+M3HMwYXfFm0xTfJaybtPuc9T1gRpTiW0i7aWke7XnVZ6Shxz4Ea5PmbcU/L3HGsUcdZ2UEeWTL1P2imvZe6Yrq3H4Q7nRBRuYZ3Reffdd3H11Vdj5MiRiIqKwmuvvdbp8SUlJYiKihI/1dXV4ewmEREReVRYBzotLS2YOHEinnzyyW6dV15ejqqqqsDPiBF2/9olIiIi+lthDV3Nnj0bs2fP7vZ5I0aMQGJiYug7RH3KtppxtJFefLp8v7yYluKshK5akwYbj+02vPQ1yhBRdf4Z1+PmUTJMpYWWtH9PxFbUuRty5GBeC5/hmLxW/Xj3c2op7u1xMpylhaXqsrv+t4+Wqq6FDcXv/qbp4hgt5KW+3wbbkBfT0IkoIhcjT5o0Cenp6bj88svx3nvvdXpsa2srGhsbXT9EREREQIQNdNLT07FmzRr8/ve/x+9//3tkZGQgLy8PO3fu7PCcoqIiJCQkBH4yMjJ6scdEREQUyaIcx3F65YmiorB+/XoUFhZ267wZM2Zg9OjR+I//+A/1/7e2tqK1tTXwuLGxERkZGWhoaEB8fLzVc1w+aE63+kRdi1Y2xdRCEm0FF4s2M+ShhVO08IlaedcIS9UqlZG1CsdaNWMze2rYURmmav9BrWg7/YcU0dYy0v1Yq4x8+GoZzvIpGVwnprrDUsP2ytCVxuZa5/17m9W1NFWXuH/X5kamgL6ZaWxZhWhTKy8btPuLG4QSRY5N7es6/f+NjY1ISEjo1t9vGxGfXj516lRs3bq1w//v9/vh9/t7sUdERETUX0RU6Eqza9cupKen93U3iIiIqB8K64xOc3Mz9u//63RyRUUFdu3aheHDh2P06NFYtmwZjh49ihdeeAEA8PjjjyMrKwvf/e538dVXX+HXv/413n77bfzP//xPOLtJREREHhXWgc6HH36ImTNnBh4vXrwYADB37lysXbsWVVVVqKysDPz/trY2LFmyBEePHkVsbCwuuOACbN682XUN6h9s10ugXq4BMVPAbdevaOtq6rLd60SS9sl1NW1xcmJTS+02U879dWfkMQ/EyhMh1xPFHXJXDd57r1yrMv4R+bpr8uTrzhRhb5lKrq17qX1ePuewt831RPL9aVh+UrSdFOfJNTnaa0zbPFi0xVqsx1HX9pxQ7i9tJ3SlzSalnYj6p7AOdPLy8tDZWue1a9e6Ht9777249957w9klIiIiGkAifo0OERERUbAiPuuK+qcmyyq42uaccds+dj2uU66lhak0ZqhKS1XXNsHUqgu3xXf97wJfvWzTXqO54aUMP+nVhrVNSc20eu33XHOdrOLc3Nwsn/Mc9+/L7CcAJN+jVEbOE0048KA7jJe2XoaptPCfFp5LLXGHwdQUdI1WQZuIBhTO6BAREZFncaBDREREnsWBDhEREXkW1+hQ902fKJrMtRzajuBaSrC2FsZn7KrdPMou/dvcVkHrR/Moue6ldoLsQ/pWuQ5FnKdsJ6G1pb56VLS1Gq9RW1ej0X6H5mu0XdsTd0hWFPfVd71lBqD04ZBc59S2Y5j7sVLR3XbN1MGb3et2xryodEuhpqGbO8dDfhFymwgi7+CMDhEREXkWBzpERETkWQxdUae0XchxQIZitGqzJi1VWQs3JX9shr3kruSa+vFdP6e2g7YZFtHO05gpz0AHFXaV36FNuEkLG2khqGD5lKrU0cZ7G5eYJftlGerzNbmP095rLUx1aI48bthe2WYqv0tWRh67aJtoO62EXhFkZWTujk4U+TijQ0RERJ7FgQ4RERF5FkNX1CktpFKlbKiZ8fxnrsdaJV4t7HLWMRkGqZx9lnGM7FebUhl53FMylGTSXs/oN5VwkFGdGQBq78x1PfbVy2tpoQw188es7KuE/pKfln3QmOFF7TVq1Zm1qsEnp7hDVVoF4qZM+d5WXSKf0wxnmRuZAnqYSqsS3ZTpfqxVTx77snyNamhJXl4cZxt+YpiKKPJxRoeIiIg8iwMdIiIi8iwOdIiIiMizuEZnABNrOyxTbIcdbRdt2poc04EbZZp4xkalom6ce92ObUr1YW3tkJJObjJTqgEAytqO5KffdzdoacoKrRKvubZDS8+3Tl0219po76N2rbGjRJvNruBaFWTtvLp89/2l7VSur8eR67ZMNru4A0DzJcpasd2y/+YXofa713CNDlHk44wOEREReRYHOkRERORZDF0NEOpUvBHyUFOjLSoEA4CvUYazTOMekiEVM7wB6KExk5Yu7WuSKc5mardWDVgLP2i/CxHqU1LQ1RCRRUjQNmyonmsRPlGP0doswpm+ehny0t5HM1TVmjS4y2MAvXSAGYI88GCsOGbs/fL1NI+S4UyNTVq9VgJA/R0SUUThjA4RERF5Fgc6RERE5Fkc6BAREZFncY1OP2ebBmtFWYOgrauwUZctx9Bt8TLVN/XVz0SbmaqupQ37lbVDIv1boe3Grm0JoLFaa9PP12xYrRVS1ibFoetUe1+9bNN2iTe3E1GvtUOuvam7UK7b0VLJq3O6Tl8fU5HY5TFE1D9wRoeIiIg8iwMdIiIi8iyGrvo561CJUsXXDBtoKbVJO+WltJ2jzfRyZ3yzPGafDC2Yab0AkFrSdTVjjU0ISgvNWFcgpo4p4SwtTGhDq9hshi9PniNLEMT+n1rR9qVy/dFX/lm0tRVc7O6D7Q7wRBTxOKNDREREnsWBDhEREXkWQ1cDhLpxpREiKP8nGWrI2Cgr1zZfJsNSw9e7w1Knj8kwlVYFV2VuUqlkg6mhhSBDUMGGqbwY8jJfU49ej/k+arQsLyXMmrTZfZyWwXfySIpVt07fKY+rH+8OhcU0yPsrrVRmiMWe6Drrsb/fE0T9HWd0iIiIyLM40CEiIiLP4kCHiIiIPItrdPo5dcdxJWVb7tktxR7Rxr1yXU36s37R1ureJBwxDXbX0iov+8w1DZa7i/f2Wggvrr0I5WsKeu2Ttp7MYJYzAIDkj2UyedUlcq2NJnGv+36NOyQrKsdW1MkTtR3NTR68T4j6E87oEBERkWdxoENERESeFdbQ1bvvvotHH30UZWVlqKqqwvr161FYWNjpOSUlJVi8eDH27NmDjIwM/PznP8e8efPC2c2IEGwoRgtTaRWOtWqzpmFHZRhJCy01j5LjY3PzxDEvKn1X0o1jtdRxM71YqbrrxbBRfxHsvapuqKqkl6vXMu4JrVRB7UQZpvI1yUu1xcm2U6JNbvzZlCkrgmub0or+K+ny2j1NROER1hmdlpYWTJw4EU8++aTV8RUVFbjqqqswc+ZM7Nq1C4sWLcIPf/hDvPXWW+HsJhEREXlUWGd0Zs+ejdmzZ1sfv2bNGmRlZeGxxx4DAIwfPx5bt27F6tWrMWvWrHB1k4iIiDwqorKuSktLkZ+f72qbNWsWFi1a1OE5ra2taG1tDTxubGwMV/fCSpuu10IE4hglE0QLU5kbeAJAdY57ej5pn8xkSdqphSRkv0RGihKm0sJsvnolH4zT+hEt2LChFqayZtwTsUoYLFZGbHH4ahluGnZU3udt8e7JbW1jWW0z27p82Y+4l4zfTw/uZ5vvAIZxiToXUYuRq6urkZqa6mpLTU1FY2MjvvxS24cYKCoqQkJCQuAnIyOjN7pKRERE/UBEDXSCsWzZMjQ0NAR+Dh8+3NddIiIioggRUaGrtLQ01NTUuNpqamoQHx+PoUOHquf4/X74/bKAHREREVFEDXRycnLwxhtvuNo2bdqEnJycPupReNim51qt21HWwlRZrkvQ1uSYyu+SfR37sgwjmruJN10oz9NSgrUquKe77BUNNGJXdWW9j5a+nrFBrrWpU+5Nk3k/A8CJqbJactpmWX6hreBi12Ot3IP6nBZrbbgeh6j7whq6am5uxq5du7Br1y4AX6eP79q1C5WVlQC+DjvdeuutgeN/9KMf4eDBg7j33nuxb98+/OpXv8Lvfvc73HPPPeHsJhEREXlUWAc6H374ISZPnozJkycDABYvXozJkydj+fLlAICqqqrAoAcAsrKy8Mc//hGbNm3CxIkT8dhjj+HXv/41U8uJiIgoKFGO4zh93YlQamxsREJCAhoaGhAfH291zuWD5oS5V6EjpvC18JYyha9N15sVjrUqsskft9j1ywxBaRWPLdOLbV4jkckMGQF62Egrc3C4wB2Cih55Uhzj2zFMtimfGZGaroSXeU/TQLSpfV2n/z+Yv982+n3WFREREVFHONAhIiIiz+JAh4iIiDwrotLLyc2m/Lu2M3Jdlqw5dHyaPDVxr/txy0h5zP9zxx7R9tnN3+q6X8q6BFtcv0DBsE3jFtuVAMjYmGQcI9fo7L1X1uvKVJYcmOvhkjbXi2Nsd3Inop7jjA4RERF5Fgc6RERE5FkMXfUBq5BUB8zdkrXdxdsmyqqu0SObRVs9Yo1j5HT9//fMVNGWeuIz0WYVblLCbNypnELFNuSpfulluUNX2k7lmetkZeTaCTGiLXm3+zhth/OkzUplZ8uK6UTUPZzRISIiIs/iQIeIiIg8i6GrPmBbzVhjboypbT5YP15u1hmzV1Z1zSh1T7FrmSZ1F8aKNqv+a1lX3MCTQqQnYR7tuNiKRNdjX738XDUp2Yyn4uT1zXBW+qr35UFa1pXyWTa/oBnKIuo+zugQERGRZ3GgQ0RERJ7FgQ4RERF5FtfohJntWgI1Pl8ndw731be5Hv95vk8ck7nujGjz1X8p2toS5bkmc00Q0MFrYlVX6kU9Watisx5O+2xU58vPwrC9Xf9bsXJlrmhLK5Wp6raVnYlCxfwu9+oaMM7oEBERkWdxoENERESexdBVmNmmkmtp1ubmgADQPMo9Nh3/yHFxzOGrZVXX5vHy+uMfkZsbmnwbPxBtNinhrPJK4WR7f6lVyJXSB2b14riXtoljMnGxaGtNkqUcbLQmDRZtspCDfE38XFEoDZR7hzM6RERE5Fkc6BAREZFnMXQVKbTNLbOmd3na3nuTRNvwHfK4tM1yqtzM9IpWKhcHa6BMiVLf6EkVZC10bG6Oe1rZgNbMeASA2AqZGVl+lzu8FHtE/nsyfas8T9/80/2YnysKJZusKy+ESzmjQ0RERJ7FgQ4RERF5Fgc6RERE5FlcoxNiajqrSUlvbSuQqatxFbKaMWDuoBwjjmhTdlRui5Nj2qSdcp2AjWBjtl6I9YaL1X0D/r5CwaaKt/Z51EotHFSqHsceMc5rktevnSgroaeWyFIRJ6dkua9dJq/Fe4Js2Hz/2n4P9bfvcs7oEBERkWdxoENERESexdBViFlN31mmrjZlmWEqucmmHt4KTo9Sdi3SFCN5arO3WU0RD0+U5ymHsXpu99j8frQNNk8q4ayzjsnrJ3/sDgnHPlItjtn3zlh5/SxZKsKsoByr3BPge0shYl3JXwn/RvL3Dmd0iIiIyLM40CEiIiLP4kCHiIiIPItrdPqCst2Dls7aFi/HoaItU6aXa5Kffl82WqYS2oiUWGwkUmPXY0e5j9G231DKENhcn+9F54L9/WjrdmIrEuWBxvu2751vi0O0tT1Nymc57tAp12Nz2xbyNpu1fLbbNtis+dO3Iem6HENH/YgUnNEhIiIiz+JAh4iIiDyLoatwU1LJo+tkReKqCXLaOnn3KdGmpaGb2hJ9Vv04bYTQIjk9MBLZ/r7UqWWL66u/e+V9NEOhfB97zvb3VXudDEulvlrvehyjVEb2NbaLtuZR8t+d5udd++443UUfA+caacI2FaKpbwVbcV4LU2mhcDOEHvfSNnmM9jdM6Uckf8dwRoeIiIg8iwMdIiIi8qxeCV09+eSTePTRR1FdXY2JEyfiiSeewNSpU9Vj165di9tuu83V5vf78dVXX/VGV0NOCyNpbRkblA39lEqpZrVks1IyAMRW1Ik2LVsjlNk6zPz5X5ahyroL3b8vrcK1+uHUsrP4u++xYLNbzKwoQGauDDsqw1Qam80/U19V3n9LDFV5g9W9alm52Px7klQns65OakshLL9jbDcJDbewz+i88sorWLx4MVasWIGdO3di4sSJmDVrFo4fl3/YvxEfH4+qqqrAz6FDh8LdTSIiIvKgsA90fvGLX2D+/Pm47bbb8J3vfAdr1qxBbGwsnn322Q7PiYqKQlpaWuAnNTW1w2NbW1vR2Njo+iEiIiICwjzQaWtrQ1lZGfLz8//6hIMGIT8/H6WlpR2e19zcjMzMTGRkZOCaa67Bnj17Ojy2qKgICQkJgZ+MjIyQvgYiIiLqv8K6Rucvf/kLzpw5I2ZkUlNTsW/fPvWccePG4dlnn8UFF1yAhoYGrFq1Crm5udizZw/OOecccfyyZcuwePHiwOPGxsZeG+zYxB+jlWqq5f8k46Ajtstr1WXLcWhaqXtNgFatVa2AqaztCOVaDq4L+V9K1Wst/TfJSPU00zwBoEZJXTZ3xgbke8v08u4LabVk47H2eR/3kFxDob3fqa9+1mUf+H4PLMG+t1rVYxvq3xjLtTeRch9GXB2dnJwc5OTkBB7n5uZi/PjxePrpp/Hggw+K4/1+P/x+f292kYiIiPqJsIauzj77bAwePBg1NTWu9pqaGqSlpVldIyYmBpMnT8b+/cwYICIiou4J64yOz+fDlClTUFxcjMLCQgBAe3s7iouLsXDhQqtrnDlzBp9++imuvPLKMPa0d417Sk7nHbx5hGhL2ifTUs2N/5oylenuko4z2v4WU8J7xnYzPS0sBSPl3EwjBjp4H202+tSqovK9DQ+LCrQZG2UJiIML5edWq6BsE25I2mlXeZufb2+yrYyslbAw08vVTWO1NiVEH8nCHrpavHgx5s6di4suughTp07F448/jpaWlkCtnFtvvRWjRo1CUVERAOCBBx7A9OnTce6556K+vh6PPvooDh06hB/+8Ifh7ioRERF5TNgHOjfccAO++OILLF++HNXV1Zg0aRI2btwYWKBcWVmJQYP+GkGrq6vD/PnzUV1djaSkJEyZMgXvv/8+vvOd74S7q0REROQxvbIYeeHChR2GqkpKSlyPV69ejdWrV/dCr4iIiMjrIi7rKhL0pGy1iINrWwIoqd4xTXKNzolrT4q2sfe727RtIrS2WG3nWsbse8Q2Nq5tAVGT536/ta0E1G07lH6Y77eaDkrhoXyuzHU12hqamAny865tFZG02Z2EcXJKlny+C+V9GPcSkze8INitSWotS1No63ZMlbPl99CYA/2rpAE39SQiIiLP4kCHiIiIPIuhK4Vt2rDGPO6gNu33YqJoO3mOnLaO3TFM9i0pyvXYV98mO6Gl/kXILrKeZxki9BnhhtoJMeKYU3Gy7axj8n4yp6QjeQp5IDBDVWp4Wfm8q//utEgvtwk/UP9k81nW/jZpoXCNWdbC1yjvy7OOBdevjvrWFzijQ0RERJ7FgQ4RERF5FkNXloINByTtk1PP2lR2TIPdmNOsZKlOWyuZXqf7WSXLvma1UaKSYXW6XGa7RI+T94D5vjWP0kKcsjKydu+YWXzaJqIUHmqY27gvtPDy2JfltczPNgC0xXf9vaBW1a6T95x2b1L/YxMO0u45LXuq9Rx3iGvYXhku9ykVu7XvNE2k3HOc0SEiIiLP4kCHiIiIPIsDHSIiIvIsrtGxZLVmA0BbwcWux9X5cufitM2DRZu2c7HGrJRqu06EusdmTZbt71k7ruonuV2eZ1ZPBux3pqc+ZJQYUHcSV6peN4+S/+4009AT98pjtIq3WpkD8zusJ2UIQnkt6h7zd62tl9Gqscc0yXtu9L+71/K0Jcrns620Hsn3AGd0iIiIyLM40CEiIiLPYujKkm215NpMd3resL3yWs2jZFvyblnJUksRNDf14waOfcd2U08tJdxM2dTCD1WXyKlmLeTRZFRZTtosuxXJ08qeY9wD6vuvfLZPxclLmWUnWkbKY5ItK0cEW2VXO4/3UwSx2FgWADI2yLC39n1i4/RY+UdMDdFGyH3CGR0iIiLyLA50iIiIyLMYulLYTt9qYQpzUzRfYw/6oayct+kDImS60OvUe0LbwLNcHpZqZErUXSjvOS2c2ZboE21mJp4Z3gQA30a70GukTDVHIuuwjpFlF3vCLsSZVirf26p/aHU9jtorN/rVvifU0ILxnFo2IN//yGK1KaZyL/nrZLavFkKNrahzPdaqc0P5PulvSyY4o0NERESexYEOEREReRYHOkRERORZXKOjsE0l1yTtdJ+rrb3QmLFSQE/9EztVM6Ye8bTKpeZ7e+Lak+IY3w65HmPY0XbRZsbQtfi5tqM5753uCfr3pa2jU1KCW5XvivRn/e5jkuT7r33HmN9DHT0nedPhAll9X6uqXZ3jrr5u7mYO6JX8Yy3v6UjBGR0iIiLyLA50iIiIyLMYurJlWfHWFFfxpWirnK1t6Cc3cNQqWZoVL+NeYvihPzJDkOnPyhTOpky7a5lVdrXUYpYc6DtaGrcWCm+LD+7fncenyTZ/nZJKbDbwnoh46jIKMxSuhIwyNmp/m2TK+Ymp7sfjH5FLKNSQlLaZdATfT5zRISIiIs/iQIeIiIg8iwMdIiIi8iyu0YGMl6uxRst0UEmWdT/rmDwq7pBM61PTRo1y/1raMPUO25IDNjHu2gkx4hBzh/OOmGXbtfLv8i6kcLFZQ2GzfQwg1+2klsh1e9r6vtiKrr/DuAVI/2SWprD9I26u5QOAzHXubwa9HIpsi3tpm2iz/j7sA5zRISIiIs/iQIeIiIg8i6GrHtBCBK1Jgzt93BGtMnJsWX1Q/aLIooUDmvJltWQbddny3yZppe6wp7qzsFadWUl7NjG80X3m79V2Sl/7Pmke5X6/tdCCVi1bOy6uwgh5GCUOqJ/Y9rH7sXJ/iVICHTBLnYx+s0Uco94nlt8LkRLO4owOEREReRYHOkRERORZDF3BbipeOya2IlFpcz8+eLPMiIhRsmn0jfm67Barm/Yh66na6RNFU3W+O0yRuU7JlFKyJJI/Fk2oneiefvZtVPolT7MKSzFM1XPqfaJkXWkhx9Ykd8hRq4Ice0T+e1Xd/NUCQ5WRT7xHWpViZUNorSJ/0j73fRJdJ0NXtpvSaiLl3uGMDhEREXlWrwx0nnzySXzrW9/CkCFDMG3aNOzYsaPT49etW4fs7GwMGTIE559/Pt54443e6CYRERF5TNgHOq+88goWL16MFStWYOfOnZg4cSJmzZqF48dl4SsAeP/993HTTTfh9ttvx0cffYTCwkIUFhZi9+7d4e4qEREReUyU4zhOOJ9g2rRpuPjii/Fv//ZvAID29nZkZGTgxz/+MZYuXSqOv+GGG9DS0oLXX3890DZ9+nRMmjQJa9asEce3traitbU18LixsREZGRloaGhAfHy8VR8vHzSny2Ns0+Rqrvt2l8do1U1tq6eax0VKDJQ6pt075f/kXnsxYrs878S1J0XbsOJhos2sqJu0U94TNqnkGq7Z6Dnb32FbwcWizVynpa2zOJUg1+OMe8pi3WGQ9wT1HvXeGTuqy/PMdXsdaRnpfmyWqgA6KFehsPle2NS+rtP/39jYiISEhG79/bYR1hmdtrY2lJWVIT8//69POGgQ8vPzUVpaqp5TWlrqOh4AZs2a1eHxRUVFSEhICPxkZGSE7gUQERFRvxbWgc5f/vIXnDlzBqmpqa721NRUVFdXq+dUV1d36/hly5ahoaEh8HP48OHQdJ6IiIj6vX6fXu73++H3+3t0DatNPS1T+FJf/cz1+OBCGcqqyZMp53GHkkSbVi2ZYYP+R3vPxj3kfnxySpY4Zvh6Wd80rkKmf1ZdcpZxjLwvxUaTgBouZXp579BCErKYgAxBaKUpnPEyxKl9x5ghc4YlI0uwVYTVlHDYha7M9PJQhqkiSVhndM4++2wMHjwYNTU1rvaamhqkpaWp56SlpXXreCIiIqKOhHWg4/P5MGXKFBQXFwfa2tvbUVxcjJycHPWcnJwc1/EAsGnTpg6PJyIiIupI2ENXixcvxty5c3HRRRdh6tSpePzxx9HS0oLbbrsNAHDrrbdi1KhRKCoqAgDcfffdmDFjBh577DFcddVVePnll/Hhhx/imWeeCVsfrabhlGn+A3fJqcb2OHcIwn9EXqp+vMySaIuLEW2+jcyK8AJ1StoIhdZOUN5/JUxh9XzKVLYWZgWzbsIi2BCEFjbw1bszbLSsK9sQp1l9XQtxstJ637ENB4k/2sqyirY4eZ72fWJuGutXQuhqOEup9i42G40gYR/o3HDDDfjiiy+wfPlyVFdXY9KkSdi4cWNgwXFlZSUGDfrrLzs3Nxcvvvgifv7zn+O+++7Deeedh9deew0TJkwId1eJiIjIY3plMfLChQuxcOFC9f+VlJSItjlz5mDOnK5r2xARERF1hntdERERkWf1+/TycNBi7Fr677iH5BoHM51crTSppY0rayhsYv39Lc1vINLeo1qjgnb6VrmmQqtueuDGoaItY6O8x0zRB47KRiXlnNVye87mM2mzbgsAmrLc7/foN7suLwAAzaNkW8YGoyK7Vl5A9or6obapzaJtmLKWy5zraE0aLI7wKfezOnCI4HIFnNEhIiIiz+JAh4iIiDyLoStLWoqdFs4yK00eLpBTgYl7u65a2pFImQoke1qYIvljdwiiLdEnjtFSRLXNP83U9FituKnFBrHUe7TPcbTyHmkbtJrSt8o2M+QFwGpDYFZL7kNKyrYacjbvE+VzPPZ+edree+UOApnrzrge+zZ+0FkPA2zvnUjBGR0iIiLyLA50iIiIyLM40CEiIiLP4hodS9p6HM3xae7H456SsczDV8s1OmZ5dgBI2lxv9ZzUD5nl0gsuFodoJdur88+ItrTN7scns5LEMdoaM5s4O9dn9FxP1r2YZSe0kgPa+r5mJeVc3hXK8/H97jNW63Eg/35kbJCnaX9P/EfkvEZs2WfuBqXkhFqGQLlPIvne4YwOEREReRYHOkRERORZDF1Zsk0vH/vyl67HB2+WYSrbasnatKX5hkXydCF9zSZE1KRVJG2Uu9xnrtOeQYazBIt7CZB9Zbpxz/Xk92VWPdYqaGthCq0yLjYkuh/zfew1NiFh9T5R2tKNcOYf3/m9OGbmP/xQtJ11TJnX0MpOmP0aO8qqX5GMMzpERETkWRzoEBERkWcxdGVL23AvM0a01RuVkGOPyEv56ttko1alVnlOhg28SauAq2VPafecWUE5eXfXm3wCdlPSvN96znYDT+07wMyo0cJUmvRnZRXcugvdmzomWWbTUM/ZhIRtwkgAcNKoon7RyrvEMb4kGfaOvuYL0VYD99KK1Fc/E8fAAxv9ckaHiIiIPIsDHSIiIvIsDnSIiIjIs7hGR6PESrX1Elr6r1kJuSZPppdrOwsnHZDdOO2B2CjpxFoIZb2MtpbLp6Sh1493P07fKs/T7iU1ddzmGK7j6BbbtGFt92pzV3t/nSwloK3b0iRt7vr7JFqpjMvvodBT7wnlOyC6TpYTML8XkpXviQM3yr8xsW+niLbkQ+71fFrJFN/G/v9554wOEREReRYHOkRERORZDF0ptBTOuIovlSN9sslI2Yw7ZJci7IXqk2RPVEpVjtE2cNRkbHRPP2uh0TjIsAiUTQS5qWfo2Yb/tDCF2aaFwltGyudsi5P/hm277tuux1oqMcNUvUNNL7fc1NMMS43YLk/TNpPW7h2t4r/ggfA1Z3SIiIjIszjQISIiIs9i6ApyGjHupW1dHgPov7yDC93Tw0n7ZGaWWckWkNkVgBoYI48wp37blGyHuEOywnGrknVlbgjbmiTvVXNzSADIqEvsqptWG39S53ry+zptbOCo3ROAXdZVaslxd0OQG71SeGi/51oj3AgAI7a7/6Y0j5LzFXXZMkyl/S2qMcOZ5j3SkX52T3BGh4iIiDyLAx0iIiLyLA50iIiIyLO4Rgd2MWgt/VtL441pcj/W4qfpW2UaKbZ93GUfyDvMNV/RxjobQK7P+JpcuVW7Osr1uOEDec+Z9yWgV/s21/uoadAeSDfta2p6sbKbeJNR6kKrjKyt+cvYINdamO+3+V4DfL97jbI+Slunl/yx/Fthruesy5afd209jk1lbI0X3mvO6BAREZFncaBDREREnsXQVQ9oYalhR91ThlpF5WitAqY2la1Mb5rT216YVhyIzPctWnuvFWol01+6p7xPZ8pDtLRkNXRhhMsYtggP2xCRKHWhbPyZ8bxFdVsAPuO91SrlJrMycliIiuPK7zn2hHz/6/LlJqvm3x0tLK2FOLUNO20qI3vhO4AzOkRERORZHOgQERGRZ3GgQ0RERJ4V1jU6J06cwI9//GNs2LABgwYNwvXXX49//dd/xbBhwzo8Jy8vD1u2bHG13XnnnVizZk04u+qixiSV9O/kxIu7vJZWej85UcZKfRs/kCf3szgo2Yse5469azF77T7U4uyHC9zbQsQekc+nrsfRnnOcXBNAoaeml2tbMhht2i73Gu0+acp0bxWh7l6uXKu/rceIRGJNnvI508pJaGttzHRyLZVcXXuj3F+iXx5Yj6MJ60DnlltuQVVVFTZt2oRTp07htttuwx133IEXX3yx0/Pmz5+PBx54IPA4NjY2nN0kIiIijwrbQGfv3r3YuHEjPvjgA1x00UUAgCeeeAJXXnklVq1ahZEjR3Z4bmxsLNLS0sLVNSIiIhogwjbQKS0tRWJiYmCQAwD5+fkYNGgQtm/fjmuvvbbDc3/729/iP//zP5GWloarr74a999/f4ezOq2trWhtbQ08bmxsDN2L6IIWDjDTv331sqJydJ1SGZnp5QOKFjYSlPdfm5LOgDtMcbhAXkpNJdYqL5ulD7R7kPdcWFiFL9Vq2ZL23RRbVh9MtygExPuoVMFuUyqVtyYNFm1nHXM/1sJbWiV/7e+OSHv36Gc7bAOd6upqjBjh/nKNjo7G8OHDUV1d3eF5N998MzIzMzFy5Eh88skn+NnPfoby8nK8+uqr6vFFRUVYuXJlSPtORERE3tDtgc7SpUvx8MMPd3rM3r17g+7QHXfcEfjv888/H+np6bjssstw4MABjB07Vhy/bNkyLF68OPC4sbERGRkZQT8/EREReUe3BzpLlizBvHnzOj1mzJgxSEtLw/Hj7o3lTp8+jRMnTnRr/c20adMAAPv371cHOn6/H36/3/p6wVJXo1tkrWjThVZhC4AhgoFOmd5WQ0mGjI1yKrsp066ShJmtY1M5lbrPtjKyWRk3aWcPvhPMe0e7vyg8zOw55W9AU+a3rS5lbuJa9Q+t4pjh6+VSD2XvV0QPkHug2wOdlJQUpKSkdHlcTk4O6uvrUVZWhilTpgAA3n77bbS3twcGLzZ27doFAEhPT+9uV4mIiGiAC1vBwPHjx6OgoADz58/Hjh078N5772HhwoW48cYbAxlXR48eRXZ2Nnbs2AEAOHDgAB588EGUlZXh888/x3//93/j1ltvxaWXXooLLrggXF0lIiIijwprZeTf/va3yM7OxmWXXYYrr7wSl1xyCZ555pnA/z916hTKy8tx8uRJAIDP58PmzZtxxRVXIDs7G0uWLMH111+PDRs2hLObRERE5FFRjuM4fd2JUGpsbERCQgIaGhoQHx9vdc7lg+aEuVduavVZJVbq1VQ/smN7n9Rc13Vs39coq6eaFVYBYPSb7jVlPVljNlBSV/ua7X1ig+9Rz9lUFz+ppJKbFc4BYOzLX4q2ytnuEgPa7uXJu0+JNptyKOraMa2Kc5C73G9qX9fp/w/m77cN7nVFREREnsWBDhEREXlWWPe66i9sptjVTfiCpG3eNlDS/KgbLMOZvkYjBXmznFbWwlvaZoChxDBIL1HuE21TT3XjYOoR9e+CRVkIreKxWhYia6hoM0NVGRuOi2O0vzE24SbbMir9DWd0iIiIyLM40CEiIiLP4kCHiIiIPItrdGC3liCkaXfKtU53fRaRes+1xbv/vaLtXNwyUl6rZaT8d05dtju2f9YxGevXdj2PNnc9V3DNTufU9RFBfjf5Nna9zpDvR89pv8Omm6aLNm2HcRvNo+RnVEsdN6llISyez6v3BGd0iIiIyLM40CEiIiLPYuiKqJ9Lfvp91+Oqn+SKY846ZnctLcRlakv0ibZobVd1lkzolmDDBqy03jtsQ4taeQcZTpafodoJMaJNSx03Hb56hDzv+c+6PG8g4YwOEREReRYHOkRERORZDF31gBcqRlLkUkMNSltbwcWuxz5lkz9tU09NW5z73z5xh2SGR1OmnGIH5CaFPptKrCHcMHAg8Grl2t4Uyuw2zeG5sgq5GUpqeF5+XjLukftrH3gwVrSlP+t3P94qM6zMSswA1O+OgYIzOkRERORZHOgQERGRZ3GgQ0RERJ7FNTpEEcp2LUFsWYXrce0EuUbA1yivr1Vdjf0/X7geHz4nWRyTsVGu2/HVt4k2Ne3ZwPUl3RPuFPFg16/0J7avx7x/tR3BtQrEWmp3zXXuz2T9Z3LN3Ok8+Xn07dB65v78aVXJvfae9RRndIiIiMizONAhIiIiz2LoiihCBTv9rE2d1+V3HUYCgJNvp7geJyqp6q0yMxaxFUqKKysj9zsMefyVGVa1LoUwfaJoMiuOt8fJTT6TP/5StFVdIsNlZphY28TX9g+7zfvthXAmZ3SIiIjIszjQISIiIs9i6IrIY7Rp5bgKOb1dly2nxWOMUJW+yaf895GWkVKb595sMPljGd4K5RQ7UbC08IxZXVgLU2nhrIOz5WdhzIvuzTlr8uRGnE1Zg0WbuqmnERLWPhun5VlB88JnjzM6RERE5Fkc6BAREZFncaBDREREnsU1OkT9nBlDV9NglfNGr3hftB14fLrrsZYGO/pNWQVZS4MddrTrHdPV1Fil2qz5ReWFdQMUeuZam57cJ+a6s6qf5IpjtDU05jo3QK7JaYuTx8Qdkp81tV/m590D6d/hxhkdIiIi8iwOdIiIiMizGLoi8hqtIrGS/q1NxY992R02akv0iWOasoaKNi1MlbTTPX2upaBr10riJoUUpGDDOlpb1Vz3RpzpW2VIte5CeX2fEroySyton6vaCTGiLbZCNAla3xnOcuOMDhEREXkWBzpERETkWRzoEBERkWdxjQ6Rx6ixeKUtfZs8zExNj1bWy2BKlmg6NEce1jzKnVKrruPZrOz+bMF6J+kQ4rqHyGfeF9o90VZwsdW1Mp7/zH0tpRSCrdqJ7vVpWnq5lqqubjuhbVdhnsf70oUzOkRERORZYRvoPPTQQ8jNzUVsbCwSExOtznEcB8uXL0d6ejqGDh2K/Px8/PnPfw5XF4mIiMjjwha6amtrw5w5c5CTk4Pf/OY3Vuc88sgj+OUvf4nnn38eWVlZuP/++zFr1iz86U9/wpAhQ8LVVSLPs9mdWXPYSLEFgOTdp0SbL65VtDWPdz9OX/WBOKbmTpniru5ybvZVSaEPd2iJ4YDIZxO+9G2U96F275w0QrRNmTL9O+OWA6Kt/v+O7rIPGq38AsOloRG2gc7KlSsBAGvXrrU63nEcPP744/j5z3+Oa665BgDwwgsvIDU1Fa+99hpuvPHGcHWViIiIPCpi1uhUVFSguroa+fn5gbaEhARMmzYNpaWlHZ7X2tqKxsZG1w8REREREEFZV9XV1QCA1NRUV3tqamrg/2mKiooCs0dEpFOrp1qcl75VTqebGSQAkPkLudFnU9Zg9+ObpotjtOwTlRGqqsuXWVdmJWYAarYZRTabrCIAaujVDP9oWYPavRNX8aVoMysVnzxHZg1mKN1qTRqstLq1xcs5htQyizLICO3GpQNFt2Z0li5diqioqE5/9u3bF66+qpYtW4aGhobAz+HDh3v1+YmIiChydWtGZ8mSJZg3b16nx4wZMyaojqSlpQEAampqkJ6eHmivqanBpEmTOjzP7/fD7/cH9ZxERETkbd0a6KSkpCAlJSUsHcnKykJaWhqKi4sDA5vGxkZs374dd911V1iek4iIiLwtbGt0KisrceLECVRWVuLMmTPYtWsXAODcc8/FsGHDAADZ2dkoKirCtddei6ioKCxatAj//M//jPPOOy+QXj5y5EgUFhaGq5tE9DfMeH+bUgVZW1cTXSdTwo/f6N6ZfOzLch1EnN2yBLGuQl2Po5k+UbZt+9juXAtM/+07NlWDbdfjaMxK3ifPkccc+O+xos0XL487MdVdkiFz3RmrPvBeCo2wDXSWL1+O559/PvB48uTJAIB33nkHeXl5AIDy8nI0NDQEjrn33nvR0tKCO+64A/X19bjkkkuwceNG1tAhIiKioIRtoLN27doua+g4juN6HBUVhQceeAAPPPBAuLpFREREA0jEpJcTUe9SU87NzTKVKrLJkJsialVdxz3lvv7hq0eIY5rHyyrL5/27TFX31xlT/UplZLXScwjDVBqGFrqnt0N92qax2uacbYm+Lq/VHifDTRnPy+tr1cTHP1LnenwyK0n2y/L3wPTy7ouYgoFEREREocaBDhEREXkWQ1dEA5RagdYICYlQFoBopYLr3n+R2VnDd7jDWVpl2cx1sgtaGEFuqCifLzbIyrKAnP7XXrfNhpH0V7a/Q6tQjJY9p4Sg6rLcmX5ahlWTcQygVyo2swvTNssuHFwow1Sj35QZiGZoV9tY1BZDVd3HGR0iIiLyLA50iIiIyLM40CEiIiLP4hodogHKJtYfre0QrayNMNNnAbkuoS1OpqAfLpDPGdMg//0V0+R+rO0Q3apUwfXXybU8qJfp62K9h1LpuS+EMpW4t9OSbdbjaLS1PSeVdVuxFfKe8xvHaetxRKmCDtrMe0xb72Nb2dvcRf10mCt2kxtndIiIiMizONAhIiIiz2LoiogCRGhBqUCsfmkoxzVd6L6WuUkiAKRvleGAytkyxJW+1R1K0lLQayeYKej65p9aFWdzU1LrVPIQhiBCmdJuVYFY6bsZYlHPs3w+jRb2NNlUKQb06sJmGQJfo7zntPtEq9DtP+KeB0jaqaSN274/5u+HYapexRkdIiIi8iwOdIiIiMizONAhIiIiz+IaHSIKsEo51nY9V9ZomCm7WjrwwZvljubZMw+Itpp97jRxrWR/xobjsq+KqkvkGp3k3ca6kCy5Q7u2xcRp5foijVtbl6Ks0Qg2HVu9vrLWRlwryPU4ah8st23Q0r1NtjuO106U72Pqq5+5HtcpJQd8TaIJaZtluQJ/nXvdjrYmyJcUunVOFD6c0SEiIiLP4kCHiIiIPIuhKyIKC7FDsxKGGf2mDD/sw1h5seyun88m3RjQ09y1Ssum2rlyp2rtWkl1ia7HaigjyNRuLQVdc3JK17u7a+EU2zRxmBWzlfICWjVjrbqwGYLSwk112V1Xy9bO1SoeH58mz4s9Iq9v9tUsQQBAfd0MU0UezugQERGRZ3GgQ0RERJ7F0BUR9ZjVdL2yQahWBTetVFap9RkbcWrnmccAQPM/tIq2qL3DRJszvrnLY846JprUyssmLYwkwnoAoISlaq9zh8u0Sr9aOMjq+komkxZu0n6vZphNCzdpmXFQQon1492vSdvUdfSbMmx0aLEj2pp3uN83LXQ1YrvsVlyFvL4aqjLYhv8YzupbnNEhIiIiz+JAh4iIiDyLAx0iIiLyLK7RIaJuCXoNgpKKG6scplYINtaXaGt0tF3Ph6+Xa1qq8+UaoLT17p6YVXEBPQVd2wndRluBrLwsV8IAqSUW1Z6V32vNnbnBXUuhrVUxKxVr63GSP5bnae9by0j3uh1tPY5WUXn4etnXtnjjPGVNkNavYHE9Tv/AGR0iIiLyLA50iIiIyLMYuiKibgk6pVZJL9eo1zJCRFpIImmfDFNpIZW0zV33oTrH7vpqNWMjrHO4QIa8tEq8WpVlwF3tWdtYVPu9Jj/9vjzOCP9pISktRNSUKTdebRnpfqyl3uubp8qQ4JgX3SE1rcK1RkurN8OXWr/UUJwWglRCguI8ppf3C5zRISIiIs/iQIeIiIg8iwMdIiIi8iyu0SGibgl2DYKWNq7RUq9Naqq0lmZ9ndxxXFvbEftIteux//+OFsdo6eXl/yS3Pkjc6/73Y+wRcYi6Hkfrl7nOpTVJPl/zKPnvVV+eXFdjbh+h7Qiu0dYm1Y93P25rsltzVDtBrn0y1wBp21xoW21oa3nM7UPUkgAW5QsAiLVPtvcv1+NEHs7oEBERkWdxoENERESexdAVEXVLuKfmYyvqZKMRljJTuAEASopw3CGZzqylUB+qTnU99ikhFl+TfMoR27UQlDt9WXs+Le29dqLsvxn+OT5N9iFjY9e7vQMyDKZVINZSwrUQV+Je9+MTU2UfMjbI97F5lAypme+RVjrANuXcvHdilWz808p5ajhLCdFS/xS2GZ2HHnoIubm5iI2NRWJiotU58+bNQ1RUlOunoKAgXF0kIiIijwvbjE5bWxvmzJmDnJwc/OY3v7E+r6CgAM8991zgsd/vD0f3iIiIaAAI20Bn5cqVAIC1a9d26zy/34+0tLQw9IiI+pJtKMAmuyVaq7KsZF1pX3Cxyrn+OndopClTubwSnhm+Q4ZZzI0r/XVnxDE33f+maPvVuqtEW2qJO0yohX5alahOdY4MQZkZSRqtcrHGzGbSfg91F8r3+1Rc19fWMurUysXbPpZtZvaUdk/Y3ofMnvKMiFuMXFJSghEjRmDcuHG46667UFtb2+nxra2taGxsdP0QERERARE20CkoKMALL7yA4uJiPPzww9iyZQtmz56NM2fkv4i+UVRUhISEhMBPRkZGL/aYiIiIIlm3BjpLly4Vi4XNn3379gXdmRtvvBHf//73cf7556OwsBCvv/46PvjgA5SUlHR4zrJly9DQ0BD4OXz4cNDPT0RERN7SrTU6S5Yswbx58zo9ZsyYMT3pj7jW2Wefjf379+Oyyy5Tj/H7/VywTNQP2K550KrUinU7FjtLd9gPbb2HQUtLTy1R0t4Ve+91L5o5799lqvfGGeeJtrQp8jkPX+1ek3PyHJnOfipO/nv1VELXVYnTt8p+aZWENWal4vK75LqXcU/J99tcC6XR1mip1bK1to0fuB5yPQ4B3RzopKSkICUlJVx9EY4cOYLa2lqkp6f32nMSERGRd4RtjU5lZSV27dqFyspKnDlzBrt27cKuXbvQ3NwcOCY7Oxvr168HADQ3N+OnP/0ptm3bhs8//xzFxcW45pprcO6552LWrFnh6iYRERF5WNjSy5cvX47nn38+8Hjy5MkAgHfeeQd5eXkAgPLycjQ0NAAABg8ejE8++QTPP/886uvrMXLkSFxxxRV48MEHGZoiGkBs0su1UINtmCL6wFHZZjw+OSVLnmi5aej4R9zp0errUcJzWkVos7KvViH4y/9XZqYO+o9k0dZsFJPWKjFrG2pq6fFm+C9jY9fHAHrFZmH6RNEUW6aUONZY3AMMUw08YRvorF27tssaOo7jBP576NCheOutt8LVHSIiIhqAIiq9nIiIiCiUuKknEUU8MyylhR96Es4yNwlVQyVKRWUtO8sMcVX9JFcckrHBrvpvdJ17400tKyr5Hke0ncySoSRzI1F1w1NlQ01N7QR3X7WKylqIUPsdCrYZdRbVsW3vCYazvI0zOkRERORZHOgQERGRZ3GgQ0RERJ7FNTpEFPHMNRQ9WWdhrscBIHfC1tb2KGtHYpXrmynn6VtblKMkc9dzAGjKGup6nLRZpqprqfBaqrqv3r2uRnu+1Fc/E23a7yvOMttbUH6HQa+PCfI8rscZeDijQ0RERJ7FgQ4RERF5FkNXRNTvWIeptONCGPLQvkBTS4zzLDYRBTqojFxW3+V5Wkq4WVEZUMJgO+XrOTxXVnrOeF6Gs0zq71mp/myTOt6TsKRNGQIaeDijQ0RERJ7FgQ4RERF5Fgc6RERE5Flco0NEA4rNGhDrtR0W2xBE225pYEHrV/LTdutjtDU5JnVrCovfhfo7tdiFXju3J+tquCaHNJzRISIiIs/iQIeIiIg8i6ErIhpQQhne0MIzNqEY613VLSpCB90vpeKxtuN4SEN9CoabKNw4o0NERESexYEOEREReRZDV0REIWQTigl3pV+bcFO0kjEWyjBSTyocE4USZ3SIiIjIszjQISIiIs/iQIeIiIg8i2t0iIgiVLBrWqzWCQVZudj6+txxnCIEZ3SIiIjIszjQISIiIs9i6ArApvZ1fd0FIiIiCgPO6BAREZFncaBDREREnsWBDhEREXkWBzpERETkWRzoEBERkWdxoENERESexYEOEREReRYHOkRERORZHOgQERGRZ3GgQ0RERJ4VtoHO559/jttvvx1ZWVkYOnQoxo4dixUrVqCtra3T87766issWLAAycnJGDZsGK6//nrU1NSEq5tERETkYWEb6Ozbtw/t7e14+umnsWfPHqxevRpr1qzBfffd1+l599xzDzZs2IB169Zhy5YtOHbsGK677rpwdZOIiIg8LMpxHKe3nuzRRx/FU089hYMHD6r/v6GhASkpKXjxxRfxd3/3dwC+HjCNHz8epaWlmD59ujintbUVra2trmuMHj0ahw8fRnx8fHheCBEREYVUY2MjMjIyUF9fj4SEhJBdt1d3L29oaMDw4cM7/P9lZWU4deoU8vPzA23Z2dkYPXp0hwOdoqIirFy5UrRnZGSEptNERETUa2pra/vnQGf//v144oknsGrVqg6Pqa6uhs/nQ2Jioqs9NTUV1dXV6jnLli3D4sWLA4/r6+uRmZmJysrKkP6iIt03I+GBNpPF183XPRDwdfN1DwTfRGQ6mxAJRrcHOkuXLsXDDz/c6TF79+5FdnZ24PHRo0dRUFCAOXPmYP78+d3vZSf8fj/8fr9oT0hIGFA3yDfi4+P5ugcQvu6Bha97YBmor3vQoNAuH+72QGfJkiWYN29ep8eMGTMm8N/Hjh3DzJkzkZubi2eeeabT89LS0tDW1ob6+nrXrE5NTQ3S0tK621UiIiIa4Lo90ElJSUFKSorVsUePHsXMmTMxZcoUPPfcc12O0qZMmYKYmBgUFxfj+uuvBwCUl5ejsrISOTk53e0qERERDXBhSy8/evQo8vLyMHr0aKxatQpffPEFqqurXWttjh49iuzsbOzYsQPA1+Gm22+/HYsXL8Y777yDsrIy3HbbbcjJyVEXImv8fj9WrFihhrO8jK+br3sg4Ovm6x4I+LpD+7rDll6+du1a3Hbbber/++YpP//8c2RlZeGdd95BXl4egK8LBi5ZsgQvvfQSWltbMWvWLPzqV79i6IqIiIi6rVfr6BARERH1Ju51RURERJ7FgQ4RERF5Fgc6RERE5Fkc6BAREZFn9fuBzueff47bb78dWVlZGDp0KMaOHYsVK1agra2t0/O++uorLFiwAMnJyRg2bBiuv/561NTU9FKvQ+Ohhx5Cbm4uYmNjxbYZHZk3bx6ioqJcPwUFBeHtaIgF87odx8Hy5cuRnp6OoUOHIj8/H3/+85/D29EQO3HiBG655RbEx8cjMTERt99+O5qbmzs9Jy8vT7zfP/rRj3qpx8F58skn8a1vfQtDhgzBtGnTAuUnOrJu3TpkZ2djyJAhOP/88/HGG2/0Uk9Dqzuve+3ateJ9HTJkSC/2NjTeffddXH311Rg5ciSioqLw2muvdXlOSUkJLrzwQvj9fpx77rlYu3Zt2PsZat193SUlJeL9joqK6nBrpEhUVFSEiy++GHFxcRgxYgQKCwtRXl7e5Xmh+Hz3+4HOvn370N7ejqeffhp79uzB6tWrsWbNGtx3332dnnfPPfdgw4YNWLduHbZs2YJjx47huuuu66Veh0ZbWxvmzJmDu+66q1vnFRQUoKqqKvDz0ksvhamH4RHM637kkUfwy1/+EmvWrMH27dtx1llnYdasWfjqq6/C2NPQuuWWW7Bnzx5s2rQJr7/+Ot59913ccccdXZ43f/581/v9yCOP9EJvg/PKK69g8eLFWLFiBXbu3ImJEydi1qxZOH78uHr8+++/j5tuugm33347PvroIxQWFqKwsBC7d+/u5Z73THdfN/D19gB/+74eOnSoF3scGi0tLZg4cSKefPJJq+MrKipw1VVXYebMmdi1axcWLVqEH/7wh3jrrbfC3NPQ6u7r/kZ5ebnrPR8xYkSYehh6W7ZswYIFC7Bt2zZs2rQJp06dwhVXXIGWlpYOzwnZ59vxoEceecTJysrq8P/X19c7MTExzrp16wJte/fudQA4paWlvdHFkHruueechIQEq2Pnzp3rXHPNNWHtT2+xfd3t7e1OWlqa8+ijjwba6uvrHb/f77z00kth7GHo/OlPf3IAOB988EGg7c0333SioqKco0ePdnjejBkznLvvvrsXehgaU6dOdRYsWBB4fObMGWfkyJFOUVGRevzf//3fO1dddZWrbdq0ac6dd94Z1n6GWndfd3c+8/0FAGf9+vWdHnPvvfc63/3ud11tN9xwgzNr1qww9iy8bF73O++84wBw6urqeqVPveH48eMOAGfLli0dHhOqz3e/n9HRNDQ0dLr7aVlZGU6dOoX8/PxAW3Z2NkaPHo3S0tLe6GKfKikpwYgRIzBu3DjcddddqK2t7esuhVVFRQWqq6td73dCQgKmTZvWb97v0tJSJCYm4qKLLgq05efnY9CgQdi+fXun5/72t7/F2WefjQkTJmDZsmU4efJkuLsblLa2NpSVlbnep0GDBiE/P7/D96m0tNR1PADMmjWr37yvQHCvGwCam5uRmZmJjIwMXHPNNdizZ09vdLdPeeH97olJkyYhPT0dl19+Od57772+7k6PNDQ0AECnf6tD9X53e6+rSLd//3488cQTWLVqVYfHVFdXw+fzifUdqamp/SrmGYyCggJcd911yMrKwoEDB3Dfffdh9uzZKC0txeDBg/u6e2HxzXuamprqau9P73d1dbWYpo6Ojsbw4cM7fQ0333wzMjMzMXLkSHzyySf42c9+hvLycrz66qvh7nK3/eUvf8GZM2fU92nfvn3qOdXV1f36fQWCe93jxo3Ds88+iwsuuAANDQ1YtWoVcnNzsWfPHpxzzjm90e0+0dH73djYiC+//BJDhw7to56FV3p6OtasWYOLLroIra2t+PWvf428vDxs374dF154YV93r9va29uxaNEifO9738OECRM6PC5Un++IndFZunSpuvjqb3/ML4GjR4+ioKAAc+bMwfz58/uo5z0TzOvujhtvvBHf//73cf7556OwsBCvv/46PvjgA5SUlITuRQQh3K87UoX7dd9xxx2YNWsWzj//fNxyyy144YUXsH79ehw4cCCEr4J6W05ODm699VZMmjQJM2bMwKuvvoqUlBQ8/fTTfd01CoNx48bhzjvvxJQpU5Cbm4tnn30Wubm5WL16dV93LSgLFizA7t278fLLL/fK80XsjM6SJUswb968To8ZM2ZM4L+PHTuGmTNnIjc3F88880yn56WlpaGtrQ319fWuWZ2ampo+31Oru6+7p8aMGYOzzz4b+/fvx2WXXRay63ZXOF/3N+9pTU0N0tPTA+01NTWYNGlSUNcMFdvXnZaWJhamnj59GidOnOjWPTtt2jQAX898jh07ttv9Daezzz4bgwcPFtmPnX0u09LSunV8JArmdZtiYmIwefJk7N+/PxxdjBgdvd/x8fGenc3pyNSpU7F169a+7ka3LVy4MJBM0dXsY6g+3xE70ElJSUFKSorVsUePHsXMmTMxZcoUPPfccxg0qPOJqilTpiAmJgbFxcW4/vrrAXy9mr2yshI5OTk97ntPdOd1h8KRI0dQW1vrGgD0hXC+7qysLKSlpaG4uDgwsGlsbMT27du7nbEWaravOycnB/X19SgrK8OUKVMAAG+//Tba29sDgxcbu3btAoA+f781Pp8PU6ZMQXFxMQoLCwF8PcVdXFyMhQsXqufk5OSguLgYixYtCrRt2rSpzz/H3RHM6zadOXMGn376Ka688sow9rTv5eTkiPTi/vZ+h8quXbsi8nPcEcdx8OMf/xjr169HSUkJsrKyujwnZJ/vYFZLR5IjR4445557rnPZZZc5R44ccaqqqgI/f3vMuHHjnO3btwfafvSjHzmjR4923n77befDDz90cnJynJycnL54CUE7dOiQ89FHHzkrV650hg0b5nz00UfORx995DQ1NQWOGTdunPPqq686juM4TU1Nzk9+8hOntLTUqaiocDZv3uxceOGFznnnned89dVXffUyuq27r9txHOdf/uVfnMTEROcPf/iD88knnzjXXHONk5WV5Xz55Zd98RKCUlBQ4EyePNnZvn27s3XrVue8885zbrrppsD/N+/z/fv3Ow888IDz4YcfOhUVFc4f/vAHZ8yYMc6ll17aVy+hSy+//LLj9/udtWvXOn/605+cO+64w0lMTHSqq6sdx3GcH/zgB87SpUsDx7/33ntOdHS0s2rVKmfv3r3OihUrnJiYGOfTTz/tq5cQlO6+7pUrVzpvvfWWc+DAAaesrMy58cYbnSFDhjh79uzpq5cQlKampsDnF4Dzi1/8wvnoo4+cQ4cOOY7jOEuXLnV+8IMfBI4/ePCgExsb6/z0pz919u7d6zz55JPO4MGDnY0bN/bVSwhKd1/36tWrnddee83585//7Hz66afO3Xff7QwaNMjZvHlzX72EbrvrrruchIQEp6SkxPV3+uTJk4FjwvX57vcDneeee84BoP58o6KiwgHgvPPOO4G2L7/80vnHf/xHJykpyYmNjXWuvfZa1+CoP5g7d676uv/2dQJwnnvuOcdxHOfkyZPOFVdc4aSkpDgxMTFOZmamM3/+/MCXaX/R3dftOF+nmN9///1Oamqq4/f7ncsuu8wpLy/v/c73QG1trXPTTTc5w4YNc+Lj453bbrvNNbgz7/PKykrn0ksvdYYPH+74/X7n3HPPdX760586DQ0NffQK7DzxxBPO6NGjHZ/P50ydOtXZtm1b4P/NmDHDmTt3ruv43/3ud863v/1tx+fzOd/97nedP/7xj73c49DozutetGhR4NjU1FTnyiuvdHbu3NkHve6Zb9KmzZ9vXuvcuXOdGTNmiHMmTZrk+Hw+Z8yYMa7PeX/R3df98MMPO2PHjnWGDBniDB8+3MnLy3Pefvvtvul8kDr6O/2371+4Pt9R/9sBIiIiIs+J2KwrIiIiop7iQIeIiIg8iwMdIiIi8iwOdIiIiMizONAhIiIiz+JAh4iIiDyLAx0iIiLyLA50iIiIyLM40CEiIiLP4kCHiIiIPIsDHSIiIvKs/x+yPq5UFwV/FwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 10_000\n",
    "\n",
    "x, _ = datasets.make_moons(n_samples=n_samples, noise=.06)\n",
    "\n",
    "scaler = preprocessing.StandardScaler()\n",
    "x = scaler.fit_transform(x)\n",
    "\n",
    "plt.hist2d(x[:, 0], x[:, 1], bins=100)\n",
    "plt.xlim(-2 ,2)\n",
    "plt.ylim(-2, 2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNFs\n",
    "\n",
    "_TODO: Harmonize notation with code and add details._\n",
    "\n",
    "The evolution of the log-density follows the instantaneous change-of-variables formula:\n",
    "$$\\frac{\\partial \\log p({z}(t))}{\\partial t}=-\\operatorname{Tr}\\left(\\frac{\\partial f}{\\partial {z}(t)}\\right)$$\n",
    "\n",
    "Get total change in log-density by integrating across time:\n",
    "$$\\log p_1\\left({z}\\left(t_1\\right)\\right)=\\log p_0\\left({z}\\left(t_0\\right)\\right)-\\int_{t_0}^{t_1} \\operatorname{Tr}\\left(\\frac{\\partial f}{\\partial {z}(t)}\\right) d t$$\n",
    "\n",
    "We can get an unbiased estimate of the trace of a matrix by taking a double product of that matrix with a noise vector.\n",
    "$$\\operatorname{Tr}(A)=E_{p({\\epsilon})}\\left[{\\epsilon}^T A {\\epsilon}\\right]$$\n",
    "\n",
    "Typically we'd also need to implement backprop-ing through an ODE with e.g. adjoints, but Diffrax will take care of this for us here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MLP(nn.Module):\n",
    "    \"\"\" A simple MLP in Flax.\n",
    "    \"\"\"\n",
    "    hidden_dim: int = 32\n",
    "    out_dim: int = 2\n",
    "    n_layers: int = 3\n",
    "\n",
    "    @nn.compact\n",
    "    def __call__(self, x):\n",
    "        for _ in range(self.n_layers):\n",
    "            x = nn.Dense(features=self.hidden_dim)(x)\n",
    "            x = nn.gelu(x)\n",
    "        x = nn.Dense(features=self.out_dim)(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Adapted from [Diffrax](https://docs.kidger.site/diffrax/examples/continuous_normalising_flow/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logp_exact(t, y, args):\n",
    "    \"\"\" Compute trace directly.\n",
    "    \"\"\"\n",
    "    y, _ = y\n",
    "    _, func = args\n",
    "    t = np.atleast_1d(t)\n",
    "\n",
    "    fn = lambda y: func(np.concatenate([y, t]))  # Augmented function\n",
    "    f, f_vjp = jax.vjp(fn, y)  # VJPs can be computed at the ~same cost as computing f through reverse-mode AD\n",
    "\n",
    "    # Compute trace\n",
    "    (size,) = y.shape\n",
    "    (dfdy,) = jax.vmap(f_vjp)(np.eye(size))\n",
    "    logp = np.trace(dfdy)\n",
    "    return f, logp\n",
    "\n",
    "def logp_approx(t, y, args):\n",
    "    \"\"\" Approx. trace using Hutchinson's trace estimator.\n",
    "    \"\"\"\n",
    "    y, _ = y\n",
    "    z, func = args\n",
    "    t = np.atleast_1d(t)\n",
    "    \n",
    "    fn = lambda y: func(np.concatenate([y, t]))  # Augmented function\n",
    "    f, f_vjp = jax.vjp(fn, y) # VJPs can be computed at the ~same cost as computing f through reverse-mode AD\n",
    "    \n",
    "    # Trace estimator\n",
    "    (z_dfdy,) = f_vjp(z)\n",
    "    logp = np.sum(z_dfdy * z)\n",
    "    return f, logp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = jax.random.PRNGKey(0)\n",
    "t = np.ones((x.shape[0], 1))\n",
    "\n",
    "f = MLP(hidden_dim=64, out_dim=2, n_layers=3)\n",
    "params = f.init(key, np.concatenate([x, t], axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([3.156627 , 2.2985032, 3.1344635, 3.863559 , 2.96205  , 2.404113 ,\n",
       "       2.4348328, 3.7143767, 3.1841269, 3.1162772, 3.385356 , 3.1589055,\n",
       "       3.0098429, 2.634412 , 2.9587815, 3.777234 , 3.6070356, 3.055339 ,\n",
       "       3.1169815, 3.5339994, 1.9420253, 2.847947 , 1.9530888, 3.5562863,\n",
       "       3.1166618, 2.3222103, 2.065152 , 2.128791 , 3.3884706, 1.9760624,\n",
       "       3.2089186, 2.6790004], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = 0.0\n",
    "t1 = 1.0\n",
    "dt0 = 1e-2\n",
    "logp = 'exact'\n",
    "\n",
    "# Runs backward-in-time to train the CNF\n",
    "def loss_fn(params, y, f):\n",
    "    if logp == 'exact':\n",
    "        term = dfx.ODETerm(logp_exact)\n",
    "    elif logp == 'approx':\n",
    "        term = dfx.ODETerm(logp_approx)\n",
    "    else:\n",
    "        raise NotImplementedError\n",
    "    solver = dfx.Heun()\n",
    "    eps = jax.random.normal(key, y.shape)\n",
    "    delta_log_likelihood = 0.0\n",
    "    y = (y, delta_log_likelihood)\n",
    "    func = lambda x: f.apply(params, x)\n",
    "    sol = dfx.diffeqsolve(term, solver, t1, t0, -dt0, y, (eps, func))\n",
    "    (z,), (delta_log_likelihood,) = sol.ys\n",
    "    log_prob = delta_log_likelihood + tfp.distributions.Normal(loc=0., scale=1.).log_prob(z).sum()\n",
    "    return - log_prob\n",
    "\n",
    "jax.vmap(loss_fn, in_axes=(None, 0, None))(params, x[:32], f)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = optax.adamw(learning_rate=3e-4, weight_decay=1e-4)\n",
    "opt_state = opt.init(params)\n",
    "\n",
    "@partial(jax.jit, static_argnums=(2,))\n",
    "def loss_fn_vmapped(params, x_batch, f):\n",
    "    loss = jax.vmap(loss_fn, in_axes=(None, 0, None))(params, x_batch, f)\n",
    "    return loss.mean()\n",
    "\n",
    "@partial(jax.jit, static_argnums=(3,))\n",
    "def train_step(params, opt_state, x_batch, f):\n",
    "    loss, grad = jax.value_and_grad(loss_fn_vmapped)(params, x_batch, f)\n",
    "    updates, opt_state = opt.update(grad, opt_state, params)\n",
    "    params = optax.apply_updates(params, updates)\n",
    "    return loss, params, opt_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2000/2000 [06:38<00:00,  5.02it/s, val=1.4153568]\n"
     ]
    }
   ],
   "source": [
    "n_steps = 2000\n",
    "n_batch = 32\n",
    "\n",
    "with trange(n_steps) as steps:\n",
    "    for step in steps:\n",
    "\n",
    "        # Draw a random batches from x\n",
    "        key, _ = jax.random.split(key)\n",
    "        idx = jax.random.choice(key, x.shape[0], shape=(n_batch,))\n",
    "\n",
    "        x_batch = x[idx]\n",
    "\n",
    "        loss, params, opt_state = train_step(params, opt_state, x_batch, f)\n",
    "        \n",
    "        # loss, grads = jax.value_and_grad(loss_fn_vmapped)(params, x_batch, f)\n",
    "        # updates, opt_state = opt.update(grads, opt_state, params)\n",
    "\n",
    "        # params = optax.apply_updates(params, updates)\n",
    "\n",
    "        steps.set_postfix(val=loss)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([ 0.60950206, -0.50103378], dtype=float64, weak_type=True)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def single_sample_fn(params, key, n_dim=2):\n",
    "    \"\"\" Produce single sample from the CNF by integrating forward.\n",
    "    \"\"\"\n",
    "    z = jax.random.normal(key, (n_dim,))\n",
    "    def func(t, x, args):\n",
    "        t = np.atleast_1d(t)\n",
    "        return f.apply(params, np.concatenate([x, t]))\n",
    "    term = dfx.ODETerm(func)\n",
    "    solver = dfx.Heun()\n",
    "    sol = dfx.diffeqsolve(term, solver, t0, t1, dt0, z)\n",
    "    (y,) = sol.ys\n",
    "    return y\n",
    "\n",
    "single_sample_fn(params, key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_fn = partial(single_sample_fn, params)\n",
    "\n",
    "n_samples = 100\n",
    "sample_key = jax.random.split(key, n_samples ** 2)\n",
    "x_sample = jax.vmap(sample_fn)(sample_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGgCAYAAABi2ofUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQXElEQVR4nO3de3RV1bnw/yfksgMk2RAIuUDQIAhWKFfB0A4ur1TwVjhajlqHgqVYfeX9iVAt+LZa29NBvZWeWo/oaQXt0draKrTY2mKEMiwRBcxQqCDQaMIlAUF2EjAXwvr90ddo8jzATNbeyV7J9zNGxpDpWnvNtffKZjKfZz4zwfM8TwAAAAKiW0d3AAAAoDUYvAAAgEBh8AIAAAKFwQsAAAgUBi8AACBQGLwAAIBAYfACAAAChcELAAAIFAYvAAAgUBi8AACAQInp4GXZsmVy0UUXSXp6uvTr109mzZolu3btOut5L7zwggwbNkxSU1NlxIgR8qc//SmW3QQAAAGSEMu9jWbMmCHXXXedXHTRRXLy5Em55557ZPv27fKPf/xDevbsaZ6zadMmmTRpkixbtkyuvPJKee655+SBBx6Qbdu2yfDhw896zVOnTsmBAwckPT1dEhISon1LAAAgBjzPk+rqasnLy5Nu3c48txLTwUtLhw8fln79+snf/vY3mTRpknnMtddeK8ePH5e1a9c2tV188cUyatQoWbFixVmvsW/fPsnPz49anwEAQPspLy+XAQMGnPGYpHbqi4iIRCIRERHJzMw87THFxcWyaNGiZm3Tp0+X1atXm8fX1dVJXV1d058/HYt9WS6XJEn22WPESlKLZ+Dk0aMd1BMAQDw4KQ3yuvxJ0tPTz3psuw1eTp06JQsXLpQvfelLZwz/VFRUSHZ2drO27OxsqaioMI9ftmyZ3H///ao9SZIlKYHBS7xK6pbSvIHPCgC6tv8XB3JJ+Wi3wcvtt98u27dvl9dffz2qr7t06dJmMzVVVVWEjQLg5JEjHd2FTiWpTx/VxnsMoLNql8HLggULZO3atbJx48azxrFycnKksrKyWVtlZaXk5OSYx4dCIQmFQlHrKwAAiG8xXSrteZ4sWLBAXnrpJXnttdekoKDgrOcUFhZKUVFRs7Z169ZJYWFhrLoJAAACJKYzL7fffrs899xzsmbNGklPT2/KWwmHw9K9e3cREbnpppukf//+smzZMhERueOOO2Ty5MnyyCOPyBVXXCHPP/+8bNmyRZ588slYdhUAAARETAcvjz/+uIiITJkypVn7ypUrZe7cuSIiUlZW1mw998SJE+W5556T7373u3LPPffIkCFDZPXq1U41XoBYSRo8SLWd3PPPDuiJP+TGAOgMYjp4cSkhs2HDBtU2e/ZsmT17dgx6BAAAgo69jQAAQKAweAEAAIHSrhV2gSAw80LiPL/Fylux7sPimgdDvgyAeMHMCwAACBQGLwAAIFAIGyFQ2iN00R6hED/34XouIR0AnRUzLwAAIFAYvAAAgEBh8AIAAAKFnBcESkflcUQ7R8XPcX7eAz/XiOZ7z7JrAH4w8wIAAAKFwQsAAAgUBi8AACBQyHlBl2bmgPQOqyZre4AOO3fwIH2u5eOIamo8L0+1Je494NY/A9sIAOgIzLwAAIBAYfACAAAChcELAAAIFHJe0Cn52v/HMT/DOjfByClxzW+xWPktVh7Mya+MU22pTlew82C8N9/VfYliHRpyYAD4wcwLAAAIFAYvAAAgUAgbIW65hn6iXX7flbns+Ohx1VY7pkC1Vecnq7b08gbVlhypU21WKCnBOM51mbWfpdxJVvjLWKKtrknYCIAPzLwAAIBAYfACAAAChcELAAAIFHJeELec8yLamHdx2nMNhydlu72epKmW43kJqi1nc71qq5iQotp6HrByY0Kq7ciF+tyU4YWqrT5D96XPjkzVllR6VLVZGjN7qrZEh+OsLx62GgDgipkXAAAQKAxeAABAoDB4AQAAgULOC2IumnkLVhn8pG2lqs2qweLq6HCdt2JxzWUR0XkrH9x4SrV51Y2qrSFd//viG3e9otp+8oerVNuJXKMn1brNzLXJ1zk+Vh0aS0NY5+Sktsih8ZPfEs1tCgAEEzMvAAAgUBi8AACAQCFshLjgGjJoGX4QEYlMHaLa0sp0mf7yS9NV2ye5OlTT/aAOBzWke6qtMbdWtd1z8zOqzfL/rbxFtV11zSbVtrb0QtX21EM6RNQ4RW8P4FXrcJX175Xkan2/H43VYa2PxupF0Bm7dVuDEXVLvjCn2Z/7r9XHWMuuXXf4BtC1xHTmZePGjXLVVVdJXl6eJCQkyOrVq894/IYNGyQhIUH9VFRUxLKbAAAgQGI6eDl+/LiMHDlSHnvssVadt2vXLjl48GDTT79+/WLUQwAAEDQxDRtddtllctlll7X6vH79+kmvXr2i3yEAABB4cZnzMmrUKKmrq5Phw4fL97//ffnSl7502mPr6uqkru6zeH9VVVV7dBGtYC1bbvng+VkSGy45pNpqC3TJe0vfrXrysT5DH9d7p855OfhVfdwvKiartsmZu1SblWvz4Ql9vyOydcj04PU6n+f43/R7XHu+zsmRGr2M2erLP//tSdV24X/dpl/PYOUHqffP2L4h0drSYfAgp2uayJcBOq24Wm2Um5srK1askN///vfy+9//XvLz82XKlCmybdu2056zbNkyCYfDTT/5+fnt2GMAANDe4mrmZejQoTJ06NCmP0+cOFH27t0ry5cvl1/96lfmOUuXLpVFixY1/bmqqooBDAAAnVhcDV4s48ePl9dff/20/z8UCkkopKfCAQBA5xT3g5eSkhLJzTXqnCMwEvceOOsxVn7L/huHqra0fbr+yIGpOsfi/KdPqLaEMZ+otuMSVm1WzsaJXF0Lxaqj8s6fdZ835w5WbXnr9eu9t1ufm1Kl+1Kfoc/taRyXs1lHhUtn6fcvIV2X/b++9H+pNlcpQ3TeWc3BFu+zUZun5+/e0C9m5UIZeTBmjRiDde7JPf90OhdA/Ijp4KWmpkb27NnT9OfS0lIpKSmRzMxMGThwoCxdulT2798vzzzzr8JeP/3pT6WgoEAuvPBCqa2tlV/84hfy2muvyV//+tdYdhMAAARITAcvW7ZskalTpzb9+dPclDlz5siqVavk4MGDUlZW1vT/6+vrZfHixbJ//37p0aOHfPGLX5RXX3212WsAAICuLcHzPD3fHGBVVVUSDodlisyUpASrPDramzVV33Ipc3W+22dl7eScXKOPu2C2Xp6895fnq7bz5r3vdJwVqrFCOq47TSdHdDl/i7XDtbW7s/X+WX22SvdbS6Vdy/5b5xas1m1HLmy+c/XJSXpZdNJGHcLrs0O/d6nGLuLWFhFmGApA3DrpNcgGWSORSEQyMoyaFZ8TV0ulAQAAzobBCwAACBQGLwAAIFDifqk02p+1bNlP+X5LaunRZn+uC+vNN60l0NYS42Nfr1Zt1pLl3lV6mbCV35K5XSfR1AzUS3HTynSZ/qw1elm49d41fGWcarNyOdLDBU7HJUf09gCJR3X/rG0TyqfpfBkrnyfF2Hmj9079eVTeppekizRv65uu+3aiSse4P7hRf2Z98vVnZuUBJYwfodoawromlPV+Wp8ZgPjBzAsAAAgUBi8AACBQGLwAAIBAIecFipnfYtRqkY91rY7aMTpHw9KyzolVkySpRrcdmKpzIPqu1bkSmdt1Hsz7c3qotvOf1rkXFiu/xWLdv5WPYtV5sd73urCuX5JqXNfagsHqi1UPJv9VnS/Ssi6LiEj/tRWqbf+VOaptwI/1v4lavvcnn0tXxxy7SDVJ6vv6btPLde2XurD+KrNyeYb8d6W+iIFtBID4xswLAAAIFAYvAAAgUBi8AACAQGFvIyhWvL8xU9c5sfIsLLvv1nU5eumthxQrD8Z1rxs798RtPyGrPorr/bvm/Fh9lt56bx9LZJSuiWPl5Lh+Po3n6Zwci1Ujxco1cfHxMP3vpoQxOofqkxqjLouRB9OQrr/GrGfMql8TLjl0um42Q84LEFvsbQQAADotBi8AACBQWCoNzVgCnWgcZoVIrCW2KUP068mu5lOC1nR+1kZjWavRNyvskbRuiz7XKBdvcrx/67ottz04Hetc1zBPeL3be2BewwiJWW2urKXXVqn+lscl6x0YpM/j3VWbtT2Ade4nufq4+gz9qR3P06HIUERvmWA+PwDiBjMvAAAgUBi8AACAQGHwAgAAAoWcly7OWhZdW6BzACxmfouRu1J1ME1ft8UyaCtPYvf8bNU25EEjH8XI2Yh87WLVFl6/W7VZJfmlTx/VZC6VNq5rLadNMHJtzHOtbRmMvlhcl3db/XMthW/dh2teUuq2Fg3GsvDDk/Tnnb7Nbbn8yUmfqLYq0c9dQrqxFUKNXnrdf5vxvltL2Y17NZ8pAFHFzAsAAAgUBi8AACBQGLwAAIBAIeelk7JyJaxYvJUXYdUqOf64rqNxslqfm7JWl3Tuu1WPkTO3Vzf7s5WzUfAdoxy7Y35G2MhFsJz8yjjVlmBsI2DWYDFyIMz8FsdzLVatFu/Nd1Wb9YvsOeatNBrnmrk2Pu6j5bNXeeNQdUzaPv2MWTlU1pYEJ6y8qhr93DWm675Z14hMHaLarG0EyG8BOgYzLwAAIFAYvAAAgEAhbNRJmctujZCLWMuMjV2La17W49xexjR/KKKXolphKHWNgToE5Rr6MUMcjuEMq2/WUnFrewBzGwHHPruyQkQW1/CF6+u5huesLxCXZ6//2gp1jPncDdDPnXVufYZeZm2Fgw5+VTXJiSurVNvxbfr5sbZlANAxmHkBAACBwuAFAAAECoMXAAAQKOS8BIxzfoeVe2G01Y4pcLpu1RC9oLbPDt2WbCwztnIZWpbqt5YEW/flWt7elfV6uli83RfzPbauYeWjOOaouC5597NE23oPXO/NYn4eDq8XLnFrs3KSsta8r9p2332+aht/3oeq7Z0/62XbvXfqfC7rdyW11PH9BBBVzLwAAIBAYfACAAACJaaDl40bN8pVV10leXl5kpCQIKtXrz7rORs2bJAxY8ZIKBSSwYMHy6pVq2LZRQAAEDAxzXk5fvy4jBw5Ur7xjW/I1VdffdbjS0tL5YorrpBbb71Vnn32WSkqKpJvfvObkpubK9OnT49lVwPDrKHheK5V8rw+I0G1HdMpAHLBcl0a3dpawMqzSJOz57OY+RkGK5/CyjGwyv5b+ThmTomVV+OYe2Id55oDYr2en/otJ60DXev/+KlX43hvqh+utXm2lTodN+j3Naqt/C39OyDDdNOhi/Tvxck0Xe3nguVOXQEQZTEdvFx22WVy2WWXOR+/YsUKKSgokEceeURERC644AJ5/fXXZfny5QxeAACAiMRZzktxcbFMmzatWdv06dOluLj4tOfU1dVJVVVVsx8AANB5xdVS6YqKCsnObl7mOzs7W6qqquSTTz6R7t27q3OWLVsm999/f3t1MS6Z4QYjbHLs69WqrX633gW6x0E9ZW4td7ZCTll7df8awiHV1jJM5BoecX1gk9ZtUW26WLyIWEvPLVZIwzXM4xhe8sP19cz+uW7D4Lgk3WWpsGv4z1oCbu2Wbb2etVN5/XC9+7S1LPrgV+v1NQ7q5xhAx4irmZe2WLp0qUQikaaf8vLyju4SAACIobiaecnJyZHKyspmbZWVlZKRkWHOuoiIhEIhCYX4FxEAAF1FXM28FBYWSlFRUbO2devWSWFhYQf1CAAAxJuYzrzU1NTInj17mv5cWloqJSUlkpmZKQMHDpSlS5fK/v375ZlnnhERkVtvvVV+/vOfy9133y3f+MY35LXXXpPf/va38vLLL8eym3HBz1JcKx/j4/xk1XbioF7q2WOITnBO2aXzYEIRvfA2rUwvPbb6YuWfNLbIZXBdJmuxHmJruwErB8J1uwHXfAxrybe59YHRF9ctA1y55hG5Hucnd6flca7L0c33yXF5u/X5ZG7Xy6etz+zQRXprgbwxB1Xb/iv1Z9t/rWpiywAgymI687JlyxYZPXq0jB49WkREFi1aJKNHj5Z7771XREQOHjwoZWVlTccXFBTIyy+/LOvWrZORI0fKI488Ir/4xS9YJg0AAJrEdOZlypQp4nnmGg8REbN67pQpU+Ttt9+OYa8AAECQxVXOCwAAwNnE1Wqjrsw5v8Vglem35K3XdVlE0lVLKNKg2syS7I55KmYOSYty9i1zYERakaNivE/tkd9ileS38jFcS/c7136x3nejVouf2i+uOVhtPdf1fXfNszFr0FhbThjv3f4b9X4Y/d7StV/K0rJV28Ar9DUOV+Xqy5Lzgi7O9TvFFTMvAAAgUBi8AACAQCFsFMdcl/tWTNShn6ohjaotY7deKp22T0+P14WNx2JMgWqylkAf/9rFqq3n795QbS23L7Bey1fpfiOMYoV+rNCCawjC3LXZWipuHeeD67JbX1sQOG6HYGnrVLCvUJXrRYznIvun+v20nuOkGv1vvcN/07+jOeXG1gI+psyjPd0OdMQzFe3XZ+YFAAAECoMXAAAQKAxeAABAoJDz0gGcl7AaS2wjRizevIYRn8/ZVO10rpVXY+UKWDkpaWX6XKtMYXKk+dYCiUaeSW1BpmpLWmfETa2+GZyXLDu9ms11KbefpfFmTo7je2C+nnXdKC+9Vsc4Llt3XSptXdN1+wbruuGSQ7ov0k+1HPu63l7jSI3OF0oLD9HXWK+v4JrjA/jhpzxCvGDmBQAABAqDFwAAECgMXgAAQKCQ8xIvrPogjuX3T+TqrJIeB62tADQrv8XabkBXiDlN3oJxXG2Lmi4iOufFYtV+MfMiDNY9WO+ndQ9+8lEsrnFj5zi0j1LzfvJv/NQlieZrWfWPrNwYK7/FtSaQpT5Dbw+Q/Xh31VYX1rWT6jPcfh+DlneATsRHbaeOwMwLAAAIFAYvAAAgUBi8AACAQCHnpQO47s0SmaprQ1j6vaVzXmoG6Bh7zUCdBxI28gI8I6ei0aiZ4VpvJHVbqT6uRXzVylFxZlwz0WgzcyV8sO7ftc6LFV/2s2eRr7oxcZJn4ZzL49o3x5wxsy/GNTK316i2o8PTVNuRKTqfK/cPKarNeh6tWkRAe4jm9097fKcw8wIAAAKFwQsAAAgUwkYx5meJbXj9btW2/8ahqs1aKt3vLb1c02RMrR83wlVWX6ywVs/ftW3q33yfjDZry4BUH9sDWJzL7/soye9aft/PkmI/S6BjvaTaKt1vfT5+tiRw3YLACt+YZQqMEGt6uEC1fVSdrNqspdIN4ZBqS3UMRQKuol3mwWn7j3YIOzPzAgAAAoXBCwAACBQGLwAAIFDIeYkiP7FFK/dk5729VFvq+/pUayuAtDK9rNPaCsCKp4eMvBI7hqlzXtqaA+EaI03dZjQa75259NrK2bByIKztEVyXXrcxRny6vph99vOcRVlbY9vWs2jlOJl5Nsb7ZH2RuebBmBy3oUgtParaMnbnqLb/9X82qbbXf3Sxaku2trWIk6Xs6DhRz1VzzOtzzU1zEe3nmJkXAAAQKAxeAABAoDB4AQAAgULOSxT5KdF+eFK2arvroj+qttU/n+r0emZOgWNtkYSILm+eaBxn1X6xWHHTRKNmhuJYQj/Bsey/1Y+TRvzWOq52jK7nYW574GfLAMfPxxfXbe+N41zzStrK9bWivY2CxbWekHVc2j5dY+m3myaotr5G7Rfz9xZdXnvkOFl1u8zv+BjnNbpi5gUAAAQKgxcAABAohI1aiPZyLtdzj+fpKeSnHrpKHzhcN6WXN6g2q/S4tQzTLMl+mj66cN4pt0V4xQytOF7TdQm4GYKxwjxGf1NdlyBa4ZYo7z5thjSsEJbFCk35KEnf5qXx7RD68bOU1FoCbX0WyUaI1WrLW69/9w5M1eGl9HLjs9U9MZ8zthHoHKL9d5D1emYpCYsVZvazLUoLqm+n6kWMXz0LMy8AACBQ2mXw8thjj8m5554rqampMmHCBHnzzTdPe+yqVaskISGh2U9qqvVvDwAA0BXFfPDym9/8RhYtWiT33XefbNu2TUaOHCnTp0+XQ4cOnfacjIwMOXjwYNPPhx9+GOtuAgCAgIh5zstPfvITmT9/vtx8880iIrJixQp5+eWX5amnnpIlS5aY5yQkJEhOji6x3R6ivZzLWnZ7dHiaaut5wFNtKVW6LRQ5qdrqwvpjtI4zc0NUi5234rS0WVpROrpF3NQ138Nc1uuYx2EtR89ao/dbcN0ywMo7sJg5Kk5nin1vVm6D65LqKOdKOG990KJ/fvJWXF6/Na/numzdYm0bYRqol6EWrG5UbdX5ycbJ+vlJWrfF7boInKgvizZ+563v81BYl4OItZb3etLT+ZunE9OZl/r6etm6datMmzbtswt26ybTpk2T4uLi055XU1Mj55xzjuTn58vMmTNlx44dpz22rq5Oqqqqmv0AAIDOK6aDl48++kgaGxslO7v5v3izs7OloqLCPGfo0KHy1FNPyZo1a+R//ud/5NSpUzJx4kTZt2+fefyyZcskHA43/eTn50f9PgAAQPyIu9VGhYWFctNNN8moUaNk8uTJ8uKLL0pWVpY88cQT5vFLly6VSCTS9FNeXt7OPQYAAO0ppjkvffv2lcTERKmsrGzWXllZ6ZzTkpycLKNHj5Y9e/aY/z8UCkkopGuaxAsrtnj80qGqrfdOXfPh42F6bHnu47qeR/XM81VbyAjFW3kbJ78yTrW51gyxYrOuWwG0PNd1i3bn3AYjj6O30dbounWB9d651ipxzFExc40c83nipX6Lq1i/fiyu6/q+W6wctPJpVn6LlrndbbsOS3u8p+gYrs+ytX1Kg59aUXEipjMvKSkpMnbsWCkqKmpqO3XqlBQVFUlhYaHTazQ2Nsq7774rubm5seomAAAIkJivNlq0aJHMmTNHxo0bJ+PHj5ef/vSncvz48abVRzfddJP0799fli1bJiIiP/jBD+Tiiy+WwYMHy7Fjx+Shhx6SDz/8UL75zW/GuqsAACAAYj54ufbaa+Xw4cNy7733SkVFhYwaNUpeeeWVpiTesrIy6dbtswmgjz/+WObPny8VFRXSu3dvGTt2rGzatEm+8IUvxLqrsWGVFK/Rh1nTyj0P6Glla3dja7mvdZy1VNp5F1vXXYaNpdLWcmw15WlMv0e7JLb1etbSbrO/1rJoH2EZc/mitZTdcQsC1zLeru+LeV0f9+sUJnQMQ7oux3dd8u+8TN91WbRxrlXOoN9buhRCzQA9GW6VVsjaqy/rZ1k54aX45md3+ZqBblsBmOUb4ngbinbZ22jBggWyYMEC8/9t2LCh2Z+XL18uy5cvb4deAQCAIIq71UYAAABnwuAFAAAESruEjboyc3nyJB077zdbF+1LvtOI2Rt5EVZsP7XUcV9xV0a837qutSzPzCtpmQNhxHRdY/N+jrOYOSCOcV7XpeKu+RPOpesdtypwvYaZA+Fj6bp6feP9NPvhmN9ift6Ox5kcc36srR+s0v31xtYU6eW6FPrBr+qSCcPu1Hu7WTltSeuiu+VCVxfE3CCr9IWVT2lvQ6GZz3ec5Lww8wIAAAKFwQsAAAgUBi8AACBQyHmJIitGWhvWWxdkP67P3TxrsGorCDeqtgbj9awaEj2NXAGztobBNV/ErEvSxnior1iyYw0aV67vkxUPtniO92bly1ift1XG28qtsrY0cP3MnHNZolgHoj3ymfzkMVj3lWy9x8Y1MrcbxZ0M3dN0zsthY/sPK18m2jli8aw97iHe3xO7fzoXKjmit5cQI+clMnWIaguv363P9VPfKoqYeQEAAIHC4AUAAAQKgxcAABAo5Ly0kZnfYtReKJ2VqM+t0WPGvlv1NVK37VFtVvzbjKcb+RM6g8auweK8h00U4+nO+ThGP1zjreY1rD2VHF8v1XE/Jue9c4x6MPrpsZl7VPnZL8oxl6UjvkD87EUU7TwGs4aPwcpdsnIRkjbq+0ip0nkwrjpDfosliPfgK9/K8TirfkuKkRNpSStz3OfOoUZVezx3zLwAAIBAYfACAAAChbCRA9cwSuo2fW7ffB3mSanyVFt9RoJqs0JEWRsrT9fN5vyENDqCj/76WepqfbauvxTW9giu51p9MbcWsLaDMJbnmlswGKXCrWfULDXvur2EFa5xCRP6eO7Mz9vxWWmPrSSsc63tOqxnIEf0M1B+abpq671TP2nJjltzdNZQUrxzfY9d/76xPkfr74fIqH6q7dhQfd30ch3aTHTcdqTl09gezxMzLwAAIFAYvAAAgEBh8AIAAAKFnBcXRj6GFW+0yitbuSwnrqxSbd42HVvM/2u1arNK0pvl4h3zMcyYvXGcJZqxcz9LBqOdB+PKep+s992SYBznp3S/3vTefi7M/JZ1W/TJrjkfbVym7nxfPmLn0T7XT36U+R1iLd03noFzn3dbBm8ub3fqnLtofkZdLffGOc/L+GytfDir9IW1JN9aAh2KuC3dd80T7IjPjZkXAAAQKAxeAABAoDB4AQAAgULOiwMrnmfV0LC2D096XtdosKT+VY8jKybqc/uvrXB6PatsebRLUfuJc7aM/3ZUboMfrvVGLK5bFTjXfHCtRVRq1G1w3OoiaZ2PfCOHkuLxlO9g3pdjno7Fujcr76lmoK7hYwlFdO6blbNg5UqcNGq/uIrmZxRPn3d78FNPyDUfLtX4Hjg8KVu1Hc/TuZhW7aBQWP8951wDKsaYeQEAAIHC4AUAAARKgud5ulZ9gFVVVUk4HJYpMlOSEvQOm2fjvKOy4fjXLlZtNQP0+NB1e4A+O+pVW52xQ6gVrnLmuHtwPGuPJZfOZeV97I7tuj2AFYKxlmi7nuvKuby5w3vl6/10vIdol2P3c67zMlnj99Eq7259r/T/1S63azjqamGdljpqKbdrGY5QxG0hfMWEFNXWe6feqdz1mYrm0nj12qfq5dWjKyUSiUhGRsYZj2XmBQAABAqDFwAAECgMXgAAQKCwVLoF12WoVtz90EU6b6Uxt1a1pb6fqtrOfV4vgba2AgiXHNLXMPIdPNflkJ0grh3NJdu+Gc+F61L7VMcliK5bEDiXkHfMv3F9Vlw+D+dtHtpjGwE/OV4+cogsjZl6qfSBqTpHLm+9zlmwlrebz5Tjc9EZvhv8iPb3irlc3nEpe8/f6XMrF05UbVY+ZXKN7p+VY5mzydiOxrFkggun7wWvwfn1mHkBAACBwuAFAAAECoMXAAAQKO2S8/LYY4/JQw89JBUVFTJy5Eh59NFHZfz48ac9/oUXXpDvfe978sEHH8iQIUPkgQcekMsvv7w9umqy4tCJxnEpQ6pU24mDaaqt5wEdl9w9X5dwHvR7I1hpsLYCaLRqhvjYMqCzape6DUbOlFVi28pxcs2DsWq6mLVfjGfA9fWsbQTMXAmHOjTt8Sw614hxvAeLax6DWO+7cV3rfc9b77ZlgFUDKtVH3pOrjqqHEmt+7sv1OLPsv2NfrK1i9l+Zo9oa9F9BIqJzXqy+uOZORnO7F1cxn3n5zW9+I4sWLZL77rtPtm3bJiNHjpTp06fLoUM68VREZNOmTXL99dfLvHnz5O2335ZZs2bJrFmzZPv27bHuKgAACICYD15+8pOfyPz58+Xmm2+WL3zhC7JixQrp0aOHPPXUU+bx//mf/ykzZsyQu+66Sy644AL54Q9/KGPGjJGf//znse4qAAAIgJiGjerr62Xr1q2ydOnSprZu3brJtGnTpLi42DynuLhYFi1a1Kxt+vTpsnr1avP4uro6qav7bDfVqiodumkNc6rQmDo7/I1C1TbgxzrMU36pDjAdN1a6DnlQl2F2XRLruruxWzHp4E0Dx1MZb+dQjRWqMPgKJTmGZqK9c7WpxXVdn0U/nJdZR3k5tmtozgpNWUtT08qMEJ51XWPa3yorb20nEuvfoXj/Tol2/1xfz7k0hyNrufPuBXqbHF3AQ0TW6ibn++iAzzKmMy8fffSRNDY2SnZ283yO7OxsqajQ8ToRkYqKilYdv2zZMgmHw00/+fn50ek8AACIS4FfbbR06VKJRCJNP+Xl5R3dJQAAEEMxDRv17dtXEhMTpbKysll7ZWWl5OTorGgRkZycnFYdHwqFJBQKRafDAAAg7sV08JKSkiJjx46VoqIimTVrloiInDp1SoqKimTBggXmOYWFhVJUVCQLFy5salu3bp0UFuock/Zixf2s8soVE9NVW+35ZnRRM3IgagbqJZJhawmrwbXsfbzEL/3oqP665hqZS4yNfIdka0m+8XmbeTDbSu1OunBdFuxj64xofkaur+Vaet01ru8n/m/+PjrmPbk+A9ZxVn6L63Yi0cx3cP3M4j03xuL6XeuaM+aaR2ZtMWLxqnXeZd56/feXHx2xVDrmdV4WLVokc+bMkXHjxsn48ePlpz/9qRw/flxuvvlmERG56aabpH///rJs2TIREbnjjjtk8uTJ8sgjj8gVV1whzz//vGzZskWefPLJWHcVAAAEQMwHL9dee60cPnxY7r33XqmoqJBRo0bJK6+80pSUW1ZWJt26fZZ6M3HiRHnuuefku9/9rtxzzz0yZMgQWb16tQwfPjzWXQUAAAHQLhV2FyxYcNow0YYNG1Tb7NmzZfbs2THuFQAACKJ2GbzEK9dYpcUqzRwZ1U+1VQ3Ra+z7bnVb5GXFq/3EEuM9dtyV2FvNb9HHOca1zfoyRg6E+QtvlYt37J9rfF7F+2OcFyPiL4/Dyk+wtgmxcpxcv1es17Nyl6zPNjlSp9pcmXVoXGuLxPg7JK7y16L9ej7yflyfgffn9FBtGbt1zkvNAH3dcInRGeO7ISnK20u0VeCXSgMAgK6FwQsAAAiULh02sjgvczSEIrroeUL6KdWWteZD1WaV8e75u7ZPxcVTiCielz9Ge5msxXy9dW7XSHAMD1ihADM0ZU3vGv1LNa5h7nbbxuWf0V5eap7r57kzwlqJVqjLONU1DGXu2Gv02TrOWiZrlVZIK9Nt1nGWnr97w+k4S4csne2g75lobyPgutv47vnZqi1J71AjM+ZsUm2vPepYeqQdwrttxcwLAAAIFAYvAAAgUBi8AACAQOnSOS9+8lusZY7V+XpZdOJBfe7hmeertqyNlfrAOM4VaY147rOfvBVXriX0reWQDWG3fbsqbhyq2nI2Vau2RNclsY5bBrgep36vfCzBjGbZ+tbwU+Le9X2yngHrS9rKtbFYz4+Vm2duL+G4ZNfKt2rrex/tHLRo58FEPa/G2iLCeFasnMhBv9cJLtYWNS+8O0a1DVvzvltf4hgzLwAAIFAYvAAAgEBh8AIAAAKlS+e8WJxzIIwaGsfz9DbjvXZ5qq0+Qx/nGttH8Lh+tmbZdiPHwKrTYW1XcXiSrgPR+ym3kvlmToo+yuSSF9BRJcZd79XKO/Dz+2h9rySMH9Hmc606QWLUfrFKyFt1Y6xtKKztIMz8mw7IaWuXXDXXZ6UdcqGsrWKsOk4ncvXfN0N+3qDa9hs5cv1/tUu1xXO+IjMvAAAgUBi8AACAQGHwAgAAAoWclxasmhxWTRdXVn6LVX9DRyptfvZ6QexFuw6EtadJyKjdUVuQqdqsHKxMK8/CuIbrniax3k/GTz0PV9HeK8k6zvy9tfY2sq7r+J2UWnpUn+y4H5NrnRdXLp+Hn9o8Zu0ki4/n2PX1fHGtrWJc18pnSq5OUW3vz+mh2i5YrnPkXPPc4iUPhpkXAAAQKAxeAABAoBA2asmYnks02hqMZWrJxnbkVojIKuHcf6/j9HucTNnBFvWS9EZIo25UP9WWVqZDP8k1ego52pyn4FtMSbuGavxM8Ud9etuYVk8wlrJb4SDri9ZPCNgq+28utXfkujWFFZqy3oOWochoh+b8bBvRYVsGOC7JN48z2qy/R3oe0AkIPQ+0fY4iXkJEFmZeAABAoDB4AQAAgcLgBQAABAo5Lw6seLC1TK3PDn3u0eFpqi2lylgYHbBlavDH9XOsXDjR8RX10tmqIY2qLWeT26tZz7xrvlWsc7XMPAZr6azrcY55B37yNkyOWxBYy6Kt7x+rz9a55pe+kS+T6JijYeUEqu0grPtyzRdyzYPpgKX8p+Oaz+T6rFilEKzf7+4HE1WblQcTbS3voz3eY2ZeAABAoDB4AQAAgcLgBQAABEqXznlxjocabce/drFqs2ptZBllt81YLyX+uzzreey/VpfxPjwp2+n1zn/6hGqzthuwmDVD/JTgj3EM3Pn3x7FcvGuOjp+cCue8CKdXOw0j5yUydYhqC6/f7dQXS1u3ZjCfxRhvSSAS/XwU1+fHej3Xz6La+J0feP5B1ZY7StcV2/vL8/V1jVpRPX/3hmprq/bINWLmBQAABAqDFwAAECidNmyUlJkpSd0+K4/uvLTOcSmlNbVnTQGGwsY2AsYyx45a0of4YT6P5pF6Ctl6HvffOFS19V/rFjZyLWXuZ9myi2hvBeAc0vETDnLYHqE1rDCz57ojtdEWPtr2vjiHXFocZx7j4z2xOH9f+nhWfD3bxv2GInqPb+vztsprrBr6P6rt1q/dptrSww2qzdyBPI7DwhZmXgAAQKAweAEAAIES08HL0aNH5YYbbpCMjAzp1auXzJs3T2pqjK2XP2fKlCmSkJDQ7OfWW2+NZTcBAECAxDTn5YYbbpCDBw/KunXrpKGhQW6++Wa55ZZb5LnnnjvjefPnz5cf/OAHTX/u0aNHq6998uhRkYTk1p9nxGYTxo9QbdZ28TUD9Fiw5++2qLZEK24KGKxYctaa91Vb7RidW5W275Rqc10iaT3z1rLbaJb9dxVPZeCjnWdh5ssYS4p1pkQrWNshWH3xsaTYZZm1n/IQfnKhXHOSzP65bulg5K3UGn9nJK3Tfz/UGGU4Dn61XrVdueVbqq3XQP07ai6D95HvGS+5mDEbvLz33nvyyiuvyFtvvSXjxo0TEZFHH31ULr/8cnn44YclL8/YO+X/6dGjh+Tk5MSqawAAIMBiFjYqLi6WXr16NQ1cRESmTZsm3bp1k82bN5/x3GeffVb69u0rw4cPl6VLl8qJE7rY1qfq6uqkqqqq2Q8AAOi8YjbzUlFRIf36NZ+iTkpKkszMTKmo0FVDP/X1r39dzjnnHMnLy5N33nlHvvOd78iuXbvkxRdfNI9ftmyZ3H///VHtOwAAiF+tHrwsWbJEHnjggTMe895777W5Q7fcckvTf48YMUJyc3Plkksukb1798p5552njl+6dKksWrSo6c9VVVWSn5+vjot2TNyq1ZKzSbclRrlOBTovM+ZsMeLQVt2Gfk/r4975s679kuaY05VqbHVh1oZwzKmIKqtmiI/fKedcCVeuORWO9XVcv8983Yfjddv63vvKPXE8Lur3anDdHiDV+P1uNH73rNovfTbo38ce1+vf+ZOSfrpunpWv968DtHrwsnjxYpk7d+4Zjxk0aJDk5OTIoUOHmrWfPHlSjh492qp8lgkTJoiIyJ49e8zBSygUklBIf7AAAKBzavXgJSsrS7Kyss56XGFhoRw7dky2bt0qY8eOFRGR1157TU6dOtU0IHFRUlIiIiK5ubmt7SoAAOiEYpbzcsEFF8iMGTNk/vz5smLFCmloaJAFCxbIdddd17TSaP/+/XLJJZfIM888I+PHj5e9e/fKc889J5dffrn06dNH3nnnHbnzzjtl0qRJ8sUvftFXf/yEZaxS3K7LCK1pQXPHXnR5LmXWRU4THjCWSh/+pd5GoHZKrb7wX3WTWT7cBzX97Djt71z233UHZOP9rC3I1McZS1hdl5f62hnaCrlFOUTia7sBx2XRJ78yrtmfzefJulfXpcjWsn1HfpYJW/xsB2HtrG2FbD8aq8seDDRe78BUvY1AWpmxsjfaYdEOENMidc8++6wMGzZMLrnkErn88svly1/+sjz55JNN/7+hoUF27drVtJooJSVFXn31Vbn00ktl2LBhsnjxYrnmmmvkj3/8Yyy7CQAAAiSmReoyMzPPWJDu3HPPFc/7bKSYn58vf/vb32LZJQAAEHDsbQQAAAIlpjMvHSkpM1OSuqU0/dlPSfG2lr8WsWPJvkp7o0txfc6snILkiJUXkKZaagamqLaQEXdPNvIMXPPBWuYPOOeoRDkOb1031SqB7vqCPpYxW1yX3fp5X8ycO9fXc1wWrZ5HxzLzrvdlPnd+nhXH/BY/eU/mlhuG8mnWtjY65+XYyzqX5YK1p6+h1oyPv/viBTMvAAAgUBi8AACAQGHwAgAAAqXT5rycPHpUJMGKHX7uGB9xaNdt0C1WzJntAeAnVyIyqp9qqxmg/21yIlfXgRj75d2qrXz5ENVm1upwrXPS4t5c8wRc7981Nu/nPY72FiN+aovUGnV9nGupuJbgNyQ4fse5vJ6f7QF8fY6O77ufHBArv8Wq6bL/Sl1tfsh/67yV9+7UhWGt3+Vo30c8/73EzAsAAAgUBi8AACBQGLwAAIBA6bQ5Ly1FO87pWo/BQp0XWPzkXvT83Ruqrf4bhcbZCaplc8lg1VYQ0U+p6x5Asebnd9l8vWjnsri+no/aIqnbjAMd99NxzW+xWHkb1vdZy/fA3I/L+q51vQc/uVCOeUCufbFeL9E4zPr9ydlUrdoOT9J7khWsblBtyZETbv0zxHMuiytmXgAAQKAweAEAAIHSZcJGziEia0ra0gmm3dC5pZfrqeYjF+qtACzWcdYUd8NXxqk2l1CSa/l0c3sNH1sLOIcv/Gwd4vrd4GNZsMX5fWnj8nYRcd4Ooc0cQ2muy6xdr+EaErWuay1bT47UqTZrKfvu+TpE1Jirzz12UG/XMeS/jaXxXQgzLwAAIFAYvAAAgEBh8AIAAAKly+S8OG8F4CPmDLQLa1mncZi1hLX/r3aptrSpeiuA+gyj9HiMWctwre0H/CyLtjjnT/jIl7G45qhE/d58LNFu63X9fIe6ntuY2VO1WTlTluSI3vbA+g2IGL8r4ZJD+kDjPbbOHfT7GtVW+yOdW9bzB47zDK7fDdHO3+oAzLwAAIBAYfACAAAChcELAAAIlC6T8+JHPMf90PW45mWlbitVbYdnnq/ajufpLQNyNtertoqJ6aotpUpnBqRbtV9a1Lgw4/A+ytb74WerD9fXM+tHuR7np0ZVlPODLNHMCYz29izWe2LlxiTuPaCvcbpOtmDViKm4Lsfp3FBE1285+YTuX22B7o31+23WDuoE+S0WZl4AAECgMHgBAACBQtjIQWeYYkPn5vo8Zq15X7Udv22oaiuflqzaeu3SISIrbGSVQW8Zvoj2js9+jvPD+XvAsay+GRJ0DDn5Odf5/XMMTbXUeJ5eimwyrmme6+N9t3Z8trg+K9ZWAL136r9aawbouYKKCXobjnMf1+UMrPfAT7izM/z9xcwLAAAIFAYvAAAgUBi8AACAQOnSOS+ucd7OutQMXY/1jJ77fIVq2z0/W7WdN0/ny3x8p47F779SLxPtv7b5n60vntoxBaotwcgncI3r+8lvcf2d9/Md4oeVe2IuNXfNtYly/9Tr+fnMfCxjtnJFrC0DXJ8VaysA63mvHlOr2hIP6mXRJmO5s7V1hut70Fkx8wIAAAKFwQsAAAgUBi8AACBQunTOi584L/ktCKKE8SN0o5FT0OOgjuNvLhmsj1tSo9oG/LhatUVG9Wv25/B6nYthlju3GPkJVm5DQ1jnGLheI55y2vxsmxD1PkdxuwFftWVc8xCN/BaL9fxYeSbW1gI5m/TzXjWkh+5L2inVdsHyw7ozDnVzEMOZlx/96EcyceJE6dGjh/Tq1cvpHM/z5N5775Xc3Fzp3r27TJs2TXbv3h2rLgIAgACK2eClvr5eZs+eLbfddpvzOQ8++KD87Gc/kxUrVsjmzZulZ8+eMn36dKmt1ZnbAACga0rwPE/X946iVatWycKFC+XYsWNnPM7zPMnLy5PFixfLt7/9bRERiUQikp2dLatWrZLrrrvO6XpVVVUSDodlisyUpARd4hzoKpyXChtLMy2HJ+nl0xaXXaqtkurWzr6WeA/ZtsdWBdEs3X+617NE9d6sHZAdQ2TO9+p4jeNfu9jpulaJ/z479A7s1fn67x7r96L3Th1KspZj+wklxfvvy+ed9Bpkg6yRSCQiGRkZZzw2bhJ2S0tLpaKiQqZNm9bUFg6HZcKECVJcXHza8+rq6qSqqqrZDwAA6LziZvBSUfGvQlnZ2c3/dZednd30/yzLli2TcDjc9JOfnx/TfgIAgI7VqsHLkiVLJCEh4Yw/O3fujFVfTUuXLpVIJNL0U15e3q7XBwAA7atVS6UXL14sc+fOPeMxgwYZsUgHOTn/WppZWVkpubm5Te2VlZUyatSo054XCoUkFHIsuwx0Ic6xbsechayN+tSMp3Us/p0/D1VtFRNSmv255wGdE5Ae1tsDWLkxVn+tJeCuOTSuor5NiJVr5Jjb4Jov4yfXxsoDCa93XP3Z8t583Jf5fhqvF5k6RLW59tc6zlo+XZ+RptqOXJii2hr0YZL/V72k2nxGrTyddtgSI2haNXjJysqSrKysmHSkoKBAcnJypKioqGmwUlVVJZs3b27ViiUAANC5xSznpaysTEpKSqSsrEwaGxulpKRESkpKpKbms6JWw4YNk5deeklERBISEmThwoXyH//xH/KHP/xB3n33XbnpppskLy9PZs2aFatuAgCAgIlZhd17771Xnn766aY/jx49WkRE1q9fL1OmTBERkV27dkkk8tn039133y3Hjx+XW265RY4dOyZf/vKX5ZVXXpHU1NRYdRMAAARMzOu8tDfqvAD+WbHzwzPPdzo3pUp/pXw8rPkkb+35uvBk6vv6HynJevcB6f+rXaqtdozOl0lat0W3WfVBLD7K4PvJR3HVHtsXdERNF+vz8VP7pbYgU7Wllh5VbVbZf8s/r9HJLFbZ/4zdiarNem4tQarLEm2BrPMCAADggsELAAAIlC69qzQAd72f0pWuSx8oVG2NubpcesuQkBUisvRfqwtUWuGrrI2V+mQjjGCFB6zdgy1meMRHOMT1y9c5XOUYcnG9D+eS9A7bS7iGfqzjrGXwnrVbtNVfI2xkhZLqwvrTaBnqFBE5mdao2nrk6thmn9XdVdv+G3UJAddQEjRmXgAAQKAweAEAAIHC4AUAAAQKOS9AF+dn2e2Q/9a5JlZeydHhzbfwSC9vUMdU5+vSBvuvzFFtVun13o65EibH3A6rXLyVe+Gce+InR8XgnFdifbY+thFwubeTXxmnj7GWshvX1Fkm9mdbG9bbxFi5LIcuSlBtvRxTTwaer5/31P+brtoS95aqtux1+j0+6XZZGJh5AQAAgcLgBQAABAqDFwAAECjkvABdXHuUI29Zh8XKi+lx90f6vNv0v692z89WbVa9GSsfx2TUKYmM6qfaev7uDbfXc+Sa32J9PlYOSWqpY60Wx5oufurLtJQcqdONRt6Kld9iSdx7QLU1GFtEhEsOqbaaATqP6qOx+sp563VuzGHReU/5Un3afkZLe2wHETTMvAAAgEBh8AIAAAKFwQsAAAgUcl4AtJlr3kZti3wEq/5G2v/Vr//BdbqGRr+3Tqk2ax+aw5N0bkxKlafaQhFdbcNqs2qLNBi1RZK26RofrjVdXFk5JOa+Tca5rvsHmbVf2libxuqH1V9LzUB9XPioztux3hPXOkGucjbrfbus/Bszr8g1h8gxv6Wr58Ew8wIAAAKFwQsAAAgUwkYAlGhPXSe1KI2eapWeN6baczbrsEy/7+uwzLFKHR7I/Kv+t1n5pToMlbNZd8XaqkCMtsztNaqtZYhMRCTVCCWZoZqjx3Wb8X5a4arU0qOqzXVbgogVmtlrlOo3Qj3mUmnHLRdaskIw4b36OOs9trYbuPVpvfx+9Zypqs0KTdUM0M9PWpkOTZmft/FZuHIN/XSlEJGFmRcAABAoDF4AAECgMHgBAACBkuB5nl47GGBVVVUSDodlisyUpAQjbg2g3ancGCO/xVrWmrNJl1638j0sRy5MUW1VQ3QZ+IzdeiGvtZy29069RLs+Q5eQz1rzvnGyvl9rKbeVQ5No5MFY2xeE1+92uq6ltiBTtVl5G21d3pxWpu/Bym85PPN81ZZe3qDaSmfpz8wq528toU8Yo/Nxzrn7E9Vmcszl6er5KG110muQDbJGIpGIZGRknPFYZl4AAECgMHgBAACBwuAFAAAECnVeAESVS+0X64un/1rdZtUpsWrE7L5b50pY2whUj9H5E1bx+nOfr1BtH1ync3J6HtApg/tvHGpct1a15f5Bn2vl81g1SKxcG4u5ZYCRQ2PVtUnVpWnEe/Nd1WbVq2mZa1IzQNfXSRs4RLVZ91U+TfctSacGSShifbY672nAj93+zW7mARk5L5Tu7xjMvAAAgEBh8AIAAAKFpdIAAsV1Z2Oz/L7BCstYjn1dL9s+cVCvqS5YrZdjW8u2rWXgrn2xWMuRjw7X/bN21g6XHFJt1tJ1V312NN992VrabLGWO1tl+vuv1WE9K8xj7TRthc1MjiEiRA9LpQEAQKfF4AUAAAQKgxcAABAoMVsq/aMf/UhefvllKSkpkZSUFDl27NhZz5k7d648/fTTzdqmT58ur7zySox6CSBeuC4vtZb/WufWGsuOrZL3xk4Ap8mL0CX5e5hLe/W/CRtz9VLp/mv1NQ5dY+SoDKlSbXaujb7f43m6f1bOi7XdwIlcfVyvXarJXN7c7/vN11lX/FkvH689X78n1tJma6sGq7+hyEnVZm1B4Lq0GfEtZoOX+vp6mT17thQWFsovf/lL5/NmzJghK1eubPpzKOS2jwkAAOgaYjZ4uf/++0VEZNWqVa06LxQKSU5O27PcAQBA5xZ3OS8bNmyQfv36ydChQ+W2226TI2dZmlZXVydVVVXNfgAAQOcVV9sDzJgxQ66++mopKCiQvXv3yj333COXXXaZFBcXS2KiXSdg2bJlTbM8AILLzEUwarqIkcdgsWp8WLU7Eo22xvPynK5h5ZRY+S3n/kr/O/G9O7NUW8Zu/XpVaTrpo0euro/f7/u69smRF3SuScvS/SIin+Tq2jQZu/VxKVV6ywXrPdj6evPtGhqN/JY+G6yUAJ1nY23BYLHymfzUZaGmS3xr1czLkiVLJCEh4Yw/O3fubHNnrrvuOvnqV78qI0aMkFmzZsnatWvlrbfekg0bNpz2nKVLl0okEmn6KS8vb/P1AQBA/GvVzMvixYtl7ty5Zzxm0CDjX0ptNGjQIOnbt6/s2bNHLrnkEvOYUChEUi8AAF1IqwYvWVlZkpWlpzpjZd++fXLkyBHJzc1tt2sCiCM+SrSbX269w6rJWnptldXP2lhpvKBespv2V727sbVbtFUKvz5Dh0j6btUT5Onl3VXbISnQ3btQN33xMr3e+eM7dZjsn8aybWuyPtnY4fnc55u/V/bSZv0+WaEfk/Fc1I7R959kbBvBjs+dQ8wSdsvKyqSkpETKysqksbFRSkpKpKSkRGpqPnvShw0bJi+99JKIiNTU1Mhdd90lb7zxhnzwwQdSVFQkM2fOlMGDB8v06dNj1U0AABAwMUvYvffee5sVnBs9erSIiKxfv16mTJkiIiK7du2SSORfI+jExER555135Omnn5Zjx45JXl6eXHrppfLDH/6QsBAAAGgSs8HLqlWrzlrj5fMbWnfv3l3+8pe/xKo7AACgk4irpdIA8HnRzkU4aeVAGMuxszbq0v1Wbky45JBqcy1dn1amr2GVs49MHaLaLNX5yarN2gpg7y/PV23HrtGv1+Ogzsnp/yudL3N4pn69lu+Vda8W18/HWsqetG6L07nWNRA8cVekDgAA4EwYvAAAgEBh8AIAAAKFnBcAnZJrvoyZZ2HUArG2EbCES4xGxy0IrDYrX8aqh5K6zeifUdemtiBTtWVt1K9nHWfl36SX63otiUdb5LhY753RN7FqsDjmqJj1Wxw/Wwu1X+IbMy8AACBQGLwAAIBAYfACAAAChZwXAIHnZ78a13NdcyXM/A5H3pvvqrbk8SOczrXyZaw9lVK3leqTjfwT67hUK0/F6kuLOi+elXvi9ErRz1Ehl6VzYOYFAAAECoMXAAAQKISNAASen1BAtMMNVpjDCulYWwHoRdF2KMlaUmyFZlKN8vgmI9Rlhs6MsJHLcmTrvqId1kPXwswLAAAIFAYvAAAgUBi8AACAQCHnBQBaiPbyXDO/Jcp5G77K41tLoI3+Wa938ivj9LnrttidPEs/rPfEddk6eTBdCzMvAAAgUBi8AACAQGHwAgAAAoWcFwBooT3qxjhvNxDla5g5JG3uiUiSkd+i6rw45rK4Ig8GzLwAAIBAYfACAAAChbARAHQAP6EPP8uMo94XYwsCa0l1rBEi6lqYeQEAAIHC4AUAAAQKgxcAABAo5LwAQJzwk6MSbc7X+Dhy1kPiJS8GnQczLwAAIFAYvAAAgEBh8AIAAAKFnBcACJh4KoXvcl3yWxBtzLwAAIBAidng5YMPPpB58+ZJQUGBdO/eXc477zy57777pL6+/ozn1dbWyu233y59+vSRtLQ0ueaaa6SysjJW3QQAAAETs7DRzp075dSpU/LEE0/I4MGDZfv27TJ//nw5fvy4PPzww6c9784775SXX35ZXnjhBQmHw7JgwQK5+uqr5e9//3usugoAgUIpfHR1CZ7nee11sYceekgef/xx+ec/7fhnJBKRrKwsee655+RrX/uaiPxrEHTBBRdIcXGxXHzxxWe9RlVVlYTDYZkiMyUpITmq/QcAALFx0muQDbJGIpGIZGRknPHYds15iUQikpmZedr/v3XrVmloaJBp06Y1tQ0bNkwGDhwoxcXF5jl1dXVSVVXV7AcAAHRe7TZ42bNnjzz66KPyrW9967THVFRUSEpKivTq1atZe3Z2tlRUVJjnLFu2TMLhcNNPfn5+NLsNAADiTKsHL0uWLJGEhIQz/uzcubPZOfv375cZM2bI7NmzZf78+VHrvIjI0qVLJRKJNP2Ul5dH9fUBoL0k9emjfgBorU7YXbx4scydO/eMxwwa9Nk+FgcOHJCpU6fKxIkT5cknnzzjeTk5OVJfXy/Hjh1rNvtSWVkpOTk55jmhUEhCoZBz/wEAQLC1evCSlZUlWVlZTsfu379fpk6dKmPHjpWVK1dKt25nnugZO3asJCcnS1FRkVxzzTUiIrJr1y4pKyuTwsLC1nYVAAB0QjHLedm/f79MmTJFBg4cKA8//LAcPnxYKioqmuWu7N+/X4YNGyZvvvmmiIiEw2GZN2+eLFq0SNavXy9bt26Vm2++WQoLC51WGgEAgM4vZnVe1q1bJ3v27JE9e/bIgAEDmv2/T1dnNzQ0yK5du+TEiRNN/2/58uXSrVs3ueaaa6Surk6mT58u//Vf/xWrbgJA3KB+C+CmXeu8tAfqvAAAEDxxW+cFAADALwYvAAAgUBi8AACAQGHwAgAAAoXBCwAACBQGLwAAIFAYvAAAgEBh8AIAAAKFwQsAAAgUBi8AACBQGLwAAIBAYfACAAAChcELAAAIlKSO7gAAAEGQ1KePajt55EgH9ATMvAAAgEDpdDMvnueJiMhJaRDxOrgzAIDO41S9ajrpNXRARzqnk/Kv9/LTv8fPpNMNXqqrq0VE5HX5Uwf3BADQqRzt6A50DdXV1RIOh894TILnMsQJkFOnTsmBAwckPT1dEhISpKqqSvLz86W8vFwyMjI6unvtpived1e8Z5Gued9d8Z5FuO+udN9d8Z49z5Pq6mrJy8uTbt3OnNXS6WZeunXrJgMGDFDtGRkZXeYB+LyueN9d8Z5FuuZ9d8V7FuG+u5Kuds9nm3H5FAm7AAAgUBi8AACAQOn0g5dQKCT33XefhEKhju5Ku+qK990V71mka953V7xnEe67K913V7zn1uh0CbsAAKBz6/QzLwAAoHNh8AIAAAKFwQsAAAgUBi8AACBQOt3g5YMPPpB58+ZJQUGBdO/eXc477zy57777pL5e70nxebW1tXL77bdLnz59JC0tTa655hqprKxsp17796Mf/UgmTpwoPXr0kF69ejmdM3fuXElISGj2M2PGjNh2NMract+e58m9994rubm50r17d5k2bZrs3r07th2NoqNHj8oNN9wgGRkZ0qtXL5k3b57U1NSc8ZwpU6aoz/rWW29tpx63zWOPPSbnnnuupKamyoQJE+TNN9884/EvvPCCDBs2TFJTU2XEiBHypz8Fc4uQ1tz3qlWr1Oeamprajr31b+PGjXLVVVdJXl6eJCQkyOrVq896zoYNG2TMmDESCoVk8ODBsmrVqpj3M9pae98bNmxQn3VCQoJUVFS0T4fjTKcbvOzcuVNOnTolTzzxhOzYsUOWL18uK1askHvuueeM5915553yxz/+UV544QX529/+JgcOHJCrr766nXrtX319vcyePVtuu+22Vp03Y8YMOXjwYNPPr3/96xj1MDbact8PPvig/OxnP5MVK1bI5s2bpWfPnjJ9+nSpra2NYU+j54YbbpAdO3bIunXrZO3atbJx40a55ZZbznre/Pnzm33WDz74YDv0tm1+85vfyKJFi+S+++6Tbdu2yciRI2X69Oly6NAh8/hNmzbJ9ddfL/PmzZO3335bZs2aJbNmzZLt27e3c8/9ae19i/yrAuvnP9cPP/ywHXvs3/Hjx2XkyJHy2GOPOR1fWloqV1xxhUydOlVKSkpk4cKF8s1vflP+8pe/xLin0dXa+/7Url27mn3e/fr1i1EP45zXBTz44INeQUHBaf//sWPHvOTkZO+FF15oanvvvfc8EfGKi4vbo4tRs3LlSi8cDjsdO2fOHG/mzJkx7U97cb3vU6dOeTk5Od5DDz3U1Hbs2DEvFAp5v/71r2PYw+j4xz/+4YmI99ZbbzW1/fnPf/YSEhK8/fv3n/a8yZMne3fccUc79DA6xo8f791+++1Nf25sbPTy8vK8ZcuWmcf/+7//u3fFFVc0a5swYYL3rW99K6b9jLbW3ndrft+DQES8l1566YzH3H333d6FF17YrO3aa6/1pk+fHsOexZbLfa9fv94TEe/jjz9ulz7Fu04382KJRCKSmZl52v+/detWaWhokGnTpjW1DRs2TAYOHCjFxcXt0cUOs2HDBunXr58MHTpUbrvtNjly5EhHdymmSktLpaKiotlnHQ6HZcKECYH4rIuLi6VXr14ybty4prZp06ZJt27dZPPmzWc899lnn5W+ffvK8OHDZenSpXLixIlYd7dN6uvrZevWrc0+o27dusm0adNO+xkVFxc3O15EZPr06YH4TD/VlvsWEampqZFzzjlH8vPzZebMmbJjx4726G6H6QyftR+jRo2S3Nxc+cpXviJ///vfO7o7HabTbczY0p49e+TRRx+Vhx9++LTHVFRUSEpKisqZyM7O7tTxxBkzZsjVV18tBQUFsnfvXrnnnnvksssuk+LiYklMTOzo7sXEp59ndnZ2s/agfNYVFRVqmjgpKUkyMzPP2P+vf/3rcs4550heXp6888478p3vfEd27dolL774Yqy73GofffSRNDY2mp/Rzp07zXMqKioC+5l+qi33PXToUHnqqafki1/8okQiEXn44Ydl4sSJsmPHDnOD2s7gdJ91VVWVfPLJJ9K9e/cO6lls5ebmyooVK2TcuHFSV1cnv/jFL2TKlCmyefNmGTNmTEd3r90FZuZlyZIlZrLS539a/oLv379fZsyYIbNnz5b58+d3UM/bri333BrXXXedfPWrX5URI0bIrFmzZO3atfLWW2/Jhg0boncTbRDr+45Hsb7nW265RaZPny4jRoyQG264QZ555hl56aWXZO/evVG8C7S3wsJCuemmm2TUqFEyefJkefHFFyUrK0ueeOKJju4aomzo0KHyrW99S8aOHSsTJ06Up556SiZOnCjLly/v6K51iMDMvCxevFjmzp17xmMGDRrU9N8HDhyQqVOnysSJE+XJJ58843k5OTlSX18vx44dazb7UllZKTk5OX667Utr79mvQYMGSd++fWXPnj1yySWXRO11WyuW9/3p51lZWSm5ublN7ZWVlTJq1Kg2vWY0uN5zTk6OSt48efKkHD16tFXP6oQJE0TkXzOT5513Xqv7G0t9+/aVxMREtdrvTL+POTk5rTo+HrXlvltKTk6W0aNHy549e2LRxbhwus86IyOj0866nM748ePl9ddf7+hudIjADF6ysrIkKyvL6dj9+/fL1KlTZezYsbJy5Urp1u3ME0xjx46V5ORkKSoqkmuuuUZE/pXRXVZWJoWFhb773latuedo2Ldvnxw5cqTZX+odIZb3XVBQIDk5OVJUVNQ0WKmqqpLNmze3eqVWNLnec2FhoRw7dky2bt0qY8eOFRGR1157TU6dOtU0IHFRUlIiItLhn7UlJSVFxo4dK0VFRTJr1iwRETl16pQUFRXJggULzHMKCwulqKhIFi5c2NS2bt26Dv39ba223HdLjY2N8u6778rll18ew552rMLCQrUMPmifdbSUlJTE5e9wu+jojOFo27dvnzd48GDvkksu8fbt2+cdPHiw6efzxwwdOtTbvHlzU9utt97qDRw40Hvttde8LVu2eIWFhV5hYWFH3EKbfPjhh97bb7/t3X///V5aWpr39ttve2+//bZXXV3ddMzQoUO9F1980fM8z6uurva+/e1ve8XFxV5paan36quvemPGjPGGDBni1dbWdtRttFpr79vzPO/HP/6x16tXL2/NmjXeO++8482cOdMrKCjwPvnkk464hVabMWOGN3r0aG/z5s3e66+/7g0ZMsS7/vrrm/5/y+d7z5493g9+8ANvy5YtXmlpqbdmzRpv0KBB3qRJkzrqFs7q+eef90KhkLdq1SrvH//4h3fLLbd4vXr18ioqKjzP87wbb7zRW7JkSdPxf//7372kpCTv4Ycf9t577z3vvvvu85KTk7133323o26hTVp73/fff7/3l7/8xdu7d6+3detW77rrrvNSU1O9HTt2dNQttFp1dXXT762IeD/5yU+8t99+2/vwww89z/O8JUuWeDfeeGPT8f/85z+9Hj16eHfddZf33nvveY899piXmJjovfLKKx11C23S2vtevny5t3r1am/37t3eu+++691xxx1et27dvFdffbWjbqFDdbrBy8qVKz0RMX8+VVpa6omIt379+qa2Tz75xPvf//t/e7179/Z69Ojh/du//VuzAU+8mzNnjnnPn79HEfFWrlzpeZ7nnThxwrv00ku9rKwsLzk52TvnnHO8+fPnN31JBkVr79vz/rVc+nvf+56XnZ3thUIh75JLLvF27drV/p1voyNHjnjXX3+9l5aW5mVkZHg333xzs8Fay+e7rKzMmzRpkpeZmemFQiFv8ODB3l133eVFIpEOugM3jz76qDdw4EAvJSXFGz9+vPfGG280/b/Jkyd7c+bMaXb8b3/7W+/888/3UlJSvAsvvNB7+eWX27nH0dGa+164cGHTsdnZ2d7ll1/ubdu2rQN63XafLgFu+fPpfc6ZM8ebPHmyOmfUqFFeSkqKN2jQoGa/30HR2vt+4IEHvPPOO89LTU31MjMzvSlTpnivvfZax3Q+DiR4nufFfn4HAAAgOgKz2ggAAECEwQsAAAgYBi8AACBQGLwAAIBAYfACAAAChcELAAAIFAYvAAAgUBi8AACAQGHwAgAAAoXBCwAACBQGLwAAIFAYvAAAgED5/wEDjeQgN3hMCgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.hist2d(x_sample[:, 0], x_sample[:, 1], bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch-mps",
   "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.9.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}