model_numpyro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import jax.numpy as jnp\n",
    "from jax import random, lax, jit, vmap\n",
    "\n",
    "import numpyro\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.infer.autoguide import AutoBNAFNormal, AutoIAFNormal\n",
    "from numpyro import optim\n",
    "from numpyro.infer import SVI, Trace_ELBO, MCMC, NUTS, init_to_value\n",
    "from numpyro.infer.reparam import NeuTraReparam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "lc = np.loadtxt(\"lc.txt\")\n",
    "interp = np.loadtxt(\"interp.txt\")\n",
    "t, F, Ferr = lc\n",
    "points, zeta_e, zeta_n, zeta_e_dot, zeta_n_dot = interp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(t, F, Ferr):\n",
    "    ln_DeltaF = numpyro.sample(\"ln_DeltaF\", dist.Normal(4.0, 4.0))\n",
    "    DeltaF = jnp.exp(ln_DeltaF)\n",
    "    ln_Fbase = numpyro.sample(\"ln_Fbase\", dist.Normal(2.0, 4.0))\n",
    "    Fbase = jnp.exp(ln_Fbase)\n",
    "    t0 = numpyro.sample(\"t0\", dist.Normal(3630.0, 50.0))\n",
    "    ln_tE = numpyro.sample(\"ln_tE\", dist.Normal(3.0, 6.0))\n",
    "    tE = jnp.exp(ln_tE)\n",
    "    u0 = numpyro.sample(\"u0\", dist.Normal(0.0, 1.0))\n",
    "    piEE = numpyro.sample(\"piEE\", dist.Normal(0.0, 0.5))\n",
    "    piEN = numpyro.sample(\"piEN\", dist.Normal(0.0, 0.5))\n",
    "\n",
    "    # Compute trajectory\n",
    "    zeta_e_t = jnp.interp(t, points, zeta_e)\n",
    "    zeta_n_t = jnp.interp(t, points, zeta_n)\n",
    "    zeta_e_t0 = jnp.interp(t0, points, zeta_e)\n",
    "    zeta_n_t0 = jnp.interp(t0, points, zeta_n)\n",
    "    zeta_e_dot_t0 = jnp.interp(t0, points, zeta_e_dot)\n",
    "    zeta_n_dot_t0 = jnp.interp(t0, points, zeta_n_dot)\n",
    "\n",
    "    delta_zeta_e = zeta_e_t - zeta_e_t0 - (t - t0) * zeta_e_dot_t0\n",
    "    delta_zeta_n = zeta_n_t - zeta_n_t0 - (t - t0) * zeta_n_dot_t0\n",
    "\n",
    "    u_per = u0 + piEN * delta_zeta_e - piEE * delta_zeta_n\n",
    "    u_par = (t - t0) / tE + piEE * delta_zeta_e + piEN * delta_zeta_n\n",
    "    u = jnp.sqrt(u_per ** 2 + u_par ** 2)\n",
    "\n",
    "    # Magnification\n",
    "    A_u = (u ** 2 + 2) / (u * jnp.sqrt(u ** 2 + 4))\n",
    "    A_u0 = (u0 ** 2 + 2) / (jnp.abs(u0) * jnp.sqrt(u0 ** 2 + 4))\n",
    "    A = (A_u - 1) / (A_u0 - 1)\n",
    "    F_pred = DeltaF * A + Fbase\n",
    "\n",
    "    ln_c = numpyro.sample(\"ln_c\", dist.Exponential(1 / 2.0))\n",
    "\n",
    "    return numpyro.sample(\"data_dist\", dist.Normal(F_pred, jnp.exp(ln_c) * Ferr), obs=F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 2500/2500 [00:11<00:00, 214.53it/s, 63 steps of size 6.80e-02. acc. prob=0.96]  \n"
     ]
    }
   ],
   "source": [
    "init_vals = {\n",
    "    \"ln_DeltaF\": np.log(250),\n",
    "    \"ln_Fbase\": np.log(188),\n",
    "    \"t0\": 3630.0,\n",
    "    \"ln_tE\": 3.1,\n",
    "    \"u0\": 0.5,\n",
    "}\n",
    "\n",
    "nuts_kernel = NUTS(\n",
    "    model,\n",
    "    dense_mass=True,\n",
    "    target_accept_prob=0.95,\n",
    "    init_strategy=init_to_value(values=init_vals),\n",
    ")\n",
    "mcmc = MCMC(nuts_kernel, num_warmup=500, num_samples=2000)\n",
    "rng_key = random.PRNGKey(0)\n",
    "mcmc.run(rng_key, t, F, Ferr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                 mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "  ln_DeltaF      5.58      0.00      5.58      5.57      5.58   2504.67      1.00\n",
      "   ln_Fbase      5.23      0.00      5.23      5.23      5.23   1794.19      1.00\n",
      "       ln_c      0.20      0.03      0.20      0.15      0.25   1565.42      1.00\n",
      "      ln_tE      4.53      0.03      4.53      4.49      4.58   1838.56      1.00\n",
      "       piEE      0.12      0.01      0.12      0.10      0.14   1680.86      1.00\n",
      "       piEN      0.30      0.04      0.30      0.22      0.37   1585.67      1.00\n",
      "         t0   3629.78      0.18   3629.78   3629.47   3630.06   1510.52      1.00\n",
      "         u0      0.49      0.02      0.49      0.45      0.53   1687.96      1.00\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "mcmc.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0, 0.5, '$\\\\pi_{EN}$'), Text(0.5, 0, '$u_0$')]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAFCCAYAAADvxZrJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZCk93nY9+/zvm+f09fc957YBbCLiyQM0rEVS+IhkLEJJpQdSrKhKFEhZAm2kpiRSLF0FcUUJaVMmzrIoGSJpCWFdklyyIigKJGJlETmBRLnAgvsYs+575m+j/d98kf3goPlHDu7s/P29Dyfqre236Pf98Fgpp/+3aKqGGOMMZtxwg7AGGNMe7NEYYwxZkuWKIwxxmzJEoUxxpgtWaIwxhizJUsUxhhjtmSJwhhj9ikR+aiIPCciz4jIX4rISOv420XkOyLyfOvfH173nje1jp8XkU+KiGz7HBtHYYwx+5OIZFR1rfX6XwCnVPX9IvIGYFZVp0TkHuArqjrauu5bwM8C3wCeBD6pql/e6jlWojDGmH3qWpJo6QK0dfxpVZ1qHT8DxEUkJiLDQEZVv67NUsLngPds9xxvl+NuG1GJaZyusMMwxuwDeZYXVLX/Vu7xIz/UpYtL/o7f953nql9R1Ydv9rki8jHgUWAV+KENLnkv8LSqVkVkFJhYd24CGN3uGR2bKOJ08WZ5a9hhGGP2ga/qn1y+1XssLvl86yuHdvw+d/jcXSLy1LpDT6jqE9d2ROSrwNAGb/2Iqn5BVT8CfEREPgw8DvzyuveeBn4deMe1QxvcZ9v2h45NFMYYs5cUCAhu5q0LqvrgpvdVfdsN3uePgS/RShQiMgb8R+BRVX21dc0EMLbuPWPAFNuwNgpjjNkViq/BjrdbISIn1u2+GzjbOp6jmTQ+rKp/+1qEqtNAXkTe0urt9Cjwhe2eYyUKY4zZBc0SxZ73Iv24iNwJBMBl4P2t448DdwC/KCK/2Dr2DlWdAz4AfAZIAF9ubVuyRGGMMbvkJquebpqqvneT478G/Nom554C7tnJcyxRGGPMLlAUv0PHpVmiMMaYXRJC1dOesERhjDG7QAHfEoUxxpitWInCGGPMphSsjcIYY8zW9rbP095piwF3IvKwiLzcmvb2Qxucz4rI/ykiz4rIGRH5qTDiNMaYzSiKfxPbfhB6iUJEXOB3gLfTHF7+bRH5oqq+uO6ynwFeVNV/JCL9wMsi8keqWgshZGOM+X4K/v743N+xdihRPAScV9ULrQ/+zwOPXHeNAunWkPMUsAQ09jZMY4w5mEIvUdCc4vbquv0J4M3XXfPbwBdpTl6VBv5r1e+fJEVEHgMeA4iTvC3BGmPMRppTeHSmdihR3Mi0tz8CPAOMAA8Avy0ime97k+oTqvqgqj4YIbb7kRpjzKYE/ya2/aAdEsUEML5uf6Npb38K+DNtOg9cBO7ao/iMMWZbCgS6820/aIdE8W3ghIgcFZEo8D6a1UzrXQHeCiAig8CdwIU9jdIYY7bRqSWK0NsoVLUhIo8DXwFc4PdV9YyIvL91/tPAR4HPiMjzNKuqfl5VF0IL2hhjrtOcwmN/fPDvVOiJAkBVnwSevO7Yp9e9nuJ7S/kZY0xbCtQShTHGmE1YicIYY8yWFMFvi2bf3WeJwhhjdolVPRljjNmUVT0ZY4zZhuCrVT0ZY4zZRHMKD0sUxhhjtmBVT8YYYzalalVPxhhjthFYicIYY8xmmr2erERhjDFmU51b9dSZ/1XGGGN2jZUojDFmF1j3WGOMMdvybQoPY4wxm7FJAY0xpo05iQSIoNUq6vuhxRF0aGO2JQpjzL7lptNIIg6uC7U6QaMRWqKw7rHGGNNGnGgUp7cHriUJ3wffD7U0oUjHtlF0ZvozxnQsN5vFGRtB+3IE2SSaiILjgB+Emiig2etpp9utEJGPishzIvKMiPyliIxcd/6QiBRE5IPrjr1JRJ4XkfMi8kkR2Ta7WaIwxuwb3uAAjA/h92doZOMEiSjqOaBKUC6HGpsq+OrseLtFv6mq96nqA8CfA7903flPAF++7tingMeAE63t4e0eYlVPxph9wevtRfty+JkEfswFB6Te+k5erxOUKyFHKHs+15Oqrq3b7aLZVNKMRuQ9wAWguO7YMJBR1a+39j8HvIfvTyav0xYlChF5WERebhWFPrTJNT/YKl6dEZG/2esYjTHh8YaHCA4PURtKU+uO0ki5+NHmx5fUfbRQBA1CjVEJpUSBiHxMRK4CP0GrRCEiXcDPA7963eWjwMS6/YnWsS2FnihExAV+B3gncAr4MRE5dd01OeB3gXer6mngH+95oMaYUHhjozQOD1IZTFLNRailnGaScASnoUixgr+0HHaYQLPX0043oE9Enlq3Pbb+niLyVRF5YYPtEQBV/YiqjgN/BDzeetuvAp9Q1cJ1IW5U5NENjr1OO1Q9PQScV9ULACLyeeAR4MV11/w48GeqegVAVef2PEpjzJ7zjhymPt5DpTdKvcsh8EAUJFCkoTjlGrqyGnaYQLPXU3BzvZ4WVPXBTe+r+rYbvM8fA18Cfhl4M/CjIvIbQA4IRKQC/Ckwtu49Y8DUdjduh0QxClxdtz9B8z9yvZNARET+GkgD/0ZVP3f9jVqZ+DGAOMnbEqwxZm94x45QO9RDuT9KLS0EkeaHsFtVJAC34uOslmi0SWkC9n4chYicUNVzrd13A2cBVPUH1l3zK0BBVX+7tZ8XkbcA3wQeBX5ru+e0Q6K4kaKQB7wJeCuQAL4uIt9Q1Vde9ybVJ4AnADLSs21xyhjTnrzjR6kc7aU8EKGaEfxY87hTBy8AtxrgrVUIZtqnckEJZWT2x0XkTiAALgPvv4H3fAD4DM3P0i+zTUM2tEeimADG1+1vVBSaoFk8KwJFEfl/gPuBVzDGdBTv2BEqx3opDkeoZptJQt1mkpAauDUlkm/gzK/SKJXCDncd2fM1s1X1vTdwza9ct/8UcM9OnhN6YzbwbeCEiBwVkSjwPuCL113zBeAHRMQTkSTNqqmX9jhOY8xt5h05TPmOPgqjEco9Qi0DjSQEra+0Th0ixYDIYpHGxGS4wV7nWolip9t+EHqJQlUbIvI48BXABX5fVc+IyPtb5z+tqi+JyF8Az9EsYv2eqr4QXtTGmN3mjY1SOdFPYSxCpUeop0A9QMENQHyIlJTocg2m2qfK6SAIPVEAqOqTwJPXHfv0dfu/CfzmXsZljNkb3kA/lbuGyY9HKPcL9TT4sWYzo1sVpAFeGWKrPpGZtbZqwF5vr6ue9kpbJApjzMHlZrPUTo2TPxShNCjUsuAnAhCQhkAAbgViawGx2TL+xcthh7whVdk3VUk7ZYnCGBMacV3800dZOxKjNChUexQ/GYCn0BCcmuDWIFKE+GID9+osjZAn/tvKboy0bkeWKIwxoZEH7mbleILisFDpU4KUD9EAAnktUXgFSCwFxKbWaMy2b9tEc81sq3oyxphd454+ycrJFMURodKv+JkGTsxHRAlqLtIQvBLE1pTEbAU9dynskLchVqIwxpjd4h0aI3+ym8KoQ3lACbrrRBJ1HEdp1F204eCVHCJ5SM75RC7O0ajVwg57S83usVaiMMaYW+Zms5ROD7N22KU8pAS9dWKpKp4bEARCverhVBwiBUgsKMkreRqT205H1BZsKVRjjNkFemSU0qBHtRv8tE80WSMerSMClWoE9R3cqjQbsFd8ZGYh7JBvyC1MCtj2LFEYY/aMd3iclbsyFMaE6oBPtKdCOlkh6vrUfBffd6DkElkTEvNK4tIKjbn5sMO+Ybe6tGm7skRhjNkTbipF8d5h1g45lId8or1lelIlEpEavjqUahHqFQ8v7xBfhNTVCv6Z/TOdW3MpVCtRGGPMTXGzWRr3HWP1mEd5NMDrq9CbKZKLl3FQVmtxKrUIFDyiy0Jq2id6bppG2IHvkFU9GWPMTaq/4Tgrx+MURxX6q/TnCgwkC8TdBoV6lGrdo1qMEllxSc4pXRfyNKZnwg57R5ptFFb1ZIwxO+becyeLJ+IUR6ExWGOgu8BIao2UV6WhDrXAJV+OI6sR4guQvlIneO5s2GHfFJvryRwIbi6Hv7ISdhimQ3j9fayc7qYwDtWhBr19BcYzK/TH8gAs17rIV+NUCjFiSw6pqYD42WkaGoQc+c7ZOApzYEg2DZYozC6p3H+EtSMO1ZE6qd4Sh7PLjMRXibkNCo0oa/UYy8UEznKE5CykX22/dSZunFU9mQNCM11hh2A6hHvvXSyciFIaDcgN5BnJrDGWWCEXKVFXl2IjxnI5SWUlTteckLlUh2f2Z5XTNZ0611Nnpj9zU9xsFj8VCzsM0wHcbJbV0znyhyAyXOJIbpmjqUWGYqsk3SrVwGOllmApnySyECE1qSRfnkcb9bBDNxuwEoV5jfT34ic93LADMfte/Q3HWT3q0BiucqJ3iWOpBQaiaySdGgU/Tr4RZ7aYprYUJz0D2fNlGhcuhR32LbFxFOZA8HtT+HHXEoW5Jd4dx5i9M05pzGdgYI3j6QUOxZZIu2Vq6lHyo8yW0yyupIjNemQv+zjPvML+a77+ftZGYTpeLRvFj3bmNyKzd/L39pM/DImhInd2z3FHco4+b615zk8wX0szmc8SzMfomoTU2WX8UinkqG+dzfVkDoRaxkOtOGFugTx4LyvHXeojVU72LnKia46RyApxqVEMYqz6Ca6WciwtpkjOOGQv1PBf2j/TdGynUxuzLVEYAJxkknqqueykMTfDzeVYPJ2iOB4wNLjK3dlZDscWyLlFikGMfJBgqpJjcjWLOxclfUWJP3tp303TsZlOHkfRFhVqIvKwiLwsIudF5ENbXPd3RMQXkR/dy/gOAmdkiHpSCKxEYW5S7U3HyR8RoqNF7u6Z4UR8hiFvhQg+pSDGQj3NxUIPawspkpNC7myexvz+mEL8RgXq7HjbD0KPUkRc4HeAdwKngB8TkVObXPfrwFf2NsKDodGfxo+HHYXZr7xDYyyfjFEeq3Osb5G7u2YYjy6SdirUcZmtZ7lS6WFiOUd0KkLuQoPgOy+EHfbu0mYbxU63/SD0RAE8BJxX1QuqWgM+DzyywXX/HPhToH1XV9/H6ukIfjTsKMx+VTvST7kf4t0VxpMrDEeWSUsFgGIQY6GeYracprwWJ7YIyYliyBHvPqXZRrHTbT9oh0QxClxdtz/ROvYaERkF/kvg03sY14FS7fYIoiAadiRmv/EG+lm+K0F5vM4dfQucTk1yKLJExqlSV4/5RoaLxT4uLPUSmYo2SxPfPRN22LdFp5Yo2qExe6Of1PUfV/8a+HlV9UU2/8GKyGPAYwBxkrsW4EFQzTgEEZBO6Mxu9lT5gcPkD0P3UJ57c1OciM3Q45YIVFgJklyq9PHqSi/FuS66r0LqudmOacBer5Mbs9shUUwA4+v2x4DrV1J/EPh8K0n0Ae8SkYaq/h/rL1LVJ4AnADLSY9+Nb5Q41FMQeODU7cdmbpw3NsrcySj1sSqn+2a4OzHFkLtGXAKWghiT9W5ezg+yNJcmPunR/XKVxqXLYYd921iiuH2+DZwQkaPAJPA+4MfXX6CqR6+9FpHPAH9+fZIwN88bH6GeAnXVEoXZkeIDoxQOK6NDy9yTnuRYZI4ep05VhTk/zbnyIBeWeolORMlcUqJPvYIfdtC3iQ24u41UtSEij9PszeQCv6+qZ0Tk/a3z1i5xm/nDPfgJwAG3anVP5sZ4Rw4zc4cHo2VOd89wKj7FkFciKsJS4HG13svZtUHyMyl6L0P3mTx+Ph922LfVfmmc3qnQEwWAqj4JPHndsQ0ThKr+N3sR00FS6Y/jJxTxwSt36vc9s9sK9w1SPKwcHljkvtRVjniLpAUqqsz4GV4uDXFhoZfEhEfPS6XO6w57Pe3cqqd26PVkQlbpcQliASi4JZvm2WzPO3KYlTs83JES93VPcWdsmn7XJyIOq4HHheoAL64OUZvqIntBkW8+H3bIt921xuy97PUkIh8VkedE5BkR+UsRGVl37j4R+bqInBGR50Uk3jr+ptb+eRH5pGzVQ6jFEoWhmhU0GiABuGuVsMMx+0DhvkGKhwKODSxwb9cEh7xVkuJS1YBJP8uLpREuzvXSdcUh9/wy6h+MkmoI3WN/U1XvU9UHgD8HfglARDzgD4H3q+pp4AeBa98CP0Wzd+iJ1vbwdg+xRGGop0G8AGkIrHZ2HbK5dd7YKCvHPbzREvfnJrkzOk2f6+AgLAYO56pDvLg8RDCVJHfex39+f69a185UdW3dbhffG1rwDuA5VX22dd1ia3jBMJBR1a+rqgKfA96z3XPaoo3ChMcbbfZ4EldxGhAsLoUdkmlzxQdGKR4OuLN/gXuTVxn3iiQlzmpQZaqR40xhlCuzPaQuC+lnZjpyzMRGbqHXU5+IPLVu/4lWV/8bIiIfAx4FVoEfah0+CaiIfAXoBz6vqr9BczDzxLq3f98A541Yojjg/NFeGl0B4ihuFYKaTR9rNucN9DN7h4c3UuD+3CQno7P0OBEAlgLhbHWEM8tDOJNxes7WO3rMxEb05hLFgqo+uNlJEfkqMLTBqY+o6hdU9SPAR0Tkw8DjwC/T/Gz/+8DfAUrA10TkO8DaBvfZtk+8JYoDrjKQIEgEeAKuNU+YbZQfOExhXLlzYIHTiQlGvAoxSVDQGlcbWZ4vjDE50033JSHx3YsHpjRxze3oHquqb7vBS/8Y+BLNRDEB/I2qLgCIyJPAG2m2W4yte89GA5y/j7VRHHDlPhcn3vxztkRhtuIkk6yciOKMlrg3O8VdsRl6nOZMkgt+wNnqCC8sDhO9GqPnxUrHTSG+HdVQej2dWLf7buBag9BXgPtEJNlq2P4HwIuqOg3kReQtrd5OjwJf2O45VqI44Ko5wYv5NGoukaKNyjab8994ksIh5ejAIvcmrzLq1khIirWgzKVGjmfy48zM5Oi9BN63XuqINbB36iarnm7Fx0XkTiAALgPXBiovi8i/ojnzhQJPquqXWu/5APAZIAF8ubVtyRLFAVfLguf51Mse0fxB/NM2N2r5ziT+SIV7ctOcjM7S7cQICJgP/GbbxOIwsStRes4UCcrlsMMNwd5P4aGq793i3B/SrGq6/vhTwD07eY4ligPM7emmEYeEG4Dv4JUsUZiNeYfHKfdDNlfiUGyJfrdCRLqoq8+8H+dSpY/FlRSpefDOTx64tolrQihR7AlLFAeYjg/TSAW4bgA1h9hyefvuD+ZAKt47TGnM5/7eWU7FJ+hzIjg4LAclztWO8NzyCEzE6TlbO3BtE9fYNOOmI1WGu6CrgesESE1wF/IH9pug2ZyTTLJy3CMxvMZ9mQmOeCsknSRVrXO1EeOZ4iEuzvSRvSjEvvlyx84Ouy1tNmh3IksUB1h5wMNNNOuS3bJDMD0bckSmHflvPElxTDnVt8Cp+BSDrouDw2pQ4pXaIZ5bHsGZiNPzUqXjZ4fdjs0eazpOpUeIxer4gYNXgqBUCjsk04ZWTiTQkQqns9Mc8RZJOTHq2uCqH+XZ0uHXShPeNw9mT6drFGujMB2omoOuaIN6wyVysL8Imk14J46TPyyMDSxzOjHBoOvj4LAQlDhXG+fZdaWJg9nTaT1buMh0GO/QGPWMEvMaFCtRugphR2TaUf50H5WRBqdyMxyPzpFtlSamfI/nSuNcmO0jc0nwvvHigS5NXGNtFKaj1A/14ad8Iq5Po+YRX7Y/c/N6TiLB6hGX7OAyd3dNM+KWiUgXy0GZV+tDPLsyChMJel6sEFRtWD9Y1ZPpMKWhOG5XFVcC/IpLfNEWLDKvFzxwkuK48obeBe6KTdHjRPA1YMYXnisd4vxsP+lLcmBHYV9P1RKF6TDlfodYvN78xS67ROeWD263RrOh1RNJGKpwKj3NiLdGTKIUtMaFeh/PrY7iTybpOVu1tol1rI3CdJRKDyTjNRrq4BYdmLCuseZ7vNER8oeEkf4VTsan6XcCApQFP+CF8hgvzw2QviREv/WKfcFYp1PbKGz22APIG+inllUSkTq1hkekIPirq2GHZdpI5dQI5WGfu3JzHIkskBSPqja41Mjx3NoY1ckuul+uH/hxEwdFWyQKEXlYRF5uLfb9oQ3O/0RrAfHnROQ/icj9YcTZKfyjw/gZn4RXp9pwiW60lIk50FaPRYkPlrira5ohr4Qrwnzg82JljJcWBui64pD47sWww2w7qrLjbT8IvepJRFzgd4C301xs49si8kVVfXHdZReBf9CaOvedwBPAm/c+2s5QGkngpGp4TkClFiGz3KHlZXNTvDuOURiFo71LHIvOkXOEugZcbaR5rjDK2lSG8fP+gZ3TaTPK/vng36l2KFE8BJxX1QuqWgM+Dzyy/gJV/U+qutza/QavX6HJ7FBx0CWerOGg1EoRkgtWy2y+p3h3H9XhOifT84x6K0RwWAkanKsOcWZxmORVl9TTk2GH2Zb0Jrb9oB0SxShwdd3+dot9/3fcwEIbZnOVXkjFawQIWvKIz9jUHeZ7Vg97ZPqL3Jmcpset4qNM+QleKI4xN5sleyGgcWUi7DDbj1rV0+200U9qw0QrIj9EM1H8/U3OPwY8BhAnuVvxdRRvdIRat9IVrVLzXdyCizMxZ/3gDQDeXScojcBd3Uscji7QJVBSn0v1fl5YGSJ6NUr2+QXr6bSZ/VJE2KF2SBQTwPi6/Q0X+xaR+4DfA96pqosb3UhVn6DZfkFGejr0f9mtCQa78RMBca9BvhbDLQuN2bmwwzJtojqaoZ7zGU6skXNKOMBK4DBR62F2LU1iAYKXrRF7M/ulhLBT7VD19G3ghIgcFZEo8D7gi+svEJFDwJ8B/0xVXwkhxo5RGk8hmTpRx6dcjxC1XrFmndWjMWIDJe5IztHrlgmAyUaGM4URStMpcq/U0YaN4t+M6s63/SD0EoWqNkTkceArgAv8vqqeEZFri4R/GvgloBf4XREBaKjqg2HFvJ8Vhl0SXQU8J6BUiZJc3v495mDwBgcojsBY9wqHowvEJSAfCBdqA7y0NEjyqkvy6Yu2uNUmbJrx20xVnwSevO7Yp9e9/mngp/c6rk5U6YNsotp8XYrSN29/9qapdvcYlUGf4+lFhrxVHGDOT/JiaYS5uSzDFwOrptyKApYozH7nHTtCtVtJxyrUfBcteiSmbJ1s07R6NEa0v8jR5Dw5p0xVhauNXl5cHSI6ESX7gjVib2e/VCXtVDu0UZg9UjnWR5BrkIrUqPge7pqLe2Um7LBMG3CzWYojwkj3KmORJSISsKZRzlcGubzQQ+oK+GeseXBbHTqQwkoUB0hhNEI0XSTq+BRqMaKrQmNuPuywTBvw7z5CZSDgcHqZAa85p8tMI8NLhSHq00m6X7b1Jra3f8ZF7JQligOkNCCkuyp4jk+hEiO+FHZEpl3kjyVx+iscTS6QdipU1OVyrY/zy30kJxwi3z1n1U43Yp+UEHbKEsUB4Y0MU+1ThuJlAhVKpRhDczbMzjQVRoT+7jzDkRXi0mAlSHCuPMj8fJaRi4HNEnsjOnjhImujOCBqJ4apd/tkolVqgUejECE5ZQvOGHBPn6Q8qIynV+j31vBxmKx380p+gMhklMwZK3resA5to7BEcUDkD8WIZqskvRqFehRv1cW7bF0dDRSP5/D76xzuWiLtVCgFUa7U+ri02ENqAvyXrBH7xslNbO3PEsUBURoUMqkyUafRbMheERpT02GHZdpAfswj011kKLpCVHxWgiTnS/2UZ7vIna+FHZ5pA5YoDgDv0BiVfqUn0ZwldqWUIL7hbFnmoPEG+ikPwnBmjX6v2Q4xU89ybrWf+LRL/LkrIUe4z3Ro1ZM1Zh8AlZND1Hsb5KIVGoFLuRile8b6sBionxyl2uczklwl7ZQpapSL1X6mF7P0XFEbib1T++SDf6csURwA+UMRYtkiqUiFtXqCYC1KcqLUqb/TZgcKhxK4PWVGYqvEnTorfhcXi30Ec3Fyr9g6JTvSwVN4WNXTAVAaEnpSJWJOg7VajMiKi3PBFp4xUBwSerNFBqJr+OowW89ycbWH5KSD8+y5sMPbd/Z69lgR+aiIPCciz4jIX4rISOt4REQ+KyLPi8hLIvLhde95U+v4eRH5pLRmWt2KJYoO5x0ao5aBdKyCI0qpHiVSAH/Jpo096Lz+PmpZ6E6U6HKqBDgs1FOsFJptWEHJShQ7tvdtFL+pqvep6gPAn9OcaRvgHwMxVb0XeBPw34vIkda5T9Fc4O1Ea3t4u4dYouhwlTuHqPfV6YsXCVRYLiaIL4QdlWkHjeOjVHt9hhJ5kk6NNT/OlXIPtbkEuVdtyo6borLz7VYep7q2breL76UeBbpExAMSQA1YE5FhIKOqX1dVBT4HvGe751ii6HBrh6PEcxWykTJlP0q5ECM1bVOLGygcTuB01xiMreESsOSnuJjvIT7rEj1jvZ1uhujOt1t+psjHROQq8BN8r0TxJ0ARmAauAP+rqi4BozRXFb1monVsS5YoOpgTi1Mchr50kYRbp9CIwUqU5GWbjsFAccghly3RHSk2R2NXu5ldzpCaUBrzVuzcsZupdmomij4ReWrd9tj624rIV0XkhQ22RwBU9SOqOg78EfB4620PAT4wAhwF/qWIHGPjEX7bpivr9dTB9N4TVPsDBpIFXFGWKgmiyw567lLYoZmQeYMDVPrgaFfhtdHYk5Uc9fk42QvVsMPbp266KmlhqxU7VfVtN3ifPwa+BPwy8OPAX6hqHZgTkb8FHgT+X2Bs3XvGgKntbmwlig62dqIL6a3SHyvgq7BUShJfgKBsczwddI1jw9R6fAbihVa32CSX893NaqcXrdrppu1xY7aInFi3+27gbOv1FeCHpakLeAtwVlWngbyIvKXV2+lR4AvbPcdKFB2sMOrQnSuSiVRYq8fJ5xMMTduMsQaKYwmcXJW+WAGAuVqGuZU06Sms2ulW7P3gpI+LyJ1AAFwG3t86/jvAHwAv0Kxu+gNVfa517gPAZ2g2cn+5tW3phhOFiPyH1svLwHeB76rqyzf6frO33LtPUh5SDnfliUmd1Xo3wUqU1JWiDbQzFIccMpkyWa9MJYgwWclRW0yQvWTVTrdkj/+4VPW9mxwv0Owiu9G5p4B7dvKcG04Uqmr3CFMAACAASURBVPpPAETkOPA/Af8bkNnJw8zeKdzZjd9XYzi5hiPKUjVBdNHFeeWqLUBzwLm5HJU+GEsWSTo18kGcyUKW2JxL7KUrWJ+4m2Qjs0FE3ioinwB+Dvg6ML5bQYjIwyLycmuk4Ic2OC+tEYTnW6MQ37hbz+5Ua0dcsj1FeiJFqoHHfDFFYgH8lZWwQzMh02Oj1HIBffESEfFZrncxu5YmOYPNKHyLwugeuxd20pj9B0AS+BvgW6q6uhsBiIhLsz7tncAp4MdE5NR1l72T740ifIzmyEKzCe+OY5SGYCSzRsqtUvBjrK0l6Jqy9gkDpfEUmq3THS3iI8xWM5SWE2SuWFnilnXo7LE7SRT/FfBRmoM4/qmI/O+7FMNDwHlVvaCqNeDzwCPXXfMI8Dlt+gaQa40wNBsonOqjMVhjKLFGxGmwXEsSLEdJXSqGHZppA8VBl0S6SsqrUQ0izJTTRBYiJF+xuefNxnaSKP5bVZ0A7lfVXwL+/S7FMApcXbe/0UjBG7nGtKwd8cj0FBmM5QnUYa6cIrbg4py9GHZopg1UeiHXVSbh1FhtJJjJp0nMQuPcq2GHZtrUzYyjGGr9+45diuFGRgre0GhCEXns2ujGOgez94Z34jjFERjNrpL1SpT8KHP5NIk58AuFsMMzIfOOHaHarfTES7gSsFTvYm01ScqqJXeFtVHAQyLyu8DdInI/uzcGY4LXN4xvNFLwRq5BVZ9Q1QdV9cEIsV0Kb3+5Vu00mlwlLg0W613kVxKkJ6z+2UBtvAc/65OLlfHVYb6SgqUo6Us2U+yu2ONJAffKth/2IpJS1YKqPigiYzSnrH0vu1f1823ghIgcBSaB99Ecfr7eF4HHReTzwJuB1dYIQ3Od1aMu2d5VhmPNvgbz1TTOUpSuVxetW6yhNBzFTVfpcmuUgyhzpRSxBQfn7GX7/bhV+6hxeqdupFTwdGvyqVdabRQTIvJSq53ilqlqQ0QeB74CuMDvq+oZEXl/6/yngSeBdwHngRLwU7vx7E7j3nMnxVHlVHaFHq9IRT1mimnic4J/5pWwwzNtoNTnkExWiTkNCo0oS4UkXbPgr+5KJ0ZzgBNFD83uq3eJyALwEvBG4I7dCkJVn6SZDNYf+/S61wr8zG49r1OVxzM0cg16YiUi0mCp0cVKKUHMhk6YlkYXdEXrOBJQrCeoFqP0Lln7xG7ZL20OO3UjieKKqv4QQKvq6W7AZpVrQyvHI6T6VhmLLxMRn/laiuJKgr6r1j5hwBsbpZqDoVhzUaKlahJZiZC6ZNO67JoO/UHeSKLIish/Bpy5VvV0m2MyN8G9+yTFMTjRvUxfpEBdXabLWdz5CMlX5mxaBkNjvI9GxicTrVJXl8VyV7N94uVL1j6xWw5wokgBHwROi0iMZtXTC6r6P9/WyMyO5E/10BiqcbhribRbZr6RZiqfITEr1j/eAFAejENXg7hbp+xHWC4mSM6Dn7eFrHbDfuruulM3kij+rqq+CiAicZrTbJy+rVGZHXGiUVaPOvT2rzHeqnaaq2VYXEkxOGn1z6ap3OcS7SoTdXzW6gmqq3EGpqysuav2SXfXndo2UVxLEq3XFVpTjN/OoMzOBA+eojQWcGd2mT4vT11dpso5dCFG+rzVP5umag6S8RoAK9U4zopH16UVq3baTR36x2YLF3WA5TsTOEMljnYtknSqzDayTBSyJGYd5Plznfq7a3bA6+2lloFstEagwkolQWxJbFncXXaQq55MG/MOj5M/LIz0rjIaXcYVfW21stykElQrYYdo2kAwNkgjExD36lT8CPlinNS8LYu76yxRmHZUuH+YykiDo+klerxmb6fJSo76fJzsBUsSpqk6mMRPBsTdBqVGhOpajMEZa5/YVR3cmH0zkwKaNiGuy8oxj8xAgaPJBZJOlRU/ydVCjvisi3fmUtghmjZR6YvgJBo4ouTrMdxVj+RVmyRy19l6FKbd6EP3UDykHOleYjjSHH49XcsxvZQlNan4S8shR2jaRSUnROLNEsRaJU5sWZBzV0KOqgN1aKKwqqd9bPnuJDJU4Y70PDm3RCmIcbXc3ax2evVgTrNuNlZPQyzSoKEOxVKM1KJNO387WNWTaSveiePkDwvDfSscii0Rd+osNVJcWushMeMSfeFy2CGaNuEN9FPvgmikQaXhUS9ESc7Z+Bpz4yxR7FNr9/VTGa1zIjtPn7eGr8JkrZuZpQypCaWxaMtamiYd6qPRpUScgHI9ipN3SU5abydz4yxR7ENeby8rxxxyA3mOJhbJuBXyQYLLxR6CuTjZ8/YhYL6n1p8kiAeIKKVahMiqg3d5NuywOpO1UZh2UX7TUUrjAQ90LzESXcYhYLae5dJqN4lpB+/5izba1rymmotA3EdVKFeiJJahMT0Tdlidx7rHmnYhrsvyySiR4RInUnP0uAUqGuFSpZf5hQyZK2qL0JjXqWYFifj4gVAvRUgsdOinWTuwEoVpB87J41T6oDddpC9SICI+K40u5itpZDVCcroWdoimzTTigrhKPXCh4hJbtfLmbbNPPvh3yhLFPrN6Xw+V0Tp3ZBcYjDRLDrONzGvVTrEXLtvaE+Z1Gklw3IBa3cXNOySm8lifp90nWNWTaQPeoTFWjzrkBvMc7/re2InL5V7m5zOkryiN+YWwwzRtxO3pbiYKR6nXPaJrgnPF2iduG6t6MmErPDBCacznDa2R2A4B8400r+b78Gaj5M7aN0XzetLfSxBtvq6XPZIr2JeJ26WDG7MtUewTXn8fK8c94kN5jqcW6PGajdhXqr1cWeqma0IIvvti2GGaNtPo7sKPKRIIVFwSSx36SdYuOvTHa1VP+0T5jUcpHlIO9y4xHlskKj6rfpJXi/1UZ5PkLjRArTxhXq+eiaARJfAdnKJLYq4edkidrUOrnkJNFCLSIyJ/JSLnWv92b3DNuIj83yLykoicEZGfDSPWMLnZLMsnI8hQmePpBXrcAr4KE7UeLq70kJhySb4wHXaYpg3VUy7qQFB3iBSE2ORa2CF1tGvrZu9k2w/CLlF8CPiaqp4Avtbav14D+JeqejfwFuBnROTUHsYYutob76Awroz2rTAeX6LLqbLid3Gx1MfSfJrMZaVx+WrYYZo21Eg6qKNo3SGyBnp5IuyQOltIJQoR+aCIqIj0rTv2YRE5LyIvi8iPrDv+JhF5vnXukyKy7ULfYSeKR4DPtl5/FnjP9Reo6rSqfrf1Og+8BIzuWYQhc+4/xfJdMYLhKsczC/R7eXwcZhsZzq/2EZmKkj1rs4CajTXizc8AqTrEViEolUKOqIPdTJLYhUQhIuPA24Er646dAt4HnAYeBn5XRNzW6U8BjwEnWtvD2z0j7EQxqKrT0EwIwMBWF4vIEeANwDc3Of+YiDwlIk/V6YxptlfuyVIYh8G+NQ4nlkg7FfJ+goulfqYXsqQmQJ96PuwwTZvyY4AKbsUhvmQD7W63kKqePgH8HK9PO48An1fVqqpeBM4DD4nIMJBR1a+rqgKfY4Mv6Ne77b2eROSrwNAGpz6yw/ukgD8F/gdV3bCiVVWfAJ4AyEjPPqn925x7392sHRHqIzWOZhfpW1eaeDXfC7NxcudtJLbZmBONEkRAAvCKQmLeflduuz3+1BGRdwOTqvrsdTVIo8A31u1PtI7VW6+vP76l254oVPVtm50TkVkRGVbV6Vamm9vkugjNJPFHqvpntynUtrNyOkt51Ke3f43xxDJpt0IpiHK53MvVxW66rgrxZ20kttmYk8uiDogveEWITK/a78ptdpMlhD4ReWrd/hOtL73Ne279ZfsXgHdsFMoGx3SL41sKexzFF4GfBD7e+vcL11/Qamj5t8BLqvqv9ja88DRLEw7xgTyHMisMRJuFqPlGmouFXuozSQYu+jTm5kOO1LStVBfqgFMTIgUIrlhDdptaUNUHNzu52ZdtEbkXOApcK02MAd8VkYdolhTG110+Bky1jo9tcHxLYbdRfBx4u4ico9kY83EAERkRkSdb1/w94J8BPywiz7S2d4UT7t5ZOZ2lNOoz2r3CoWSzbaLYmq7j8mIPyUmH9AuWJMzmNBlDBdwqxFaVoGZVT7fdHjZmq+rzqjqgqkdU9QjNJPBGVZ2h+SX8fSISE5GjNButv9VqC86LyFtaX8IfZYMv6NcLtUShqovAWzc4PgW8q/X6/2Pj4lLHcu4/xdphh+hAgUOpZQaieQAW6s3SRGU2yfClgMb5CyFHatpZEG/O3eFWIb5klU63XRsNoFPVMyLyH4AXaQ4x+BlVvdab4QPAZ4AE8OXWtqWwq57MBlZPZymPBIx3rzEYy5N0apSCGBOVbi4vNgfYZV9cssWJzJY06iIBOGWILVba5TOsYwnhfqNtlSrW738M+NgG1z0F3LOTe4dd9WSu4x07Qqlf0Gyd3niJpFPDkYC8H2e2kqaSj5FYAD1/OexQTZsLPEECcGvgLtv4iT3RoVN4WImizRTuHaQ4qvT0FhhOrBJ36lSDCFPVLJdWuolMRcm9UiUo27rYZmvqOjgNiBSUYNKmeNkL+2VKjp2yRNFGvDvvYPWoiz9QZTi1Rsar4COsNpJcLvawOp+m+ypEn37Vqp3MttQVnDrE1gIbkb1XLFGY223tdC/FESXbU2QgkSfu1FltJJmuZLmy3E10KkL3uSr+ykrYoZp9wq0qsWWbMXbPWKIwt5N7712sHXbxB2oMpApkvAr1wGW2muZKIUdhvoveCYg9c8kGTZkbIqp4FSWyVLYS6F7YR7PB7pQlijaxdleO0rCSzpXoixeJOAFrjQQTxRxTK1ni0x65V6o0FhfDDtXsE+IrkWKAzC+HHcrBYYnC3C7y4L2sHXLw+2v0pQqkvSr1wGGmkmFyJUt1NknvFWubMDsjjYBI0acxbWtk7xUrUZjbwk2lWDmZojykJLNlcrEKjgSs1JNMFTMUF7pITrt0v1K2tgmzI07Nx6n4nfoltz116A/bEkXI/PvuID8q+L01urvKJL0aZT/KTDnN7HKa+LRH+oriPX3OShNmR5xKAylWrE1rD1mJwuw6b3CAxTsSVAaVeKZKJlYBYKWWYCafxp9LkL0KubMF/IItTmR2RkpVdNHaJ/bMPhpAt1OWKEJUPTVOYUTwe2r0JitEHZ9SI8pcKcXqUorktEPulTL67TNhh2r2o5JVV+45SxRmN3nHjrB8NEq1T4mmayQjNQIVVqtxFte6cOcipK8okafP42sQdrhmH9K1fNghHCiCVT2ZXVY62U9pSGjk6qTjNVwnoFCPsljsoraQIDMl5M6u4eftj93cHKuuDIElCrNb3HvuJH8oQqVXiaRqxL0G9cClUI2RX0kQm3PJXvQJvvNC2KEaY3ZAtDMzhc0eu8ecRILi8RzlfvAzDaLRZp+UYi3K6loSdyFKagJSZ2xRImNMe7ASxR7T03dQGHGpdStuooHnBlQaHqVyFH85SmpayJ2v2KJExuw31uvJ7AZvcIDVI12U+8BP+USjPn7gUK07VFdjxOddMlcCIs9dtDETxuxD1phtblkw2k+l28FPANEABHzfoVaO4BY8YsuQnKniL1nfd2P2JUsU5lZ4/X3kD6Wo9kAjFSBegAZCre6iqxESc0Lmso/3zHkrTRizT1mJwtySxh2jlAZdahnQWIAINKouWowQW3TpmlJSLy9bd1hj9rMOTRSh9noSkR4R+SsROdf6t3uLa10ReVpE/nwvY9wN3pHDFMcSVHPgJxV1FfUdtOThrbgkZiF7oYL/0ithh2qMuVmt9Sh2uu0HYXeP/RDwNVU9AXyttb+ZnwVe2pOodpG4LrUjvZT7HOpdEHja/NZRdnHXXBLzkLnawHv+YtihGmNuld7Etg+EnSgeAT7bev1Z4D0bXSQiY8B/AfzeHsW1a5yTxykNRamlwY83v0JI1cEtOMSWhK7pgK6zizYnjzH73LUpPDqxRBF2G8Wgqk4DqOq0iAxsct2/Bn4OSO9ZZLvAzeWojGeodDs0kqAuiC84NSG6IiRnlcyrRRrnXg07VGPMbujQkdm3PVGIyFeBoQ1OfeQG3/8PgTlV/Y6I/OA21z4GPAYQJ7nDSHefHh2h0hehngT1QAJwyg5eERILkLlSQ54/v19Kn8aYbeyXEsJO3fZEoapv2+yciMyKyHCrNDEMzG1w2d8D3i0i7wLiQEZE/lBV/+kGz3oCeAIgIz2h/i/zRkcoDXVRzQhBrHnMqQluBeJLkJpsEDs7TaNcDjNMY8xu2UdtDjsVdhvFF4GfbL3+SeAL11+gqh9W1TFVPQK8D/i/NkoSbUUc/LE+qt0efhzUAWmAV4LYMnRN+3SdX6YxORV2pMaYXSTBzrf9IOxE8XHg7SJyDnh7ax8RGRGRJ0ON7BZ4xw5T7YtTTwqB2/xlcCsQXYPkfEDqUoHG2XNhh2mM2W0d2usp1MZsVV0E3rrB8SngXRsc/2vgr297YLfATaVoDGSopR38aPOYUwWvAvGlgNSVMpyxxmtjOpG1UZgbImPD1HJRGjEBAbcOTh1iq0rXZBXv/CSNaiXsMI0xu02xXk9me15/H42eLhpJt9UVFpwaRAsByZk60QtzNOYXwg7TGHObWInCbEv7eqinPPxos13CqyheWUks1IldXqQxMRl2iMYYs2OWKHaJNzxEI5sgiDmIglsDtxIQW24Qm1ihceFS2CEaY263Di1RhN3rqXN0JQliLirNrrBuJSBSDIislNHp2bCjM8bcZmFO4SEiHxQRFZG+1v7bReQ7IvJ8698fXnftm1rHz4vIJ0VEtru/JYpd4PX2EqTjBBGnWeVUDYitNojNFJELk/iFQtghGmNuN9Wb226RiIzTHF5wZd3hBeAfqeq9NMeo/bt15z5FcwaLE63t4e2eYYliN2TTaNRDVHGrAZHVOtG5InJ5Cn91NezojDF7JKQSxSdozoX32t1U9enWMAOAM0BcRGKtGTAyqvp1VVXgc2wyGet61kZxi9yebjQRRUVwagFSD3DzZZicsxlhjTlobu6Dv09Enlq3/0RrOqJtici7gUlVfXaLGqT3Ak+ralVERoGJdecmgNHtnmOJ4hZJMtlMEnUfqnWcQhWdW7CShDEH0E2WEBZU9cFN77n1xKq/ALxji/eeBn593TUbZZNto7ZEcQvcdBo8D3xFSjWkVCaYXyQolcIOzRiz1xQIdr/b02YTq4rIvcBR4FppYgz4rog8pKozrXV8/iPwqKpemw5ionXdNWPAtpPOWaK4BRKNQuAjpTLU6gSLSwS1WthhGWPCsofdY1X1eeC1NXxE5BLwoKouiEgO+BLwYVX923XvmRaRvIi8Bfgm8CjwW9s9yxqzb5K4LqoBVKroap7G9IwlCWMOuDZa4e5x4A7gF0XkmdZ2Lal8gOZqoeeBV4Evb3czK1HcAq1U0UbREoQxpinEuZ5aSzFce/1rwK9tct1TwD07ubclipukvo9aW4QxZh2b68kYY8zm9tH6EjtlicIYY3ZBcwqPzswUliiMMWa37JOlTXfKEoUxxuwSK1EYY4zZXAe3Udg4CmOMMVuyEoUxxuyK3Zk2vB2FWqIQkR4R+SsROdf6t3uT63Ii8iciclZEXhKRv7vXsRpjzHbaaGT2rgq76ulDwNdU9QTwtdb+Rv4N8BeqehdwP/DSHsVnjDE3LoSFi/ZC2IniEeCzrdefZYMFNEQkA/znwL8FUNWaqtpCD8aY9qIgwc63/SDsRDGoqtPQnNWQdTMhrnMMmAf+QESeFpHfE5GuvQzSGGNuiJUobo6IfFVEXthge+QGb+EBbwQ+papvAIpsUkUlIo+JyFMi8lSd6i79FxhjzA3Sm9j2gdve62mzRTcARGRWRIZbc6QPA3MbXDYBTKjqN1v7f8ImiaK1fOATABnp2Sf/C4wxnaJTB9yFXfX0ReAnW69/EvjC9Reo6gxwVUTubB16K/Di3oRnjDE7YFVPt8XHgbeLyDng7a19RGRERJ5cd90/B/5IRJ4DHgD+lz2P1BhjtqI053ra6bYPhDrgTlUXaZYQrj8+Bbxr3f4zwKaLjxtjTNgE7diqJxuZbYwxu8UShTHGmC1ZojDGGLOpa20UHcgShTHG7BJrozDGGLM1SxTGGGM2t3/GRexU2OMojDHGtDkrURhjzG5QOrZEYYnCGGN2i/V6MsYYsxXr9WSMMWZrliiMMcZsSoHAEoUxxphNdW73WEsUxhizWyxRGGOM2VKHJgobcGeMMbvhWhvFTrddICIfFBEVkb7rjh8SkYKIfHDdsTeJyPMicl5EPikist39LVEYY8yuUNBg59stEpFxmiuEXtng9CeAL1937FPAY8CJ1vbwds+wRGGMMbslnDWzPwH8HM0yzWtE5D3ABeDMumPDQEZVv66qCnwOeM92D7BEYYwxuyGEqicReTcwqarPXne8C/h54Feve8soMLFuf6J1bEvWmG2MMbvl5koIfSLy1Lr9J1T1iWs7IvJVYGiD930E+AXgHRuc+1XgE6pauK4JYqP2iG2DtkRhjDG75eYSxYKqPrj5LfVtGx0XkXuBo8CzrWQwBnxXRB4C3gz8qIj8BpADAhGpAH/auu6aMWBquwBDTRQi0gP8e+AIcAn4J6q6vMF1/yPw0zQz3/PAT6lqZe8iNcaY7eztgDtVfR4YuLYvIpeAB1V1AfiBdcd/BSio6m+39vMi8hbgm8CjwG9t96yw2yg+BHxNVU8AX2vtv46IjAL/guYP4B7ABd63p1EaY8x2FAiCnW977wPA7wHngf+/vbuLsaOs4zj+/fnS4Cpi11pYWhCMxPhWIyGkVhJriqYtwYqJCcaXGo2ECxK50KSmiTThBo16YRQTRZOqRKJYpSFFQYUQJW2oDW0h1ZYKaunaQsWiiQaDfy/mWTuePWfOnN152TP7+yQnO2fmec78dnb2/PPMzJlzlNlXRc3S9qGnTcDaNL0deIDsBEyvlwAvk/RvYIISQyUzs8UkIi4aMH9bz/O9wFtGee22RxTnRsQ0QPq5vLdBRDwFfInsGuFp4HRE3NtoSjOzMtq5PLZ2tY8ohpyxL9N/KdnI42Lgb8CPJH0kIr7fp+11ZB8k4Swm5pzZzGxOxuSNf1S1F4pBZ+wBJJ2QNBUR0+mDICf7NLsSeCIink59dgBrgFmFIl1S9k2AV2qym38xM1ugqrslx0LT9qGnncDmNL0ZuKtPmz8BqyVNpHuSrAMONZTPzKycgIj/jPwYB20XiluA90g6QnavklsAJJ0vaRdAROwB7gT2kV0a+yLSqMHMbEFp6aaAdWv1qqeIOEU2QuidfxzYmHt+E3BTg9HMzEbncxRmZjZQRFufi6idC4WZWVU8ojAzsyLhEYWZmQ02Ph+gG5ULhZlZFWa+j6KDXCjMzKoyJp+LGJULhZlZBQIIjyjMzGygCI8ozMysmEcUZmZWrKMjCkVHL+eS9DTwxyHNlgHPNBBnvsYh5zhkBOesWldyvjYiXjOfFUj6WVrPqJ6JiPXzWXfdOlsoypC0t+hLzReKccg5DhnBOavmnItD23ePNTOzBc6FwszMCi32QjEu32sxDjnHISM4Z9WccxFY1OcozMxsuMU+ojAzsyE6XSgkTUq6T9KR9HNpnzZvkPRI7vGcpBvTsm2Snsot2zh7Lc3kTO2elHQwZdk7av8mckq6QNL9kg5JekzSp3PLat2ektZL+r2kxyVt6bNckr6alh+QdGnZvg3n/HDKd0DSQ5LellvWdx9oIeNaSadzf8vPl+3bcM7P5jI+KukFSZNpWSPbshMiorMP4IvAljS9BfjCkPYvBv5Cdk01wDbgMwslJ/AksGy+v2edOYEp4NI0fTZwGHhT3dsz/e2OAq8DlgD7Z9aba7MRuAcQsBrYU7ZvwznXAEvT9IaZnEX7QAsZ1wJ3z6Vvkzl72l8N/KrJbdmVR6dHFMAmYHua3g68f0j7dcDRiBj2Qb2qjZqz6v6VrScipiNiX5r+O3AIWFFTnrzLgccj4g8R8TxwR8qbtwn4bmR2A6+SNFWyb2M5I+KhiHg2Pd0NrKwpy5wz1tS37pwfAn5QU5ZO63qhODcipiF7AwOWD2l/LbN3pBvSIYDv1HVIh/I5A7hX0m8lXTeH/k3lBEDSRcDbgT252XVtzxXAn3PPjzG7QA1qU6ZvVUZd1yfJRkEzBu0DVSqb8R2S9ku6R9KbR+xbhdLrkjQBrAd+nJvdxLbshLG/15OkXwDn9Vm0dcTXWQK8D/hcbvY3gJvJdqibgS8Dn2gx5zsj4rik5cB9kn4XEQ/OJc8gFW7PV5D9U94YEc+l2ZVtz36r7DOv95K+QW3K9K1K6XVJejdZobgiN7v2faBkxn1kh2j/kc41/RS4pGTfqoyyrquB30TEX3PzmtiWnTD2hSIirhy0TNIJSVMRMZ0OMZwseKkNwL6IOJF77f9NS/oWcHebOSPiePp5UtJPyIbeDwKj/J6155T0UrIicXtE7Mi9dmXbs49jwAW55yuB4yXbLCnRtyplciJpFXAbsCEiTs3ML9gHGs2YK/5ExC5Jt0paVqZvkzlzZh0taGhbdkLXDz3tBDan6c3AXQVtZx2/TG+GM64BHq003RlDc0p6uaSzZ6aB9+byjPJ71p1TwLeBQxHxlZ5ldW7Ph4FLJF2cRofXprx5O4GPpaufVgOn0yG0Mn0byynpQmAH8NGIOJybX7QPNJ3xvPS3RtLlZO8lp8r0bTJnyncO8C5y+2uD27Ib2j6bXucDeDXwS+BI+jmZ5p8P7Mq1myDbyc/p6f894CBwgGwHnGorJ9mVHfvT4zFg67D+LeW8gmz4fwB4JD02NrE9ya5qOkx2JczWNO964Po0LeDraflB4LKivjXul8Ny3gY8m9t+e4ftAy1kvCFl2E92wn3NQtyW6fnHgTt6+jW2Lbvw8CezzcysUNcPPZmZ2Ty5UJiZWSEXCjMzK+RCYWZmhVwozMyskAuFmZkVcqEwM7NCLhTWaZJ2p5sTImmF9+qg/AAAANlJREFUv3fAbHQuFNZZ6RYTFwIzt41fRfaJbDMbwdjfFNCswOuBJ+LM7QdWAQfTvX1uBZ4HHoiI29sKaDYOPKKwLnsr/z+CuIzsPlMfAO6MiE+R3VrezAq4UFiXTQL/BJD0RuAqssKxkjNfePNCO9HMxocLhXXZz4F1kn4IfBA4Fdl3YhzjzNeL+n/AbAjfPdYWnXSO4mvAv4Bf+xyFWTEXCjMzK+Rht5mZFXKhMDOzQi4UZmZWyIXCzMwKuVCYmVkhFwozMyvkQmFmZoVcKMzMrJALhZmZFfovqukBza5HoVMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def eval_ll_surface(param=\"piEN\"):\n",
    "    ln_DeltaF = float(mcmc.get_samples()[\"ln_DeltaF\"][-1])\n",
    "    DeltaF = jnp.exp(ln_DeltaF)\n",
    "    ln_Fbase = float(mcmc.get_samples()[\"ln_Fbase\"][-1])\n",
    "    Fbase = jnp.exp(ln_Fbase)\n",
    "    t0 = float(mcmc.get_samples()[\"t0\"][-1])\n",
    "    ln_tE = float(mcmc.get_samples()[\"ln_tE\"][-1])\n",
    "    tE = jnp.exp(ln_tE)\n",
    "    piEE = float(mcmc.get_samples()[\"piEE\"][-1])\n",
    "\n",
    "    piEN_ = jnp.linspace(-0.9, 0.9, 800)\n",
    "    u0_ = jnp.linspace(-0.8, 0.8, 800)\n",
    "    u0_grid, piEN_grid = jnp.meshgrid(u0_, piEN_)\n",
    "\n",
    "    # Compute trajectory\n",
    "    zeta_e_t = jnp.interp(t, points, zeta_e)\n",
    "    zeta_n_t = jnp.interp(t, points, zeta_n)\n",
    "    zeta_e_t0 = jnp.interp(t0, points, zeta_e)\n",
    "    zeta_n_t0 = jnp.interp(t0, points, zeta_n)\n",
    "    zeta_e_dot_t0 = jnp.interp(t0, points, zeta_e_dot)\n",
    "    zeta_n_dot_t0 = jnp.interp(t0, points, zeta_n_dot)\n",
    "\n",
    "    delta_zeta_e = zeta_e_t - zeta_e_t0 - (t - t0) * zeta_e_dot_t0\n",
    "    delta_zeta_n = zeta_n_t - zeta_n_t0 - (t - t0) * zeta_n_dot_t0\n",
    "    ln_c = float(mcmc.get_samples()[\"ln_c\"][-1])\n",
    "\n",
    "    def eval_ll(u0, piEx):\n",
    "        u_per = u0 + piEx * delta_zeta_e - piEE * delta_zeta_n\n",
    "        u_par = (t - t0) / tE + piEE * delta_zeta_e + piEx * delta_zeta_n\n",
    "\n",
    "        u = jnp.sqrt(u_per ** 2 + u_par ** 2)\n",
    "\n",
    "        # Magnification\n",
    "        A_u = (u ** 2 + 2) / (u * jnp.sqrt(u ** 2 + 4))\n",
    "        A_u0 = (u0 ** 2 + 2) / (jnp.abs(u0) * jnp.sqrt(u0 ** 2 + 4))\n",
    "        A = (A_u - 1) / (A_u0 - 1)\n",
    "        F_pred = DeltaF * A + Fbase\n",
    "\n",
    "        return jnp.sum(-0.5 * (F_pred - F) ** 2 / (jnp.exp(ln_c) * Ferr) ** 2)\n",
    "\n",
    "    eval_ll_vectorized = jit(vmap(vmap(eval_ll)))\n",
    "\n",
    "    return u0_grid, piEN_grid, eval_ll_vectorized(u0_grid, piEN_grid)\n",
    "\n",
    "\n",
    "u0_grid, piEN_grid, ll = eval_ll_surface()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "im = ax.pcolormesh(u0_grid, piEN_grid, ll, vmin=-450)\n",
    "fig.colorbar(im, orientation=\"vertical\")\n",
    "\n",
    "ax.set(xlabel=r\"$u_0$\", ylabel=r\"$\\pi_{EN}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AutoBNAF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize guide\n",
    "# Hidden layer i has ``hidden_factors[i]`` hidden units per input dimension (H in the paper)\n",
    "guide = AutoBNAFNormal(model, num_flows=5, hidden_factors=[25, 25, 25, 25])\n",
    "\n",
    "schedule = optim.optimizers.polynomial_decay(0.01, 10, 0.001)\n",
    "svi = SVI(model, guide, optim.Adam(schedule), Trace_ELBO(num_particles=10))\n",
    "svi_state = svi.init(random.PRNGKey(1), t, F, Ferr)\n",
    "\n",
    "last_state, losses = lax.scan(\n",
    "    lambda state, i: svi.update(state, t, F, Ferr), svi_state, jnp.zeros(40000)\n",
    ")\n",
    "params = svi.get_params(last_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbf0lEQVR4nO3deZwU9Z3/8fdnhuEUUQJRF9yAghp+SbxmSdQk3goa1kc2Mer+cu26cTfG3ZhsohiPNWsSyc9kE/OLUTG65PCKVzxAUfHAgwUGBEG5EcMAcggznDPMTH/2j64pe4aeoWe6a6qr+/V8PHhMd3XVtz6WM/Oe+n6rvmXuLgAAJKki7gIAAMWDUAAAhAgFAECIUAAAhAgFAECoV9wF5GPIkCE+YsSIuMsAgESZN2/eFncfmu2zRIfCiBEjVFNTE3cZAJAoZvZuR5/RfQQACBEKAIAQoQAACBEKAIBQ0Qw0m1mFpJskHSipxt1/F3NJAFB2Ij1TMLN7zGyTmS1ut3ycmS0zs5VmNjFYfIGkYZKaJNVGWRcAILuou4+mSBqXucDMKiXdJmm8pDGSLjGzMZKOljTL3b8r6ZsR1wUAyCLSUHD3mZK2tls8VtJKd1/t7nslPaD0WUKtpG3BOi1R1rVt115NfXNDlLsAgESKY6B5mKS1Ge9rg2WPSjrXzP6/pJkdbWxml5lZjZnVbN68uVsF/N/fzta37puv9+oburU9AJSqOELBsixzd9/t7pe6+7+6+20dbezuk9292t2rhw7Nepf2fr29YbskadbqLd3aHgBKVRyhUCvp8Iz3wyWtj6EO7d4baS8VACROHKEwV9JoMxtpZr0lXSzpia40YGYTzGxyfX19XoWkUjyKFAAyRX1J6v2SZkk62sxqzexSd2+WdIWk6ZKWSPqTu7/VlXbd/Ul3v2zQoEF51ddCKABAG5HevObul3SwfJqkaVHuOxdkAgC0lchpLgrVfUQmAEBbiQyFQnUfuRMLAJApkaFQKIwpAEBbZR0KZAIAtJXIUCjYJal0HwFAG4kMhUKNKfSqyHZzNQCUr0SGQqF8fHh+oQIApaasQ4FrUgGgrbIOhUXr8huTAIBSk8hQKNRA8/S33itQRQBQGhIZCoUaaAYAtJXIUCiUjdsb4y4BAIpKWYfCuro9cZcAAEWlrEMBANBWIkOhUAPNAIC2EhkKDDQDQDQSGQoAgGgQCgCAEKEAAAgRCgCAEKEAAAglMhS4JBUAopHIUOCSVACIRiJDAQAQDUIBABAiFAAAIUIBABAiFAAAIUIBABAiFAAAoUSGAjevAUA0EhkK3LwGANFIZCgAAKJBKAAAQoQCACBEKAAAQoQCACBEKAAAQoQCACBEKAAAQoQCACBEKAAAQokMBeY+AoBoJDIUmPsIAKKRyFAAAESDUAAAhAgFAECIUAAAhAgFAECIUAAAhAgFAECIUAAAhAgFAECIUAAAhAgFAECo7EPhjb9s059q1sZdBgAUhV5xFxC3z//mdUnSl6oPj7kSAIhf2Z8pAAA+QCgAAEKEAgAgVDShYGanmdkrZnaHmZ0Wdz0AUI4iDQUzu8fMNpnZ4nbLx5nZMjNbaWYTg8UuaaekvpJqo6wLAJBd1GcKUySNy1xgZpWSbpM0XtIYSZeY2RhJr7j7eElXS/phxHUBALKINBTcfaakre0Wj5W00t1Xu/teSQ9IusDdU8Hn2yT16ahNM7vMzGrMrGbz5s2R1A0A5SqOMYVhkjLvFquVNMzM/s7M7pT0B0m/7mhjd5/s7tXuXj106NCISwWA8hJHKFiWZe7uj7r7P7v7Re7+UpQF/POpR0TZPAAkVhyhUCsp8/bh4ZLWd6UBM5tgZpPr6+u7VcBHBg/o1nYAUOriCIW5kkab2Ugz6y3pYklPdKUBd3/S3S8bNGhQtwqwbOcqAIDIL0m9X9IsSUebWa2ZXeruzZKukDRd0hJJf3L3t6KsY5+6enJnAJAgkU6I5+6XdLB8mqRpUe67M5wpAEB2RXNHc1fkO6YAAMgukaGQ95gCHUgAkFUiQyFvZAIAZFWeoQAAyCqRoZDvmAInCgCQXSJDIf/7FIgFAMgmkaEAAIgGoQAACJVlKNB5BADZJTIU8h5ozjMV3F1NLan9rwgACZPIUMh3oDlfU15fo9HXPq3NOxpj2T8ARCWRoRC3x95YJ0laX7cn5koAoLDKMhS4IhUAsivPUGCoGQCySmQoMEsqAEQjkaEQ90BzWEeseweAwktkKOQr3zEFOp8AlKqcQsHMvm1mB1ra3WY238zOibo4AEDPyvVM4R/dfbukcyQNlfQPkiZFVhUAIBa5hkJrj8l5kv7b3ReKXhQAKDm5hsI8M3tW6VCYbmYDJSV2nodCTZ3tzlAzgNLSK8f1LpV0nKTV7r7bzAYr3YUUCzObIGnCqFGjurd9/gXk2wIAFKVczxROkrTM3evM7MuSrpMU200CxXJJKgCUmlxD4XZJu83sWElXSXpX0u8jqwoAEItcQ6HZ0x3oF0i61d1vlTQwurKiRe8PAGSX65jCDjO7RtJXJH3GzColVUVXVrQKNfcRw8wASk2uZwoXSWpU+n6F9yQNk3RLZFUVOU40AJSqnEIhCIJ7JQ0ys89JanD3xI4p0H0EANnlOs3FlyTNkXShpC9Jmm1mX4yysCidcuSQuEsAgKKU65jCtZL+xt03SZKZDZX0vKSHoyosSv37VMZdAgAUpVzHFCpaAyHwfhe2LbhieZ4CNzQDKDW5/mJ/xsymm9nXzezrkqZKmhZdWZ3L9+a1fIcUGJMAUKpy6j5y9++b2RcknaL079TJ7v5YpJUBAHpcrmMKcvdHJD0SYS0AgJh1GgpmtkPZ79EySe7uB0ZSFQAgFp2GgrsndiqLzhRq6mzuaQZQasryGc35YpwZQKkiFAAAobIMBf7SB4DsyjIUAADZEQp54I5mAKWGUOiGwl29BADFJZGhkO/cR/xOB4DsEhkK+c59BADILpGhEDdnMAFAiSIUAAChsgyFfAeKGWgGUKrKMhQAANkRCgCAEKEQSKUYPAYAQiFw5YMLurwNMQKg1BAKgScWrs95XYaZAZQqQgEAECIUMnzpzllxlwAAsSIUMsx5Z2vcJQBArAiFPDDbBYBSQyi009DUst91uKEZQKkiFNo55vpn4i4BAGJDKJSA+t1NenrRhrjLAFACCIUScMX98/XNe+drXd2euEsBkHBFFQpmNsDM5pnZ5+KuJRfF8lyFtVt3S5L2NqdirgRA0kUaCmZ2j5ltMrPF7ZaPM7NlZrbSzCZmfHS1pD9FWVMuFq6t6/Rz455mACUq6jOFKZLGZS4ws0pJt0kaL2mMpEvMbIyZnSXpbUkbI65pvy647TUmyANQliINBXefKan9HWFjJa1099XuvlfSA5IukHS6pE9J+ntJ3zCzrLWZ2WVmVmNmNZs3b46s9k//9AWt3LQzsvYBoBjFMaYwTNLajPe1koa5+7XufqWk+yTd5e5ZO8jdfbK7V7t79dChQyMrcn19g8679ZXI2geAYhRHKGTrkA/7atx9irs/1YP1dGhvS0rurttfWqX1Wa7soYMJQKmJIxRqJR2e8X64pNznrZZkZhPMbHJ9fX1BC8tm8brt+ukzS3XZH2oyCoh8twAQizhCYa6k0WY20sx6S7pY0hNdacDdn3T3ywYNGhRJgZkm/PpVSdKuxv1PfwEASdcrysbN7H5Jp0kaYma1kv7D3e82syskTZdUKeked38ryjoK4Z0tu3Tz00voMwJQ0iINBXe/pIPl0yRNi3LfUbjz5dWSpLEjB8dcCQBEo6juaM5VT44pZNP63IU4b2h2d93x8ipt3bU3viIAlJxEhkJPjil0JhVjKtS8u02Tnl6qqx5+M7YaAJSeRIZCsbh/zl9i23dTMM/Rrsbm2GoAUHoIhTxMXbRBLUyHAaCEJDIU4h5TaOUuHfmDaWps5nJVAKUhkaFQLGMKrer3NMVdAgAURCJDoRDOHnNIwdpKxfAYAzqtAEShbEPhrq9W6xcXHVuQtuK8CkmSWorkYT8Aki/Sm9eiYmYTJE0YNWpUXu18/vjh+vzxw9ssS6VcLy/frH+YMjfndk6e9ILWTDo/r1q6a9bq9wva3o6GJvWrqlSvyrL9ewEoa4n8yY9yTKGiwnT6MR/Wmknnh/9m/Pup+91uxMSpBa8lDh+/8Vl976GFcZcBICaJDIWeduTQA8KAWP2T8zpc786XV3V7H6+s2Kw//s+73d6+kP68oEuT1gIoIYRCF1VUWIfhcPPTS7VpR0O32v3K3XN03Z8X73/FQBTDCFxFBYBQ6KbWcGhv7I9nxFBN9yxYW6cRE6dqdjAu0djE/RZAuSMU8pQtGBaurYuhkvQkeV3x2sotkqSXlkf3rGsAyZLIUCiWO5pbtQ+GC257LZL9HP+fz+rqKCfA44lyQNlLZCgU2x3NkjT7B2e2eX/5vfMKvo9tu5v0YM3agrcLAK0SGQrF6JAD+7Z5P23Re13uzukKL+A9za1lGqcKQNkjFApoVbsrkq647w3taIjmip5C5I3FlAEPzv2LXl2xJZ6dA+gUoVBAlRVtf8tOXbRBJ9z0XKfbzFiyUe9s2dXlfTW1xDDhUoFc/cgiffnu2XGXASALQqHAlt40rs37phbXb19Z3WFX0qW/q9HpP3sp62c/mbZED87N/iCfS39Xk1ed2cR15gCgeBAKBda3qnKfZT+auqTN3cpNLSmNvGaqHp5X22lbk2eu1tWPLCp4je0VcnwCQLIlMhSK7ZLU9p7/7r5zJV3/+Fu68I7X9fC8Wm3Z2Sh3ZZ1j6NIuTMSXr/YDy5woAEhkKBTjJamZRn34gKzL567Zpu89tFAn3fxCh9vOWLpJX/5t2/729XV7ul1L3e69unTKXG3dtbfbbQAoH4kMhSR48opPd3vbV1duaTPr6smTXtCvZqwI3+9qbM65rd/Pelczlm7S5Jmr2yx397DbaGdDs2584i01Nid38BpAYSTyeQpJ8PHhhT2L+a/nloevb3rq7Zy3a+0SuuPlVZo4/phw+chrpoWv752dfTC7M5t2NKhfVaUG9q3q8rYAihehEKGlN43TMdc/U/B2H5ib/a7mM37+sr552pFqbklp6psbtL6+7Yyt97z6joYd3E9rOrgEdsrra3KuYeyPZ2jIAb1Vc93ZOW8DoPgRChHqW1WpRy8/WX/3m9d7bJ+3v9TxMx3+swtnGLnYspNxCqDUMKYQsRP++mCtmXS+Ft14joYd1C/ucrqlfk+TUikuWwXKAWcKPWRg3yq9NvGM8P2y93Zo5JAB2rO3RYP6V6m5JaWFtXX6wu2zYqxyX+/vbNSJP3pe3z5ztL5z9lGdrvuX93dryMDe6t+bbysgqRJ5plDs9ynk4uhDB6p3rwoN6p8eqO1VWaETPzJYayadr3v/6ZMxV/eBzTsbJUlPL96w33U/e8uL+vo9PXefBYDCS2QoFPt9Cvk6ZdQQXXf+R2Ot4fJ752ld3Z7wBrdcJ+Cbs2ZrhFUBiFoiQ6Ec/NNnjtCiG8+Jbf/TFr2nm6ctUescf2u37c46f9Mvn1/e4dxNAJKHUChiA/tW6ewxh8S2/wqzcJK8hqaU7puz7/0Mv3x+RbdmeQVQnAiFInfXV6t16lFDY9l3+izhgxmRrn1scSx1AOg5hEIC/O4fx+rWi4/r8f1mnil018vLN2su4wxAYhAKCVFVmf5f9cUTh/fYPs1MFe1SoatdRV+7Z44uvKO4LrMF0DFCISHOGXOIvnX6kbr+/DGa84MzNf/66KeXqLB9p9N+qCb7FBsASgOhkBC9Kiv0/XOP0aD+VfrwgX01eEDvyLuUNu5o3Kf76DedTKMBIPkIhQS74LhhkbY/c/nmfR7EA6C0EQoJd/lpR+53nYf/5aTI6+joGdQAkoVQSLirxh2jV646XV84oeMB6OoRg/XzC4/dZ/ljl5+83/ZzvfqITABKA6FQAg4f3F8/u/ATWv6j8R2uM3bk4H2WjfjQgP22ne8lqQCSJZHTWZrZBEkTRo0aFXcpRcPM1LuXadY1Z2R9BvThg/trzaTztXBtnbbsbNTRhw7UwQN667a/P0Hfum9+p+3mghMFoDQkMhTc/UlJT1ZXV38j7lqKzWGD+mnNpPMlSS8s3agD+rR9XOaxhx/U5v34jx3aaXu5nigwpgCUhkSGAnJzxjH7nzdpfycC7W9ek6R5727rbkkAihxjCuhUttD4wu37Pl6U8wSgNBAKZS5zzGBA78p9P8+xnWy9Rz94bFE3qwIQF0IBoWx/7e9paslp29nvvL/Psvtm7zvVNoDiRihA3znrKI34UH/94qJ9p8049ZaXcmrjK3fPCV+v2LhDU157p8t1pFKuPXtzC6FCenLhev1+1poe3y9QjAgF6NtnjdZL3z9d5/6fQzX/+rO16ifn6fnvfrbb7Z39i5m68cm32yxraGrRzsZmvbB0Y7hs8bp6Pb5gnZpaUnp1xRZd+eACffSGZzRi4lSd+4uZkqS9zSn9YdYataQ+OI85ZdILGjFxqqYt2qC7Zq7W2q27JUm3vbhSqzfv1Lq6PWpoatFdM1dryYbtkqSN2xv09KIN2rSjQc0tKZ16y4uasWSjmltS+tf739ANj7+lrbv2hvtoakmpoalF7q5UsO+l723XjoamcB13167GZu1sbA7fP/ZGrfY2p7f979fe0SdunK55725VY3OL6vekt31kXq3W1+3Rrc+v0E1PpY9T7bbdqmk3xfiOhqawppaU62fTl4VXeTW3pPR+8Pzs9sf58QXrwpr+/MY6TX2z7fO1W49B+yvGmlpSmrXq/bCdzBlxN25v0Lq6PfvsrzuaWlJqbkl1aZtVm3cWZN9d0dyS0s7GZrm7Nu1oiHx/G+r3hN9rcbIkX0pYXV3tNTU1cZdRdty9zRjC7qYWzX93m655dFHBfnGgtFRWWJtgL5TDBvXVzoZm7QhC8PDB/bR2a/p7sH/vSu3u4MyzV4WpOaOeQf2qwtCWpH87c7R+NWNF1m3//eyj9PPnluuK00fp1y+u1OABvdv8QdHqtKOHau47W7UrqGH4wf1Uu63tz8d3zjpKv31ldVi/JA05oI+2ZAn89lb+eLx6VXbv73ozm+fu1Vk/IxRQaCMmTo27BKDkHf/XB+mxy0/p1radhQL3KaDgWm+ea+Xuak65mlpS2tnQrD1NLdrbnFLdnqbw/UH9qvTG2jod1L9Khx7YVy8t26wD+vbS2+u3a0CfSm3fk/5rcMeeJh08oLc+NKC3tuxs1IK1dUq51K+qcr+D4n16VaixuWvdFkCxuuWLn4ikXUIBkTMzVVWaqior1L93x99yJ48aEr4+86P7v/EOQOEx0AwACBEKAIAQoQAACBEKAIAQoQAACBEKAIAQoQAACBEKAIBQoqe5MLPNkt7t5uZDJG0pYDmFQl1dQ11dU6x1ScVbWynW9RF3H5rtg0SHQj7MrKajuT/iRF1dQ11dU6x1ScVbW7nVRfcRACBEKAAAQuUcCpPjLqAD1NU11NU1xVqXVLy1lVVdZTumAADYVzmfKQAA2iEUAAChsgwFMxtnZsvMbKWZTeyB/a0xs0VmtsDMaoJlg83sOTNbEXw9OGP9a4LalpnZuRnLTwzaWWlmvzIz60Yt95jZJjNbnLGsYLWYWR8zezBYPtvMRuRR141mti44bgvM7LyerMvMDjezF81siZm9ZWbfLobj1UldcR+vvmY2x8wWBnX9sEiOV0d1xXq8MtqsNLM3zOypYjhewUPYy+efpEpJqyQdIam3pIWSxkS8zzWShrRb9v8kTQxeT5T00+D1mKCmPpJGBrVWBp/NkXSSJJP0tKTx3ajls5JOkLQ4ilokXS7pjuD1xZIezKOuGyV9L8u6PVKXpMMknRC8HihpebDvWI9XJ3XFfbxM0gHB6ypJsyV9qgiOV0d1xXq8Mvb3XUn3SXqqGH4eY/vlHNe/4MBNz3h/jaRrIt7nGu0bCsskHRa8PkzSsmz1SJoe1HyYpKUZyy+RdGc36xmhtr98C1ZL6zrB615K33Fp3ayrox/aHq0ro73HJZ1dLMcrS11Fc7wk9Zc0X9Ini+l4tasr9uMlabikGZLO0AehEOvxKsfuo2GS1ma8rw2WRcklPWtm88zssmDZIe6+QZKCrx/eT33DgtftlxdCIWsJt3H3Zkn1kj6UR21XmNmblu5eaj2N7vG6gtPu45X+K7Nojle7uqSYj1fQFbJA0iZJz7l7URyvDuqS4v/++qWkqySlMpbFerzKMRSy9cNHfV3uKe5+gqTxkr5lZp/tZN2O6ouj7u7UUsg6b5d0pKTjJG2Q9PM46jKzAyQ9IulKd9/e2aox1xX78XL3Fnc/Tum/gMea2cc6+0+Iua5Yj5eZfU7SJneft7/6e7KucgyFWkmHZ7wfLml9lDt09/XB102SHpM0VtJGMztMkoKvm/ZTX23wOoq6C1lLuI2Z9ZI0SNLW7hTl7huDH+aUpLuUPm49WpeZVSn9i/ded380WBz78cpWVzEcr1buXifpJUnjVATHK1tdRXC8TpH0t2a2RtIDks4wsz8q5uNVjqEwV9JoMxtpZr2VHnx5IqqdmdkAMxvY+lrSOZIWB/v8WrDa15TuF1aw/OLgqoGRkkZLmhOcRu4ws08FVxZ8NWObfBWylsy2vijpBQ86NLuq9Qcj8Hmlj1uP1RW0cbekJe7+XxkfxXq8OqqrCI7XUDM7KHjdT9JZkpYWwfHKWlfcx8vdr3H34e4+QunfQy+4+5fjPl5dHqQshX+SzlP6io1Vkq6NeF9HKH3FwEJJb7XuT+l+vRmSVgRfB2dsc21Q2zJlXGEkqVrpb9xVkn6tLg5IBm3cr/SpcpPSf0VcWshaJPWV9JCklUpfEXFEHnX9QdIiSW8G39yH9WRdkj6t9Kn2m5IWBP/Oi/t4dVJX3MfrE5LeCPa/WNINhf5eL3BdsR6vdjWepg8GmmM9XkxzAQAIlWP3EQCgA4QCACBEKAAAQoQCACBEKAAAQoQC0E1mdqWZ9Y+7DqCQuCQV6KbgTtRqd98Sdy1AoXCmAOQguDN9qqXn5F9sZv8h6a8kvWhmLwbrnGNms8xsvpk9FMxN1Po8jZ9aek7/OWY2Klh+YdDWQjObGd9/HfABQgHIzThJ6939WHf/mNKzW66XdLq7n25mQyRdJ+ksT09+WKP0PPmttrv7WKXvNv1lsOwGSee6+7GS/ran/kOAzhAKQG4WSTor+Iv/M+5e3+7zTyn9EJTXgimavybpIxmf35/x9aTg9WuSppjZN5R++BMQu15xFwAkgbsvN7MTlZ5j6GYze7bdKqb0PP2XdNRE+9fu/i9m9klJ50taYGbHufv7ha4d6ArOFIAcmNlfSdrt7n+U9DOlHx26Q+nHYUrS/0g6JWO8oL+ZHZXRxEUZX2cF6xzp7rPd/Qaln4iVOS0yEAvOFIDcfFzSLWaWUnom128q3Q30tJltCMYVvi7pfjPrE2xzndKz8UpSHzObrfQfYq1nE7eY2WilzzJmKD2TLhArLkkFIsalq0gSuo8AACHOFAAAIc4UAAAhQgEAECIUAAAhQgEAECIUAACh/wU5aY6CKb30EgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(losses)\n",
    "plt.ylabel(\"loss\")\n",
    "plt.xlabel(\"steps\")\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAFDCAYAAAAd/vg3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAei0lEQVR4nO3df5Td9V3n8ed7fk8yJDMhkx8kFLAlEegWLUhbtW4iYgE9Ul3dg7UVu3ZzONpud1211B717PHsbl33uNXTdlkOreLqymYrp8UaxZYmVU+lAi0CKQYiFAgJhPzOJDOZH/e9f9w702E6ydwkn7l3bvJ8nHPPfL/3+7n3vnInM6/5fr/3+/1GZiJJ0tlqa3YASdK5wUKRJBVhoUiSirBQJElFWCiSpCIsFElSEQuiUCLixojYERE7I+KOWZYvjYg/j4h/jIjtEfHeZuSUJJ1cNPs4lIhoB54GbgB2AQ8DP52Z35g25teApZn5oYgYBHYAqzJztBmZJUnfbiGsoVwH7MzMZ2sFcS9wy4wxCVwQEQH0AQeA8cbGlCSdykIolDXAi9Pmd9Xum+7jwBXAbuAJ4IOZWWlMPElSPTqaHQCIWe6buR3uHcBjwA8Crwe+EBF/m5lHXvNEEZuATQA9PT3XvO51r5uHuGVVKhXa2hZCr5+aOcsyZzmtkBFaJ+fTTz+9LzMHz+jBmdnUG/A24IFp8x8GPjxjzF8Ab582/yXgulM977p167IVbN26tdkR6mLOssxZTitkzGydnMAjeYa/zxdCXT4MXB4Rl0VEF3ArcP+MMS8A1wNExEpgPfBsQ1NKkk6p6Zu8MnM8It4PPAC0A5/OzO0RcXtt+Z3AbwF/GBFPUN1E9qHM3Ne00JKkb9P0QgHIzC3Alhn33Tltejfww43OJUmq30LY5CVJOgdYKJKkIiwUSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUUsiEKJiBsjYkdE7IyIO04yZkNEPBYR2yPiy43OKEk6tY5mB4iIduATwA3ALuDhiLg/M78xbUw/8Engxsx8ISJWNCetJOlkFsIaynXAzsx8NjNHgXuBW2aMeRdwX2a+AJCZexucUZI0h4VQKGuAF6fN76rdN906YCAitkXEoxHxsw1LJ0mqS2RmcwNE/BTwjsx8X23+PcB1mfmBaWM+DlwLXA/0An8P/EhmPj3juTYBmwAGBwev2bx5c2P+EWdhaGiIvr6+ZseYkznLMmc5rZARWifnxo0bH83Ma8/ksU3fh0J1jeTiafNrgd2zjNmXmceAYxHxN8DVwGsKJTPvAu4CWL9+fW7YsGG+Mhezbds2zFmOOctqhZytkBFaJ+fZWAibvB4GLo+IyyKiC7gVuH/GmM8Bb4+IjohYBLwFeKrBOSVJp9D0NZTMHI+I9wMPAO3ApzNze0TcXlt+Z2Y+FRF/BTwOVIC7M/PJ5qWWJM3U9EIByMwtwJYZ9905Y/53gN9pZC5JUv0WwiYvSdI5wEKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqYkEUSkTcGBE7ImJnRNxxinHfExETEfGTjcwnSZpb0wslItqBTwA3AVcCPx0RV55k3G8DDzQ2oSSpHk0vFOA6YGdmPpuZo8C9wC2zjPsA8GfA3kaGkyTVZyEUyhrgxWnzu2r3TYmINcCPA3c2MJek81hmMjExwfj4OBMTE2RmXcvOZx3NDgDELPfN/O58DPhQZk5EzDa89kQRm4BNAIODg2zbtq1UxnkzNDRkzoLMWVYr5JyvjJMlMfk7JzOJiNeUx8xlzci5kCyEQtkFXDxtfi2we8aYa4F7a9+w5cDNETGemZ+dPigz7wLuAli/fn1u2LBhvjIXs23bNsxZjjnLaoWcJTNmJuPj41QqldeshbS3t1OpVGhvb58qlba2Njo7O6fmI4L29vaG5FyoFkKhPAxcHhGXAS8BtwLvmj4gMy+bnI6IPwQ+P7NMJOlsVSqVqYKY3JwVEZw4cYLh4WEyk87OTjo7O6fKZXqpnO+aXiiZOR4R76f66a124NOZuT0ibq8td7+JpHmVmVQqFcbGxogIxsbGGB0dBaols2/fPiYmJmhvb6e3t5fx8XG6u7sB6Oio/hqda5PX+aDphQKQmVuALTPum7VIMvPnGpFJ0vlhcm1kshCGh4cZHR1lbGyM9vZ2Dh8+zNGjR5mYmKCjo4ORkRF6e3un1mYmN4mdanPX+WJBFIokNUulUgFgYmKC4eFhxsbGmJiYoFKpMDw8zJEjR6b2q7S3tzMyMsKJEycYHR2lq6uLsbExurq6XEPBQpF0nqtUKoyPjzM6Osrw8DDAVGEMDw9TqVQYGRlheHiYiOCCCy6go6Njah/K5HO0tS2EozCay0KRdF6b3Hdy/PhxTpw4MXUDOH78OE8//TQRwaJFi+jq6uLYsWP09fXR3d1Ne3s7bW1tU5vDzne+A5LOa5VKhRMnTnDo0CEmJiYYGRlh7969HDhwAICRkRG6uroYGhqis7OTjo4Ouru76e3tndoxryoLRdJ5KzMZGxtjaGiII0eO8Pzzz3PixAn279/PoUOHGB8f58ILL2R4eJje3t6pjwd3dnbS3d09Ne8O+SoLRdJ5a3JT1+7du/nmN7/JgQMHeOWVV3jppZc4evQoPT09VCoVent7GR4eZvny5Rw5coShoaGpgx+7uroslBoLRdJ5KTM5duwYIyMjvPrqqzz55JPs27ePXbt2MTIyAlSPLTl+/DiXXHIJ/f39TExM0NPTQ2dnJydOnGDx4sVTBzbKQpF0npqYmGBoaIjdu3ezY8cOXnzxRfbv38/evXs5duwYvb29wLc+wTU2NsaSJUvo7+9nYGDgNZu8VGWhSDrvVCoVjh07xoEDB9i9ezf79u3j8OHDvPDCC4yOjvLKK6+QmQwMDDA6OkpmsnLlSnp6ehgbG5s6oj4zLZRpLBRJ55XM5MSJExw/fpzh4WG2b9/Oc889x8svv8yePXs4cuQIExMTAHR3d7No0aKpxx4/fpxKpUKlUpka4+aub/FIHEnnlYmJCU6cOMHBgwd5/vnn2bNnDy+99BL79u3j4MGDU0UBsH//fsbGxujs7JxaOzlx4sTUkfER4QGN0/hOSDqvjI+PMzw8zKuvvsqePXvYt28fL7/8Mi+99NKs4ydPCjk6Okp3d/fUUfJ9fX3ukJ/BQpF0XhkfH+fo0aNT+0wOHjxIR0cHY2Njs44/duwYo6OjLF68mMWLF9Pd3U1XVxednZ2unczguyHpvDIyMsK+ffvYvXs3Bw4cYHh4mOPHj5/yMZNF0t7ezqJFi+ju7vb8XbPw3ZB03hgbG2Pfvn1kJkNDQ+zatYu9e/dy6NChWcd3dHSwbNky+vr6aGtrY8WKFfT19dHT0/Oak0Oqyk95STovZCZHjx4FqpuxhoaGaGtrO+kaSk9PD8uWLWPt2rUMDg6yZMkSFi1axODgIN3d3a6dzMJ3RNI5b/o5u44dO8ahQ4c4duwY7e3tLF26dNbHLFmyhGXLlnHhhReyZMkSVq1axbJly7jgggtoa2uzUGbhGoqkc97ExMTUZX0PHz7M0NAQw8PDU+fg6uzsnDqtSk9PD93d3axcuZLLLruMK6+8kvXr13P55ZezYsUKent73dx1EhaKpHPexMQEExMTdHV1cfDgQQ4fPkxXVxdHjx5lfHycwcFBxsfHaWtrY/HixfT397Ny5UquuuoqrrnmGi699FIGBgYYGBjwuien4Dsj6bwwed347u5uDh06xIEDB7jgggum9p8cP36c5cuXc8kll0yVxxvf+EZe//rX09/fP3WBLZ2chSLpnNfe3s74+DgHDx4kIujv72f37t1TBygCDA4O0t/fz+rVq7nooou44ooruOqqqxgcHJy6qJanqT81C0XSOW+yCEZGRhgdHZ066n3yFCqLFy+e+gTXypUrWb16NRdeeOHURbUmC8X9JqdmoUg6L2Tm1PEjw8PDU5f4nTz5Y2dnJ11dXSxatIgLL7yQnp4eurq6pi7166e65mahSDrnVSoVgKlT0i9btowVK1awd+9eoFomfX19rFixglWrVrF06VL6+/unriHvpq76WCiSznmTp5yfnF6+fDnj4+NT1zKZPNHj8uXLp67I2N/f70eET9OCWIeLiBsjYkdE7IyIO2ZZ/jMR8Xjt9pWIuLoZOSW1nsxkYmKCiGBgYGCqSAYGBnjDG97AmjVrptZKVq9ezSWXXMKqVavo6+vzU12nqelrKBHRDnwCuAHYBTwcEfdn5jemDXsO+JeZeTAibgLuAt7S+LSSWk2lUqG9vZ2uri7Gx8dZvXo1lUqFV199lcxk6dKlRATt7e2sWrWKgYEBli5dSltbG5np2slpaHqhANcBOzPzWYCIuBe4BZgqlMz8yrTxDwFrG5pQUsvKzKkDFkdGRli0aNHUGsjhw4cZHx+nvb2d1atX093dTU9PD729vXR0dEw9VvVZCIWyBnhx2vwuTr328fPAX85rIknnjMlrv3d0dDAwMMDx48fp7Oykt7eXNWvWTO1f6ezsZNGiRbS1tU19nLijo8M1lNMQkzulmhYg4qeAd2Tm+2rz7wGuy8wPzDJ2I/BJ4Pszc/8syzcBmwAGBwev2bx587xmL2FoaGjqwKqFzJxlmbOcejJO33Q1uXN++jL41rXhJy/tW3pzVyu8lwAbN258NDOvPZPHLoQ1lF3AxdPm1wK7Zw6KiDcBdwM3zVYmAJl5F9X9K6xfvz43bNhQPGxp27Ztw5zlmLOsVshZT8bMpFKpTJXH5I768fFxOjo6pjZvTV7ud/JswiULpRXey7O1EArlYeDyiLgMeAm4FXjX9AER8TrgPuA9mfl04yNKamWTO92n6+zsfE3RtLW1uYnrLDW9UDJzPCLeDzwAtAOfzsztEXF7bfmdwG8AFwKfrH2zx890lUySJs1WNDpzTS8UgMzcAmyZcd+d06bfB7yv0bkkSfXz83CSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSiuiod2BEbK5NPg98DfhaZu6Yl1SSpJZTd6Fk5r8GiIjXA78E/C9gyTzlkiS1mLo3eUXE9RHxP4BfBf4euLhUiIi4MSJ2RMTOiLhjluUREb9fW/54RLy51GtLkso4nX0ofwAsAr4M/ENmHi4RICLagU8ANwFXAj8dEVfOGHYTcHnttgn4nyVeW5JUzukUyk8AvwUcA94dEX9aKMN1wM7MfDYzR4F7gVtmjLkF+KOsegjoj4jVhV5fklTA6RTKv8nMXcDVmfkbwP8tlGEN8OK0+V21+053jCSpiereKT/NqtrXHwY+WyBDzHJfnsEYImIT1U1iDA4Osm3btrMON9+GhobMWZA5y2qFnK2QEVon59k4nUK5LiI+CVwREVef5mNPZRev3cG/Fth9BmPIzLuAuwDWr1+fGzZsKBRx/mzbtg1zlmPOslohZytkhNbJeTbmLIWI6MvMocy8NiLWAtcA/4pym5weBi6PiMuAl4BbgXfNGHM/8P6IuBd4C3A4M/cUen1JUgH1rGV8PSJuAZ6u7UPZFRFP1fajnLXMHI+I9wMPAO3ApzNze0TcXlt+J7AFuBnYCRwH3lvitSVJ5dRTKMuofqz3OyNiH/AU8GbgDaVCZOYWqqUx/b47p00n8IulXk+SVF49hfJCZm4EqG3yugIYntdUkqSWU0+hLI2I7wW2T27ymudMkqQWVE+h9AG/DFwVEd1UN3k9mZm/Mq/JJEktpZ5CeVtm/jNARPRQPT3KVfOaSpLUcuYslMkyqU2PUDt1/XyGkiS1Hi+wJUkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFWCiSpCIsFElSERaKJKkIC0WSVISFIkkqwkKRJBVhoUiSirBQJElFNLVQImJZRHwhIp6pfR2YZczFEbE1Ip6KiO0R8cFmZJUknVqz11DuAB7MzMuBB2vzM40D/zEzrwDeCvxiRFzZwIySpDo0u1BuAe6pTd8DvHPmgMzck5lfq00fBZ4C1jQsoSSpLs0ulJWZuQeqxQGsONXgiLgU+G7gq/OeTJJ0WiIz5/cFIr4IrJpl0UeAezKzf9rYg5n5bftRasv6gC8D/zkz7zvJmE3AJoDBwcFrNm/efLbx593Q0BB9fX3NjjEnc5ZlznJaISO0Ts6NGzc+mpnXntGDM7NpN2AHsLo2vRrYcZJxncADwC/V+9zr1q3LVrB169ZmR6iLOcsyZzmtkDGzdXICj+QZ/k5v9iav+4HbatO3AZ+bOSAiAvgU8FRm/m4Ds0mSTkOzC+WjwA0R8QxwQ22eiLgoIrbUxnwf8B7gByPisdrt5ubElSSdTEczXzwz9wPXz3L/buDm2vTfAdHgaJKk09TsNRRJ0jnCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIRFookqQgLRZJUhIUiSSqiqYUSEcsi4gsR8Uzt68ApxrZHxNcj4vONzChJqk+z11DuAB7MzMuBB2vzJ/NB4KmGpJIknbZmF8otwD216XuAd842KCLWAj8C3N2gXJKk09TsQlmZmXsAal9XnGTcx4BfBSqNCiZJOj2RmfP7AhFfBFbNsugjwD2Z2T9t7MHMfM1+lIj4UeDmzPyFiNgA/HJm/uhJXmsTsAlgcHDwms2bNxf6V8yfoaEh+vr6mh1jTuYsy5zltEJGaJ2cGzdufDQzrz2jB2dm027ADmB1bXo1sGOWMf8V2AV8E3gZOA788VzPvW7dumwFW7dubXaEupizLHOW0woZM1snJ/BInuHv9GZv8rofuK02fRvwuZkDMvPDmbk2My8FbgW+lJnvblxESVI9ml0oHwVuiIhngBtq80TERRGxpanJJEmnpaOZL56Z+4HrZ7l/N3DzLPdvA7bNezBJ0mlr9hqKJOkcYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklSEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklSEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklSEhSJJKsJCkSQV0dRCiYhlEfGFiHim9nXgJOP6I+IzEfFPEfFURLyt0VklSafW7DWUO4AHM/Ny4MHa/Gx+D/irzPxO4GrgqQblkyTVqdmFcgtwT236HuCdMwdExBLgB4BPAWTmaGYealhCSVJdml0oKzNzD0Dt64pZxnwH8CrwBxHx9Yi4OyIWNzKkJGlukZnz+wIRXwRWzbLoI8A9mdk/bezBzHzNfpSIuBZ4CPi+zPxqRPwecCQzf32W19oEbAIYHBy8ZvPmzQX/JfNjaGiIvr6+ZseYkznLMmc5rZARWifnxo0bH83Ma8/owZnZtBuwA1hdm14N7JhlzCrgm9Pm3w78xVzPvW7dumwFW7dubXaEupizLHOW0woZM1snJ/BInuHv9GZv8rofuK02fRvwuZkDMvNl4MWIWF+763rgG42JJ0mqV7ML5aPADRHxDHBDbZ6IuCgitkwb9wHgTyLiceC7gP/S8KSSpFPqaOaLZ+Z+qmscM+/fDdw8bf4x4My26UmSGqLZayiSpHOEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklSEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklSEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBRJUhEWiiSpCAtFklREUwslIpZFxBci4pna14GTjPsPEbE9Ip6MiD+NiJ5GZ5UknVqz11DuAB7MzMuBB2vzrxERa4B/B1ybmW8E2oFbG5pSkjSnZhfKLcA9tel7gHeeZFwH0BsRHcAiYHcDskmSTkOzC2VlZu4BqH1dMXNAZr4E/HfgBWAPcDgz/7qhKSVJc+qY7xeIiC8Cq2ZZ9JE6Hz9AdU3mMuAQ8P8i4t2Z+cezjN0EbKrNnoiIJ88sdUMtB/Y1O0QdzFmWOctphYzQOjnXn+kD571QMvOHTrYsIl6JiNWZuSciVgN7Zxn2Q8Bzmflq7TH3Ad8LfFuhZOZdwF21cY9k5rUl/g3zyZxlmbOsVsjZChmhtXKe6WObvcnrfuC22vRtwOdmGfMC8NaIWBQRAVwPPNWgfJKkOjW7UD4K3BARzwA31OaJiIsiYgtAZn4V+AzwNeAJqpnvak5cSdLJzPsmr1PJzP1U1zhm3r8buHna/G8Cv3maT98qpWPOssxZVivkbIWMcB7kjMwsGUSSdJ5q9iYvSdI54pwplFY5jctp5OyPiM9ExD9FxFMR8baFmLM2tj0ivh4Rn29kxtprz5kzIi6OiK2193F7RHywQdlujIgdEbEzImY7C0RExO/Xlj8eEW9uRK4zyPkztXyPR8RXIuLqhZhz2rjviYiJiPjJRuab9vpz5oyIDRHxWO3/45cbnbGWYa7v+9KI+POI+MdazvfO+aSZeU7cgP8G3FGbvgP47VnGrAGeA3pr85uBn1toOWvL7gHeV5vuAvoXYs7a8l8C/g/w+QX6fV8NvLk2fQHwNHDlPOdqB/4Z+I7a9+8fZ74m1f2EfwkE8Fbgq014/+rJ+b3AQG36poWac9q4LwFbgJ9ciDmBfuAbwOtq8ysWaM5fm/x5AgaBA0DXqZ73nFlDoXVO4zJnzohYAvwA8CmAzBzNzEMNS1hV1/sZEWuBHwHublCumebMmZl7MvNrtemjVD92vmaec10H7MzMZzNzFLi3lnW6W4A/yqqHgP7a8ViNNGfOzPxKZh6szT4ErG1wRqjv/QT4APBnzH5MWyPUk/NdwH2Z+QJAZjYjaz05E7igdrhGH9VCGT/Vk55LhdIqp3GZMyfVvxpeBf6gtinp7ohY3MiQ1JcT4GPArwKVRgWbod6cAETEpcB3A1+d51xrgBenze/i20usnjHz7XQz/DzVtapGmzNn7USyPw7c2cBcM9Xzfq4DBiJiW0Q8GhE/27B031JPzo8DV1D9o/sJ4IOZecqf86Z+bPh0NfI0LmfjbHNS/b68GfhAZn41In6P6uacXy8UESjyfv4osDczH42IDSWzzXids30/J5+nj+pfr/8+M4+UyHaql5vlvpkfqaxnzHyrO0NEbKRaKN8/r4lmV0/OjwEfysyJ6h/VTVFPzg7gGqqHTPQCfx8RD2Xm0/Mdbpp6cr4DeAz4QeD1wBci4m9P9bPTUoWSDTyNS5Nz7gJ2ZfWgTqge2HnSnZBNzPl9wI9FxM1AD7AkIv44M9+9wHISEZ1Uy+RPMvO+kvlOYhdw8bT5tXz75tV6xsy3ujJExJuobta8KavHjzVaPTmvBe6tlcly4OaIGM/MzzYmIlD/931fZh4DjkXE3wBXU9231yj15Hwv8NGs7kTZGRHPAd8J/MPJnvRc2uTVKqdxmTNnZr4MvBgRkydpu57qTrxGqifnhzNzbWZeSvUaNV8qXSZ1mDNn7Xv9KeCpzPzdBuV6GLg8Ii6LiC6q78/9M8bcD/xs7dNeb6W6CXZPg/LVnTMiXgfcB7ynwX9FTzdnzsy8LDMvrf1//AzwCw0uk7pyUv0/+vaI6IiIRcBbaPzvoXpyvkDtwPOIWEn1pJHPnvJZG/3pgvm6ARdSvUjXM7Wvy2r3XwRsmTbuPwH/BDwJ/G+ge4Hm/C7gEeBx4LPUPmWz0HJOG7+B5nzKa86cVDfRZO29fKx2u7kB2W6m+lfnPwMfqd13O3B7bTqAT9SWP0H1InINff/qzHk3cHDae/fIQsw5Y+wf0oRPedWbE/gVqn8kPkl1E+yCy1n7Gfrr2v/NJ4F3z/WcHikvSSriXNrkJUlqIgtFklSEhSJJKsJCkSQVYaFIkoqwUCRJRVgokqQiLBSpASLiodqJKYmINRHxSHMTSeVZKNI8q5365XXA87W73kT16GPpnNJSJ4eUWtQbqJ6UdPK0FG8CnqhdkuCTwCiwLTP/pFkBpRJcQ5Hm37/gtWsk11I9r9hPAJ/JzH8L/FgzgkklWSjS/FsGDANExBVUr3D5BNVThk9e5GiiOdGkciwUaf49AFwfEZuBnwL2Z+YrVK9JMXk5XX8W1fI827DUJLV9KB8HRoC/cx+KWp2FIkkqwtVsSVIRFookqQgLRZJUhIUiSSrCQpEkFWGhSJKKsFAkSUVYKJKkIiwUSVIR/x8/3MKdsE1IzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get samples from guide\n",
    "guide_samples = guide.sample_posterior(\n",
    "    random.PRNGKey(2), params, sample_shape=(500,)\n",
    ").copy()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "\n",
    "ax.scatter(guide_samples[\"u0\"], guide_samples[\"piEN\"], color=\"black\", alpha=0.01)\n",
    "ax.set(xlim=(-0.8, 0.8), ylim=(-0.8, 0.8))\n",
    "ax.set(xlabel=r\"$u_0$\", ylabel=r\"$\\pi_{EN}$\")\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:caustic]",
   "language": "python",
   "name": "conda-env-caustic-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}