demo.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "extending-jax-demo.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Vgw63bl7SQTf"
      },
      "source": [
        "# Demo: Extending JAX with custom C++ and CUDA code\n",
        "\n",
        "This demo is adapted from an example in [the exoplanet project's documentation](https://docs.exoplanet.codes/en/stable/tutorials/intro-to-pymc3/#A-more-realistic-example:-radial-velocity-exoplanets) to work with [numpyro](http://num.pyro.ai/) and the custom C++ op defined in the [Extending JAX with custom C++ and CUDA code](https://github.com/dfm/extending-jax/) tutorial. See those tutorial for all the details."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CQ1VJptcACy0",
        "outputId": "3f3726e5-8c41-4ffa-b933-04acbc26133b"
      },
      "source": [
        "%matplotlib inline\n",
        "!python -m pip install -q numpyro\n",
        "!python -m pip install -q git+https://github.com/dfm/extending-jax.git"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\u001b[K     |████████████████████████████████| 184kB 4.5MB/s \n",
            "\u001b[K     |████████████████████████████████| 481kB 7.6MB/s \n",
            "\u001b[K     |████████████████████████████████| 32.1MB 150kB/s \n",
            "\u001b[?25h  Building wheel for jax (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Installing build dependencies ... \u001b[?25l\u001b[?25hdone\n",
            "  Getting requirements to build wheel ... \u001b[?25l\u001b[?25hdone\n",
            "    Preparing wheel metadata ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for kepler-jax (PEP 517) ... \u001b[?25l\u001b[?25hdone\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w02puAcQSyG-"
      },
      "source": [
        "Download some data (see [here](https://docs.exoplanet.codes/en/stable/tutorials/intro-to-pymc3/#A-more-realistic-example:-radial-velocity-exoplanets) for more info):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "jCWrLqr3AzMS",
        "outputId": "b9ddf8b1-1cf1-4dfa-da9e-aec11674970d"
      },
      "source": [
        "import requests\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Download the dataset from the Exoplanet Archive:\n",
        "url = \"https://exoplanetarchive.ipac.caltech.edu/data/ExoData/0113/0113357/data/UID_0113357_RVC_001.tbl\"\n",
        "r = requests.get(url)\n",
        "if r.status_code != requests.codes.ok:\n",
        "    r.raise_for_status()\n",
        "data = np.array(\n",
        "    [\n",
        "        l.split()\n",
        "        for l in r.text.splitlines()\n",
        "        if not l.startswith(\"\\\\\") and not l.startswith(\"|\")\n",
        "    ],\n",
        "    dtype=float,\n",
        ")\n",
        "t, rv, rv_err = data.T\n",
        "t -= np.mean(t)\n",
        "\n",
        "# Plot the observations \"folded\" on the published period:\n",
        "# Butler et al. (2006) https://arxiv.org/abs/astro-ph/0607493\n",
        "lit_period = 4.230785\n",
        "plt.errorbar(\n",
        "    (t % lit_period) / lit_period, rv, yerr=rv_err, fmt=\".k\", capsize=0\n",
        ")\n",
        "plt.xlim(0, 1)\n",
        "plt.ylim(-110, 110)\n",
        "plt.annotate(\n",
        "    \"period = {0:.6f} days\".format(lit_period),\n",
        "    xy=(1, 0),\n",
        "    xycoords=\"axes fraction\",\n",
        "    xytext=(-5, 5),\n",
        "    textcoords=\"offset points\",\n",
        "    ha=\"right\",\n",
        "    va=\"bottom\",\n",
        "    fontsize=12,\n",
        ")\n",
        "plt.ylabel(\"radial velocity [m/s]\")\n",
        "_ = plt.xlabel(\"phase\")"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de3hU9bX3v2sykwsDDMglJAIhIEREQmKjNUp9Q/G8rfVU7bGHaqFQTnuo11aP72vLsRfbnoqtWmmltaX1BtJa3qLVKvXGMVo80TbKJAURCUTCJQkYYCCBTGYn6/1jZm9+e2fPZCeZWybr8zzzZGbfZmXPzF573YmZIQiCIAjxwJVqAQRBEITMQZSKIAiCEDdEqQiCIAhxQ5SKIAiCEDdEqQiCIAhxw51qARLJ+PHjedq0aakWQxAEYUjxzjvvfMTMEwayb0YrlWnTpqG2tjbVYgiCIAwpiGjfQPcV95cgCIIQN0SpCIIgCHFDlIogCIIQN0SpCIIgCHFDlIogCIIQN1KqVIjoUSI6TETblWVnEdErRLQ78ndsZDkR0c+JqIGI6onogtRJLgiCINiRakvlcQCftiz7FoAtzDwTwJbIawC4AsDMyGMFgIeTJKMgCILgkJQqFWZ+A8BRy+KrATwRef4EgGuU5es4zFsAxhBRQXIkFQRBEJyQakvFjnxmbo48bwGQH3l+NoD9ynYHIstMENEKIqolotojR44kVlJBEATBRDoqFQMOTxDr1xQxZl7LzBXMXDFhwoC6DAiCIAgDJB2VSqvu1or8PRxZfhDAFGW7yZFlgiAIQpqQjkrlOQDLIs+XAXhWWb40kgV2MYCA4iYTBEEQ0oCUNpQkot8DqAIwnogOAPgegHsBbCSirwDYB2BRZPPNAD4DoAHAKQDLky6wIAiCEJOUKhVmvj7KqoU22zKAmxMrkSAIgjAY0tH9JQiCIAxRRKkIgiAIcUOUiiAIghA3RKkIgiAIcUOUiiAIghA3RKkIgiAIcUOUiiAIghA3RKkIgiAIcUOUiiAIghA3RKkIgiAIcUOUiiAIghA3RKkIgiAIcUOUiiAIghA3RKkIgiAIcSOjlcquXbtQVVWVajEEQRCGDRmtVFJNVVWVKDVBEIYVKR3SFQ0iKgHwB2XRdADfBTAGwL8DOBJZ/p/MvDnJ4gmCIAhRSEulwsy7AJQBABFlATgI4BmERwg/yMz3p1A8QRAEIQpDwf21EMAeZt6XakEEQRCE2KSlpWLhOgC/V17fQkRLAdQCuIOZj6kbE9EKACsAICcnJ2lC2uH3+1P6/oIgCMkmrS0VIsoGcBWA/xdZ9DCAGQi7xpoBPGDdh5nXMnMFM1d4PJ6kyRoLCdgLgjBcSGulAuAKAO8ycysAMHMrM3czcw+A3wC4KKXS9YGmaejs7EQgEEi1KIIgCEkh3ZXK9VBcX0RUoKz7HIDtSZfIITU1Nejo6EAwGER9ff2gFYtYO4IgDAXSVqkQkRfAPwF4Wln8EyL6BxHVA1gA4PZYx+ju7kZTUxNqamoSKKk91dXVxvOenp64WyuiZARBSEfSNlDPzB0AxlmWfak/xzh9+jQaGxuxcOFCbNmyBZWVlXGVMRbqBd/lcsHn8yX8vVRFJgiCkArS1lKJJ11dXUm/4FZWViI3NxdZWVkoLCxEIBBIicUkCIKQTIaFUsnOzk66q6impgadnZ3o7u7GgQMHDItJFIsgCJlMRiuVvLw8FBcXJ931Bdi7ouwspnjGRiTOIghCqknbmEo8yMrKwtSpU5OuUADYXtx1iykeMZBAIGC41PQiy7KysgEfTxCE4UMi47AZrVRSSWVlJbxeLzRNw9SpU6FpGjZs2BAXBVdTU4P6+nr09PRg4cKFcLlccLvloxQEoX8kQrnIlciC3Unu74nXt3e73XC73SgsLASAuFlM1dXV6OnpARB2qenvIwiCkGoyOqaSTJIZz6iqqoLLFf7osrOzRaEIgpA2yNWoHzixWKqqquD3+1FWVtbvGEes4+sxlPLycvh8PpSWliIQCGDDhg34p3/6J6MdTCLrYQRByAz064nP54v7NUMsFQuBQCBlVfhAuLOx1eLRYyiNjY1Gyxefz4epU6cCQFzbwQiCkNmo15O6ujp88MEHcb3eDXulorqt1JOt1pToisbpBVvTNDQ1NWHVqlVxCYCpMRRryxdrO5h9+/ZJLYwgCFFd8ur1hJnR3Nwc1xo6cX8pWAPgixcvhs/nMzKtXC4X8vLyUFVVZass/H4/NE1DR0eHqT1MtG2rqqpMqcHR0GMougy6uer3+3spumPHjqWkLY0gCEMD9Xqio9fQxeOakdGWSklJSZ+WgurusgbAfT4fAoGAyUroq5W9pmnG877awwQCAZNlFO24lZWVKC0tRXFxMUpLS00+UJ/Ph9zcXBCR4/cVBGH4ol9PCgoKjOtGPLuOZLRS6Yvy8nL4/X7jog7AuHhv2bLFCGLpigYIdz6ur6/vZVkEAgF0dnaaLu59fVCqwurq6kJra2tUpaXGUHTXWllZmfG+zOz4fQVBGN74fD7MmjUL8+bNM3UdiUcW67BWKq2trcZz3d3V2NhoVOHrrqkZM2Zg7NixxrY9PT1Yt26d8VqPxQSDQXR2diI3NzdmexhdEbjdbkNhud1utLS0IBgMoq6uDmvXrrXdz2rZWBVQVlaWuL4EQXCUdKTfrMbzejFslUpNTQ1aWlqM12632+RWUoP2e/bswfjx401WyGOPPWZ8WGosBggrKJ/PZ/qg9DsAVQHt2bMHM2bMQHFxMZYvX25YG8yMW265xdYaUi0buxRil8uFm266SSwVQRiGWK8zqWhkO2yVSnV1tclltHz5cvh8PpSVlaG6urpXxlVTUxMmTZpkbK9pmhG3qKqqMimcnp4e1NXV2X6Q1uO2t7dj6tSpWLp0qekY3d3dveIiqitOj/n4fD54vV54PB4AQCgUktRiQRjmWJOOrNcS/RqXCNJWqRDRh5Epj34iqo0sO4uIXiGi3ZG/Y/s6TjTUoLzL5cLSpUujBu2BsCWTn59vvFbjFpWVlSaFA4StDbsPzaqAWlpaEAgEUFlZiXPOOcdYnpOT08va0Ise1ZiPLpsqayImTQqCMHSwJh1F81wkQrmkrVKJsICZy5i5IvL6WwC2MPNMAFsirweEmlG1detWADCZi8CZoH1ZWRkqKioMq8AuXqIqHAAgol4fpN/vx8qVK00KiJkNBVBYWBj1+DrRfKBqq5ZET5oUBCG9Ua9vTmKsaoBe73o+UNJdqVi5GsATkedPALhmMAdTL9B25qI14yoQCMDtdtte1FU3VEFBAd58882oH6SqgKwKINrxY6ErPa/Xi5ycnF5px4IgDA9Ub8tAgvCBQACnTp0CgKkDlSGdlQoDeJmI3iGiFZFl+czcHHneAiDfuhMRrSCiWiKqPXLkSMw3UE2/aOaimnFVX1+P4uLiqOai2+3GiBEjMGvWrF4fpJ7xpQfXdYtEVQB+vx/t7e29jqv3E7NirfR3u93Izc0FgJS2mhEEIflYg/P9dYEHAgHU1dUhFAoBwISByhG1op6IznKwfw8zHx/om/fBfGY+SEQTAbxCRO+rK5mZiYitOzHzWgBrAaCioqLX+mjo5qLepFG3XlatWmVc0IkIixYt6vc/os4/qa+vR2lpqWGR9BddoanHdLlcKC0tNepW1Fkrkl4sCMMDq7elPw1mA4EA9u3bZ0peGiix2rQcijwoxjZZGISZFAtmPhj5e5iIngFwEYBWIipg5mYiKgBwOJ7vqWdTWVOB9ZYGfQW8olkUixcvdtS7S9M0uN1urFq1qtcx9I7HeosYu35g1g4A8Wy9IAhCemO9VjkdCqjeoMaDWO6vncw8nZmLoz0AtMVFCgtE5CWiUfpzAP8bwHYAzwFYFtlsGYBnE/H+Kv0NeNmh3i1EC6LX1NTg9OnTCAaDjvLKVXddXl4eNmzYgOrqamzYsMHYRirrBWHoMNhq9oFeq6x1dqNGjQKA2LGDGMSyVJxIlKhb4HwAz0RSb90AfsfMLxLR3wFsJKKvANgHoP++qBhEi5XYWTDR9rf7UugxFE3TMHv2bFulolozTiwMO3edvlx/L3F9CcLQpL/TZnWcXqus76U2rJ0xYwb8fn9Tv95YIapSYeZOACCiGQAOMHOQiKoAlAJYx8zH9W3iDTPvBTDPZnkbgIWJeM9Eo6f8WoNn+pdHVTRWC0MdqNPY2Gi4wKJ9gfTxwqJQBCEzUJWMncKJh4UTr6FdTrK/NgHoJqJzEA6ATwHwu0G96xAjHgVCakt8uywyNSOspKQEK1euBNB7QJfaBTnWe3V2dhpFnOICE4TMYyC/7Wj76OnH8ShFcDJPpYeZNSL6HICHmPkhIto26HfOUKIpH1UZRKt4t8sIswbku7q6YqYL19TUoKOjAwCwcOFClJSUSM2KIAxRoikNPSEo2sjyMWPGAACOH09Ucm50nFgqISK6HuHA+PORZZ7EiZR5qBlZwJnAur5OrzXR+46pWNvFhEKhmHno6v56WqEgCEOL/k6b1dG9Kpqm4cSJEygvLx/Q/oPBiaWyHMANAH7EzI1EVAxg/aDeNckMNOgVD9R0PSLCpEmTsGnTJlRWVtqus1ogur+ztbUVzc3NxvJgMGibh25NKxQrRRDSn6qqKtTW1mLixIlYu3Ztrxo0AKYJsZqmQdM0Q+lYJ8jq3gp99pOTJKN4Eav4cS2AvwB4lZm/ri9n5kYAP46bBBmOdR50bm6u8QFb1+mzoq0uK73+RCUrK8tWYViDbrt378bEiRMdfbEEQUgNgUDAiLnefPPNJpe3ekM5f/58ZGdno7MznCNVX19vUkILFy7EsmXLjOPqs5/0376qfBJ1PYjl/noE4QyszUS0hYi+SUS9MrKE2Fi7IauKwOraAqK7rKwTKNesWRPVCtGDbh0dHcYXNdkzFQRBcI76m9c9FzofffSRaV2kjYrxetOmTaZyBCv67KdkzViJqlSY+W1mvpuZP4FwPUgTgDsiregfJaK41ohkKrHmy0ebFW2nLPS290QEr9eLFStW9NpGJRAIoKGhwXgdDAZlbr0gpCnqbz4nJ8d4eL1eZGdnG+tcLheysrJMr6+99lpT30JrHCUUCvXqwmE3YyVeOImp6PUhv488QEQfA/DphEiUgcTK/dbX5efnG4WMejqx3bYul8uoeYn1pQgEAqY+PllZWZJaLAhphvqb9Hg8GD9+PDZt2oQrrrgCAHDy5EnDbc7MmDp1Knbv3m3sM3bsWMydO9dUCG29LjCz8T5OWk4Nlj6VChGNAbAUwDR1ezXOIgwetZDR6vdUOymPHDkyahqhTnV1NWpqajB//nzDlF6zZo3EVAQhDVGbwLa2thrLNU0zLIvOzk54vd5eXczb2tqMOKza5l5XHgBwzjnnGMvtunDEGyeWymYAbwH4B4D4dBwTDFSFAZizxey6DGuaZtSp9NXGZevWrVi8eDF8Pl+f7jJBEFKDWnLQ09ODBQsWGM1lVU6fPm1kdalYOxLHqpAfSBuX/uJEqeQy838kTIJhhBMfpp3fU83cUAPvffX2qqysHFB7fUEQkofu1tZ/926321AqwWDQ2C5aF2G7jsTxaLcyUJwolfVE9O8IFz4a/yEzH02YVBmGU2UChC2VaH7PRYsWGZW00tZeEDIDPQln586dRt++kSNHori4OOpoX73xo6ZpMV1Zuiu9vLw8aUrGiVLpAnAfgLsQnsaIyN/piRIq3iQjNzteROs+DDif7aIiGV+CkJ6oBY8+nw8XX3wxgHALFk3TsG/fvl77EBHcbjfmzJljcnfZocZq9CLKbdsS32HLiVK5A8A5zPxRn1umIX3FKNKRaH7PWApHEIShherOVivnT506hVAoZIqf6B038vPz0djY6MjqsMZqktWyyUnvrwYApxItSKJIVm52stALG0WhCMLQxlrw2Nraivr6elNxo868efMwa9Ysxy4sfWBftMLrROJEqXQA8BPRr4no5/ojUQIR0RQieo2I3iOiHUT0jcjyu4noYKT40k9En3FyPLVqXSYhCoKQKPrbil69yOfl5eHqq6+2Dcar1fX9QS28njFjhqk3WCJx4v76U+SRLDQAdzDzu5GRwu8Q0SuRdQ8y8/39OVimuYzSwdJKZYNOQUh3nP4+9BlKEydONLqW/+pXv+q1HTOjvr7eiIkMRHElMwTQp1Jh5icS9u7279cMoDny/CQR7QRw9mCOmYzc7ExElIcgJI5AIABN00zXJo/HY+v+GkxMRI2tJCNrNKr7K9KlOCZOthkMRDQNQDmAtyOLbiGi+kjvsbFR9llBRLVEVHvkyJFEiicIgjAg9ASiYDAIv99v9OtS+3oRkeH6Umcw6X28nMw+scZWkhECiGWpXENEsWbQE4AFcZbnzMGJRiI8yvg2Zj5BRA8D+CHC6cw/BPAAgH+z7sfMaxEee4yKigq2rhcSi1g3gtA3agIRcGZWit7SHgi3V/F6vYN23Sc7BBBLqfxfB/v/NV6CqBCRB2GFsoGZnwYAZm5V1v8GZ6ZQCjYk8uI+lOp+BCFZ6L8LJ9Xsas2Znpll/a3qrrF4uO6TGQKIqlSSHUvRobC99wiAncz8U2V5QSTeAgCfA7A9FfIlg3S5y7dTHmrdz/z585NWUCUI6Yz6uwCAgoKCmDdddv25VEUDIGVtVgaLk5TiZHMpgC8B+KQlffgnRPQPIqpH2O12e0qlzHCiDfRRzXY9eNjfVEpByDSs7ix9imtNTU3U34dec2ZtBJmTk4OysrK43qzFY/a8UxzNU0kmzLwV4XiNlc3JlmU4E62xpfrjSGZBlSCkkr7cyVYrA3BebK16BHw+H3Jzc4f076pPS4WI5iZDECG+BAIBo0X+QIhWNFpZWQmv14ucnJxekyzj8b6CMBSwWh+6leHxeIxlTjKt9P5cukcgWa1UEokT99cviehvRHQTEQ1d9TmMiMcsarUat6SkxDSN0u12295NJWsGtiCkG1VVVWhsbMSIESPg9XptfzdWqqursWjRIpNHYNGiRTh+/HjaxFUHgpPix08Q0UyE03ffIaK/AXiMmV/pY9e0YSh/QAMh1kyW/hAti0XTNGiaZgQZdfN93bp1SS2yEoTBEI8MSTtLxO12m+YYRcuWrKqqQiAQSMqI32TidEb9biL6NoBaAD8HUB7J0vpPPeVXSB8G0iLfKTU1NUb31Lq6Opx11lk4evQomBmPPvqosd1A31fqXIRMQm0/b9ciRZ+lEq2GZCj+DpzEVEqJ6EEAOwF8EsBnmXl25PmDCZZPGACq6yrefX7ULzkzo62tDczhGtPu7m54PB4QEUpKSsRKETIGp7FCTdPQ2dlpxEbsWqRYybTO404slYcA/BZhq+S0vpCZD0WsFyENiUexk/4DKC8vN8x3uywXnezsbGPdUM5eEQQVu5lMdoWOxcXFxnb19fWYMWMGOjs7QURg5oxxb/WFk0D9M8y8XlUoejt6Zl6fMMmEtMAafAeA0tJSFBQUmLbzeDzYsmULOjs70d3d3euuTmpZhHTDqfVhjVGuW7fO+E3U19fbWiU9PT1oaGhAc3O4XrugoKCX10B//0zI+FJxYqksBbDasuzLAH4Wd2mEtMMu6K/fnZ04cQKaphlztSsrK41tdSU0FCZtCsMPq/VRUlJi2yoF6B2jBGBSHvv27UNRURF8Pp+xndvthqZpAMJu4tzcXNPvQH3/vLw8bNmyJfH/dJKI1aX4eiL6M4BiInpOebwG4GjyRBQGQrwqaPsz5KympsaIrwBAMBgckoFGIfOx3izFshbUSveSkhIsXbrU+E0AwLFjx1BXV4fW1lbMmDEDxcXF+MUvfhFz6mKmTaRViWWp/A/Cc03GI9wRWOckgPpECiWkD3YdTqurq1FTU4NLLrkEQFh5eDwe/OQnPzHtS0SGEgoEAti5cyfKy8ulV5iQcqzWR18xQLXSfeXKlcjLy0N2djaOHTsGIGyNNDc3G7PmV6xYgYcffjhqVlciMzRTTVRLhZn3MXM1M1cy8+vK411m1pIppJBa7LJTrHdWoVAIf/7zn03LPvvZz6KystI0O6K+vl6KIoWUY82Q7G9iidvtRlFRkcliAczDtGJldSUyQzPVxHJ/bY38PUlEJ5THSSI6kTwRhXTE7s5KdX0REe68804AvZtQxjL1pc2LkCwGm8qr15gUFBQYw7T60w8v01KJdWK1vp8f+TsqeeIIQwW9B5heCAkAOTk5mDx5MpqamjB79mzjx2KdHRHN1LdL3cy0H5yQnvQ1I6i6uhrl5eV46623jMQU4Ezqfn5+vuNZKpmOk+LHi4lolPJ6FBF9PLFiCelEtKB/RUUFvF4vPB6PkTJZWFiIiy++2BQ3UQOdpaWlURVFJgcvhfTFSc861YXb0dFhZHbpWNvYA30nyySzHX0ycVKn8jCAduV1R2SZMMwJBALo6OhAKBRCa2trzG3VQKdas2J9nsxZ2oIAOLuZsc5Lyc7OxqpVq5Il4pDCSZ0KseIsZ+YeIkq7OSxC8lm0aBH8fj+AvmMlQLiFRVNTU1QXQbJnaQvDG/37WlNTEzMTy9r4EQinEav1LdZjDmecWCp7iejrROSJPL4BYG+iBYsGEX2aiHYRUQMRfStVcghmyyIvLy+mZbFq1SqcPn26VxWylUwNXgrpi5NMLD0on5WVZSzTW9WLIjHjRKncAOASAAcjj48DWJFIoaJBRFkAfgHgCgDnAbieiM5LhSyC/Y8xmp/YbgwxINleQnrg5GbG5/P1ewjXcMTJPJXDAK5LgixOuAhAAzPvBQAiegrA1QDeS6lUwxinjSutGWA+n88226u/SKt8YTDo3x/djRuLQCCAzs5OAOGU+dWrV4tFbUOfSoWIJiPcqfjSyKK/AvgGMx9IpGBROBvAfuX1AYQtJwMiWoGIJaUOyhFSixovcbvdUYd6xUJVIHKHKCSKqqoq+P1+lJWVobq62kg31hUKcGbsg47c1JzBifvrMQDPASiMPP4cWZaWMPNaZq5g5ooJEyakWhxBQbdq9uzZg8bGRjz66KNG0Vh2djY2btwY9Y5R/6FbnwvCQHDSNdvv96O8vNxIN25paTHWxaq3Gu44USoTmPkxZtYij8cBpOpqfRDAFOX15MgyIUX0N9debQ/e3d2NSZMm9WqV4ff75QcrJJWysjKUlZX1Wq5+X10uFwoKCvqstxruOFEqbUS0hIiyIo8lANr63Csx/B3ATCIqJqJshGM9z6VIFmEA6O3BgbB1kp+fL9leQtqhT3B0u929vq96vZVgjxOl8m8AFgFoQbhr8ecBLE+kUNGINLK8BcBLCI833sjMO1IhizAw9NTMaI38ot0xOkWGgQn9Qc8+PHTokJGFeOjQIXR0dCAYDGLPnj1GO3v9+6rHWgR7nGR/7QNwVRJkcQQzbwawOdVyCP1H/yFWVVUNaNSxWjypowdRy8vL5e5RsMUuQzAQCKC1tRUtLS2mRqhVVVXo6uoyXvf09KC9vR2zZs0Sa9ohUZUKET0EgKOtZ+avJ0QiYdjQn3Rg/e6xsbERRASXy4Xt27fj6NGjYGZjjoUoFqEv1FR2K6FQqNeylpYW5OfnJ0O0jCCWpVKbNCkEAWesDn0Uq94xtqamBg0NDcZ2zIzu7m5TSqdeUClKRegLax8vFY/Hg1AoZLJemNko1hW3V9/Ean3/hPqaiEYw86nEiyQMR+zuHvWCyOrqatOP3I7+zLEQhjdqIS4QViTTpk2DpmnYsGEDbrrpJuzYscOwWvLy8rBhw4ZUijykcNL6vpKI3gPwfuT1PCL6ZcIlE4YVdnePekGk2mPMjoKCAnF9CVGxtgLSxwEXFxfD6/VixIgRKCwsNLIQfT4fLrnkEpSVlWXkZMZE46Tb8GoAn0IkdZeZ64josoRKJWQ0asBex3r3CJzpraRX4+/cudMYjtTV1YXx48cjPz/fVN8CYFDZY0LmoHcXtmsF5Ha7cfToUZw+fRoej8fWdeq0BZFgxlELe2ber1c+R+hOjDjCcMI6bU9t46K7IvQftD6PRcftdmPWrFmmY3V2doKIjLtSuRgIavGibvmq35Wenh4Eg0HU19ejtLQ0xdJmBk6Uyn4iugQAE5EHwDcQrhERhAFj10zSOmfFqVJQ70YBGBP8xG0h6MW2+qyUcePG2WZ+EREWLVoEwN6SFpzjtPX9zQg3czwIoCzyWhAGzEBHB+uVzuo4V/VuVEfGEWce/S1s1S1htXixra3NNvNL2tjHD6eTHxcnXBJhWKHGUPQf9EsvvRRzn+LiYuMuMxgMGn5w9W5Ux3qRkBb56UmiPhfVEtZrmHSr1fpdidbGXr4rA8OJpfImEb1MRF8hojEJl0gYFjiZtqdSXV2NRYsWGReDrKwsw13h8/kwY8YMY1uZdZE59GWdRFsfbSic/r0bO3assa21jb0wOPpUKsw8C8C3AcwB8C4RPR9pKikIg8I6ba+vjsdqarFuiej7qO4wuUgI6nfF5XLh8OHDhgJqbGxEUVGRab24vuKHE0sFzPw3Zv4PhCcvHgXwRB+7CMKAiaZcYlk3aoBfLhKZibXepKqqCrW1tWhqakJ5ebnpM1e/K1u3bkVFRYXpWGpj061bt4pVG0ecFD+OJqJlRPQXAP+DcKfiixIumSDYYDdLXL/79Hq9A5p1IZ2N0xN1GJseI9Ez+2pqahAIBIx+cPX19YaLSyfa3Hm9y7CTufRC/3ESqK8D8CcAP2DmmgTLIwwjVGskHgFbt9sNt9vtqLJeAvdDg/b2dvj9fttsQVWJ9PT0YN++fVHrk/RMMDVtXT77xODE/TWdmW8XhSIIQqpQYyTMjI0bNxrdFXSOHTtmWDHAGaWhjgS2s2iE+OIkUB+7k58gpBhrexYZRzx08Pv9xucXC7Vfl175vmfPnl7bqf3i9O+AWsekZoIJicFRoD5ZENF9RPQ+EdUT0TN6CjMRTSOi00Tkjzx+lWpZhfTh5MmTOHHiBFatWoXjx49L768hiF1cS1cGgUAAtbW1Rr+uxsZG24JXIOwCtR5HHWFNROjs7DSsGSH+pJVSAfAKgPOZuRTABwBWKuv2MHNZ5HFDasQTUo01M6ympgY9PT1gZpProy+smURCelFTU4O6ujpj5EFHR4eRNg0AK6MAACAASURBVK53VbD0IwQALF++vFdMRc/0KigoAAA0Nzf367si9I+0mvzIzC8rL98C8Pl4v4eQWagKJhgMGgFctVGlFbXaev78+cjLy+uVciokB03ToGlary7BdjN0Tp8+jezsbHR2dqKjowNEhHHjxpmmfy5dutS0jxqgz83NNY6pu8kk8yv+pPPkx38D8AfldTERbQNwAsC3mfmvdjsR0QoAKwBg6tSpCRdSiA99KQKgd4O/6upqjBs3znjd09OD48eP92pUWVlZaTq+tdpaLZwUkkdNTQ06OjoAoFeXYLuYWE9PDzo7O43XzIzRo0djypQpCAQCpq7W1hn0LpcLM2bM6NUaSIg/jic/xgsiehXAJJtVdzHzs5Ft7gKgAdDHrTUDmMrMbUT0MQB/IqI5zHzCRu61ANYCQEVFhSQZDAHsOhY7vYNUK+ddLhd+85vf2DaqVC0T9eLicrmMLCJJM3ZOPM6Vuq81JbiyshIzZ87E7t27o+6vT/u0zj2xmyJKRFi+fDk2btzYSwEJ8aXPOhUimgDgmwDOA2AMtGDmTw7kDZn58j7e78sA/hnAQj3zjJmDAIKR5+8Q0R4As5B6a0qIA3Y1CLF+8KrVYbVUvF4vjh07BsDcysVqmeizW3w+H3bv3o2mpqZerfeFxGK1FPSU4JKSEvh8PhQWFqKhocF2lLTH48GcOXNsPy+7KaJq01IZvJVYnBQ/bkDYDXUlwm3wlwE4kghhiOjTAO4E8L+Y+ZSyfAKAo8zcTUTTAcwEsDcRMgjJx65jcTSsk/yWLVtmrHO5XPB4PPB6vZg4caLpblS1TNS7W7UqW+9mK6SOrq4u2ymMVkKhEH75y1/aKgf1+0REmDRpEjZt2iSKJEk4USrjmPkRIvoGM78O4HUi+nuC5FkDIAfAK5HMjrcimV6XAfgBEYUA9AC4gZmPJkgGIcmoUx9juSUCgQD27dtnsmpUcnJy4PP5cPToUaP9hq6gVMtEvWBZq7KlhiF52LnOsrOz4Xa7DcsxJyfHFEdRWbdune13Rf0+yUjg5OMkpTgU+dtMRFcSUTmAsxIhDDOfw8xTrKnDzLyJmedEll3AzH9OxPsLqaOvPky6n1x3bQHhC9DSpUuNnl/6/PHOzs6ofaAAoKmpyVivVmXrVkwspE9YmHikZFvPY0FBAVavXo09e/agsbERu3btMrWo7w/65y3uzOTjxFL5LyLyAbgDwEMARgO4PaFSCYIFq5987NixeOGFF1BZWWlKB9ZdY7t27ep1wVNdZ3l5eVi9ejVuvPFGAOFA7owZM+Dz+SRg3weDSaxQqayshNfrRUdHBzweD/Lz802TGa2WqJXy8vIByS8klj6VCjM/H3kaALAgseIIgj2qn9zlcqGoqMh0IfP7/Vi8eHHMEcVqFXZXVxc2bdpkvGZmtLe3o6mpCYcPHzZVZqszy/WWInoCwHCkv4kVTgiFQqivrzeUPBC2RPPz89Ha2oqsrCxommYK2t92222YO3eu7XtH+2yG62eWTKK6v4jozsjfh4jo59ZH8kQUBPN8jNLSUlu3htqOwy7gb11/7bXXGq+BcKV1Y2OjqXrbKVVVVRgzZsywcI3ZDUsbKOp57unpQVtbm8md6fP5QEQYMWIE5s2bZ9rX7sZBSD2xYio7I39rAbxj8xCEpKL7ybdt22Z7MVEHL9m5ZKzrV6xYgdLSUuTk5MDj8Zi2lYLI6DgdBe0k/mSNadltP3LkSJSVlZluCoDBKzQhMcQqfvxz5K9MeRQSzmCL6PSLizXbxzpHw7peb9+haRpCoZBxzFAohLfffhs5OTm93k/TNCNIPVyziuKVVeV2uw135oYN4VpnvcreGq+pqqoCEZm2H67nP52J1fvrz4jd++uqhEgkCHFCDSjrNSjRsoHcbje8Xi86OzvR3d1ttAQJhUIm5REMBtHZ2YnGxkbMnz9/2I2ijbdlUFZWBr/fb2T+TZ8+3VhnjdeosS0Aw+q8DyViBervj/z9F4TbqjwZeX09gNZECiUIdsSyZqJN9rPO0YimVPTGhtZK7O7ubuOO+dChQ6aaiZ6eHmlKOABiZdepn4+4t4YmsdxfrwMAET3AzGoL1z8TkbRHEdIG1SLJy8sz6lWsGWM+n8/2QqZpmuFysaOrqwvr1q1DQ0ODaTkRmQZBnTp1Ch988EGfbrFMTlnWM+SczrRRXZc6Ho8H48ePN6rgM/l8ZSJOih+9kdYoAAAiKgbgTZxIgtA/7FJcAWcZY9XV1Zg4cWLM42dlZQFArx5U55xzDiorK43ZH6FQCM3NzViwYMGwn9Xh9/tRW1trWyBpLZzUZ+ToNwehUAitreIMGao4USq3A6gmomoieh3AawBuS6xYguCcWCmufWWMAcCGDRtMWUVWiAjl5eW9Mo80TTNa6asKJ9Wpromq+lfrdPp6P936a2xsNA3E0hWHdTlgvjkgIrFMhihOZtS/iHADx28A+DqAEmZ+KdGCCYJTnKa4Otl/5syZICJTirGmaWhrazPSj2fOnAlN04wL47hx40xTCJ3EAvx+/5Coa7EqDD3zza5Hmqp01JRsVclGsyr1/Z3Uv1infwrphdNxwjMBlACYB+ALRLS0j+0FIan01TvM6f6FhYVwuVzIzs421mVnZ2Pjxo1obGw00o/1C+Pp06dx7733Yt68efB4PCgoKMDs2bOxcuXKaG+VcgZaqKl2dPb7/Ub8yI5oSlZVHMxskiHazYGMfh5a9KlUiOh7CPf8egjhNi0/ASDpxELGMnLkSFRUVMDlcoGIjMpuHbUIT22lP2LECMyaNctRE0N9znoiuiIn6iJslbW5uRmXXHIJamvP5O1omob29nYjS46IsHr1akNB6IojJycHpaWltgWq6s1BLHeZkJ44sVQ+D2AhgBZmXo6wtSKtP4WMpaysLKZ7Ra3Mj1X7Eg39jj8YDKK+vr7XhXIwMRH1Ijx//vy4Nl2M9n92dXWhqakJhw4dQkdHB7q7u411zGyazgkA27Ztw+zZs41Ba7GI5S4T0hMnXYpPM3MPEWlENBrAYQBTEiyXIPSLwV5srPvX1NQYF7OFCxdi8uTJ6OzsxOzZs4002IG2VVczm3p6eqLOBXGKmnJrV5tj3U7H7/f32RhTnbLp8/ng9XrR1dXVq/tAY2Oj7f5q2rVOrC7HVlnUtHCpWxkaOLFUaoloDIDfINzz610ACbNBiehuIjpIRP7I4zPKupVE1EBEu4joU4mSQcgcBhrUVfcJBoNoaGhAMBjs1VJfdzWtWrUKx48fR3V1dUz3U01NDVpaWkzLfv3rX6O8vNxU8zJQ95V60bWbD+P3+20zuOyOU15ebnI9HTp0CJqmITs722j6aO2ZZkVPu1bpj/Ux2CQMIfnEtFQoHG1bxczHAfyKiF4EMJqZ6xMs14PMfL+6gIjOA3AdgDkACgG8SkSzmLnb7gCCMBjUO+SsrCzjzly9CFpHG+tFl7FmjVjTj4Gwi0iv9h/srBJ9RommaZg9e7atNaVpmjGrJFZMRx0VoCtWZkYwGITX64Xb7Y4588TlcqGwsLDX8v5aHzK9cWgR01Lh8Ld/s/L6wyQolGhcDeApZg4ycyOABgAXpUgWIcNR75DXrFljm+pqnc9idT/Z3YXbXUBVi8Ju/4HGWDo6OkwWj9/vx8mTJ9HR0YFQKIRQKAS/34+1a9f2ktHv9/dKSFCV4enTp43jAMDo0aNNGV96WvaqVat6ySXWR2bjxP31LhFdmHBJzNxCRPVE9CgR6fNEzwawX9nmQGSZCSJaQUS1RFR75MiRZMgqZCh6JpLeIl+9CFZXV5uKJnVl01ethW5JEBFyc3ONLChdqTit1YhGTU2NkQSwe/duU9aUXW8zALjllltsXW1qQkJxcbFpnfU4J0+exDnnnAMiAhEZFk20jK1t27Zh7969olAyECdK5eMAaohoT+RC/w8iGpS1QkSvEtF2m8fVAB4GMANAGYBmAA/059jMvJaZK5i5YsKECYMRUxAM7Opg7O64nd6FZ2dn49xzz0Vubi58Pp8RRwGAvLw8Y0jVypUr+4yBqDEYu/iE3rssWn+z7u5u2/38fj8aGxuN+h2v1xt1ZjwzQ9M0Q6Go7x2PJArJ+ho6OMn+intAnJkvd7IdEf0GgD7O+CDMWWeTI8sEIWXY+ftjxQB0SwIIx16ys7PxwQcfoKWlBcyMhQsXwuVyITc319g/1vwWawxm9erVvd5TLeS0Iycnx7bNiqZppiFabrcbRUVFhttPVSAulwtut7uXBSMZW8MPJ21a9tk9EiUQERUoLz8HYHvk+XMAriOinEhTy5kA/pYoOQRBvUOO192yegx9Zktzc7Nxce7q6sLEiRONLr9qFbudK8kag1HH8c6cOdOwmJYu7d0Ew+PxoLi4GCUlJaYOAGodTUdHhymYr7rE5s2bB6/Xi+LiYmzdurXXtMyxY8dKzGQY4rRNSzL5ieJiW4BwQ0sw8w4AGwG8B+BFADdL5pcw1FDv2tXAtk52djbcbrdhmagX9GiBf2sMxu12Izc3F4WFhTFb12RnZ2Pq1Km9MsSsHYJ37tyJ2tpaowOA7gr0+Xxwu93Ge6gxJn06oyiU4YcT91dSYeYvxVj3IwA/SqI4ghCT/lovasrvz3/+c9x4442GK2nSpEm4++67jWV6s0qdaIH/0tJS42IfreeYVU6Px2NybenY1dFMmDABBw4cABB22UXrImCVZaDFocLQJu2UiiAMdfpSNG63G263GytWrMDDDz9sugi3tbWZakMOHgyHDa09tFTUC7ieNqyn86rz3XVcLhdef/11WwVkV0ejywCYJ2jq6c7RZBGGJ+no/hKEjMGuxuT48eM4fvw4AJhcSfr2OllZWcYF3q6Hlo5aJa+nDVvTeSsrK1FWVmZq5KjHidTsMdWdpqMqmby8PFEaQkxEqQhCmuD3+7Fy5Uoj0L5mzRpjncvlcpRFFW2OCRBWYHoKs461C/BNN92EvLw8U/sVIoLb7UZBQUGvjs3REhgkDXj4Iu4vQUgCatNHp8ydOxe5ubkIhUL45S9/2WfQOxAImDoEO0nntWaPBQIBwz2naRpGjx6NlpYWaJpmBPBj/Q+iSARRKoKQQPQuv9FiDepFuKqqykjnBYAFCxYgGAwCAG677TbMnTvXVrGUlZUhEAigrq7O5KrSYzCqYrG29bf24fL5fDh69CgAGFlk+jF7enqwePFi7N27d8DnQ8h8xP0lCAlCdS3V19f3OZArEAhg374zJWBqs8ZYlen6ftYAe7QYjIq1A4BV8an9v/TXghALsVQEIUHEmm1iRa2M18nOzjYslWiuLLv9APs5JtGUkl0HAL34EoCkCQv9QiwVQUgQaiYVEaGzszPqjBRVAQHhzK/XXnsNZWVlMfuIWfdzuVzweDyYN2/egAsPdReZ/rBmqAlCLESpCEKC0F1LBQXhzkPNzc1Ru/ZaU3n13l92jSxj7ZeXl4cRI0aYFEB/hn7F2nYww8OE4YMoFUFIIHoar9rby84NpcY2cnNzoWmao4u3vl9OTo4xOEsNxltThmMdM9a2+kAyJ8cRhjeiVAQhwajB7lhpvnrMoquryyhe1IdcxUoN1hWXXdsVJ6N7VVeX3bbV1dVYtGiR4xHAwvBGAvWCkGD0zr6BQAAbNmyIGeuwmybpFGu6MNC/0b2xtu3vCGBh+CKWiiAkEGuwO5ZCiTZN0sl7FBcX28Y7+jO6N9a2MgJYcIpYKoKQRqidfvuyanSsg7qsF/1YQ8OsxNq2P8cRhi9iqQhCmuHEqlFxEjcRhGQhloogJIHBXOj1Vi9244SB5MU7RFkJTkgrpUJEfwBQEnk5BsBxZi4jomkAdgLYFVn3FjPfkHwJBSHxqBfvvlxbQN8uM1EGQjJJK6XCzF/QnxPRAwDUvhZ7mLms916CkLnYubYk3iGkM2mlVHQoPLx7EYBPploWQUglyU7lFatGGCzpGqj/BIBWZt6tLCsmom1E9DoRfSLajkS0gohqiaj2yJEjiZdUEBKIpPIKQ42kWypE9CqASTar7mLmZyPPrwfwe2VdM4CpzNxGRB8D8CcimsPMJ6wHYea1ANYCQEVFBVvXC8JQQ1xbwlAi6UqFmS+PtZ6I3AD+BcDHlH2CAIKR5+8Q0R4AswDUJlBUQRgyiNtKSBfS0f11OYD3mfmAvoCIJhBRVuT5dAAzAcj4OUEQhDQjHQP118Hs+gKAywD8gIhCAHoA3MDMR5MumSAIghCTtFMqzPxlm2WbAGxKvjSCIAhCf0g7pSIIghmJlwhDiXSMqQiCIAhDFFEqgiAIQtwQpSIIgiDEDVEqgiAIQtwQpSIIgiDEDVEqgiAIQtwQpSIIgiDEDVEqgiAIQtwQpSIIgiDEDVEqgiAIQtwQpSIIgiDEDVEqgiCYmDNnzoD7jRERGhoa4iuQgGnTpuHVV19NtRiOEKUiCIKJHTt2oKqqKtVi9GL37t3Izc3FkiVLom5z33334fzzz8eoUaNQXFyM++67z7R+wYIFmDBhAkaPHo158+bh2WefNa3/3e9+h6KiIni9XlxzzTU4evTMhI2RI0eaHllZWbj11luN9Rs3bsTs2bMxatQonHfeefjTn/5krHv88ceRlZVl2j9TG4WKUhEEAQCgaVqqRYjJzTffjAsvvDDmNsyMdevW4dixY3jxxRexZs0aPPXUU8b6n/3sZ2hubsaJEyewdu1aLFmyBM3NzQDCyvRrX/sa1q9fj9bWVowYMQI33XSTsW97e7vxaGlpQV5eHv71X/8VAHDw4EEsWbIEP/3pT3HixAncd999+OIXv4jDhw8b+1dWVpqOkY6KOx6IUhGEIc60adOwatUqnHfeeRg7diyWL1+Ozs5OY/3zzz+PsrIyjBkzBpdccgnq6+tN+/74xz9GaWkpvF4vNE0zuVqCwSBuu+02FBYWorCwELfddhuCwaCx/3333YeCggIUFhbi0UcfTdj/+NRTT2HMmDFYuHBhzO3uvPNOXHDBBXC73SgpKcHVV1+NN99801hfWloKtzs88YOIEAqFsH//fgDAhg0b8NnPfhaXXXYZRo4ciR/+8Id4+umncfLkyV7vs2nTJkycOBGf+MQnAAAHDhzAmDFjcMUVV4CIcOWVV8Lr9WLPnj0D+n/Xr1+PoqIijBs3Dj/60Y9M6/72t7+hsrISY8aMQUFBAW655RZ0dXUBCCveO+64w7T9VVddhQcffBAA8OMf/xhnn302Ro0ahZKSEmzZsmVA8sWEmZP+APCvAHYgPMWxwrJuJYAGALsAfEpZ/unIsgYA33LyPh/72MdYEDKdoqIinjNnDjc1NXFbWxtfcsklfNdddzEz87vvvssTJkzgt956izVN48cff5yLioq4s7PT2HfevHnc1NTEp06dMpa98sorzMz8ne98hz/+8Y9za2srHz58mCsrK/nb3/42MzP/5S9/4YkTJ/I//vEPbm9v5+uvv54B8O7du23lvPHGG9nn89k+5s6dG/X/CwQCPHPmTN6/fz9/73vf48WLFzs6Lz09PVxWVsYPP/ywafmVV17JOTk5DIA/9alPcXd3NzMzX3XVVXzvvfeatvV6vVxbW9vr2AsWLODvfe97xmtN0/iyyy7jZ599ljVN42eeeYbPPvtsbm9vZ2bmxx57jEeMGMHjxo3jmTNn8g9+8AMOhUK2cu/YsYO9Xi+//vrr3NnZybfffjtnZWUZn0ltbS3X1NRwKBTixsZGPvfcc/nBBx9kZua3336bCwoKjP/pyJEjnJeXxy0tLfz+++/z5MmT+eDBg8zM3NjYyA0NDbYyAKjlgV7fB7rjYB4AZgMoAVCtKhUA5wGoA5ADoBjAHgBZkcceANMBZEe2Oa+v9xGlIgwHioqKTBfOF154gadPn87MzDfccIOhBHRmzZrF1dXVxr6PPPJIr+PpF7Dp06fzCy+8YKx78cUXuaioiJmZly9fzt/85jeNdbt27YqpVAbK17/+deNi3x+l8t3vfpdLS0sNBarS1dXFmzdv5gceeMBY9slPfrKXAiosLOTXXnvNtOzDDz9kl8vFe/fuNS3/7W9/y16vl7OysjgvL4+ff/55Y92ePXt479693N3dzfX19Tx79my+5557bOX+/ve/z1/4wheM1+3t7ezxeIzPxMqDDz7I11xzjfH63HPP5ZdffpmZmR966CG+4oormJl59+7dPGHCBH7llVe4q6vL9lg6g1EqKXF/MfNOZt5ls+pqAE8xc5CZGxG2Si6KPBqYeS8zdwF4KrKtIAgApkyZYjwvKirCoUOHAAD79u3DAw88gDFjxhiP/fv3G+ut+1o5dOgQioqKbI996NChXu8bb/x+P1599VXcfvvt/dpvzZo1WLduHV544QXk5OT0Wu/xeHDFFVfg5ZdfxnPPPQcgHIg/ceKEabsTJ05g1KhRpmXr16/H/PnzUVxcbCx79dVXceedd6K6uhpdXV14/fXX8dWvfhV+vx8AMH36dBQXF8PlcmHu3Ln47ne/iz/+8Y+2slvPq9frxbhx44zXH3zwAf75n/8ZkyZNwujRo/Gf//mf+Oijj4z1y5Ytw5NPPgkAePLJJ/GlL30JAHDOOedg9erVuPvuuzFx4kRcd911pu9BvEi3mMrZAPYrrw9ElkVb3gsiWkFEtURUe+TIkYQJKgjphB4XAICmpiYUFhYCCCuMu+66C8ePHzcep06dwvXXX29sT0RRj1tYWIh9+/bZHrugoKDX+8bihhtu6JVBpT/mzJlju091dTU+/PBDTJ06FZMmTcL999+PTZs24YILLoj6Po8++ijuvfdebNmyBZMnT44pk6ZpRtxjzpw5qKurM9bt3bsXwWAQs2bNMu2zbt06LFu2zLTM7/fjsssuQ0VFBVwuFy688EJ8/OMfj5oGTES6d6YX1vN66tQptLW1Ga9vvPFGnHvuudi9ezdOnDiBe+65x3SsJUuW4Nlnn0VdXR127tyJa665xlj3xS9+EVu3bsW+fftARPjmN78Z8/wMiIGaOH09ALwKYLvN42plm2qY3V9rACxRXj8C4PORx2+V5V8CsKYvGcT9JQwHioqK+Pzzz+f9+/dzW1sbX3rppbxy5UpmZv773//OkydP5rfeeot7enq4vb2dn3/+eT5x4oSxr9Wtoi676667uLKykg8fPsxHjhzhSy+91IjXbN68mfPz83nHjh3c0dHBixcvjrv7q6Ojg5ubm43HHXfcwddeey0fPnzYdvsnn3yS8/Pz+b333uu1bufOnbx582Y+deoUd3V18fr169nj8fA777zDzMzbt2/nUaNG8RtvvMHt7e28ePFikxuKmfnNN9/kESNGGOdPp7q6mseNG8fbtm1j5nAs66yzzuKXXnqJmcPnqqWlxZBjzpw5fPfdd9v+D9u3b2ev18t//etfORgM8h133GGKqVx44YX8/e9/n3t6enjnzp08a9YsvvTSS03HuPzyy3nu3Lm8fPlyY9n777/PW7Zs4c7OTg4Gg7x8+XJeunSprQwYajEV4817K5WVAFYqr18CUBl5vBRtu2gPUSrCcKCoqIjvuecenj17Nvt8Pl66dCl3dHQY6//yl79wRUUF+3w+njRpEn/+8593rFROnz7Nt956K0+aNIknTZrEt956K58+fdrYdtWqVZyfn88FBQX8yCOPJCSmomKNqbzxxhvs9XqN19OmTWO3281er9d4fO1rX2Nm5vfee48vuugiHjlyJPt8Pq6oqOCnn37adPwNGzbwlClTeMSIEXzVVVdxW1ubaf2KFSt4yZIltrI99NBDPGPGDB45ciQXFxfz/fffb6y74447eOLEiTxixAguLi7m73znOzHjGo8//jhPmTKFzzrrLP6v//ov02fy+uuvc0lJCXu9Xp4/fz5/5zvf6aVU1q9fzwD4v//7v41ldXV1fOGFF/LIkSN57NixfOWVVxpBeyuDUSrEUUywZEBE1QD+DzPXRl7PAfA7hGMohQC2AJgJgAB8AGAhgIMA/g7gi8y8I9bxKyoquLa2NmHyC0I6MG3aNPz2t7/F5ZdfnmpRhDThjTfewJIlSww3V38honeYuWIg7+0eyE6DhYg+B+AhABMAvEBEfmb+FDPvIKKNAN4DoAG4mZm7I/vcgrDlkgXg0b4UiiAIwnAkFArhZz/7Gb761a8OSKEMlpQoFWZ+BsAzUdb9CMCPbJZvBrA5waIJgiAMWXbu3ImKigrMmzcPjz32WEpkSIlSEQQhfnz44YepFkFIE2bPno2Ojo6UypBuKcWCIAjCEEaUiiAIghA3RKkIgiAIcSOlKcWJhohOItyEUgDGA/ioz62GB3IuziDn4gxyLs5Qwsyj+t6sN5keqN810FzrTIOIauVchJFzcQY5F2eQc3EGIhpwgZ+4vwRBEIS4IUpFEARBiBuZrlTWplqANELOxRnkXJxBzsUZ5FycYcDnIqMD9YIgCEJyyXRLRRAEQUgiolQEQRCEuJERSoWIPk1Eu4iogYi+ZbM+h4j+EFn/NhFNS76UycHBufgPInqPiOqJaAsRxX8GbJrQ17lQtruWiJiIMjad1Mm5IKJFke/GDiL6XbJlTBYOfiNTieg1ItoW+Z18JhVyJhoiepSIDhPR9ijriYh+HjlP9UQUfdymykAHsaTLA+FW+HsATAeQDaAOwHmWbW4C8KvI8+sA/CHVcqfwXCwAMCLy/MbhfC4i240C8AaAt6AMjMukh8PvxUwA2wCMjbyemGq5U3gu1gK4MfL8PAAfplruBJ2LywBcAGB7lPWfAfAXhOdZXQzgbSfHzQRL5SIADcy8l5m7ADwF4GrLNlcDeCLy/I8AFlIqBg0knj7PBTO/xsynIi/fAhB7iPfQxcn3AgB+CODHADqTKVyScXIu/h3AL5j5GAAw8+Eky5gsnJwLBjA68twH4FAS5UsazPwGgKMxNrkawDoO8xaAMURU0NdxM0GpnA1gv/L6QGSZ7TbMrAEIABiXFOmSi5NzofIVhO9ERZOsqAAAA/ZJREFUMpE+z0XEnJ/CzC8kU7AU4OR7MQvALCJ6k4jeIqJPJ0265OLkXNwNYAkRHUB4htOtyREt7ejv9QRA5rdpEaJAREsAVAD4X6mWJRUQkQvATwF8OcWipAtuhF1gVQhbr28Q0VxmPp5SqVLD9QAeZ+YHiKgSwHoiOp+Ze1It2FAgEyyVgwCmKK8nR5bZbkNEboRN2rakSJdcnJwLENHlAO4CcBUzB5MkW7Lp61yMAnA+gGoi+hBhn/FzGRqsd/K9OADgOWYOMXMjgA8QVjKZhpNz8RUAGwGAmWsA5CLcbHK44eh6YiUTlMrfAcwkomIiykY4EP+cZZvnACyLPP88gP/mSCQqw+jzXBBROYBfI6xQMtVvDvRxLpg5wMzjmXkaM09DOL50FTMPuJFeGuPkN/InhK0UENF4hN1he5MpZJJwci6aACwEACKajbBSOZJUKdOD5wAsjWSBXQwgwMzNfe005N1fzKwR0S0AXkI4s+NRZt5BRD8AUMvMzwF4BGETtgHhwNR1qZM4cTg8F/cBGAng/0VyFZqY+aqUCZ0gHJ6LYYHDc/ESgP9NRO8B6Abwf5k546x5h+fiDgC/IaLbEQ7afzkTb0KJ6PcI30iMj8SPvgfAAwDM/CuE40mfAdAA4BSA5Y6Om4HnShAEQUgRmeD+EgRBENIEUSqCIAhC3BClIgiCIMQNUSqCIAhC3BClIgiCIMQNUSqCkCCI6MNIzYcgDBtEqQiCIAhxQ5SKIAwSIppGRO8T0QYi2klEfySiEZHVtxLRu0T0DyI6N7L9RURUE5nX8T9EVBJZPoeI/kZE/sj8ipmR5UuU5b8moqwU/auC0CeiVAQhPpQA+CUzzwZwAuEZPgDwETNfAOBhAP8nsux9AJ9g5nIA3wVwT2T5DQB+xsxlCDf7PBBpE/IFAJdGlncDWJyMf0gQBsKQb9MiCGnCfmZ+M/L8SQBfjzx/OvL3HQD/EnnuA/BExBJhRFpjAKgBcBcRTQbwNDPvJqKFAD4G4O+Rtjp5ADK5Z5swxBGlIgjxwdrvSH+td4Huxpnf2w8BvMbMn6PwaOtqAGDm3xHR2wCuBLCZiL6G8NS9J5h5ZeJEF4T4Ie4vQYgPUyOzNwDgiwC2xtjWhzMtxL+sLySi6QD2MvPPATwLoBTAFgCfJ6KJkW3OIqKiOMsuCHFDlIogxIddAG4mop0AxiIcQ4nGTwCsIqJtMHsLFgHYTkR+hGe9rGPm9wB8G8DLRFQP4BUAfY50FYRUIV2KBWGQRFxYzzPz+SkWRRBSjlgqgiAIQtwQS0UQBEGIG2KpCIIgCHFDlIogCIIQN0SpCIIgCHFDlIogCIIQN0SpCIIgCHHj/wOCjk0CqC82NwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nasxBsFpS63I"
      },
      "source": [
        "Fit these data using numpyro:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 334
        },
        "id": "4Vn6SYstAUdV",
        "outputId": "820bf5e5-dd46-480a-f5b0-26c2cb5bf2c7"
      },
      "source": [
        "from jax.config import config\n",
        "\n",
        "config.update(\"jax_enable_x64\", True)\n",
        "\n",
        "from jax import random\n",
        "import jax.numpy as jnp\n",
        "\n",
        "import numpyro\n",
        "import numpyro.distributions as dist\n",
        "from numpyro.infer import MCMC, NUTS, init_to_value\n",
        "\n",
        "from kepler_jax import kepler\n",
        "\n",
        "\n",
        "def model(t, rv_err, rv=None):\n",
        "    # Parameters\n",
        "    K = numpyro.sample(\"K\", dist.Uniform(10.0, 100.0))\n",
        "    P = numpyro.sample(\"P\", dist.LogNormal(np.log(4.23), 5.0))\n",
        "    ecc = numpyro.sample(\"ecc\", dist.Uniform(0.0, 1.0))\n",
        "\n",
        "    # Handle wrapping of angles appropriately\n",
        "    phi_angle = numpyro.sample(\"phi_angle\", dist.Normal(0.0, 1.0), sample_shape=(2,))\n",
        "    phi = numpyro.deterministic(\"phi\", jnp.arctan2(phi_angle[0], phi_angle[1]))\n",
        "    w_angle = numpyro.sample(\"w_angle\", dist.Normal(0.0, 1.0), sample_shape=(2,))\n",
        "    norm = jnp.sqrt(jnp.sum(w_angle ** 2))\n",
        "    sinw = w_angle[0] / norm\n",
        "    cosw = w_angle[1] / norm\n",
        "    omega = numpyro.deterministic(\"omega\", jnp.arctan2(sinw, cosw))\n",
        "\n",
        "    # RV trend parameters\n",
        "    rv0 = numpyro.sample(\"rv0\", dist.Normal(0.0, 10.0))\n",
        "    rv_trend = numpyro.sample(\"rv_trend\", dist.Normal(0.0, 10.0))\n",
        "\n",
        "    # Deterministic transformations\n",
        "    bkg = numpyro.deterministic(\"bkg\", rv0 + rv_trend * t / 365.25)\n",
        "    mean_anom = 2 * np.pi * t / P - (phi + omega)\n",
        "\n",
        "    # Solve Kepler's equation\n",
        "    sinE, cosE = kepler(mean_anom, ecc)\n",
        "\n",
        "    # MAGIC: convert to true anomaly\n",
        "    tanf2 = jnp.sqrt((1 + ecc) / (1 - ecc)) * sinE / (1 + cosE)\n",
        "    sinf = 2 * tanf2 / (1 + tanf2 ** 2)\n",
        "    cosf = (1 - tanf2 ** 2) / (1 + tanf2 ** 2)\n",
        "\n",
        "    # Evaluate the RV model\n",
        "    rv_model = numpyro.deterministic(\n",
        "        \"rv_model\", bkg + K * (cosw * (cosf + ecc) - sinw * sinf)\n",
        "    )\n",
        "\n",
        "    # Condition on the observations\n",
        "    numpyro.sample(\"obs\", dist.Normal(bkg + rv_model, rv_err), obs=rv)\n",
        "\n",
        "\n",
        "# It's often useful to initialize well\n",
        "init_values = {\n",
        "    \"K\": 56.0,\n",
        "    \"P\": 4.230785,\n",
        "    \"ecc\": 0.01,\n",
        "    \"phi_angle\": np.array([0.85, 0.5]),\n",
        "    \"w_angle\": np.array([0.0, 1.0]),\n",
        "    \"rv0\": -1.8,\n",
        "    \"rv_trend\": -1.6,\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_values),\n",
        ")\n",
        "mcmc = MCMC(nuts_kernel, num_warmup=1000, num_samples=1000, num_chains=1)\n",
        "rng_key = random.PRNGKey(42307)\n",
        "%time mcmc.run(rng_key, t, rv_err, rv=rv)\n",
        "\n",
        "samples = mcmc.get_samples()\n",
        "plt.hist(samples[\"K\"], 20, histtype=\"step\")\n",
        "plt.yticks([])\n",
        "plt.xlabel(\"RV semiamplitude [m/s]\");"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "sample: 100%|██████████| 2000/2000 [02:16<00:00, 14.65it/s, 1023 steps of size 2.19e-03. acc. prob=0.97]\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "CPU times: user 2min 17s, sys: 1.31 s, total: 2min 18s\n",
            "Wall time: 2min 18s\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEGCAYAAABbzE8LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAOgElEQVR4nO3dfYxldX3H8c9HtrBTtQ52ocVBGDFitdgCZenDSgurpSursQ+kVWsr2sTYKv6hfZjUhJDSJpOSxrSx2qgxJF2DkPqIi8U1y2K78WGXx+XJAtshOKkKtaMYx1Lh2z9+vwnX2Tt3Zu7DfM+5834lk7lz7vnd8z2/mXzm3N+553ccEQIAbLxnZBcAAJsVAQwASQhgAEhCAANAEgIYAJJsWc/K27Zti+np6RGVAgDj6dZbb30sIk5avnxdATw9Pa3Dhw8PryoA2ARsP9xtOUMQAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkGRdV8IB42jH7H7NLyz23X5qckIHZ3YOsSJsFgQwNr35hUXNze7uu/30zN4hVoPNhCEIAEhCAANAEgIYAJIwBoyxMMiJtKnJiSFXA6wNAYyxMOiJNCADQxAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIMmW7AKAJTtm92t+YbGvtlOTE0OuBhg9AhiNMb+wqLnZ3dllABuGIQgASEIAA0AShiCAAU1NTmh6Zm/fbQ/O7BxyRWgLAhgY0CAB2m9wYzwwBAEASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIS7IgOJuKX95kYAA4m4pf3mxhAEACQhgAEgCQEMAEkIYABIwkk4DNWO2f2aX1jsq+3U5MSQqwGajQDGUM0vLGpudnd2GUArMAQBAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQJIt2QWgWXbM7tf8wmLf7acmJ4ZYDTDeCGD8iPmFRc3N7s4uA9gUGIIAgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCXNBAC01NTmh6Zm9fbc9OLNzyBVhvQhgoKUGCdB+gxvDxRAEACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkXIgxhga5tTy3lQc2DgE8hri1PNAODEEAQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkW7ILALDxpiYnND2zt++2B2d2DrmizYkABjahQQK03+DGsRiCAIAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASpqMEsC6DzCW81J75hAsCGMC6DBqezCf8NIYgACAJAQwASRiCaKgds/s1v7DYV9upyYkhVwNgFAjghppfWNTc7O7sMgCMEEMQAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJGEuiBFiQh0AvRDAI8SEOsCxBrmjxrjdTYMABrChBgnQcbubBmPAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIwsfQAGwKg14YNYrPH499AA/S6YPiajagOQa5MGpUnz/esADO+u/D1WjA+Bj0Krqm2bAAbuJ/HwDtMk6XIUuchAOANAQwACQhgAEgiSNi7Svbj0p6eHTl9LRN0mNJ2+5X22puW71S+2puW71S+2puYr2nR8RJyxeuK4Az2T4cEedl17Eebau5bfVK7au5bfVK7au5TfUyBAEASQhgAEjSpgD+YHYBfWhbzW2rV2pfzW2rV2pfza2ptzVjwAAwbtp0BAwAY4UABoAkjQhg23O2j9i+w/bhZc+923bY3rZC2zfZfqB+vWljKh645idruztsfyarXttX2p7vqOWSFdrusv012w/antmIeodQ84q/n42sty6/3Pb9tu+x/bcrtG1MH6+j5kb0se3rOv4e5mzfsULblD7uKSLSvyTNSdrWZfnzJd2kcvFHt+efK+lo/X5ifXxik2uu63yvCX0s6UpJf7pKu+MkPSTpDEnHS7pT0kubXHOv309CvRdJ+oKkE+rPJ7egj1etuUl9vOz5v5N0RZP6uNdXI46Ae3ivpD+XtNKZwt+QtC8ivh0R/yNpn6RdG1XcClaruW3Ol/RgRByNiCckfUzSa5NrapM/ljQbEf8rSRHxrS7rNK2P11Jz49i2pN+VdG2Xp5vWx5IaMgShElaft32r7bdKku3XSpqPiDt7tJuS9EjHz1+vyzZCvzVL0lbbh21/2fZvjrzS4ph6q3fYvsv2R2yf2KVdo/q4Wq3mXm1Hqds2z5R0ge2v2L7F9vYu7ZrWx2upeaW2o9ZrmxdI+mZEPNClXWYfr6gpd8R4eUTM2z5Z0j7b90v6S0kXJ9fVyyA1n17bniFpv+0jEfHQSKvtXu8HJF2l8kd9lcrbt7eMuI71GKTmY9pGxBcT6t2iMkT2S5K2S7re9hlR3xc3wCA1N6KPO7b5enU/+m2sRhwBR8R8/f4tSZ+U9GuSXiDpTttzkk6VdJvtn17WdF5lzHXJqXXZyA1Qc2fbo5IOSDonod7zI+KbEfFkRDwl6UMqb9OWa1Ifr7Xmrm0z6lU50vpEFF+V9JTKZDGdGtXHWlvNTepj2d4i6bclXbdC07Q+7iU9gG0/0/azlx6rHEEeioiTI2I6IqZV/iDOjYhvLGt+k6SLbZ9Y34peXJc1tuZa6wn18TZJOyTdm1Dv3bZP6VjttyTd3aX5IUkvsv0C28dLep2kkX9yY5CaV2qbUa+kT6mc1JLtM1VOAC2fqatRfbyWmhvWx5L0Skn3R8TXV2ie0seryj4LqHJW8s76dY+k93RZZ071zKek8yR9uOO5t0h6sH69uek1S/oVSUdq2yOS/iirXkn/XGu4S+WP8ZS6/HmSbuxof4mk/1A5i3zMvjat5rX8fjaw3uMl7VEJitsk7WxBH69ac5P6uD53jaS3LVs/vY9X++JSZABIkj4EAQCbFQEMAEkIYABIQgADQBICGACSEMBjzE/Puna37RtsT7rMHnftsvW22X506fPJI67pr2y/csTbOGD7vPr4xrrfk7b/pI/XutD2Z9e5/nds39jHtk6x/fkez99s+3tL+4b2I4DH22JEnB0RZ0n6tqS3q1w99Ou2f7xjvUsl3RB18pVRiogrIuILo95Ox/YuiYgFSZOS1h3Affq3iOg6TeYqdqnHhUQRcZGkDZn2ERuDAN48viRpKiK+K+kWSa/peO516nINve132r63TnzzsbrsmXUSnK/avt1lAiLZvsz2p2zvc5mT9R2231XX+bLt59b1rrF9aX18he1D9Qj9g7Zdlx+w/V6XCYvus73d9idc5nz+67rOtMt8tR+t6/zLsn8qS/swV684nJX0wvqO4OrlR7a232f7svp4V33t21Qub11ap+u+91K3c4vtT9s+anvW9u/X1zhi+4Udq++S9Ll6JPzFjncvF6y2HbQTAbwJ2D5O0iv09KWX16qErmw/T2X2q/1dms5IOicifk7S2+qy90jaHxHnq1yuenW9LFSSzlIJrO2S/kbS9yPiHJXw/8Mur/++iNhej9AnJL2647knIuI8Sf8k6dMqR+9nSbrM9k/WdV4s6f0R8RJJ31XvI9wZSQ/VdwR/ttJKtreqzDHxGkm/IKlzLo9e+97Lz6v030sk/YGkM+trfFjS5XW7x0l6cUTcK+kNkm6KiLNr264TjKP9CODxNuFyd4BvSPoplfmSJWmvpB22f0Jl/tSPR8STXdrfJemjtt8o6Yd12cWSZurrHpC0VdJp9bmbI+LxiHhU0nck3VCXH5E03eX1L3KZ8vCIpJ2SfrbjuaV/Fkck3RMR/1WHSI7q6UlVHomIg/XxHkkv79kba/Mzkv4zIh6Icpnono7neu17L4c66n9I0tI4b2e//KKkryytL+nNtq+U9LKIeLz/3UGTEcDjbbEeRZ0uySpHkYqIRUn/qjKZTdfhh2q3pH+UdK6kQy4zTlnS79QjybMj4rSIuK+u3zmG/FTHz09p2dSn9Ujz/ZIujYiXqRx1bu1YpbPt8tddeq3l19Gv57r6H+pH//63rrRih1773sta+uVVKr8TRZle8VdVZuu6xna3dw8YAwTwJhAR35f0TknvriEqldB9l8qR8ZeWt7H9DEnPj4ibJf2FpOdIepbKSaLLO8Zr+51KcynwHrP9LJUTget1mu1fro/fIOnfe6z7uKRnd/z8sKSX2j7B9qTKEI0k3S9pumNs9vUdbYa17928QuU2QLJ9usrE4h9SGaY4d4jbQYMQwJtERNyuMqSwFCj7VGaLui66z8h0nKQ9dXjgdkn/UD9NcJWkH5N0l+176s/91LOgctR7t0qwHerjZb4m6e2271O5J+AHemzvvyUdrCe1ro6IRyRdX7d/vco+KiJ+IOmtkvbWk3Cdt+MZyr4vZ/skST/oGGq4UGVe6dsl/Z6kvx/GdtA8zIaGVrI9Lemz9QReY9i+UOWmoa9ebd2ONm+UdGpEzK5h3QP19fk42hhoyi2JgHHxhKSzbN+41s8CR8Se1dcqF2KozIn7fwPUhwbhCBgAkjAGDABJCGAASEIAA0ASAhgAkhDAAJDk/wERYADW8r3OBgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JYFjEFNbTB5Z"
      },
      "source": [
        "This works, but samples a lot slower than the _exoplanet_ example so I'm not totally sure what's up with that!"
      ]
    }
  ]
}