Merge pull request #302 from sblunt/add_e2e_tests

add 1 planet HR 8799 e fit to end to end tests.

tests/end-to-end-tests/hr8799e_1epochgravity.csv

@@ -0,0 +1,13 @@
+epoch,object,raoff,raoff_err,decoff,decoff_err,radec_corr,sep,sep_err,pa,pa_err
+55042,1,-306,7,-211,7,,,,,
+55136,1,-310,9,-187,9,,,,,
+55390,1,-323,6,-166,6,,,,,
+55499,1,-341,16,-143,16,,,,,
+55763,1,-352,8,-130,8,,,,,
+56130,1,-373,8,-84,8,,,,,
+56226,1,-370,9,-76,9,,,,,
+56581,1,-373,13,-17,13,,,,,
+56855,1,-387,11,3,11,,,,,
+56609,1,,,,,,382.6,2.1,265.13,0.24
+57650,1,,,,,,384.8,1.7,281.68,0.25
+58358.24,1,-357.640,0.07,163.340,0.18,-0.530,,,,

tests/end-to-end-tests/hr8799e_only.py

@@ -0,0 +1,74 @@
+import os
+import matplotlib.pyplot as plt
+
+import orbitize
+from orbitize import system, read_input, priors, sampler
+
+"""
+Attempts to reproduce the 1 planet orbit fit in GRAVITY Collaboration et al. (2019)
+"""
+
+# End to end example with beta Pic astrometry
+# Based on driver.py
+
+import numpy as np
+from orbitize import read_input, system, sampler, priors
+import multiprocessing as mp
+import orbitize
+
+
+# System parameters
+datafile='hr8799e_1epochgravity.csv'
+num_secondary_bodies=1
+system_mass=1.52 # Msol
+plx=24.2175 #mas
+mass_err=0.15 # Msol
+plx_err=0.0881 #mas
+
+# Sampler parameters
+likelihood_func_name='chi2_lnlike'
+n_temps=20
+n_walkers=1000
+n_threads=mp.cpu_count()
+total_orbits=10000000 # n_walkers x num_steps_per_walker
+burn_steps=50000
+
+tau_ref_epoch = 50000
+
+
+# Read in data
+data_table = read_input.read_file(datafile)
+
+# Initialize System object which stores data & sets priors
+my_system = system.System(
+    num_secondary_bodies, data_table, system_mass,
+    plx, mass_err=mass_err, plx_err=plx_err, tau_ref_epoch=tau_ref_epoch
+)
+
+my_sampler = sampler.MCMC(my_system, n_temps, n_walkers, n_threads)
+
+# Run the sampler to compute some orbits, yeah!
+# Results stored in bP_sampler.chain and bP_sampler.lnlikes
+my_sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=10)
+
+my_sampler.results.save_results("hr8799e_gravity_chains.hdf5")
+
+#import orbitize.results as results
+#my_sampler.results = results.Results()
+#my_sampler.results.load_results("hr8799e_gravity_chains.hdf5")
+
+# make corner plot
+fig = my_sampler.results.plot_corner()
+plt.savefig('corner_hr8799e_gravity.png', dpi=250)
+
+# print SMA, ecc, inc
+labels = ["sma", "ecc", "inc"]
+paper_vals = ["16.4 (+2.1/-1.1)", "0.15 +/- 0.08", "25 +/- 8"]
+for i in range(len(labels)):
+    med_val = np.median(my_sampler.results.post[:,i])
+    ci_vals = np.percentile(my_sampler.results.post[:,i], [84, 16]) - med_val
+    if labels[i] == 'inc':
+        med_val = np.degrees(med_val)
+        ci_vals = np.degrees(ci_vals)
+    print("{0}: paper value is {1}".format(labels[i], paper_vals[i]))
+    print("{3}: this fit obtained {0:.2f} (+{1:.2f}/-{2:.2f})".format(med_val, ci_vals[0], ci_vals[1], labels[i]))