{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Modifying Priors\n", "\n", "by Sarah Blunt (2018)\n", "\n", "Most often, you will use the `Driver` class to interact with `orbitize`. This class automatically reads your input file, creates all of the `orbitize` objects you need to run an orbit fit, and allows you to run the orbit fit. See the introductory OFTI and MCMC tutorials for examples of working with this class.\n", "\n", "However, sometimes you will want to work with the underlying methods directly. Doing this gives you control over the functionality `Driver` executes automatically, and allows you more flexibility.\n", "\n", "Modifying priors is an example of something you might want to use the underlying API for. This tutorial walks you through how to do that. \n", "\n", "**Goals of this tutorial**:\n", "- Learn to modify priors in `orbitize`\n", "- Learn how to fix a parameter at a specific value\n", "- Learn about the structure of the `orbitize` code base" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "import orbitize\n", "from orbitize import read_input, system, priors, sampler" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in Data\n", "\n", "First, let's read in our data table. This is accomplished with `orbitize.read_input`:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " epoch object quant1 quant1_err quant2 quant2_err quant12_corr quant_type instrument\n", "-------------- ------ ------ ---------- ------ ---------- ------------ ---------- ----------\n", " 55645.95 1 2479.0 16.0 327.94 0.39 nan seppa defsp\n", " 55702.89 1 2483.0 8.0 327.45 0.19 nan seppa defsp\n", " 55785.015 1 2481.0 33.0 326.84 0.94 nan seppa defsp\n", " 55787.935 1 2448.0 24.0 325.82 0.66 nan seppa defsp\n", "55985.19400184 1 2483.0 15.0 326.46 0.36 nan seppa defsp\n", "56029.11400323 1 2487.0 8.0 326.54 0.18 nan seppa defsp\n", "56072.30200459 1 2499.0 26.0 326.14 0.61 nan seppa defsp\n" ] } ], "source": [ "data_table = read_input.read_file('{}/GJ504.csv'.format(orbitize.DATADIR))\n", "\n", "print(data_table)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialize `System` Object\n", "\n", "Next, we initialize an `orbitize.system.System` object. This object stores information about the system you're fitting, such as your data, the total mass, and the parallax." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# number of secondary bodies in system\n", "num_planets = 1\n", "\n", "# total mass & error [msol]\n", "total_mass = 1.22\n", "mass_err = 0.08 \n", "\n", "# parallax & error[mas]\n", "plx = 56.95\n", "plx_err = 0\n", "\n", "sys = system.System(\n", " num_planets, data_table, total_mass, \n", " plx, mass_err=mass_err, plx_err=plx_err\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `System` object has a few handy attributes to help you keep track of your fitting parameters. `System.labels` is a list of the names of your fit parameters, and `System.sys_priors` is a list of the priors on each parameter. Notice that the \"prior\" on parallax (`plx`) is just a float. That's because we fixed this parameter at the printed value by specifying that `plx_err`=0.\n", "\n", "Finally, `System.param_idx` is a dictionary that maps the parameter names from `System.labels` to their indices in `System.sys_priors`." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['sma1', 'ecc1', 'inc1', 'aop1', 'pan1', 'tau1', 'plx', 'mtot']\n", "[Log Uniform, Uniform, Sine, Uniform, Uniform, Uniform, 56.95, Gaussian]\n", "{'sma1': 0, 'ecc1': 1, 'inc1': 2, 'aop1': 3, 'pan1': 4, 'tau1': 5, 'plx': 6, 'mtot': 7}\n" ] } ], "source": [ "print(sys.labels)\n", "print(sys.sys_priors)\n", "print(sys.param_idx)\n", "\n", "# alias for convenience\n", "lab = sys.param_idx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Explore & Modify Priors\n", "\n", "Priors in `orbitize` are Python objects. You can view an exhaustive list [here](https://orbitize.readthedocs.io/en/latest/priors.html). Let's print out the attributes of some of our priors:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'minval': 0.0, 'maxval': 1.0}\n", "{'minval': 0.001, 'maxval': 10000.0, 'logmin': -6.907755278982137, 'logmax': 9.210340371976184}\n" ] } ], "source": [ "print(vars(sys.sys_priors[lab['ecc1']]))\n", "print(vars(sys.sys_priors[lab['sma1']]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check out the priors documentation (linked above) for more info about the attributes of each of these priors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we understand how priors are represented and where they are stored, we can modify them! Here's an example of changing the prior on eccentricity from the current uniform prior to a Gaussian prior:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['sma1', 'ecc1', 'inc1', 'aop1', 'pan1', 'tau1', 'plx', 'mtot']\n", "[Log Uniform, Gaussian, Sine, Uniform, Uniform, Uniform, 56.95, Gaussian]\n", "{'mu': 0.2, 'sigma': 0.05, 'no_negatives': True}\n" ] } ], "source": [ "mu = 0.2\n", "sigma = 0.05\n", "\n", "sys.sys_priors[lab['ecc1']] = priors.GaussianPrior(mu, sigma)\n", "\n", "print(sys.labels)\n", "print(sys.sys_priors)\n", "print(vars(sys.sys_priors[lab['ecc1']]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's do one more example. Say we want to fix the inclination to a particular value (i.e. not allow it to vary in the fit at all), perhaps the known inclination value of a disk in the system.\n", "We can do that as follows:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['sma1', 'ecc1', 'inc1', 'aop1', 'pan1', 'tau1', 'plx', 'mtot']\n", "[Log Uniform, Gaussian, 2.5, Uniform, Uniform, Uniform, 56.95, Gaussian]\n", "Inclination \"prior:\" 2.5\n", "Eccentricity prior: Gaussian\n" ] } ], "source": [ "sys.sys_priors[lab['inc1']] = 2.5\n", "\n", "print(sys.labels)\n", "print(sys.sys_priors)\n", "print('Inclination \"prior:\" {}'.format(sys.sys_priors[sys.param_idx['inc1']]))\n", "print('Eccentricity prior: {}'.format(sys.sys_priors[sys.param_idx['ecc1']]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run OFTI\n", "\n", "All right! We're in business. To finish up, I'll demonstrate how to run an orbit fit with our modified `System` object, first with OFTI, then with MCMC." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [] }, { "data": { "text/plain": [ "Text(0, 0.5, 'number of orbits')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAS6klEQVR4nO3df7BndV3H8eerRUURZVcWh/jhLtNqiWnCqpWGFFYiKEximmGrUjuSCkVNLFlDOU1hmoUNoZtkq2P+yExoUItIMsdEd1dgQSJ+igsEqyiiJAi9++N79nhd7r179t57vuf+eD5mztxzzvece96f+e7uaz/nc36kqpAkCeAHhi5AkjR/GAqSpJahIElqGQqSpJahIElq7TV0AbOx//7716pVq4YuQ5IWlC1btny1qlZO9tmCDoVVq1axefPmocuQpAUlyZen+szTR5KklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKk1oK+o1mar1ZtuHiwY99yznGDHVsLnz0FSVLLUJAktQwFSVLLUJAktXoLhSR/k+SuJFdPWLciySVJrm9+Lm/WJ8k7ktyQ5KokR/RVlyRpan32FP4WeOEu6zYAl1bVGuDSZhngWGBNM60Hzu+xLknSFHoLhar6NHD3LqtPADY185uAEyesf2+NfA7YL8mBfdUmSZrcuMcUnlhVdwA0Pw9o1h8EfGXCdtubdQ+TZH2SzUk279ixo9diJWmpmS8DzZlkXU22YVVtrKq1VbV25cpJXzEqSZqhcd/RfGeSA6vqjub00F3N+u3AIRO2Oxi4fcy1aREa8s5iaSEad0/hImBdM78OuHDC+l9prkL6ceCenaeZJEnj01tPIckHgKOB/ZNsB84GzgE+nOQU4FbgZc3mHwdeBNwA3Ae8pq+6JElT6y0UquqXpvjomEm2LeD1fdUiSepmvgw0S5LmAUNBktTyfQoaC68CkhYGewqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpNYgoZDkN5Nck+TqJB9IsneS1UkuT3J9kg8leeQQtUnSUjb2UEhyEHAasLaqngYsA14BvAX486paA3wdOGXctUnSUjfU6aO9gEcn2Qt4DHAH8DPAR5rPNwEnDlSbJC1ZYw+FqroNeBtwK6MwuAfYAnyjqh5sNtsOHDTu2iRpqRvi9NFy4ARgNfCDwD7AsZNsWlPsvz7J5iSbd+zY0V+hkrQEDXH66AXAzVW1o6q+C3wU+Elgv+Z0EsDBwO2T7VxVG6tqbVWtXbly5XgqlqQlYohQuBX48SSPSRLgGOBLwKeAk5pt1gEXDlCbJC1pQ4wpXM5oQHkrsK2pYSNwJnBGkhuAJwAXjLs2SVrq9tr9JnOvqs4Gzt5l9U3AswcoR5LU8I5mSVLLUJAktQwFSVLLUJAktQwFSVJrt6GQ5LlJ9mnmT07y9iRP6r80SdK4dekpnA/cl+QZwO8AXwbe22tVkqRBdAmFB6uqGD2v6NyqOhfYt9+yJElD6HLz2r1JzgJOBo5Ksgx4RL9lSZKG0KWn8HLgfuCUqvofRo+0fmuvVUmSBtGlp/CbVXXmzoWqujXJ4T3WJEkaSJeews9Osm6y9x9Ikha4KXsKSU4Ffh04LMlVEz7aF/hs34VJmplVGy4e5Li3nHPcIMfV3Jru9NHfAZ8A/gTYMGH9vVV1d69VSZIGMV0oVFXdkuT1u36QZIXBIEmLz+56CscDWxi9LzkTPivgsB7rkiQNYMpQqKrjm5+rx1eOJGlInd68luQXgOcx6iH8R1V9rNeqJEmD6PJAvL8CXsfofcpXA69Lcl7fhUmSxq9LT+H5wNOa5x+RZBOjgJAkLTJdbl67Djh0wvIhwFVTbCtJWsCmu3ntnxiNITweuDbJ55vl5+DNa5K0KE13+uhtY6tCkjQvTHdJ6r83j8n+56p6wRhrkiQNZNoxhap6iNFb1x4/pnokSQPqcvXRd4BtSS4Bvr1zZVWd1ltVkqRBdAmFi5tJkrTI7TYUqmpTkkcCT25WXVdV3+23LEnSEHYbCkmOBjYBtzB6KN4hSdZV1af7LU2SNG5dTh/9GfBzVXUdQJInAx8AjuyzMEnS+HW5o/kROwMBoKr+G3hEfyVJkobSJRQ2J7kgydHN9NeM3rEwY0n2S/KRJP+V5NokP5FkRZJLklzf/Fw+m2NIkvZcl1A4FbgGOA04HfgSo6emzsa5wCer6oeBZwDXMnrl56VVtQa4lO9/BagkaQy6XH10P/D2Zpq1JI8DjgJe3fz+B4AHkpwAHN1stgm4DDhzLo4pSeqmS09hrh0G7ADek+SLSd6dZB/giVV1B0Dz84DJdk6yPsnmJJt37NgxvqolaQkYIhT2Ao4Azq+qZzK6S7rzqaKq2lhVa6tq7cqVK/uqUZKWpOkenf2+qnpVktOr6tw5POZ2YHtVXd4sf4RRKNyZ5MCquiPJgcBdc3hMAas2eGO6pOlN11M4MsmTgNcmWd5cHdROMz1gVf0P8JUkT2lWHcNo8PoiYF2zbh1w4UyPIUmamekGmt8JfJLRGMAWRncz71TN+pl6I/D+5vEZNwGvYRRQH05yCnAr8LJZ/H5J0gxM9z6FdwDvSHJ+VZ06lwetqiuAtZN8dMxcHkeStGe6XJJ6apJnAD/VrPp0VfmOZklahHZ79VGS04D3M7pE9ABGp33e2HdhkqTx6/JAvF8FnlNV3wZI8hbgP4G/7LMwSdL4dblPIcBDE5Yf4vsHnSVJi0SXnsJ7gMuT/GOzfCJwQX8lSZKG0mWg+e1JLgOex6iH8Jqq+mLfhUmSxq9LT4Gq2gps7bkWSdLAhnj2kSRpnjIUJEmtaUMhybIk/zquYiRJw5o2FKrqIeC+JI8fUz2SpAF1GWj+DrAtySWM3n0AQFWd1ltVkqRBdAmFi5tJkrTIdblPYVOSRwOHVtV1Y6hJkjSQLg/EezFwBaN3K5Dkx5Jc1HdhkqTx63JJ6h8Azwa+Ae27EFb3WJMkaSBdQuHBqrpnl3XVRzGSpGF1GWi+OskrgWVJ1gCnAZ/ttyxJ0hC69BTeCBwO3A98APgm8Bt9FiVJGkaXq4/uA97UvFynqure/suSJA2hy9VHz0qyDbiK0U1sVyY5sv/SJEnj1mVM4QLg16vqPwCSPI/Ri3ee3mdhkqTx6zKmcO/OQACoqs8AnkKSpEVoyp5CkiOa2c8neRejQeYCXg5c1n9pkqRxm+700Z/tsnz2hHnvU5CkRWjKUKiqnx5nIZKk4e12oDnJfsCvAKsmbu+jsyVp8ely9dHHgc8B24D/67ccSdKQuoTC3lV1Ru+VSJIG1+WS1Pcl+bUkByZZsXPqvTJJ0th16Sk8ALwVeBPfu+qogMP6KkqSNIwuoXAG8ENV9dW5PHCSZcBm4LaqOj7JauCDwApgK/CqqnpgLo8pSZpel9NH1wD39XDs04FrJyy/BfjzqloDfB04pYdjSpKm0SUUHgKuSPKuJO/YOc3moEkOBo4D3t0sB/gZ4CPNJpuAE2dzDEnSnuty+uhjzTSX/gL4HWDfZvkJwDeq6sFmeTtw0GQ7JlkPrAc49NBD57gsSVraurxPYdNcHjDJ8cBdVbUlydE7V0926Cnq2QhsBFi7dq2P25CkOdTljuabmeQf6Kqa6dVHzwVekuRFwN7A4xj1HPZLslfTWzgYuH2Gv1+SNENdTh+tnTC/N/AyRlcIzUhVnQWcBdD0FH67qn45yd8DJzG6AmkdcOFMjyFJmpndDjRX1dcmTLdV1V8wGhSea2cCZyS5gdEYwwU9HEOSNI0up4+OmLD4A4x6DvtOsfkeqarLaN7NUFU3Ac+ei98rSZqZLqePJr5X4UHgFuAXe6lGkjSoLlcf+V4FSVoiupw+ehTwUh7+PoU391eWJGkIXU4fXQjcA2wB7u+3HEnSkLqEwsFV9cLeK5EkDa7Ls48+m+RHe69EkjS4Lj2F5wGvbu5svp/RIymqqp7ea2WSpLHrEgrH9l6FJGle6HJJ6pfHUYgkaXhdxhQkSUuEoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqRWl5fsSNJurdpw8WDHvuWc4wY79mJjT0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEmtsYdCkkOSfCrJtUmuSXJ6s35FkkuSXN/8XD7u2iRpqRviPoUHgd+qqq1J9gW2JLkEeDVwaVWdk2QDsAE4c4D6ejfk9dySNJ2x9xSq6o6q2trM3wtcCxwEnABsajbbBJw47tokaakbdEwhySrgmcDlwBOr6g4YBQdwwBT7rE+yOcnmHTt2jKtUSVoSBguFJI8F/gH4jar6Ztf9qmpjVa2tqrUrV67sr0BJWoIGCYUkj2AUCO+vqo82q+9McmDz+YHAXUPUJklL2RBXHwW4ALi2qt4+4aOLgHXN/DrgwnHXJklL3RBXHz0XeBWwLckVzbrfBc4BPpzkFOBW4GUD1CZJS9rYQ6GqPgNkio+PGWctkqTv5x3NkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJag3xlNR5wfckS9LD2VOQJLUMBUlSy1CQJLWW7JiCpMVjqDHCW845bpDj9smegiSpZShIklqGgiSpZShIklqGgiSpZShIklqGgiSp5X0KkjRDQz5Dra97JOwpSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqTWvQiHJC5Ncl+SGJBuGrkeSlpp5EwpJlgHnAccCTwV+KclTh61KkpaWeRMKwLOBG6rqpqp6APggcMLANUnSkjKf7mg+CPjKhOXtwHN23SjJemB9s/itJNeNobZx2R/46tBF9MB2LTyLtW2Lpl15y8NW7UnbnjTVB/MpFDLJunrYiqqNwMb+yxm/JJurau3Qdcw127XwLNa2LdZ2wdy1bT6dPtoOHDJh+WDg9oFqkaQlaT6FwheANUlWJ3kk8ArgooFrkqQlZd6cPqqqB5O8AfhnYBnwN1V1zcBljduiPC2G7VqIFmvbFmu7YI7alqqHnbaXJC1R8+n0kSRpYIaCJKllKIzB7h7fkeSoJFuTPJjkpF0+W5fk+mZaN76qd2+W7XooyRXNNO8uKOjQtjOSfCnJVUkuTfKkCZ8t5O9sunYt9O/sdUm2NfV/ZuITE5Kc1ex3XZKfH2/l05tpu5KsSvK/E76zd3Y6YFU59TgxGjS/ETgMeCRwJfDUXbZZBTwdeC9w0oT1K4Cbmp/Lm/nlQ7dptu1qPvvW0G2YZdt+GnhMM38q8KFF8p1N2q5F8p09bsL8S4BPNvNPbbZ/FLC6+T3Lhm7THLRrFXD1nh7TnkL/dvv4jqq6paquAv5vl31/Hrikqu6uqq8DlwAvHEfRHcymXfNdl7Z9qqruaxY/x+i+Glj439lU7ZrvurTtmxMW9+F7N8eeAHywqu6vqpuBG5rfNx/Mpl0zYij0b7LHdxw0hn37Ntva9k6yOcnnkpw4t6XN2p627RTgEzPcd5xm0y5YBN9ZktcnuRH4U+C0Pdl3ILNpF8DqJF9M8u9JfqrLAefNfQqLWKfHd/Swb99mW9uhVXV7ksOAf0uyrapunKPaZqtz25KcDKwFnr+n+w5gNu2CRfCdVdV5wHlJXgn8HrCu674DmU277mD0nX0tyZHAx5IcvkvP4mHsKfRvNo/vmM+P/phVbVV1e/PzJuAy4JlzWdwsdWpbkhcAbwJeUlX378m+A5lNuxbFdzbBB4GdvZ0F/51N0LarOR32tWZ+C6OxiSfv9ohDD6Qs9olRb+wmRgNYOweKDp9i27/l4QPNNzMasFzezK8Yuk1z0K7lwKOa+f2B69ll8Gy+t43RP4g3Amt2Wb+gv7Np2rUYvrM1E+ZfDGxu5g/n+weab2L+DDTPpl0rd7aD0UD1bV3+LA7e6KUwAS8C/rv5y/amZt2bGf1PDOBZjP5H8G3ga8A1E/Z9LaOBrxuA1wzdlrloF/CTwLbmD/g24JSh2zKDtv0rcCdwRTNdtEi+s0nbtUi+s3OBa5p2fWriP66MekY3AtcBxw7dlrloF/DSZv2VwFbgxV2O52MuJEktxxQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZqhJCcn+XzzrPp3JVnWPPt+a5Irk1zabPfYJO9pnnl/VZKXDl27NBUfiCfNQJIfAV4OPLeqvpvkr4CTgT8Cjqqqm5OsaDb/feCeqvrRZt/lgxQtdWAoSDNzDHAk8IUkAI8GngN8ukbP5Keq7m62fQHwip071ug9C9K85OkjaWYCbKqqH2umpwB/yOSPXM4U66V5x1CQZuZS4KQkBwA0p4quBJ6fZPWEdQD/Arxh546ePtJ85gPxpBlK8nLgLEb/ufou8HpGj5j+42bdXVX1s0keC5zH6HTTQ8AfVtVHh6lamp6hIElqefpIktQyFCRJLUNBktQyFCRJLUNBktQyFCRJLUNBktT6f5a3Gdn/2M4SAAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAURUlEQVR4nO3df7RlZX3f8fdHQLCRMCADpQzjQMVWTf2BI7IKSRFsIsR2aAuF/BBCaWdFiWDpWpWka9U27VqVlQZ/pAohYjO4iMAiKhNFUwoSTRFwUALoiE4JwpSpgAKCROyM3/5xntk5Dnfu7Jm5+5y5975fa5119n7Oc+75PuvO3M95nrP3PqkqJEkCeMG0C5Ak7TkMBUlSx1CQJHUMBUlSx1CQJHX2nnYBu+Pggw+uFStWTLsMSZpX7rrrrseraulMj83rUFixYgXr1q2bdhmSNK8k+fb2HnP5SJLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSZ1BQyHJg0nuTXJ3knWt7aAkNyX5Vrs/sLUnyQeTbEhyT5JjhqxNkvR8k5gpvKmqXltVK9v+xcDNVXU0cHPbBzgFOLrdVgOXTaA2SdKYaSwfrQLWtO01wGlj7VfVyO3AkiSHTaE+SVq0hj6juYD/kaSA36+qK4BDq2oTQFVtSnJI63s48PDYcze2tk3jPzDJakYzCZYvXz5w+dKuWXHxZ6b22g++9xen9tqa/4YOheOr6pH2h/+mJN+YpW9maHve18K1YLkCYOXKlX5tnCTNoUGXj6rqkXb/KPBJ4FjgO1uXhdr9o637RuCIsacvAx4Zsj5J0k8aLBSS/FSS/bduAz8P3AesBc5p3c4Bbmjba4Gz21FIxwFPbV1mkiRNxpDLR4cCn0yy9XX+qKo+l+TLwHVJzgMeAs5o/W8ETgU2AM8C5w5YmyRpBoOFQlU9ALxmhvbvAifP0F7A+UPVI0naMc9oliR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUmfwUEiyV5KvJvl02z8yyR1JvpXk2iQvbO37tv0N7fEVQ9cmSfpJk5gpXAisH9u/BHhfVR0NPAGc19rPA56oqpcB72v9JEkTNGgoJFkG/CLwkbYf4CTg+tZlDXBa217V9mmPn9z6S5ImZOiZwvuBfwv8uO2/BHiyqja3/Y3A4W37cOBhgPb4U63/T0iyOsm6JOsee+yxIWuXpEVnsFBI8lbg0aq6a7x5hq7V47G/bqi6oqpWVtXKpUuXzkGlkqSt9h7wZx8P/OMkpwL7AT/NaOawJMnebTawDHik9d8IHAFsTLI3cADwvQHrkyRtY7CZQlX9ZlUtq6oVwFnALVX1K8DngdNbt3OAG9r22rZPe/yWqnreTEGSNJxpnKfwbuCiJBsYfWZwZWu/EnhJa78IuHgKtUnSojbk8lGnqm4Fbm3bDwDHztDnh8AZk6hHkjQzz2iWJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHV2GApJjk/yU237V5NcmuSlw5cmSZq0PjOFy4Bnk7yG0beofRu4atCqJElT0ScUNrfvNVgFfKCqPgDsP2xZkqRp6HPp7KeT/Cbwq8DPJdkL2GfYsiRJ09BnpnAm8BxwXlX9X+Bw4HcGrUqSNBV9Zgr/uqrevXWnqh5K8qoBa5IkTUmfmcI/nKHtlLkuRJI0fdudKSR5O/AO4Kgk94w9tD9w29CFSZImb7bloz8CPgv8F+Disfanq+p7g1YlSZqK2UKhqurBJOdv+0CSgwwGSVp4djRTeCtwF1BAxh4r4KgB65IkTcF2Q6Gq3truj5xcOZKkaepzSCpJ/ilwAqMZwher6lODViVJmoo+F8T7MPDrwL3AfcCvJ/nQ0IVJkiavz0zhHwA/065/RJI1jAJCkrTA9Dl57X5g+dj+EcA92+krSZrHZjt57U8YfYZwALA+yZ1t/4148pokLUizLR/914lVIUnaI8x2SOqftctk/2lVvXmCNUmSpmTWzxSqagujb107YEL1SJKmqM/RRz8E7k1yE/CDrY1VdcFgVUmSpqJPKHym3XZKkv2ALwD7tte5vqrek+RI4BrgIOArwNuq6kdJ9mX03c+vB74LnFlVD+7s60qSdt0OQ6Gq1iR5IfDy1nR/Vf2/Hj/7OeCkqnomyT7Anyf5LHAR8L6quibJ5cB5wGXt/omqelmSs4BLGH3rmyRpQvqc0Xwi8C3gQ8CHgW8m+bkdPa9Gnmm7+7RbAScB17f2NcBpbXtV26c9fnKS8YvwSZIG1mf56HeBn6+q+wGSvBz4OKNlnlm1o5fuAl7GKFT+N/BkVW1uXTYy+s5n2v3DAFW1OclTwEuAx3uPRpK0W/qc0bzP1kAAqKpvMnrXv0NVtaWqXgssA44FXjFTt3Y/06ygtm1IsjrJuiTrHnvssT5lSJJ66hMK65JcmeTEdvsDRu/+e6uqJ4FbgeOAJUm2zlCWAY+07Y2MLqFBe/wA4Hlf5FNVV1TVyqpauXTp0p0pQ5K0A31C4e3A14ALgAuBrzO6auqskixNsqRtvwh4M7Ae+Dxweut2DnBD217b9mmP37L1InySpMnoc/TRc8Cl7bYzDgPWtM8VXgBcV1WfTvJ14Jok/xn4KnBl638l8LEkGxjNEM7aydeTJO2mXl+ysyuq6h7gdTO0P8Do84Vt238InDFUPZKkHeuzfCRJWiS2GwpJPtbuL5xcOZKkaZptpvD6JC8F/kWSA5McNH6bVIGSpMmZ7TOFy4HPAUcxOgR1/DyCau2SpAVkuzOFqvpgVb0C+GhVHVVVR47dDARJWoD6HJL69iSvAX62NX2hHVkkSVpg+lwQ7wLgauCQdrs6yTuHLkySNHl9zlP4l8Abq+oHAEkuAb4E/N6QhUmSJq/PeQoBtoztb2Hmi9dJkua5PjOF/w7ckeSTbf80/vrSFJKkBaTPB82XJrkVOIHRDOHcqvrq0IVJkiav17WPquorjL5PWZK0gHntI0lSx1CQJHVmDYUkeyX5n5MqRpI0XbOGQlVtAZ5NcsCE6pEkTVGfD5p/CNyb5CbgB1sbq+qCwaqSJE1Fn1D4TLtJkha4PucprEnyImB5Vd0/gZokSVPS54J4/wi4m9F3K5DktUnWDl2YJGny+hyS+h+AY4EnAarqbuDIAWuSJE1Jn1DYXFVPbdNWQxQjSZquPh8035fkl4G9khwNXADcNmxZkqRp6DNTeCfwKuA54OPA94F3DVmUJGk6+hx99Czw79qX61RVPT18WZKkaehz9NEbktwL3MPoJLa/SPL64UuTJE1an88UrgTeUVVfBEhyAqMv3nn1kIVJkiavz2cKT28NBICq+nPAJSRJWoC2O1NIckzbvDPJ7zP6kLmAM4Fbhy9NkjRpsy0f/e42++8Z2/Y8BUlagLYbClX1pkkWIkmavh1+0JxkCXA2sGK8v5fOlqSFp8/RRzcCtwP3Aj8ethxJ0jT1CYX9quqinf3BSY4ArgL+JqMwuaKqPpDkIOBaRjOPB4F/XlVPJAnwAeBU4Fng16rqKzv7upKkXdfnkNSPJflXSQ5LctDWW4/nbQb+TVW9AjgOOD/JK4GLgZur6mjg5rYPcApwdLutBi7b2cFIknZPn1D4EfA7wJeAu9pt3Y6eVFWbtr7Tb5fGWA8cDqwC1rRua4DT2vYq4KoauR1YkuSwnRiLJGk39Vk+ugh4WVU9vqsvkmQF8DrgDuDQqtoEo+BIckjrdjjw8NjTNra2Tdv8rNWMZhIsX758V0uSJM2gz0zha4zW+HdJkhcDfwy8q6q+P1vXGdqedz5EVV1RVSurauXSpUt3tSxJ0gz6zBS2AHcn+Tyjy2cD/Q5JTbIPo0C4uqo+0Zq/k+SwNks4DHi0tW8Ejhh7+jLgkR71SZLmSJ9Q+FS77ZR2NNGVwPqqunTsobXAOcB72/0NY+2/keQa4I3AU1uXmSRJk9Hn+xTW7KjPdhwPvI3R5bbvbm2/xSgMrktyHvAQcEZ77EZGh6NuYLRcde4uvq4kaRf1OaP5L5l5bf+o2Z7XrqY60+cEACfP0L+A83dUjyRpOH2Wj1aObe/H6J19n/MUJEnzzA6PPqqq747d/k9VvR84aQK1SZImrM/y0TFjuy9gNHPYf7CKJElT02f5aPx7FTbTrlc0SDWSpKnqc/SR36sgSYtEn+WjfYF/xvO/T+G3hytLkjQNfZaPbgCeYnQhvOd20FeSNI/1CYVlVfWWwSuRJE1dnwvi3Zbk7w1eiSRp6vrMFE4Afq2d2fwco7OUq6pePWhlkqSJ6xMKpwxehSRpj9DnkNRvT6IQSdL09flMQZK0SBgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqTOYKGQ5KNJHk1y31jbQUluSvKtdn9ga0+SDybZkOSeJMcMVZckafuGnCn8IfCWbdouBm6uqqOBm9s+wCnA0e22GrhswLokSdsxWChU1ReA723TvApY07bXAKeNtV9VI7cDS5IcNlRtkqSZTfozhUOrahNAuz+ktR8OPDzWb2Nre54kq5OsS7LuscceG7RYSVps9pQPmjNDW83UsaquqKqVVbVy6dKlA5clSYvLpEPhO1uXhdr9o619I3DEWL9lwCMTrk2SFr1Jh8Ja4Jy2fQ5ww1j72e0opOOAp7YuM0mSJmfvoX5wko8DJwIHJ9kIvAd4L3BdkvOAh4AzWvcbgVOBDcCzwLlD1SVJ2r7BQqGqfmk7D508Q98Czh+qFklSP3vKB82SpD2AoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqTOHhUKSd6S5P4kG5JcPO16JGmx2WNCIclewIeAU4BXAr+U5JXTrUqSFpc9JhSAY4ENVfVAVf0IuAZYNeWaJGlR2XvaBYw5HHh4bH8j8MZtOyVZDaxuu88kuX8XX+9g4PFdfO585ZgXgVyy+MbMIvw9s3tjfun2HtiTQiEztNXzGqquAK7Y7RdL1lXVyt39OfOJY14cHPPiMNSY96Tlo43AEWP7y4BHplSLJC1Ke1IofBk4OsmRSV4InAWsnXJNkrSo7DHLR1W1OclvAH8K7AV8tKq+NuBL7vYS1DzkmBcHx7w4DDLmVD1v2V6StEjtSctHkqQpMxQkSZ0FHQpJjkjy+STrk3wtyYUz9EmSD7ZLa9yT5Jhp1DpXeo75V9pY70lyW5LXTKPWudJnzGN935BkS5LTJ1njXOs75iQnJrm79fmzSdc5l3r+2z4gyZ8k+YvW59xp1DpXkuyX5M6x8fzHGfrsm+Ta9jfsjiQrdutFq2rB3oDDgGPa9v7AN4FXbtPnVOCzjM6TOA64Y9p1T2DMfx84sG2fshjG3B7bC7gFuBE4fdp1T+D3vAT4OrC87R8y7bonMObfAi5p20uB7wEvnHbtuzHmAC9u2/sAdwDHbdPnHcDlbfss4Nrdec0FPVOoqk1V9ZW2/TSwntGZ0+NWAVfVyO3AkiSHTbjUOdNnzFV1W1U90XZvZ3ROyLzV8/cM8E7gj4FHJ1jeIHqO+ZeBT1TVQ63fvB53zzEXsH+SAC9mFAqbJ1roHGp/l55pu/u027ZHB60C1rTt64GT2/h3yYIOhXFtSvU6Rkk7bqbLa8z0B2XemWXM485jNFNaELY35iSHA/8EuHzyVQ1rlt/zy4EDk9ya5K4kZ0+6tqHMMub/BryC0Ymv9wIXVtWPJ1rcHEuyV5K7Gb2Zuamqtvs3rKo2A08BL9nV19tjzlMYUpIXM3qH+K6q+v62D8/wlHl/nO4Oxry1z5sYhcIJk6xtKDsY8/uBd1fVlt14E7XH2cGY9wZeD5wMvAj4UpLbq+qbEy5zTu1gzL8A3A2cBPxt4KYkX9ze/4H5oKq2AK9NsgT4ZJKfqar7xrrM6d+wBT9TSLIPo39AV1fVJ2bosuAur9FjzCR5NfARYFVVfXeS9Q2hx5hXAtckeRA4HfhwktMmWOKc6/lv+3NV9YOqehz4AjDfDyrY0ZjPZbRkVlW1AfhL4O9OssahVNWTwK3AW7Z5qPsblmRv4ABGy2a7ZEGHQltXuxJYX1WXbqfbWuDsdhTSccBTVbVpYkXOsT5jTrIc+ATwtvn+rhH6jbmqjqyqFVW1gtG66zuq6lMTLHNO9fy3fQPws0n2TvI3GF11eP2kapxrPcf8EKOZEUkOBf4O8MBkKpx7SZa2GQJJXgS8GfjGNt3WAue07dOBW6p96rwrFvry0fHA24B725ocjI5OWA5QVZczOhLlVGAD8CyjdxrzWZ8x/3tGa44fbkspm2t+X2Gyz5gXmh2OuarWJ/kccA/wY+Aj2yw7zDd9fs//CfjDJPcyWlZ5d5slzVeHAWsy+hKyFwDXVdWnk/w2sK6q1jIKyo8l2cBohnDW7rygl7mQJHUW9PKRJGnnGAqSpI6hIEnqGAqSpI6hIEnqGArSbkhy27RrkOaSh6RKkjrOFKTdkOSZdn9iu/Dc9Um+keTqrVeqbN/hcFu7Jv6dSfafbtXS9i30M5qlSXod8CpG1876X8DxSe4ErgXOrKovJ/lp4K+mWKM0K0NBmjt3VtVGgHYZhhWMLmO8qaq+DDCfr9apxcHlI2nuPDe2vYXRm66wAC7FrsXDUJCG9Q3gbyV5A0CS/dvljaU9kv84pQFV1Y+SnAn8Xrv08V8xuvzxM7M/U5oOD0mVJHVcPpIkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdf4/hybY4/pdp3UAAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ofti_sampler = sampler.OFTI(sys)\n", "\n", "# number of orbits to accept\n", "n_orbs = 500\n", "\n", "_ = ofti_sampler.run_sampler(n_orbs)\n", "\n", "plt.figure()\n", "accepted_eccentricities = ofti_sampler.results.post[:, lab['ecc1']]\n", "plt.hist(accepted_eccentricities)\n", "plt.xlabel('ecc'); plt.ylabel('number of orbits')\n", "\n", "plt.figure()\n", "accepted_inclinations = ofti_sampler.results.post[:, lab['inc1']]\n", "plt.hist(accepted_inclinations)\n", "plt.xlabel('inc'); plt.ylabel('number of orbits')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run MCMC" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting Burn in\n", "\n", "Burn in complete. Sampling posterior now.\n", "10/10 steps completed\n", "Run complete\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'number of orbits')" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG0CAYAAADU2ObLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuT0lEQVR4nO3df1xUdaL/8fcgv8xgEE1+uCh4s1BTJElCMy0pNW9p2S27rJHraj+0UrSSzZ+tBXlb9WouXs3Stsy1Ld02i0eKJncNydB+qdekKKkEtgwQDUQ594++zt35osXIDDN8fD0fj/N4OJ9z5viez2MezPtx5pwzNsuyLAEAABjKz9sBAAAAPImyAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACM5tWyk5+fr5tvvlnR0dGy2WzatGnTObe97777ZLPZtGTJEqfxo0ePKi0tTaGhoQoLC9OECRNUU1Pj2eAAAKDV8Pfmf378+HElJCToN7/5jW677bZzbrdx40bt2rVL0dHRjdalpaXpyJEj2rJli+rr6zV+/HhNmjRJ69ata3KOhoYGffvttwoJCZHNZjuv1wIAAFqWZVk6duyYoqOj5ef3M8dvLB8hydq4cWOj8a+//trq3Lmz9emnn1pdu3a1Fi9e7Fi3f/9+S5K1e/dux9jbb79t2Ww265tvvmny/11aWmpJYmFhYWFhYWmFS2lp6c9+znv1yM4vaWho0Lhx4/TII4+oV69ejdYXFBQoLCxMSUlJjrHU1FT5+fmpsLBQt95661n3W1dXp7q6Osdj6//98HtpaalCQ0Pd/CoAAIAnVFdXKyYmRiEhIT+7nU+Xnaefflr+/v566KGHzrq+rKxMnTp1chrz9/dXeHi4ysrKzrnfrKwszZ8/v9F4aGgoZQcAgFbml05B8dmrsYqKivSf//mfWrNmjdvPo8nMzFRVVZVjKS0tdev+AQCA7/DZsvPf//3fqqioUJcuXeTv7y9/f3999dVXmj59umJjYyVJkZGRqqiocHreqVOndPToUUVGRp5z30FBQY6jOBzNAQDAbD77Nda4ceOUmprqNDZs2DCNGzdO48ePlySlpKSosrJSRUVF6tevnyRp27ZtamhoUHJycotnBgAAvserZaempkbFxcWOxyUlJfrwww8VHh6uLl26qEOHDk7bBwQEKDIyUpdffrkkqUePHho+fLgmTpyoFStWqL6+XlOmTNHYsWPPepk6AAC48Hj1a6wPPvhAiYmJSkxMlCRlZGQoMTFRc+bMafI+Xn75ZcXHx2vo0KG66aabdM0112jlypWeigwAAFoZm3XmuusLWHV1tex2u6qqqjh/BwCAVqKpn98+e4IyAACAO1B2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACj+ewPgQKAu8TO3OztCC77MnuktyMAxuDIDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEbz93YAAEBjsTM3ezuCy77MHuntCMBZcWQHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEbzatnJz8/XzTffrOjoaNlsNm3atMmxrr6+Xo899ph69+6tdu3aKTo6Wnfffbe+/fZbp30cPXpUaWlpCg0NVVhYmCZMmKCampoWfiUAAMBXebXsHD9+XAkJCVq+fHmjdSdOnNCePXs0e/Zs7dmzR6+//roOHjyoW265xWm7tLQ07du3T1u2bNGbb76p/Px8TZo0qaVeAgAA8HE2y7Isb4eQJJvNpo0bN2r06NHn3Gb37t3q37+/vvrqK3Xp0kUHDhxQz549tXv3biUlJUmScnNzddNNN+nrr79WdHT0WfdTV1enuro6x+Pq6mrFxMSoqqpKoaGhbn1dALyvNf6CeGvEr56jpVVXV8tut//i53erOmenqqpKNptNYWFhkqSCggKFhYU5io4kpaamys/PT4WFhefcT1ZWlux2u2OJiYnxdHQAAOAlrabs1NbW6rHHHtNdd93laG9lZWXq1KmT03b+/v4KDw9XWVnZOfeVmZmpqqoqx1JaWurR7AAAwHv8vR2gKerr63XHHXfIsizl5OQ0e39BQUEKCgpyQzIAAODrfL7snCk6X331lbZt2+b0nVxkZKQqKiqctj916pSOHj2qyMjIlo4KAAB8kE9/jXWm6Bw6dEhbt25Vhw4dnNanpKSosrJSRUVFjrFt27apoaFBycnJLR0XAAD4IK8e2ampqVFxcbHjcUlJiT788EOFh4crKipKt99+u/bs2aM333xTp0+fdpyHEx4ersDAQPXo0UPDhw/XxIkTtWLFCtXX12vKlCkaO3bsOa/EAgAAFxavlp0PPvhA1113neNxRkaGJCk9PV3z5s3TG2+8IUnq27ev0/O2b9+uIUOGSJJefvllTZkyRUOHDpWfn5/GjBmjpUuXtkh+AADg+7xadoYMGaKfu81PU24BFB4ernXr1rkzFgAAMIhPn7MDAADQXJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEaj7AAAAKNRdgAAgNEoOwAAwGiUHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEaj7AAAAKNRdgAAgNEoOwAAwGheLTv5+fm6+eabFR0dLZvNpk2bNjmttyxLc+bMUVRUlNq2bavU1FQdOnTIaZujR48qLS1NoaGhCgsL04QJE1RTU9OCrwIAAPgyr5ad48ePKyEhQcuXLz/r+oULF2rp0qVasWKFCgsL1a5dOw0bNky1tbWObdLS0rRv3z5t2bJFb775pvLz8zVp0qSWegkAAMDH+XvzPx8xYoRGjBhx1nWWZWnJkiWaNWuWRo0aJUl68cUXFRERoU2bNmns2LE6cOCAcnNztXv3biUlJUmSli1bpptuuknPPPOMoqOjW+y1AAAA3+Sz5+yUlJSorKxMqampjjG73a7k5GQVFBRIkgoKChQWFuYoOpKUmpoqPz8/FRYWnnPfdXV1qq6udloAAICZfLbslJWVSZIiIiKcxiMiIhzrysrK1KlTJ6f1/v7+Cg8Pd2xzNllZWbLb7Y4lJibGzekBAICv8Nmy40mZmZmqqqpyLKWlpd6OBAAAPMRny05kZKQkqby83Gm8vLzcsS4yMlIVFRVO60+dOqWjR486tjmboKAghYaGOi0AAMBMPlt24uLiFBkZqby8PMdYdXW1CgsLlZKSIklKSUlRZWWlioqKHNts27ZNDQ0NSk5ObvHMAADA93j1aqyamhoVFxc7HpeUlOjDDz9UeHi4unTpoqlTp2rBggXq3r274uLiNHv2bEVHR2v06NGSpB49emj48OGaOHGiVqxYofr6ek2ZMkVjx47lSiwAACDJy2Xngw8+0HXXXed4nJGRIUlKT0/XmjVr9Oijj+r48eOaNGmSKisrdc011yg3N1fBwcGO57z88suaMmWKhg4dKj8/P40ZM0ZLly5t8dcCAAB8k82yLMvbIbyturpadrtdVVVVnL8DGCh25mZvR7ggfJk90tsRcIFp6ue3z56zAwAA4A6UHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGa3bZqa6u1qZNm3TgwAF35AEAAHArl8vOHXfcoWeffVaS9OOPPyopKUl33HGH+vTpo9dee83tAQEAAJrD5bKTn5+vQYMGSZI2btwoy7JUWVmppUuXasGCBW4PCAAA0Bwul52qqiqFh4dLknJzczVmzBhddNFFGjlypA4dOuT2gAAAAM3hctmJiYlRQUGBjh8/rtzcXN14442SpB9++EHBwcFuDwgAANAc/q4+YerUqUpLS9PFF1+srl27asiQIZJ++nqrd+/e7s4HAADQLC6XnQceeEDJyck6fPiwbrjhBvn5/XRwqFu3bnryySfdHhAAAKA5XP4a64knnlCPHj1066236uKLL3aMX3/99dq6datbwwEAADSXy2Vn/vz5qqmpaTR+4sQJzZ8/3y2hAAAA3MXlsmNZlmw2W6Pxjz76yHGVFgAAgK9o8jk77du3l81mk81m02WXXeZUeE6fPq2amhrdd999HgkJAABwvppcdpYsWSLLsvSb3/xG8+fPl91ud6wLDAxUbGysUlJSPBISAADgfDW57KSnp0uS4uLiNGDAAAUEBHgsFAAAgLs0qexUV1crNDRUkpSYmKgff/xRP/7441m3PbMdAACAL2hS2Wnfvr2OHDmiTp06KSws7KwnKJ85cfn06dNuDwkAAHC+mlR2tm3b5rjSavv27R4NBAAA4E5NKjuDBw8+678BAAB8ncs/FyH99KOfq1ev1oEDByRJPXv21Pjx47nPDgAA8Dku31QwPz9fsbGxWrp0qX744Qf98MMPWrp0qeLi4pSfn++JjAAAAOfN5SM7kydP1p133qmcnBy1adNG0k83FXzggQc0efJkffLJJ24PCQAAcL5cPrJTXFys6dOnO4qOJLVp00YZGRkqLi52azgAAIDmcrnsXHnllY5zdf7ZgQMHlJCQ4JZQAAAA7tKkr7E+/vhjx78feughPfzwwyouLtbVV18tSdq1a5eWL1+u7Oxsz6QEAAA4TzbLsqxf2sjPz082m02/tGlrvalgdXW17Ha7qqqquAM0YKDYmZu9HeGC8GX2SG9HwAWmqZ/fTTqyU1JS4rZgAAAALalJZadr166SpPr6et17772aPXu24uLiPBoMAADAHVw6QTkgIECvvfaap7IAAAC4nctXY40ePVqbNm3yQBQAAAD3c/mmgt27d9cTTzyhnTt3ql+/fmrXrp3T+oceesht4U6fPq158+bppZdeUllZmaKjo3XPPfdo1qxZjl9etyxLc+fO1apVq1RZWamBAwcqJydH3bt3d1sOAADQerlcdlavXq2wsDAVFRWpqKjIaZ3NZnNr2Xn66aeVk5OjtWvXqlevXvrggw80fvx42e12x/+zcOFCLV26VGvXrlVcXJxmz56tYcOGaf/+/QoODnZbFgAA0Dq5XHZa8sqs9957T6NGjdLIkT9dzhgbG6tXXnlF77//vqSfjuosWbJEs2bN0qhRoyRJL774oiIiIrRp0yaNHTu2xbICAADf5PI5O//MsqxfvPdOcwwYMEB5eXn67LPPJEkfffSR/v73v2vEiBGSfipeZWVlSk1NdTzHbrcrOTlZBQUF59xvXV2dqqurnRYAAGCm8yo7L774onr37q22bduqbdu26tOnj/70pz+5O5tmzpypsWPHKj4+XgEBAUpMTNTUqVOVlpYmSSorK5MkRUREOD0vIiLCse5ssrKyZLfbHUtMTIzbswMAAN/gctlZtGiR7r//ft10003asGGDNmzYoOHDh+u+++7T4sWL3Rpuw4YNevnll7Vu3Trt2bNHa9eu1TPPPKO1a9c2a7+ZmZmqqqpyLKWlpW5KDAAAfI3L5+wsW7ZMOTk5uvvuux1jt9xyi3r16qV58+Zp2rRpbgv3yCOPOI7uSFLv3r311VdfKSsrS+np6YqMjJQklZeXKyoqyvG88vJy9e3b95z7DQoKUlBQkNtyAgAA3+XykZ0jR45owIABjcYHDBigI0eOuCXUGSdOnJCfn3PENm3aqKGhQZIUFxenyMhI5eXlOdZXV1ersLBQKSkpbs0CAABaJ5fLzqWXXqoNGzY0Gv/zn//s9nvb3HzzzXryySe1efNmffnll9q4caMWLVqkW2+9VdJPl7pPnTpVCxYs0BtvvKFPPvlEd999t6KjozV69Gi3ZgEAAK2Ty19jzZ8/X3feeafy8/M1cOBASdLOnTuVl5d31hLUHMuWLdPs2bP1wAMPqKKiQtHR0br33ns1Z84cxzaPPvqojh8/rkmTJqmyslLXXHONcnNzuccOAACQJNms87h2vKioSIsXL9aBAwckST169ND06dOVmJjo9oAtoak/EQ+gdYqdudnbES4IX2aP9HYEXGCa+vnt8pEdSerXr59eeuml8w4HAADQUpp1U0EAAABfR9kBAABGo+wAAACjNansfPzxx4572wAAALQmTSo7iYmJ+u677yRJ3bp10/fff+/RUAAAAO7SpKuxwsLCVFJSok6dOunLL7/kKA/gBq31cmguLwbQ2jSp7IwZM0aDBw9WVFSUbDabkpKS1KZNm7Nu+8UXX7g1IAAAQHM0qeysXLlSt912m4qLi/XQQw9p4sSJCgkJ8XQ2AACAZmvyTQWHDx8u6ae7Jz/88MOUHQAA0Cq4fAflF154wfHvr7/+WpL0q1/9yn2JAAAA3MjlstPQ0KAFCxboD3/4g2pqaiRJISEhmj59uh5//HH5+XHrHgC4ELXGk+454f7C4HLZefzxx7V69WplZ2c7fvX873//u+bNm6fa2lo9+eSTbg8JAABwvlwuO2vXrtVzzz2nW265xTHWp08fde7cWQ888ABlBwAA+BSXv3M6evSo4uPjG43Hx8fr6NGjbgkFAADgLi6XnYSEBD377LONxp999lklJCS4JRQAAIC7uPw11sKFCzVy5Eht3bpVKSkpkqSCggKVlpbqrbfecntAAACA5nD5yM7gwYP12Wef6dZbb1VlZaUqKyt122236eDBgxo0aJAnMgIAAJw3l4/sSFJ0dDQnIgMAgFaBm+IAAACjUXYAAIDRKDsAAMBoLpUdy7J0+PBh1dbWeioPAACAW7lcdi699FKVlpZ6Kg8AAIBbuVR2/Pz81L17d33//feeygMAAOBWLp+zk52drUceeUSffvqpJ/IAAAC4lcv32bn77rt14sQJJSQkKDAwUG3btnVaz+9jAQAAX+Jy2VmyZIkHYgAAAHiGy2UnPT3dEzkAAAA84rzus/P5559r1qxZuuuuu1RRUSFJevvtt7Vv3z63hgMAAGgul8vOjh071Lt3bxUWFur1119XTU2NJOmjjz7S3Llz3R4QAACgOVwuOzNnztSCBQu0ZcsWBQYGOsavv/567dq1y63hAAAAmsvlsvPJJ5/o1ltvbTTeqVMnfffdd24JBQAA4C4ul52wsDAdOXKk0fjevXvVuXNnt4QCAABwF5fLztixY/XYY4+prKxMNptNDQ0N2rlzp2bMmKG7777bExkBAADOm8tl56mnnlJ8fLxiYmJUU1Ojnj176tprr9WAAQM0a9YsT2QEAAA4by7fZycwMFCrVq3S7Nmz9emnn6qmpkaJiYnq3r27J/IBAAA0i8tl54wuXbooJiZGkmSz2dwWCAAAwJ3O66aCq1ev1hVXXKHg4GAFBwfriiuu0HPPPefubAAAAM3m8pGdOXPmaNGiRXrwwQeVkpIiSSooKNC0adN0+PBhPfHEE24PCQAAcL5cLjs5OTlatWqV7rrrLsfYLbfcoj59+ujBBx+k7AAAAJ/i8tdY9fX1SkpKajTer18/nTp1yi2h/tk333yjX//61+rQoYPatm2r3r1764MPPnCstyxLc+bMUVRUlNq2bavU1FQdOnTI7TkAAEDr5HLZGTdunHJychqNr1y5UmlpaW4JdcYPP/yggQMHKiAgQG+//bb279+vP/zhD2rfvr1jm4ULF2rp0qVasWKFCgsL1a5dOw0bNky1tbVuzQIAAFqnJn2NlZGR4fi3zWbTc889p3feeUdXX321JKmwsFCHDx92+00Fn376acXExOiFF15wjMXFxTn+bVmWlixZolmzZmnUqFGSpBdffFERERHatGmTxo4d69Y8AACg9WlS2dm7d6/T4379+kmSPv/8c0lSx44d1bFjR+3bt8+t4d544w0NGzZM//Zv/6YdO3aoc+fOeuCBBzRx4kRJUklJicrKypSamup4jt1uV3JysgoKCs5Zdurq6lRXV+d4XF1d7dbcAADAdzSp7Gzfvt3TOc7qiy++UE5OjjIyMvS73/1Ou3fv1kMPPaTAwEClp6errKxMkhQREeH0vIiICMe6s8nKytL8+fM9mh0wVezMzd6OAAAuOa/77LSUhoYGXXnllXrqqaeUmJioSZMmaeLEiVqxYkWz9puZmamqqirHUlpa6qbEAADA17h86Xltba2WLVum7du3q6KiQg0NDU7r9+zZ47ZwUVFR6tmzp9NYjx499Nprr0mSIiMjJUnl5eWKiopybFNeXq6+ffuec79BQUEKCgpyW04AAOC7XC47EyZM0DvvvKPbb79d/fv39+hPRQwcOFAHDx50Gvvss8/UtWtXST+drBwZGam8vDxHuamurlZhYaHuv/9+j+UCAACth8tl580339Rbb72lgQMHeiKPk2nTpmnAgAF66qmndMcdd+j999/XypUrtXLlSkk/XRk2depULViwQN27d1dcXJxmz56t6OhojR492uP5AACA73O57HTu3FkhISGeyNLIVVddpY0bNyozM1NPPPGE4uLitGTJEqf7+Tz66KM6fvy4Jk2apMrKSl1zzTXKzc1VcHBwi2QEAAC+zWZZluXKE95++23HTfzOfJ3U2lVXV8tut6uqqkqhoaHejoMLBFc1Ad73ZfZIb0dAMzT189vlIztJSUmqra1Vt27ddNFFFykgIMBp/dGjR11PCwAA4CEul5277rpL33zzjZ566ilFRER49ARlAACA5nK57Lz33nsqKChQQkKCJ/IAAAC4lcs3FYyPj9ePP/7oiSwAAABu53LZyc7O1vTp0/Xuu+/q+++/V3V1tdMCAADgS1z+Gmv48OGSpKFDhzqNW5Ylm82m06dPuycZAACAG7hcdrz1o6AAAADnw+WyM3jwYE/kAAAA8AiXy05+fv7Prr/22mvPOwwAAIC7uVx2hgwZ0mjsn++1wzk7AADAl7h8NdYPP/zgtFRUVCg3N1dXXXWV3nnnHU9kBAAAOG8uH9mx2+2Nxm644QYFBgYqIyNDRUVFbgkGAADgDi4f2TmXiIgIHTx40F27AwAAcAuXj+x8/PHHTo8ty9KRI0eUnZ2tvn37uisXAACAW7hcdvr27SubzSbLspzGr776aj3//PNuCwYAAOAOLpedkpISp8d+fn665JJLFBwc7LZQAAAA7uJy2enatasncgAAAHiEy2VHkvLy8pSXl6eKigo1NDQ4reOrLAAA4EtcLjvz58/XE088oaSkJEVFRTndUBAAAMDXuFx2VqxYoTVr1mjcuHGeyAMAAOBWLt9n5+TJkxowYIAnsgAAALidy2Xnt7/9rdatW+eJLAAAAG7n8tdYtbW1WrlypbZu3ao+ffooICDAaf2iRYvcFg4AAKC5zusOymfulPzpp586reNkZQAA4GtcLjvbt2/3RA4AAACPcNsPgQIAAPgiyg4AADAaZQcAABiNsgMAAIxG2QEAAEaj7AAAAKNRdgAAgNEoOwAAwGiUHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgtFZVdrKzs2Wz2TR16lTHWG1trSZPnqwOHTro4osv1pgxY1ReXu69kAAAwKe0mrKze/du/dd//Zf69OnjND5t2jT97W9/06uvvqodO3bo22+/1W233eallAAAwNe0irJTU1OjtLQ0rVq1Su3bt3eMV1VVafXq1Vq0aJGuv/569evXTy+88ILee+897dq1y4uJAQCAr2gVZWfy5MkaOXKkUlNTncaLiopUX1/vNB4fH68uXbqooKDgnPurq6tTdXW10wIAAMzk7+0Av2T9+vXas2ePdu/e3WhdWVmZAgMDFRYW5jQeERGhsrKyc+4zKytL8+fPd3dUAADgg3z6yE5paakefvhhvfzyywoODnbbfjMzM1VVVeVYSktL3bZvAADgW3y67BQVFamiokJXXnml/P395e/vrx07dmjp0qXy9/dXRESETp48qcrKSqfnlZeXKzIy8pz7DQoKUmhoqNMCAADM5NNfYw0dOlSffPKJ09j48eMVHx+vxx57TDExMQoICFBeXp7GjBkjSTp48KAOHz6slJQUb0QGAAA+xqfLTkhIiK644gqnsXbt2qlDhw6O8QkTJigjI0Ph4eEKDQ3Vgw8+qJSUFF199dXeiAwAAHyMT5edpli8eLH8/Pw0ZswY1dXVadiwYfrjH//o7VgAAMBH2CzLsrwdwtuqq6tlt9tVVVXF+TtoMbEzN3s7AnDB+zJ7pLcjoBma+vnt0ycoAwAANBdlBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEaj7AAAAKNRdgAAgNEoOwAAwGiUHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADCaT5edrKwsXXXVVQoJCVGnTp00evRoHTx40Gmb2tpaTZ48WR06dNDFF1+sMWPGqLy83EuJAQCAr/HpsrNjxw5NnjxZu3bt0pYtW1RfX68bb7xRx48fd2wzbdo0/e1vf9Orr76qHTt26Ntvv9Vtt93mxdQAAMCX+Hs7wM/Jzc11erxmzRp16tRJRUVFuvbaa1VVVaXVq1dr3bp1uv766yVJL7zwgnr06KFdu3bp6quv9kZsAADgQ3z6yM7/r6qqSpIUHh4uSSoqKlJ9fb1SU1Md28THx6tLly4qKCg4537q6upUXV3ttAAAADO1mrLT0NCgqVOnauDAgbriiiskSWVlZQoMDFRYWJjTthERESorKzvnvrKysmS32x1LTEyMJ6MDAAAvajVlZ/Lkyfr000+1fv36Zu8rMzNTVVVVjqW0tNQNCQEAgC/y6XN2zpgyZYrefPNN5efn61e/+pVjPDIyUidPnlRlZaXT0Z3y8nJFRkaec39BQUEKCgryZGQAAOAjfPrIjmVZmjJlijZu3Kht27YpLi7OaX2/fv0UEBCgvLw8x9jBgwd1+PBhpaSktHRcAADgg3z6yM7kyZO1bt06/fWvf1VISIjjPBy73a62bdvKbrdrwoQJysjIUHh4uEJDQ/Xggw8qJSWFK7EAAIAkHy87OTk5kqQhQ4Y4jb/wwgu65557JEmLFy+Wn5+fxowZo7q6Og0bNkx//OMfWzgpAADwVT5ddizL+sVtgoODtXz5ci1fvrwFEgEAgNbGp8/ZAQAAaC7KDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBolB0AAGA0f28HANwhduZmb0cA0Aq1xr8dX2aP9HaEVocjOwAAwGiUHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAwGmUHAAAYjbIDAACMRtkBAABGo+wAAACjUXYAAIDRKDsAAMBo/t4OYLrYmZu9HcFlX2aP9HYEAMA58LniOo7sAAAAo1F2AACA0Sg7AADAaJQdAABgNE5QRiOt8eQ3AADOhSM7AADAaMaUneXLlys2NlbBwcFKTk7W+++/7+1IAADABxhRdv785z8rIyNDc+fO1Z49e5SQkKBhw4apoqLC29EAAICXGVF2Fi1apIkTJ2r8+PHq2bOnVqxYoYsuukjPP/+8t6MBAAAva/UnKJ88eVJFRUXKzMx0jPn5+Sk1NVUFBQVnfU5dXZ3q6uocj6uqqiRJ1dXVbs/XUHfC7fsEAKA18cTn6z/v17Ksn92u1Zed7777TqdPn1ZERITTeEREhP7nf/7nrM/JysrS/PnzG43HxMR4JCMAABcy+xLP7v/YsWOy2+3nXN/qy875yMzMVEZGhuNxZWWlunbtqsOHD//sZOH8VFdXKyYmRqWlpQoNDfV2HOMwv57F/HoW8+s5F8LcWpalY8eOKTo6+me3a/Vlp2PHjmrTpo3Ky8udxsvLyxUZGXnW5wQFBSkoKKjRuN1uN/YN4QtCQ0OZXw9ifj2L+fUs5tdzTJ/bphykaPUnKAcGBqpfv37Ky8tzjDU0NCgvL08pKSleTAYAAHxBqz+yI0kZGRlKT09XUlKS+vfvryVLluj48eMaP368t6MBAAAvM6Ls3HnnnfrHP/6hOXPmqKysTH379lVubm6jk5bPJSgoSHPnzj3rV1toPubXs5hfz2J+PYv59Rzm9v/YrF+6XgsAAKAVa/Xn7AAAAPwcyg4AADAaZQcAABiNsgMAAIxmbNlZvny5YmNjFRwcrOTkZL3//vs/u/2rr76q+Ph4BQcHq3fv3nrrrbec1t9zzz2y2WxOy/Dhwz35EnyaK/O7b98+jRkzRrGxsbLZbFqyZEmz92k6d8/vvHnzGr1/4+PjPfgKfJcrc7tq1SoNGjRI7du3V/v27ZWamtpoe8uyNGfOHEVFRalt27ZKTU3VoUOHPP0yfJa755e/vc5cmd/XX39dSUlJCgsLU7t27dS3b1/96U9/ctrmgnn/WgZav369FRgYaD3//PPWvn37rIkTJ1phYWFWeXn5WbffuXOn1aZNG2vhwoXW/v37rVmzZlkBAQHWJ5984tgmPT3dGj58uHXkyBHHcvTo0ZZ6ST7F1fl9//33rRkzZlivvPKKFRkZaS1evLjZ+zSZJ+Z37ty5Vq9evZzev//4xz88/Ep8j6tz++///u/W8uXLrb1791oHDhyw7rnnHstut1tff/21Y5vs7GzLbrdbmzZtsj766CPrlltuseLi4qwff/yxpV6Wz/DE/PK39/+4Or/bt2+3Xn/9dWv//v1WcXGxtWTJEqtNmzZWbm6uY5sL5f1rZNnp37+/NXnyZMfj06dPW9HR0VZWVtZZt7/jjjuskSNHOo0lJydb9957r+Nxenq6NWrUKI/kbW1cnd9/1rVr17N+GDdnn6bxxPzOnTvXSkhIcGPK1qm577NTp05ZISEh1tq1ay3LsqyGhgYrMjLS+o//+A/HNpWVlVZQUJD1yiuvuDd8K+Du+bUs/vb+M3f8nUxMTLRmzZplWdaF9f417muskydPqqioSKmpqY4xPz8/paamqqCg4KzPKSgocNpekoYNG9Zo+3fffVedOnXS5Zdfrvvvv1/ff/+9+1+Ajzuf+fXGPlsrT87FoUOHFB0drW7duiktLU2HDx9ubtxWxR1ze+LECdXX1ys8PFySVFJSorKyMqd92u12JScn895V8+f3DP72Nn9+LctSXl6eDh48qGuvvVbShfX+Na7sfPfddzp9+nSjuydHRESorKzsrM8pKyv7xe2HDx+uF198UXl5eXr66ae1Y8cOjRgxQqdPn3b/i/Bh5zO/3thna+WpuUhOTtaaNWuUm5urnJwclZSUaNCgQTp27FhzI7ca7pjbxx57TNHR0Y4PhzPP473rmfmV+Nt7xvnOb1VVlS6++GIFBgZq5MiRWrZsmW644QZJF9b714ifi2gJY8eOdfy7d+/e6tOnj/7lX/5F7777roYOHerFZMAvGzFihOPfffr0UXJysrp27aoNGzZowoQJXkzWemRnZ2v9+vV69913FRwc7O04xjnX/PK3t3lCQkL04YcfqqamRnl5ecrIyFC3bt00ZMgQb0drUcYd2enYsaPatGmj8vJyp/Hy8nJFRkae9TmRkZEubS9J3bp1U8eOHVVcXNz80K3I+cyvN/bZWrXUXISFhemyyy67oN6/zZnbZ555RtnZ2XrnnXfUp08fx/iZ5/He9cz8ng1/e12bXz8/P1166aXq27evpk+frttvv11ZWVmSLqz3r3FlJzAwUP369VNeXp5jrKGhQXl5eUpJSTnrc1JSUpy2l6QtW7acc3tJ+vrrr/X9998rKirKPcFbifOZX2/ss7VqqbmoqanR559/fkG9f893bhcuXKjf//73ys3NVVJSktO6uLg4RUZGOu2zurpahYWFvHfV/Pk9G/72Nu9vQ0NDg+rq6iRdYO9fb58h7Qnr16+3goKCrDVr1lj79++3Jk2aZIWFhVllZWWWZVnWuHHjrJkzZzq237lzp+Xv728988wz1oEDB6y5c+c6XXp+7Ngxa8aMGVZBQYFVUlJibd261bryyiut7t27W7W1tV55jd7k6vzW1dVZe/futfbu3WtFRUVZM2bMsPbu3WsdOnSoyfu8kHhifqdPn269++67VklJibVz504rNTXV6tixo1VRUdHir8+bXJ3b7OxsKzAw0PrLX/7idOnzsWPHnLYJCwuz/vrXv1off/yxNWrUKCMv3W0Kd88vf3uduTq/Tz31lPXOO+9Yn3/+ubV//37rmWeesfz9/a1Vq1Y5trlQ3r9Glh3Lsqxly5ZZXbp0sQIDA63+/ftbu3btcqwbPHiwlZ6e7rT9hg0brMsuu8wKDAy0evXqZW3evNmx7sSJE9aNN95oXXLJJVZAQIDVtWtXa+LEiRfkB/EZrsxvSUmJJanRMnjw4Cbv80Lj7vm98847raioKCswMNDq3Lmzdeedd1rFxcUt+Ip8hytz27Vr17PO7dy5cx3bNDQ0WLNnz7YiIiKsoKAga+jQodbBgwdb8BX5FnfOL397G3Nlfh9//HHr0ksvtYKDg6327dtbKSkp1vr16532d6G8f22WZVkteywJAACg5Rh3zg4AAMA/o+wAAACjUXYAAIDRKDsAAMBolB0AAGA0yg4AADAaZQcAABiNsgMAAIxG2QEAAEaj7AAAAKNRdgAAgNEoOwBarYaGBmVlZSkuLk5t27ZVQkKC/vKXvzjW79u3T//6r/+q0NBQhYSEaNCgQfr8888d659//nn16tVLQUFBioqK0pQpU7zxMgB4mL+3AwDA+crKytJLL72kFStWqHv37srPz9evf/1rXXLJJbr00kt17bXXasiQIdq2bZtCQ0O1c+dOnTp1SpKUk5OjjIwMZWdna8SIEaqqqtLOnTu9/IoAeAK/eg6gVaqrq1N4eLi2bt2qlJQUx/hvf/tbnThxQrGxsVq/fr0OHjyogICARs/v3Lmzxo8frwULFrRkbABewJEdAK1ScXGxTpw4oRtuuMFp/OTJk0pMTFRlZaUGDRp01qJTUVGhb7/9VkOHDm2puAC8iLIDoFWqqamRJG3evFmdO3d2WhcUFKSpU6ee87lt27b1ZDQAPoayA6BV6tmzp4KCgnT48GENHjy40fo+ffpo7dq1qq+vb3R0JyQkRLGxscrLy9N1113XUpEBeAllB0CrFBISohkzZmjatGlqaGjQNddc4zjJODQ0VFOmTNGyZcs0duxYZWZmym63a9euXerfv78uv/xyzZs3T/fdd586deqkESNG6NixY9q5c6cefPBBb780AG5G2QHQav3+97/XJZdcoqysLH3xxRcKCwvTlVdeqd/97nfq0KGDtm3bpkceeUSDBw9WmzZt1LdvXw0cOFCSlJ6ertraWi1evFgzZsxQx44ddfvtt3v5FQHwBK7GAgAARuOmggAAwGiUHQAAYDTKDgAAMBplBwAAGI2yAwAAjEbZAQAARqPsAAAAo1F2AACA0Sg7AADAaJQdAABgNMoOAAAw2v8CWFmBPMSoDSAAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# number of temperatures & walkers for MCMC\n", "num_temps = 3\n", "num_walkers = 50\n", "\n", "# number of steps to take\n", "n_orbs = 500\n", "\n", "mcmc_sampler = sampler.MCMC(sys, num_temps, num_walkers)\n", "\n", "# number of orbits to accept\n", "n_orbs = 500\n", "\n", "_ = mcmc_sampler.run_sampler(n_orbs)\n", "\n", "accepted_eccentricities = mcmc_sampler.results.post[:, lab['ecc1']]\n", "plt.hist(accepted_eccentricities)\n", "plt.xlabel('ecc'); plt.ylabel('number of orbits')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The MCMC sampler will automatically take the priors from the system used to instantiate it. However, updates to the system priors made after the sampler instantiation will not be propogated to the sampler, and thus will not be used when running MCMC. The recommended practice when modifying priors is to first instantiate and modify system priors, then instantiate the sampler. When modifying priors we also recommend instantiating the system and sampler objects individually, rather than instantiating a driver object." ] } ], "metadata": { "interpreter": { "hash": "e899b22145868d3cd465733d82c36c2ae3ac0d3591d6a0807ec2e5e577a9cf5c" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 2 }