{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%run notebook_setup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Getting started with The Joker\n", "\n", "*The Joker* (pronounced Yo-kurr) is a highly specialized Monte Carlo (MC) sampler that is designed to generate converged posterior samplings for Keplerian orbital parameters, even when your data are sparse, non-uniform, or very noisy. This is *not* a general MC sampler, and this is *not* a Markov Chain MC sampler like emcee, or pymc3: This is fundamentally a [rejection sampler](https://en.wikipedia.org/wiki/Rejection_sampling) with some tricks that help improve performance for the two-body problem.\n", "\n", "*The Joker* shines over more conventional MCMC sampling methods when your radial velocity data is imprecise, non-uniform, sparse, or has a short baseline: In these cases, your likelihood function will have many, approximately equal-height modes that are often spaced widely, all properties that make conventional MCMC bork when applied to this problem. In this tutorial, we will not go through the math behind the sampler (most of that is covered [in the original paper](https://arxiv.org/abs/1610.07602)). However, some terminology is important to know for the tutorial below or for reading the documentation. Most relevant, the parameters in the two-body problem (Kepler orbital parameters) split into two sets: nonlinear and linear parameters. The nonlinear parameters are always the same in each run of The Joker: period $P$, eccentricity $e$, argument of pericenter $\\omega$, and a phase $M_0$. The default linear parameters are the velocity semi-ampltude $K$, and a systemtic velocity $v_0$. However, there are ways to add additional linear parameters into the model (as described in other tutorials).\n", "\n", "For this tutorial, we will set up an inference problem that is common to binary star or exoplanet studies, show how to generate posterior orbit samples from the data, and then demonstrate how to visualize the samples. Other tutorials demonstrate more advanced or specialized functionality included in *The Joker*, like:\n", "- [fully customizing the parameter prior distributions](2-Customize-prior.ipynb), \n", "- [allowing for a long-term velocity trend in the data](3-Polynomial-velocity-trend.ipynb), \n", "- [continuing sampling with standard MCMC methods](4-Continue-sampling-mcmc.ipynb) when *The Joker* returns one or few samples,\n", "- [simultaneously inferring constant offsets between data sources](5-Calibration-offsets.ipynb) (i.e. when using data from multiple instruments that may have calibration offsets)\n", "\n", "But let's start here with the most basic functionality!\n", "\n", "First, imports we will need later:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] } ], "source": [ "import astropy.table as at\n", "from astropy.time import Time\n", "import astropy.units as u\n", "from astropy.visualization.units import quantity_support\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "%matplotlib inline\n", "\n", "import thejoker as tj" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# set up a random generator to ensure reproducibility\n", "rnd = np.random.default_rng(seed=42)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading radial velocity data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, we need some radial velocity data to play with. Our ultimate goal is to construct or read in a thejoker.RVData instance, which is the main data container object used in *The Joker*. For this tutorial, we will use a simulated RV curve that was generated using a separate script and saved to a CSV file, and we will create an RVData instance manually. \n", "\n", "Because we previously saved this data as an Astropy [ECSV](http://docs.astropy.org/en/latest/api/astropy.io.ascii.Ecsv.html#astropy.io.ascii.Ecsv) file, the units are provided with the column data and read in automatically using the [astropy.table read/write interface](http://docs.astropy.org/en/latest/table/index.html):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "QTable length=2\n", "\n", "\n", "\n", "\n", "\n", "\n", "
bjdrvrv_err
km / skm / s
float64float64float64
2458512.359806455-12.0713970989103073.0870942699366886
2458512.7746735318-12.6683283530919330.22630471074546044
" ], "text/plain": [ "\n", " bjd rv rv_err \n", " km / s km / s \n", " float64 float64 float64 \n", "------------------ ------------------- -------------------\n", " 2458512.359806455 -12.071397098910307 3.0870942699366886\n", "2458512.7746735318 -12.668328353091933 0.22630471074546044" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_tbl = at.QTable.read('data.ecsv')\n", "data_tbl[:2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The full simulated data table has many rows (256), so let's randomly grab 4 rows to work with:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sub_tbl = data_tbl[rnd.choice(len(data_tbl), size=4, replace=False)]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "QTable length=4\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
bjdrvrv_err
km / skm / s
float64float64float64
2458580.047609328-15.4455319515741251.4677548574751984
2458526.5749810245-18.1531452641809070.11717712773471835
2458609.74259345-4.4730160900183230.4059567786902223
2458626.8563795867-20.4494721676787852.4368445579891596
" ], "text/plain": [ "\n", " bjd rv rv_err \n", " km / s km / s \n", " float64 float64 float64 \n", "------------------ ------------------- -------------------\n", " 2458580.047609328 -15.445531951574125 1.4677548574751984\n", "2458526.5749810245 -18.153145264180907 0.11717712773471835\n", " 2458609.74259345 -4.473016090018323 0.4059567786902223\n", "2458626.8563795867 -20.449472167678785 2.4368445579891596" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sub_tbl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like the time column is given in Barycentric Julian Date (BJD), so in order to create an RVData instance, we will need to create an astropy.time.Time object from this column:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "t = Time(sub_tbl['bjd'], format='jd', scale='tcb') \n", "data = tj.RVData(t=t, rv=sub_tbl['rv'], rv_err=sub_tbl['rv_err'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have an RVData object, so we could continue on with the tutorial. But as a quick aside, there is an alternate, more automatic (automagical?) way to create an RVData instance from tabular data: RVData.guess_from_table. This classmethod attempts to guess the time format and radial velocity column names from the columns in the data table. It is very much an experimental feature, so if you think it can be improved, please open an issue in the [GitHub repo for The Joker](https://github.com/adrn/thejoker/issues). In any case, here it successfully works:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "data = tj.RVData.guess_from_table(sub_tbl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the handy features of RVData is the .plot() method, which generates a quick view of the data:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "findfont: Font family ['cursive'] not found. Falling back to DejaVu Sans.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEGCAYAAABGnrPVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAal0lEQVR4nO3dfZRddX3v8fcnCXAZRAMyCBIyh6eKgBr1FEuLFjEqpoEA8hDvWOG2dujtTZcP91phTau0XbNWRSntvSpyQNbleo+C5fKQBBCJxWptwU4gwERAA0xCAsUBQbBpUcj3/rF/EzbDmTPnTPach5nPa62zZu/f/u19vr+cyXzO3vucvRURmJmZFWVeuwswM7PZxcFiZmaFcrCYmVmhHCxmZlYoB4uZmRVqQbsL6AT77bdflEqldpdhZtY11q9f/2RE9NZa5mABSqUSw8PD7S7DzKxrSNo82TIfCjMzs0I5WMzMrFAOFjMzK5SDxczMCuVgMTOzQjlYzMysUA4WMzMrlIPFzMwK5WAxM2vQmjVr2l1CV3CwmJk1aO3ate0uoSs4WMzMplCtVimVSlQqFUqlEtVqtd0ldTRfK8zMrI5qtcrAwADbt28HYPPmzQwMDADQ39/fztI6lvdYzMzqGBwc3Bkq47Zv387g4GCbKup8DhYzszq2bNnSVLs5WMzM6lq8eHFT7dYlwSLpQknbJG1Ij2WT9DtJ0oOSNkk6v9V1mtnsMzQ0RE9Pz8vaenp6GBoaalNFna8rgiW5JCKWpMfNExdKmg98CfgAcBTwIUlHtbpIM5td+vv7qVQq9PX1AdDX10elUvGJ+zq6KVimciywKSIejohfAlcDK9pck5nNAv39/YyOjrJ69WpGR0cdKlPopmBZJeleSVdK2qfG8oOAR3PzW1NbTZIGJA1LGh4bGyu6VjObhU4++eR2l9AVOiZYJK2TNFLjsQK4FDgMWAI8Dly8q88XEZWIKEdEube3d1c3Z2ZmScd8QTIiljbST9LlQK3rKmwDDs7NL0ptZmbWQh2zx1KPpANzs6cBIzW6/QtwhKRDJO0OrARWt6I+MzN7ScfssUzhIklLgABGgfMAJL0euCIilkXEC5JWAbcC84ErI2Jjuwo2M5uruiJYIuJ3J2l/DFiWm78ZeMVHkc3MrHW64lCYmZl1DweLmZkVysFiZmaFcrCYmVmhHCxmZlYoB4uZmRXKwWJmZoVysJiZWaEcLGZmVigHi5mZFcrBYmZmhXKwmJlZoRwsZmZWKAeLmZkVysFiZmaFcrCYmVmhHCxmZlaorriDpKRrgDek2YXAMxGxpEa/UeA54EXghYgot6xIMzMDuiRYIuLs8WlJFwM/r9P93RHx5MxXZWZmtXRFsIyTJOAs4MR212JmZrV12zmWdwJPRMRPJlkewLclrZc0UG9DkgYkDUsaHhsbK7xQM7O5qmP2WCStAw6osWgwIm5M0x8CvlFnM8dHxDZJ+wO3SXogIr5Xq2NEVIAKQLlcjl0o3czMcjomWCJiab3lkhYApwNvr7ONbennTyVdDxwL1AwWMzObGd10KGwp8EBEbK21UNJekvYenwbeB4y0sD4zM6O7gmUlEw6DSXq9pJvT7OuAf5R0D/BD4KaI+FaLazQzm/M65lDYVCLi3BptjwHL0vTDwFtaXJaZmU3QTXssZmbWBRwsZmZWKAeLmZkVysFiZmaFcrCYmVmhHCxmZlYoB4uZmRXKwWJmZoVysJiZWaEcLGZmVigHi5mZFcrBYmZmhXKwmJlZoRwsZmZWKAeLmZkVysFiZmaFcrCYmVmhOipYJJ0paaOkHZLKE5ZdIGmTpAclvX+S9Q+RdGfqd42k3VtTuZmZjeuoYAFGgNOB7+UbJR1Fds/7o4GTgC9Lml9j/c8Bl0TE4cDTwO/PbLlmZjZRRwVLRNwfEQ/WWLQCuDoino+IR4BNwLH5DpIEnAhcm5quAk6dyXrNzOyVOipY6jgIeDQ3vzW15b0WeCYiXqjTx8zMZtiCVj+hpHXAATUWDUbEjS2sYwAYAFi8eHGrntbMbNZrebBExNJprLYNODg3vyi15T0FLJS0IO211OqTr6MCVADK5XJMoyYzM6uhWw6FrQZWStpD0iHAEcAP8x0iIoDbgTNS0zlAy/aAzMws01HBIuk0SVuB44CbJN0KEBEbgW8CPwK+Bfy3iHgxrXOzpNenTXwa+KSkTWTnXL7a6jGYmc11yt7oz23lcjmGh4fbXYaZWdeQtD4iyrWWddQei5mZdT8Hi5mZFcrBYmZmhXKwmJlZoRr6HoukfRvotiMintnFeszMrMs1+gXJx9JDdfrMB/wVdjOzOa7RYLk/It5ar4Okuwuox8zMulyj51iOK6iPmZnNcg0FS0T8x/i0pE9P1cfMzOauKQ+FSfpmfhZYQnZDLTMzs1do5BzLsxHx0fEZSZfOYD1mZtblGjkUNjRhfnAmCjEzs9lhymBJtwJG0n5p/mczXZSZmXWvZr55f+WMVWFmZrNGM8FS78uRZmZmQHPB4hu3mJnZlLzHYmZmhWomWC6YsSrMzGzWaDhYImJE0pmS9gaQ9KeSrpP0tqKKSdvfKGmHpHKu/b2S1ku6L/08cZL1L5S0TdKG9FhWVG1mZtaYZu/H8mcR8Zyk44GlwFeBIr8wOQKcDnxvQvuTwMkR8SbgHOBrdbZxSUQsSY+bC6zNzMwa0GywvJh+/g5QiYibgN2LKiYi7o+IB2u03x0Rj6XZjcCekvYo6nnNzKw4zQbLNkmXAWcDN6c/7q2+C+UHgbsi4vlJlq+SdK+kKyXtM9lGJA1IGpY0PDY2NjOVmpnNQc2GwlnArcD7090i9wU+1cwGJK2TNFLjsaKBdY8muwDmeZN0uRQ4jOxCmY8DF0+2rYioREQ5Isq9vb3NDMHMzOpo9EZfAETEduC63PzjZH/Am9nG0mb6j5O0CLge+EhEPDTJtp/I9b8cWDud5zIzs+lr9WGsaZG0ELgJOD8iflCn34G52dPIPgxgZmYt1FHBIuk0SVvJ7kZ5k6Rb06JVwOHAZ3IfJd4/rXNF7qPJF6WPJN8LvBv4RKvHYGY21ymi8Su1pD/gg0Af2WE0ARERb56Z8lqjXC7H8PBwu8swM+saktZHRLnWsqbOsQBVspP19wE7drUwMzObfZoNlrGIWD0jlZiZ2azQbLB8VtIVwHeAnd8jiYjrJl/FzMzmkmaD5b8ARwK78dKhsCD3EWQzM5vbmg2WX4+IN8xIJWZmNis0+3Hjf5J01IxUYmZms0Kzeyy/AWyQ9AjZOZZZ8XFjMzMrTrPBctKMVGFmZrNGs8GyX0SszzdIWg5sLq4kMzPrZs2eY7lc0jHjM5I+BPxZsSWZWbPWrFnT7hLMdmo2WM4A/o+kIyX9AfBHwPuKL8vMmrF2rS/kbZ2j2cvmPyxpJXADsAV4X0T8+4xUZmZmXamhYJF0H9kXIcftC8wH7pSEPxVmZmbjGt1jWT6jVZiZ2azRULBEhD/1ZWZmDWno5L2ku4roY2Zms1+jh8LemO7KOBkBrymgHjMz63KNBsuRDfR5cVcKAZB0JnAh8Ebg2IgYTu0l4H7gwdT1joj4wxrr7wtcA5SAUeCsiHh6V+syM7PGddo5lhHgdOCyGsseioglU6x/PvCdiPgrSeen+U8XXKOZmdXR7BckZ1RE3B8RD07dc1IrgKvS9FXAqbtelZmZNaOjgmUKh0i6W9I/SHrnJH1eFxGPp+l/BV432cYkDUgaljQ8NjZWeLFmZnNVo58K+5Kk3yriCSWtkzRS47GizmqPA4sj4q3AJ4GvS3p1veeJiODlX+qcuLwSEeWIKPf29k5rLGZm9kqN7rH8GPiCpFFJF0l663SfMCKWRsQxNR431lnn+Yh4Kk2vBx4Cfq1G1yckHQiQfv50unWadYNqtUqpVKJSqVAqlahWq+0uyayxYImIv42I44DfBp4CrpT0gKTPSqr1B75QknolzU/ThwJHAA/X6LoaOCdNnwNMGlZm3a5arTIwMMDmzdlnazZv3szAwIDDxdpO2RGjaayY7bVcCbw5IuYXUox0GvC/gF7gGWBDRLxf0geBvwB+BewAPhsRa9I6VwBfiYhhSa8FvgksJrtHzFkR8bOpnrdcLsfw8HARQzBrmVKptDNU8vr6+hgdHW19QTanSFofEeWay5oJFkkLgA8AK4H3AN8FvlHvMFY3cLBYN5o3bx61/v9KYseOHW2oyOaSesHS6Mn790q6EtgG/AFwE3BYRKzs9lAx61aLFy9uqt2sVRo9eX8B8E/AkRFxSkR8PSL+TdI8Sf0zWJ+ZTWJoaIienp6XtfX09DA0NNSmiswyjQbLacD+wF+mvRdJ+mOyE+hnzVh1Zjap/v5+KpUKfX19QHZupVKp0N/v93rWXg2dY5F0I/Az4A6ycyv7k1148mMRsWFGK2wBn2Oxbnfeeedx2WW1roRkNjPqnWNp9CKUh0bEirSxK3jpC4v/UVCNZmY2SzR6KOxX4xMR8SKw1aFiZma1NLrH8hZJz6ZpAXumeZFdPaXu5VXMzGzuaPSy+YV8AdLMzGa/brq6sZmZdQEHi5mZFcrBYmZmhXKwmJlZoRwsZmZWKAeLmZkVysFiNgssX7683SWY7eRgMZsFTj755HaXYLaTg8XMzArVUcEi6UxJGyXtkFTOtfdL2pB77JC0pMb6F0raluu3rLUjMDOzRq8V1iojwOnAy67/HRFVoAog6U3ADXUu139JRHxhRqs0M7NJdVSwRMT9kN2zu44PAVe3pCAzM2taRx0Ka9DZwDfqLF8l6V5JV0raZ7JOkgYkDUsaHhsbK75KM7M5quXBImmdpJEajxUNrPsOYHtEjEzS5VLgMGAJ2c3ILp5sWxFRiYhyRJR7e3unMxQzM6uh5YfCImLpLqy+kjp7KxHxxPi0pMuBtbvwXGZmNg1dcyhM0jzgLOqcX5F0YG72NLIPA5iZWQt1VLBIOk3SVuA44CZJt+YWvwt4NCIenrDOFbmPJl8k6T5J9wLvBj4xU7VWq1VKpRLz5s2jVCpRrVZn6qnMzLqKIqLdNbRduVyO4eHhhvtXq1UGBgbYvn37zraenh4qlQr9/f0zUaKZWUeRtD4iyrWWddQeS7cYHBx8WagAbN++ncHBwTZVZGbWvDVr1szIdh0s07Bly5am2s3MOtHatTPz+SYHyzQsXry4qXYzs7nEwTINQ0ND9PT0vKytp6eHoaGhNlVkZtY5HCzT0N/fT6VSoa+vD0n09fX5xL2ZWdJR1wrrJv39/Q4SM7MavMdiZmaFcrCYmVmhHCxmZlYoB4uZmRXKwWJmZoVysJiZWaEcLGZmVigHi5mZFcrBYmZmhXKwmJlZoRwsZmZWqI4KFkmfl/SApHslXS9pYW7ZBZI2SXpQ0vsnWf8QSXemftdI2r111ZuZGXRYsAC3AcdExJuBHwMXAEg6ClgJHA2cBHxZ0vwa638OuCQiDgeeBn6/JVWbmdlOHRUsEfHtiHghzd4BLErTK4CrI+L5iHgE2AQcm19XkoATgWtT01XAqTNftZmZ5XVUsEzwe8Atafog4NHcsq2pLe+1wDO5YKrVx8zMZljL78ciaR1wQI1FgxFxY+ozCLwAVGewjgFgAHxLYTOzIrV8jyUilkbEMTUe46FyLrAc6I+ISKttAw7ObWZRast7ClgoaUGdPvk6KhFRjohyb29vASMzM+sO1WqVUqlEpVKhVCpRrRb7Hr6jDoVJOgn4E+CUiNieW7QaWClpD0mHAEcAP8yvm0LoduCM1HQOcOPMV21m1j2q1SoDAwNs3rwZgM2bNzMwMFBouOilnYL2k7QJ2INs7wPgjoj4w7RskOy8ywvAxyPiltR+M/DRiHhM0qHA1cC+wN3AhyPi+amet1wux/DwcOHjMTPrNKVSaWeo5PX19TE6OtrwdiStj4hyzWWdFCzt4mAxs7li3rx51Pq7L4kdO3Y0vJ16wdJRh8LMzGxmTfZhpSI/xORgMTObQ4aGhujp6XlZW09PD0NDQ4U9h4PFzGwO6e/vp1Kp0NfXB2TnViqVCv39/YU9h8+x4HMsZjY3nXfeeVx22WXTWtfnWMzMrGUcLGZmVigHi5mZFcrBYmZmhXKwmJlZoRwsZmZWKAeLmZkVysFiZmaFcrCYmVmhHCxmZlYoB4uZmRXKwWJmZoVysJiZWaEcLGZmVqiOChZJn5f0gKR7JV0vaWFqf6+k9ZLuSz9PnGT9CyVtk7QhPZa1dgRmZtZRwQLcBhwTEW8GfgxckNqfBE6OiDcB5wBfq7ONSyJiSXrcPLPlmpnZRB0VLBHx7Yh4Ic3eASxK7XdHxGOpfSOwp6Q92lGjmZnV11HBMsHvAbfUaP8gcFdEPD/JeqvSobQrJe0z2cYlDUgaljQ8NjZWRL1mZkYbgkXSOkkjNR4rcn0GgReA6oR1jwY+B5w3yeYvBQ4DlgCPAxdPVkdEVCKiHBHl3t7eXRyVmZmNW9DqJ4yIpfWWSzoXWA68JyIi174IuB74SEQ8NMm2n8j1vxxYW0TNZmbWuI46FCbpJOBPgFMiYnuufSFwE3B+RPygzvoH5mZPA0ZmqlYzM6uto4IF+CKwN3Bb+rjwV1L7KuBw4DO5jxLvDyDpCknl1O+i9JHke4F3A59o9QDMzLrF8uXLZ2S7yh1tmrPK5XIMDw+3uwwzs64haX1ElGst67Q9FjMz63IOFjMzK5SDxczMCuVgMTOzQjlYzMysUA4WMzMrlIPFzMwK5WAxM7NC+QuSgKQxYHO765hh+5Hd12YumWtjnmvjBY+5nfoiouYVfB0sc4Sk4cm+JTtbzbUxz7XxgsfcqXwozMzMCuVgMTOzQjlY5o5Kuwtog7k25rk2XvCYO5LPsZiZWaG8x2JmZoVysJiZWaEcLF1G0mi6S+YGScOpbYmkO8bbJB2b2k+Q9PPcXTc/M2Fb8yXdLWltru0QSXdK2iTpGkm7t3aEr1TUmCUtlHStpAck3S/puNS+r6TbJP0k/dynPSN9SYFj/oSkjZJGJH1D0n9K7R31Ojcz3rTshNS+UdI/5NpPkvRgGtf5ufaOGm+qaZfHLOlgSbdL+lFq/1iuf/t+ryPCjy56AKPAfhPavg18IE0vA76bpk8A1tbZ1ieBr+f7AN8EVqbprwD/dbaMGbgK+Gia3h1YmKYvAs5P0+cDn5sNYwYOAh4B9sy9tud24uvc5HgXAj8CFqf5/dPP+cBDwKHp9b0HOKoTx1vgmA8E3pam9wZ+nBtz236vvccyOwTw6jT9GuCxqVaQtAj4HeCKXJuAE4FrU9NVwKmFVlqcpsYs6TXAu4CvAkTELyPimbR4BdlYYRaNOVkA7ClpAdADPNZFr/Nk4/3PwHURsQUgIn6a2o8FNkXEwxHxS+BqYEUXjReaHHNEPB4Rd6Xp54D7yd5QQDt/r9ud2n409yB7B3oXsB4YSG1vBLYAjwLbyC61ANk72afI3rndAhyd2861wNvJvdslu1TEplyfg4GR2TBmYAnwQ+B/A3eTBepeadkzuedSfr6bx5yWfQz4BTAGVDv1dW5yvH8DfAn4bur/kdR+BnBFbpu/C3yxE8db1JgnbK+U1n11u3+v2/oP68c0XjA4KP3cP/0heRfwP4EPpvazgHVp+tXAq9L0MuAnaXo58OU0fQKdHyxFjLkMvAC8I83/LfCXafqZCc/39CwZ8z7A3wO9wG7ADcCHO/F1bnK8XwTuAPZKY/kJ8Gt0X7Ds8phz23oVWeCcnmtr2++1D4V1mYjYln7+FLiebPf/HOC61OXvUhsR8WxE/CJN3wzsJmk/4LeAUySNkh0uOFHS/yV717swHTYBWET2rqmtChrzVmBrRNyZ1rkWeFuafkLSgQDp5/ihlbYpaMxLgUciYiwifpXW/U068HVuZrxkr+WtEfFvEfEk8D3gLWRjODi32fFxddx4obAxI2k34P+R7ZFe99IztO/32sHSRSTtJWnv8WngfcAI2XHY307dTiR7N4OkA9LxZdKnS+YBT0XEBRGxKCJKwErg7yPiw5G9rbmd7J0fZL/kN7ZkcJMocMz/Cjwq6Q1pnfeQnQwFWE02VphFYyY7LPIbknrS8vcA93fa69zseMlqPV7SAkk9wDvIzi38C3BE+gTY7mS/26s7bbxQ3JjT6/pVstf1ryc8Tft+r9u9O+hH4w+yT7vckx4bgcHUfjzZbvA9wJ3A21P7qtTvHrLd6N+ssc0TePmnwg4lOxexiewd0x6zZcxk51mGgXvJDgvtk9pfC3yH7D/xOmDfWTTmPwceIPuj9bXx17OTXudmx5uWfYrsjcEI8PFc+zKyT0Y9NL6dThtvkWNO/SP9Tm9Ij2Xt/r32JV3MzKxQPhRmZmaFcrCYmVmhHCxmZlYoB4uZmRXKwWJmZoVysJg1QNmVkf8oN/96SdfWW2eaz3OhpG2S/iLNnytpLHdV22vT9xjG+4akw3Prfzy1ldP8aPqyJJJezG3nHkn/XdK8tOyd6Qq5I0WPyeYeB4tZYxYCO4MlIh6LiDPq9N8Vl0RE/hYH10TEkog4GvglcHZu2X1kXwQcdybZ9yJq+ffcdt4LfAD4LEBEfJ/sOyBmu8zBYtaYvwIOS+/4Py+pNP7uPu1V3JDueTEqaZWkTyq7180dkvZN/Q6T9C1J6yV9X9KRzRSQLkmyF/B0rvkGsqvYIukw4OfAk1NtK7LLiAwAq8a/tW9WFAeLWWPOBx5K7/g/VWP5McDpwK8DQ8D2iHgr8M/AR1KfCvDHEfF24H8AX27wuc+WtIHs+lb7Amtyy54lu1TNMWR7Ltc0OqCIeJjsHib7N7qOWSMcLGbFuD0inouIMbK9hvE//vcBJUmvIrsA5N+lkLiM7CZNjbgmIpYAB6TtTQy2q8lC5VSyixmatZWDxawYz+emd+Tmd5DdbGse2WXMl+Qeb2zmCSK7/tIassur560lu0T8loh4ttHtSToUeJEOuJqzzS4OFrPGPEd269dpSX/wH5F0JmR365T0lmls6niyCyzmt70d+DTZIbiGSOolu0XvF8MXDLSCLZi6i5lFxFOSfpBO2N9Cdje/ZvUDl0r6U7Ibb11NdhXbqZwt6XiyN4JbgXNr1Hf1JOsu4KW9pz3TYbjdyG569jVg4qXWzXaZr25s1kEkXQj8IiK+UMC2eoENEXHQlJ2z/iWyWygcs6vPbXObD4WZdZZfAAPjX5CcLkmnAN8HLmiw/zvJzt9M+VFls6l4j8XMzArlPRYzMyuUg8XMzArlYDEzs0I5WMzMrFAOFjMzK9T/B//A4uwueoVyAAAAAElFTkSuQmCC\n", "text/plain": [ "