Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $orvara$ tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Please make sure you are running this notebook from within the root directory!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34mbuild\u001b[0m/                    HGCA_vEDR3.fits  \u001b[01;34morvara\u001b[0m/           README.rst\r\n",
      "CHANGES.md                LICENCE          \u001b[01;34morvara.egg-info\u001b[0m/  setup.py\r\n",
      "HGCA_vDR2_corrected.fits  \u001b[01;34mmisc\u001b[0m/            pytest.ini        Tutorial.ipynb\r\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The above verifies that you are in the root directory.\n",
    "You should see something like:\n",
    "\n",
    "\n",
    "    build/                    HGCA_vEDR3.fits  orvara/           README.rst\n",
    "    CHANGES.md                LICENCE          orvara.egg-info/  setup.py\n",
    "    HGCA_vDR2_corrected.fits  misc/            pytest.ini        Tutorial.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now that you are in the root directory, lets install orvara.\n",
    "\n",
    "### Warning about path names:\n",
    "DO NOT run orbit fits from paths that contain files with spaces. This can cause python to fail at finding various input files. So if you install orvara and run from:\n",
    "\n",
    "    /user/your name/Research Projects/orvara/\n",
    "That may cause failures. But:\n",
    "\n",
    "    /user/yourname/ResearchProjects/orvara/\n",
    "\n",
    "is OK. If you are stuck with path names with spaces: You can always try first to see if it works with your directory paths, if you have spaces involved. But if you get a wierd error (especially when you the orbit fit has yet to actually start and is only trying to load in data), your path names containing spaces are probably the issue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Note: there is no need to run this following cell if you already have orvara installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -e ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you see: \"Successfully installed orvara-X.X.X\"\n",
    "\n",
    "With X.X.X being the version number,\n",
    "Great!\n",
    "\n",
    "## if you installed orvara with the above cell, then restart the Kernal in this jupyter notebook, then continue with the next cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run pytest to make sure that the tests pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will run a small suite of tests. This should take < 30 seconds. It will complete with a few warnings -- that is OK.\n",
    "\n",
    "**If this cell hangs for a long time (more than 3 minutes)-- please skip this cell and continue.**\n",
    "\n",
    "If any of the tests fail, something is wrong with the install. Try running \"pip install -e .\" again If the issue persists, please submit an issue ticket!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1m============================= test session starts ==============================\u001b[0m\n",
      "platform linux -- Python 3.8.10, pytest-6.2.3, py-1.10.0, pluggy-0.13.1 -- /home/gmbrandt/envs/orvara_new/bin/python\n",
      "cachedir: .pytest_cache\n",
      "rootdir: /home/gmbrandt/Documents/Repositories/orvara, configfile: pytest.ini\n",
      "collected 10 items / 4 deselected / 6 selected                                 \u001b[0m\u001b[1m\n",
      "\n",
      "orvara/tests/test_config.py::test_example_config_file \u001b[32mPASSED\u001b[0m\n",
      "orvara/tests/test_main.py::test_set_initial_parameters \u001b[32mPASSED\u001b[0m\n",
      "orvara/tests/test_main.py::test_get_priors \u001b[32mPASSED\u001b[0m\n",
      "orvara/tests/test_main.py::test_run Unable to load relative RV data from file \n",
      "Loadied 630 RV data points from file orvara/tests/data/Gl758_RV.dat\n",
      "Loaded data from 3 RV instruments.\n",
      "Loading astrometric data for 1 planets\n",
      "Loaded 4 relative astrometric data points from file orvara/tests/data/Gl758_relAST.txt\n",
      "Loading absolute astrometry data for Hip 95319\n",
      "Recognized a 5-parameter fit in Gaia for Hip 95319\n",
      "Not using companion proper motion from Gaia.\n",
      "Using ptemcee.\n",
      "Running MCMC ... \n",
      "[###############################]100%50%\n",
      "Total Time: 2 seconds\n",
      "Mean acceptance fraction (cold chain): 0.174100\n",
      "Writing output to /tmp/tmp5pzv4df9/HIP95319_chain000.fits\n",
      "\u001b[32mPASSED\u001b[0m\n",
      "orvara/tests/test_main.py::test_run_with_hip2javatoolIAD Unable to load relative RV data from file \n",
      "Loadied 49 RV data points from file orvara/tests/data/HD4747_RV.dat\n",
      "Loaded data from 1 RV instruments.\n",
      "Loading astrometric data for 1 planets\n",
      "Loaded 8 relative astrometric data points from file orvara/tests/data/HD4747_relAST.txt\n",
      "Loading absolute astrometry data for Hip 3850\n",
      "Recognized a 5-parameter fit in Gaia for Hip 3850\n",
      "Not using companion proper motion from Gaia.\n",
      "Using ptemcee.\n",
      "Running MCMC ... \n",
      "[###############################]100%50%\n",
      "Total Time: 2 seconds\n",
      "Mean acceptance fraction (cold chain): 0.432900\n",
      "Writing output to /tmp/tmp42g2h08w/HIP3850_chain000.fits\n",
      "\u001b[32mPASSED\u001b[0m\n",
      "orvara/tests/test_main.py::test_run_with_secondary_priors Unable to load relative RV data from file \n",
      "Loadied 630 RV data points from file orvara/tests/data/Gl758_RV.dat\n",
      "Loaded data from 3 RV instruments.\n",
      "Loading astrometric data for 1 planets\n",
      "Loaded 4 relative astrometric data points from file orvara/tests/data/Gl758_relAST.txt\n",
      "Loading absolute astrometry data for Hip 95319\n",
      "Recognized a 5-parameter fit in Gaia for Hip 95319\n",
      "Not using companion proper motion from Gaia.\n",
      "Using ptemcee.\n",
      "Running MCMC ... \n",
      "[###############################]100%50%\n",
      "Total Time: 2 seconds\n",
      "Mean acceptance fraction (cold chain): 0.217600\n",
      "Writing output to /tmp/tmphuicseu3/HIP95319_chain000.fits\n",
      "\u001b[32mPASSED\u001b[0m\n",
      "\n",
      "\u001b[33m=============================== warnings summary ===============================\u001b[0m\n",
      "orvara/tests/test_main.py::test_run\n",
      "orvara/tests/test_main.py::test_run_with_secondary_priors\n",
      "  /home/gmbrandt/envs/orvara_new/lib/python3.8/site-packages/htof/parse.py:511: UserWarning: For source 095319. The DVD IAD does not indicate which observation epochs were rejected for the final solution. htof will attempt to find which epochs to reject in order to reproduce the catalog parameters. However, if this source also has some corrupted residuals (see Brandt et al. 2021, Section 4), then this will fail. We recommend you switch to using the IAD from the Java tool, since that version of the IAD indicates rejected epochs with negative uncertainties.\n",
      "    warnings.warn(f\"For source {self.meta['star_id']}. The DVD IAD does not indicate which observation epochs were \"\n",
      "\n",
      "-- Docs: https://docs.pytest.org/en/stable/warnings.html\n",
      "\u001b[33m================= \u001b[32m6 passed\u001b[0m, \u001b[33m\u001b[1m4 deselected\u001b[0m, \u001b[33m\u001b[1m2 warnings\u001b[0m\u001b[33m in 7.48s\u001b[0m\u001b[33m ==================\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!pytest -sv -m \"not e2e\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the usage of fit_orbit and plot_orbit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fit_orbit and plot_orbit are the key commands to operate this package. Run --h to see their functionality. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "both commands require one optional argument --output-dir immediately followed by the desired directory to write output to. fit_orbit outputs the MCMC chain from the fit and plot_orbit outputs a set of publication-quality plots. If --output-dir flag is omitted, the output directory will be automatically set to the current working directory by default. \n",
    "\n",
    "Following the output directory is the configuration file tailored to the system you are fitting. Example configuration files can be found in /orvara/tests/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: fit_orbit [-h] [--output-dir OUTPUT_DIR] configfile\r\n",
      "\r\n",
      "Fit a Keplerian orbit or orbits using orvara.\r\n",
      "\r\n",
      "positional arguments:\r\n",
      "  configfile            Path to the planet-specific configuration file.\r\n",
      "\r\n",
      "optional arguments:\r\n",
      "  -h, --help            show this help message and exit\r\n",
      "  --output-dir OUTPUT_DIR\r\n",
      "                        Directory within which to save the MCMC results,\r\n",
      "                        default current directory.\r\n"
     ]
    }
   ],
   "source": [
    "!fit_orbit --h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: plot_orbit [-h] [--output-dir OUTPUT_DIR] configfile\r\n",
      "\r\n",
      "Plot astrometry or RV orbits and other relavant plots.\r\n",
      "\r\n",
      "positional arguments:\r\n",
      "  configfile            Path to the planet-specific configuration file.\r\n",
      "\r\n",
      "optional arguments:\r\n",
      "  -h, --help            show this help message and exit\r\n",
      "  --output-dir OUTPUT_DIR\r\n",
      "                        Directory within which to save the plots, default\r\n",
      "                        current directory.\r\n"
     ]
    }
   ],
   "source": [
    "!plot_orbit --h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running an example of fit_orbit, we will fit orbit to the White Dwarf - Main Sequence (WD-MS) binary system HD159062. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The customizable configuration file for HD 159062 is located in orvara/tests/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's load the config.ini file for HD 159062 here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34mchains\u001b[0m/\r\n",
      "config_Gl758.ini\r\n",
      "config_HD159062.ini\r\n",
      "config_HD4747.ini\r\n",
      "config_hip21.ini\r\n",
      "config.ini\r\n",
      "config_test_relativeRV.ini\r\n",
      "config_with_secondary_priors.ini\r\n",
      "\u001b[01;34mdata\u001b[0m/\r\n",
      "diagnostic_config_fake_fulldr3_7pfit.ini\r\n",
      "diagnostic_config_fake_fulldr3_9pfit.ini\r\n",
      "diagnostic_config.ini\r\n",
      "\u001b[01;34mplots\u001b[0m/\r\n",
      "\u001b[01;34m__pycache__\u001b[0m/\r\n",
      "start_Gl758.dat\r\n",
      "start_hip3850_after7pfit.dat\r\n",
      "start_hip3850.dat\r\n",
      "test_config.py\r\n",
      "test_main.py\r\n",
      "utils.py\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_file(fpath):\n",
    "    with open(fpath, 'r') as f:\n",
    "        lines = f.readlines()\n",
    "        for line in lines:\n",
    "            line = line.rstrip('\\n')\n",
    "            print(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is what an orvara configuration file looks like:\n",
      "[data_paths]\n",
      "# Hipparcos ID of the star, used for fetching absolute astrometry.\n",
      "# If the Hipparcos ID is invalid the fit will not use HGCA astrometry. \n",
      "# In that case you must supply a parallax and parallax_error in [priors].\n",
      "HipID = 85653\n",
      "# The file containing the Hipparcos Gaia Catalog of Accelerations.\n",
      "HGCAFile = HGCA_vEDR3.fits\n",
      "# The file containing the radial velocity time series for the star.\n",
      "RVFile = orvara/tests/data/HD159062_RV.dat\n",
      "# The file containing the relative astrometry for the star.\n",
      "AstrometryFile = orvara/tests/data/HD159062_relAST.txt\n",
      "# The path to the Gaia/GOST epochs and scan angles as a csv file.\n",
      "GaiaDataDir = orvara/tests/data/gaia\n",
      "# The path to all the Hipparcos (original reduction) intermediate data. Subfolders (as used on the original CD) are automatically handled.\n",
      "Hip1DataDir = orvara/tests/data/hip1\n",
      "# The path to all the Hipparcos (Floor van Leeuwen/second reduction) intermediate data. Subfolders (as used on the original DVD) are automatically handled.\n",
      "Hip2DataDir = orvara/tests/data/hip2\n",
      "# The file path to the initial parameter guesses (starting conditions for the walkders) to the orbit. Set to None for default guess.\n",
      "start_file = None\n",
      "\n",
      "[mcmc_settings]\n",
      "# Number of temperatures to use in the parallel tempering chain.\n",
      "ntemps = 10\n",
      "# Number of walkers. Each walker will have ntemps number of chains.\n",
      "nwalkers = 100\n",
      "# Number of planets to fit.\n",
      "nplanets = 1\n",
      "# Number of steps contained in each chain.\n",
      "nstep = 20000\n",
      "# How much to thin the chain (save every thin-th step).\n",
      "thin = 50\n",
      "# Number of threads to use with emcee. Note this built-in parallelization is poor.\n",
      "nthreads = 1\n",
      "# True if you want to use the epoch astrometry in GaiaDataDir, Hip1DataDir,\n",
      "# Hip2DataDir. Hip1 and Hip2 data is combined in a 60/40 mix as described in\n",
      "# the HGCA paper (Brandt 2018 and 2021). False if not. \n",
      "use_epoch_astrometry = True\n",
      "# Use a separate jitter for each RV instrument?\n",
      "jit_per_inst = False\n",
      "\n",
      "[priors_settings]\n",
      "# priors on primary mass (solar), if prior is not specified, mpri should be inf\n",
      "mpri = 0.8\n",
      "mpri_sig = 0.05\n",
      "minjitter = 1e-5\n",
      "maxjitter = 1e3\n",
      "\n",
      "[secondary_gaia]\n",
      "# If the secondary star is in Gaia, set companion_ID to a nonnegative number\n",
      "# matching the ID of the companion in the relative astrometry file.  Setting\n",
      "# companion_ID to a negative number ignores the rest of this.\n",
      "companion_ID = -1\n",
      "# The rest of these should be from Gaia in units of mas.\n",
      "pmra = 0\n",
      "pmdec = 0\n",
      "epmra = 100\n",
      "epmdec = 100\n",
      "corr_pmra_pmdec = 0\n",
      "\n",
      "[plotting]\n",
      "# Path to mcmc chains. This is what's produced by the orbit_fit command.\n",
      "McmcDataFile = orvara/tests/chains/HD159062B_chain000.fits\n",
      "# Diagnostic plots to check for convergence.\n",
      "check_convergence = True\n",
      "# Define burnin for chains.\n",
      "burnin = 50\n",
      "# Which companion to plot?\n",
      "iplanet = 0\n",
      "\n",
      "# Name of the target\n",
      "target = HD159062B\n",
      "# This is a customized range of epochs you want to plot.\n",
      "start_epoch = 1990.\n",
      "end_epoch = 2015.\n",
      "# Number of random orbits drawn from the posterior distribution you want to plot.\n",
      "num_orbits = 50\n",
      "# Define step size for plotting.\n",
      "num_steps = 1500\n",
      "# Plot random predicted epoch positions on the Astrometry plot.\n",
      "predicted_years = 1990,2000,2010,2020,2030\n",
      "# Prediction of the location of the companion in a future epoch.\n",
      "position_predict = 2010\n",
      "\n",
      "# Select which plot you want to generate\n",
      "Astrometry_orbits_plot = True\n",
      "Astrometric_prediction_plot = True\n",
      "RV_orbits_plot = True\n",
      "RV_plot = True\n",
      "RV_Instrument = All\n",
      "Relative_separation_plot = True\n",
      "Position_angle_plot = True\n",
      "Proper_motion_plot = True\n",
      "Proper_motion_separate_plots = False\n",
      "Corner_plot = True\n",
      "\n",
      "############# Advanced plotting settings #############\n",
      "# 1. Set axes limits. Set the upper and lower limit for x and y axes.\n",
      "set_limit = False\n",
      "xlim = 1980, 2025\n",
      "ylim = -2.8,2.8\n",
      "# Choose the color of the marker for plotting the observed data points.\n",
      "marker_color = blue\n",
      "\n",
      "# 2. Turn on/off colorbar. Choose a matplotlib colormap, and choose to color\n",
      "# code by the secondary mass (msec_solar or msec_jup) or eccentricity (ecc).\n",
      "use_colorbar = True\n",
      "colormap = viridis\n",
      "reference = msec_solar\n",
      "\n",
      "# 3. Turn on/off the title of the plot? Additionally, if user wants to add\n",
      "# text somewhere on the plot, enter the text and its x and y position (with\n",
      "# (0, 0) being the bottom left and (1, 1) being the top right).\n",
      "show_title = False\n",
      "add_text = True\n",
      "text_name = HD 159062B\n",
      "x_text = 0.5\n",
      "y_text = 0.9\n",
      "\n",
      "\n",
      "[save_results]\n",
      "# quantiles: median and 1sigma = 0.16, 0.5, 0.84, \n",
      "# median and 2sigma = 0.025, 0.5, 0.975, or arbitrary\n",
      "save_params = True\n",
      "err_margin = 0.16, 0.5, 0.84\n"
     ]
    }
   ],
   "source": [
    "print('This is what an orvara configuration file looks like:')\n",
    "display_file('orvara/tests/config_HD159062.ini')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take note of the setting above in mcmc_settings:\n",
    "    use_epoch_astrometry = True\n",
    "\n",
    "This uses a precise treatment for the Hipparcos and Gaia absolute astrometry. If it is set to True, then the directory paths in the config file: GaiaDataDir, Hip1DataDir, and Hip2DataDir, have to exist. Orvara will not make them for you. Orvara will, however, download the appropriate Intermediate Data Files for your source into those folders. Note that for this tutorial the directories have already been created and the files have already been downloaded.\n",
    "\n",
    "For a new fit, it is possible that downloading of the intermediate data files does not work. In this case ,you can just set use_epoch_astrometry = False and proceed with your fits. In most cases, setting it to False is entirely fine and will result in nearly identical fits to the True case. You can troubleshoot (or download the files manually and place them in the appropriate directories) and set this parameter to True later. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's load the data sets needed for orbit fitting. There are three data types in general: RV, relative astrometry, or/and absolute astrometry. HD 159062 has all. Example data files are located in orvara/tests/data. Let's display them here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34mgaia\u001b[0m/                HD159062_RV.dat             HGCA_vDR3_FAKE_7p.fits  \u001b[01;34mhip21\u001b[0m/\r\n",
      "Gl758_relAST.txt     HD4747_FAKE_relativeRV.dat  HGCA_vDR3_FAKE_9p.fits\r\n",
      "Gl758_RV.dat         HD4747_relAST.txt           \u001b[01;34mhip1\u001b[0m/\r\n",
      "HD159062_relAST.txt  HD4747_RV.dat               \u001b[01;34mhip2\u001b[0m/\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the RV data, it has four columns: Epoch[JD], RV[m/s], RV_err[m/s] and InstrumentID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is what an orvara RV data input file looks like:\n",
      "# Instrument: HIRES on Keck\n",
      "# Published in Hirsch, L. A., Ciardi, D. R., Howard, A. W., et al. 2019, ApJ, 878, 50\n",
      "# File format:\n",
      "# Epoch[JD]   RV[m/s] RV_err[m/s] InstrumentID\n",
      "2452832.91668\t89.22\t1.62\t0\n",
      "2453074.09155\t83.23\t1.68\t0\n",
      "2453196.81631\t76.88\t1.58\t0\n",
      "2453430.15779\t65.68\t0.93\t0\n",
      "2453547.94513\t60.98\t0.84\t0\n",
      "2453604.84622\t57.70\t1.03\t0\n",
      "2453807.13748\t50.30\t1.00\t0\n",
      "2453926.45063\t43.33\t0.57\t0\n",
      "2453927.88530\t44.96\t0.83\t0\n",
      "2454248.03230\t31.71\t1.44\t0\n",
      "2454249.95458\t32.00\t1.49\t0\n",
      "2454252.04380\t32.81\t1.02\t0\n",
      "2454255.93510\t32.28\t1.01\t0\n",
      "2454277.86337\t31.99\t1.33\t0\n",
      "2454278.90748\t32.22\t1.22\t0\n",
      "2454279.94547\t31.42\t1.29\t0\n",
      "2454285.90800\t30.58\t1.34\t0\n",
      "2454294.90713\t30.60\t1.30\t0\n",
      "2454634.43562\t19.74\t1.03\t0\n",
      "2454635.94255\t19.00\t1.43\t0\n",
      "2454636.93316\t17.04\t1.53\t0\n",
      "2454638.33247\t18.76\t0.99\t0\n",
      "2454640.38633\t14.57\t1.06\t0\n",
      "2454641.94950\t19.55\t1.47\t0\n",
      "2454644.08325\t15.11\t1.52\t0\n",
      "2454688.88580\t17.38\t1.47\t0\n",
      "2454689.89920\t9.37\t1.58\t0\n",
      "2454957.13845\t3.77\t1.67\t0\n",
      "2455401.78282\t-16.12\t1.40\t0\n",
      "2455722.86450\t-22.42\t1.56\t0\n",
      "2455782.89029\t-26.95\t1.55\t0\n",
      "2456164.72278\t-40.66\t1.54\t0\n",
      "2456451.97582\t-50.77\t1.03\t0\n",
      "2456475.78721\t-49.97\t1.43\t0\n",
      "2456709.08945\t-54.30\t1.43\t0\n",
      "2456883.74892\t-65.35\t1.49\t0\n",
      "2456910.83550\t-66.37\t1.62\t0\n",
      "2457061.17682\t-66.74\t1.43\t0\n",
      "2457211.95513\t-76.17\t1.40\t0\n",
      "2457671.70374\t-91.12\t1.72\t0\n",
      "2457831.16397\t-94.97\t1.59\t0\n",
      "2457853.14530\t-91.51\t1.92\t0\n",
      "2457970.94501\t-103.25\t1.86\t0\n",
      "2458263.03029\t-108.99\t1.22\t0\n",
      "2458349.72648\t-110.69\t1.58\t0\n"
     ]
    }
   ],
   "source": [
    "print('This is what an orvara RV data input file looks like:')\n",
    "display_file('orvara/tests/data/HD159062_RV.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the relative astrometry file for HD 159062. It has seven columns: Date      , Sep (arcsec) , Err_sep  ,    PA (degrees),  Err_PA  Corr_Sep_PA, and  PlanetID. Planet ID depends on how many planets are in the system, 0 denotes the innermost planet, the second innermost is denoted 1 and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is what an orvara relative astrometry data input file looks like:\n",
      "# Instrument: NIRC2 on Keck\n",
      "# Published in Hirsch, L. A., Ciardi, D. R., Howard, A. W., et al. 2019, ApJ, 878, 50\n",
      "##\n",
      "# File format\n",
      "#   Date       Sep  Err_sep      PA  Err_PA  Corr_Sep_PA  PlanetID  \n",
      " 2012.48     2.594    0.014   298.6     0.1            0         0\n",
      " 2014.78     2.637    0.014   299.5     0.1            0         0\n",
      " 2018.32     2.673    0.003  301.27    0.04            0         0\n"
     ]
    }
   ],
   "source": [
    "print('This is what an orvara relative astrometry data input file looks like:')\n",
    "display_file('orvara/tests/data/HD159062_relAST.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For absolute astrometry from the HGCA, make sure that the HGCA file (DR2 or eDR3) is specified. We are using Gaia eDR3 version of the HGCA.\n",
    "\n",
    "To use intermediate astrometry, make sure that the parameter use_epoch_astrometry in the config file is set to True.\n",
    "\n",
    "We first recommend running fits with use_epoch_astrometry set to False, then turning it on later.\n",
    "\n",
    "Make sure that the directories for the Gaia and Hipparcos intermediate data, set in the .ini config file, are valid directories. You will need one directory each for Hip1, Hip2, and Gaia.\n",
    "\n",
    "orvara will download any file you need automatically (and save them to the directories you specified), but:\n",
    "\n",
    "You may download the files by hand if you want, and simply put them in the data directory, with filenames containing the hip ID. The Gaia eDR3 observations are publicly accessible via the Gaia Observation Forecast Tool (https : //gaia.esac.esa.int/gost/index.jsp). Gaia files are .csv files queried from GOST. The Hip1 and Hip2 data are (if not, soon) avaliable publically from ESA.\n",
    "\n",
    "Here are those files, with filenames containing HIP85653 are the ones we want for HD 159062. Remember, orvara should download these files automatically in normal use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HIP003850.csv  HIP085653.csv  HIP095319.csv\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/data/gaia/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "003850.txt  085653.txt  95319.txt\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/data/hip1/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HIP003850.d  HIP085653.d  HIP095319.d\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/data/hip2/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, let's run fit_orbit and plot_orbit commands in the root directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to run an example fit to HD 159062. Verifty that the RV data, relative astrometry data and absolute astrometry data are loaded in correctly to ensure a suceesful fit. Progress bar shows the progress of MCMC. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unable to load relative RV data from file \n",
      "Loadied 45 RV data points from file orvara/tests/data/HD159062_RV.dat\n",
      "Loaded data from 1 RV instruments.\n",
      "Loading astrometric data for 1 planets\n",
      "Loaded 3 relative astrometric data points from file orvara/tests/data/HD159062_relAST.txt\n",
      "Loading absolute astrometry data for Hip 85653\n",
      "Recognized a 5-parameter fit in Gaia for Hip 85653\n",
      "Not using companion proper motion from Gaia.\n",
      "/home/gmbrandt/envs/orvara_new/lib/python3.8/site-packages/htof/parse.py:511: UserWarning: For source 085653. The DVD IAD does not indicate which observation epochs were rejected for the final solution. htof will attempt to find which epochs to reject in order to reproduce the catalog parameters. However, if this source also has some corrupted residuals (see Brandt et al. 2021, Section 4), then this will fail. We recommend you switch to using the IAD from the Java tool, since that version of the IAD indicates rejected epochs with negative uncertainties.\n",
      "  warnings.warn(f\"For source {self.meta['star_id']}. The DVD IAD does not indicate which observation epochs were \"\n",
      "Using ptemcee.\n",
      "Running MCMC ... \n",
      "[###############################]100% 1%  2%  3%  4%  5%  6%  7%  8%  9% 10% 11% 12% 13% 14% 15% 16% 17% 18% 19% 20% 21% 22% 23% 24% 25% 26% 27% 28% 29% 30% 31% 32% 33% 34% 35% 36% 37% 38% 39% 40% 41% 42% 43% 44% 45% 46% 47% 48% 49% 50% 51% 52% 53% 54% 55% 56% 57% 58% 59% 60% 61% 62% 63% 64% 65% 66% 67% 68% 69% 70% 71% 72% 73% 74% 75% 76% 77% 78% 79% 80% 81% 82% 83% 84% 85% 86% 87% 88% 89% 90% 91% 92% 93% 94% 95% 96% 97% 98% 99%\n",
      "Total Time: 541 seconds\n",
      "Mean acceptance fraction (cold chain): 0.119032\n",
      "Writing output to orvara/tests/chains/HD159062B_chain000.fits\n",
      "[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7fc557731340>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x7fc5573396d0>]\n"
     ]
    }
   ],
   "source": [
    "!fit_orbit --output-dir orvara/tests/chains orvara/tests/config_HD159062.ini"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have the MCMC chain output by fit_orbit, feed it to McmcDataFile= in the config file. You can set a burn-in phase for the chain. After that, you are ready to go on to plot the orbits. plot_orbit will generate a suite of 7 plots relavant to RV and Astrometry, and a plot of the MCMC chain to breifly check convergence, and two files containting the beststep and posterior parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating plots for target HD159062B\n",
      "Generating diagnostic plots to check convergence.\n",
      "Plotting Astrometry orbits, your plot is generated at orvara/tests/plots/\n",
      "WARNING: You will want to use more than {self.num_orbits} orbits to predict the position. Ideally you should use 500 or more. Set num_orbits=500 or more in the config file.\n",
      "simulating 50 orbits for astrometric prediction. Will take roughly 4 seconds per 100 orbits.\n",
      "The mean RA and Dec position of the companion at 2010.0 is:\n",
      "RA: -2263.138 +- 11.5629 mas, Dec: 1169.0681 +- 8.4776 mas\n",
      "Plotting RV orbits, your plot is generated at orvara/tests/plots/\n",
      "Plotting RV, your plot is generated at orvara/tests/plots/\n",
      "Plotting Separation, your plot is generated at orvara/tests/plots/\n",
      "Plotting Position Angle, your plot is generated at orvara/tests/plots/\n",
      "Plotting Proper Motions, your plot is generated at orvara/tests/plots/\n",
      "Plotting Corner plot, your plot is generated at orvara/tests/plots/\n",
      "Saving beststep parameters to orvara/tests/plots/\n",
      "Saving posterior parameters to orvara/tests/plots/\n"
     ]
    }
   ],
   "source": [
    "!plot_orbit --output-dir orvara/tests/plots/ orvara/tests/config_HD159062.ini"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Great! Your plots are generated in the default or specified output directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see the whole sets of plots in the output directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "astrometric_orbit_HD159062B.pdf       posterior_params_HD159062B.txt\r\n",
      "astrometric_prediction_HD159062B.pdf  Proper_Motions_HD159062B.pdf\r\n",
      "beststep_params_HD159062B.txt         relsep_OC_HD159062B.pdf\r\n",
      "Corner_HD159062B.pdf                  RV_OC_HD159062B_InstAll.pdf\r\n",
      "\u001b[0m\u001b[01;35mDiagnostic_HD159062B.png\u001b[0m              RV_orbit_HD159062B.pdf\r\n",
      "PA_OC_HD159062B.pdf\r\n"
     ]
    }
   ],
   "source": [
    "ls orvara/tests/plots/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For instance, here is the visual orbit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gmbrandt/envs/orvara_new/lib/python3.8/site-packages/IPython/core/display.py:717: UserWarning: Consider using IPython.display.IFrame instead\n",
      "  warnings.warn(\"Consider using IPython.display.IFrame instead\")\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<iframe src=./orvara/tests/plots/astrometric_orbit_HD159062B.pdf width=700 height=350></iframe>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<iframe src=./orvara/tests/plots/astrometric_orbit_HD159062B.pdf width=700 height=350></iframe>')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We can also inspect the posteriors and chisquareds from the fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are the chisquared values of the model with respect to the data (for all the input data types):\n",
      "plx_ML  0.04618635408\n",
      "pmra_ML  0.1698076766\n",
      "pmdec_ML  0.08777529253\n",
      "chisq_sep  0.4189472422\n",
      "chisq_PA  2.122972047\n",
      "chisq_H  4.00100902\n",
      "chisq_HG  0.1445642573\n",
      "chisq_G  0.2667943578\n",
      "chisq_relRV  0\n",
      "RV_ZP_0_ML  -1456.099472\n"
     ]
    }
   ],
   "source": [
    "print('Here are the chisquared values of the model with respect to the data (for all the input data types):')\n",
    "display_file('orvara/tests/plots/beststep_params_HD159062B.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `chisq_relRV` is equal to 0. This is because this fit *did not* use any relative RV data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here are the posteriors with +- 1 sigma confidence intervals:\n",
      "Jit (m/s)         ${1.27}_{-0.30}^{+0.33}$\n",
      "Mpri (solar)         ${0.797}_{-0.049}^{+0.049}$\n",
      "Msec (solar)         ${0.6028}_{-0.0071}^{+0.0080}$\n",
      "a (AU)         ${55.6}_{-7.4}^{+5.3}$\n",
      "sqesinw         ${-0.19}_{-0.10}^{+0.21}$\n",
      "sqecosw         ${0.32}_{-0.18}^{+0.20}$\n",
      "inclination (deg)         ${61.0}_{-3.8}^{+2.0}$\n",
      "ascending node (deg)         ${134.9}_{-1.4}^{+3.1}$\n",
      "mean longitude (deg)         ${143.9}_{-10}^{+6.1}$\n",
      "parallax (mas)         ${46.1855}_{-0.0042}^{+0.0042}$\n",
      "period (yrs)         ${350}_{-68}^{+53}$\n",
      "argument of periastron (deg)         ${328}_{-253}^{+14}$\n",
      "eccentricity         ${0.142}_{-0.092}^{+0.18}$\n",
      "semimajor axis (mas)         ${2568}_{-343}^{+244}$\n",
      "T0 (JD)         ${2518039}_{-9521}^{+25485}$\n",
      "mass ratio         ${0.757}_{-0.044}^{+0.050}$\n"
     ]
    }
   ],
   "source": [
    "print('Here are the posteriors with +- 1 sigma confidence intervals:')\n",
    "display_file('orvara/tests/plots/posterior_params_HD159062B.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}