OFTI Introduction

by Isabel Angelo and Sarah Blunt (2018)

OFTI (Orbits For The Impatient) is an orbit-generating algorithm designed specifically to handle data covering short fractions of long-period exoplanets (Blunt et al. 2017). Here we go through steps of using OFTI within orbitize!

import orbitize

Basic Orbit Generating

Orbits are generated in OFTI through a Driver class within orbitize. Once we have imported this class:

import orbitize.driver

we can initialize a Driver object specific to our data:

myDriver = orbitize.driver.Driver('{}/GJ504.csv'.format(orbitize.DATADIR), # path to data file
                                  'OFTI', # name of algorithm for orbit-fitting
                                  1, # number of secondary bodies in system
                                  1.22, # total mass [M_sun]
                                  56.95, # total parallax of system [mas]
                                  mass_err=0.08, # mass error [M_sun]
                                  plx_err=0.26) # parallax error [mas]

Because OFTI is an object class within orbitize, we can assign all of the OFTI attributes onto a variable (s). We can then generate orbits for s using a function called run_sampler, a method of the OFTI class. The run_sampler method takes in the desired number of accepted orbits as an input.

Here we use run OFTI to randomly generate orbits until 1000 are accepted:

s = myDriver.sampler
orbits = s.run_sampler(1000)

We have now generated 1000 possible orbits for our system. Here, orbits is a (1000 x 8) array, where each of the 1000 elements corresponds to a single orbit. An orbit is represented by 8 orbital elements.

Here is an example of what an accepted orbit looks like from orbitize:

orbits[0]

array([4.93916907e+01, 8.90197501e-03, 2.63925411e+00, 2.44962990e+00,
       9.31508665e-01, 1.20302112e-01, 5.74242058e+01, 1.22728974e+00])

To further inspect what each of the 8 elements in your orbit represents, you can view the system.param_idx variable. This is a dictionary that tells you the indices of your orbit that correspond to semi-major axis (a), eccentricity (e), inclination (i), argument of periastron (aop), position angle of nodes (pan), and epoch of periastron passage (epp). The last two indices are the parallax and system mass, and the number following the parameter name indicates the number of the body in the system.

s.system.param_idx

{'sma1': 0,
 'ecc1': 1,
 'inc1': 2,
 'aop1': 3,
 'pan1': 4,
 'tau1': 5,
 'plx': 6,
 'mtot': 7}

Plotting

Now that we can generate possible orbits for our system, we want to plot the data to interpret our results. Here we will go through a brief overview on ways to visualize your data within orbitize. For a more detailed guide on data visualization capabilities within orbitize, see the Orbitize plotting tutorial.

Histogram

One way to visualize our results is through histograms of our computed orbital parameters. Our orbits are outputted from run_sampler as an array of orbits, where each orbit is represented by a set of orbital elements:

print(orbits.shape)
orbits[:5]

(1000, 8)

array([[4.93916907e+01, 8.90197501e-03, 2.63925411e+00, 2.44962990e+00,
        9.31508665e-01, 1.20302112e-01, 5.74242058e+01, 1.22728974e+00],
       [4.69543031e+01, 1.31571508e-01, 2.52917998e+00, 1.34963602e+00,
        4.18692436e+00, 4.17659289e-01, 5.73207900e+01, 1.23162413e+00],
       [5.15848551e+01, 1.18074455e-01, 2.26110475e+00, 2.98346893e+00,
        2.31713931e+00, 3.94202277e-03, 5.69191065e+01, 1.15389146e+00],
       [3.89558225e+01, 3.71357464e-01, 2.94801131e+00, 3.08542398e+00,
        4.77645562e+00, 7.15043369e-01, 5.72827098e+01, 1.05721164e+00],
       [8.19463988e+01, 1.45955646e-02, 2.11512811e+00, 4.52064036e+00,
        4.44802306e+00, 9.32004660e-01, 5.72429430e+01, 1.35323242e+00]])

We can effectively view outputs from run_sampler by creating a histogram of a given orbit element to see its distribution of possible values. Our system.param_idx dictionary is useful here. We can use it to determine the index of a given orbit that corresponds to the orbital element we are interested in:

s.system.param_idx

{'sma1': 0,
 'ecc1': 1,
 'inc1': 2,
 'aop1': 3,
 'pan1': 4,
 'tau1': 5,
 'plx': 6,
 'mtot': 7}

If we want to plot the distribution of orbital semi-major axes (a) in our generated orbits, we would use the index dictionary s.system.param_idx to index the semi-major axis element from each orbit:

sma = [x[s.system.param_idx['sma1']] for x in orbits]

%matplotlib inline
import matplotlib.pyplot as plt
plt.hist(sma, bins=30)
plt.xlabel('orbital semi-major axis [AU]')
plt.ylabel('occurrence')
plt.show()

You can use this method to create histograms of any orbital element you are interested in:

ecc = [x[s.system.param_idx['ecc1']] for x in orbits]
i = [x[s.system.param_idx['inc1']] for x in orbits]

plt.figure(figsize=(10,3))
plt.subplot(131)
plt.hist(sma, bins=30)
plt.xlabel('orbital semi-major axis [AU]')
plt.ylabel('occurrence')

plt.subplot(132)
plt.hist(ecc, bins=30)
plt.xlabel('eccentricity [0,1]')
plt.ylabel('occurrence')

plt.subplot(133)
plt.hist(i, bins=30)
plt.xlabel('inclination angle [rad]')
plt.ylabel('occurrence')

plt.show()

In addition to our orbits array, Orbitize also creates a Results class that contains built-in plotting capabilities for two types of plots: corner plots and orbit plots.

Corner Plot

After generating the samples, the run_sampler method also creates a Results object that can be accessed with s.results:

myResults = s.results

We can now create a corner plot using the function plot_corner within the Results class. This function requires an input list of the parameters, in string format, that you wish to include in your corner plot. We can even plot all of the orbital parameters at once! As shown below:

corner_figure = myResults.plot_corner(param_list=['sma1', 'ecc1', 'inc1', 'aop1', 'pan1','tau1'])

A Note about Convergence

Those of you with experience looking at corner plots will note that the result here does not look converged (i.e. we need more samples for our results to be statistically significant). Because this is a tutorial, we didn’t want you to have to wait around for a while for the OFTI results to converge.

It’s safe to say that OFTI should accept a minimum of 10,000 orbit for convergence. For pretty plots to go in publications, we recommend at least 1,000,000 accepted orbits.

Orbit Plot

What about if we want to see how the orbits look in the sky? Don’t worry, the Results class has a command for that too! It’s called plot_orbits. We can create a simple orbit plot by running the command as follows:

epochs = myDriver.system.data_table['epoch']

orbit_figure = myResults.plot_orbits(
    start_mjd=epochs[0] # Minimum MJD for colorbar (here we choose first data epoch)
)

WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 8 of "dubious year (Note 6)" [astropy._erfa.core]

<Figure size 1008x432 with 0 Axes>

Advanced OFTI and API Interaction

We’ve seen how the run_sampler command is the fastest way to generate orbits within OFTI. For users interested in what’s going on under-the-hood, this part of the tutorial takes us each step of run_sampler. Understanding the intermediate stages of orbit-fitting can allow for more customization that goes beyond Orbitize’s default parameters.

We begin again by intializing a sampler object on which we can run OFTI:

myDriver = orbitize.driver.Driver('{}/GJ504.csv'.format(orbitize.DATADIR), # path to data file
                                  'OFTI', # name of algorith for orbit-fitting
                                  1, # number of secondary bodies in system
                                  1.22, # total mass [M_sun]
                                  56.95, # total parallax of system [mas]
                                  mass_err=0.08, # mass error [M_sun]
                                  plx_err=0.26) # parallax error [mas]

s = myDriver.sampler

In orbitize, the first thing that OFTI does is prepare an initial set of possible orbits for our object through a function called prepare_samples, which takes in the number of orbits to generate as an input. For example, we can generate 100,000 orbits as follows:

samples = s.prepare_samples(100000)

Here, samples is an array of randomly generated orbits that have been scaled-and-rotated to fit our astrometric observations. The first and second dimension of this array are the number of orbital elements and total orbits generated, respectively. In other words, each element in samples represents the value of a particular orbital element for each generated orbit:

print('samples: ', samples.shape)
print('first element of samples: ', samples[0].shape)

samples:  (8, 100000)
first element of samples:  (100000,)

Once our initial set of orbits is generated, the orbits are vetted for likelihood in a function called reject. This function computes the probability of an orbit based on its associated chi squared. It then rejects orbits with lower likelihoods and accepts the orbits that are more probable. The output of this function is an array of possible orbits for our input system.

orbits, lnlikes = s.reject(samples)

Our orbits array represents the final orbits that are output by OFTI. Each element in this array contains the 8 orbital elements that are computed by orbitize:

orbits.shape

(1, 8)

We can synthesize this sequence with the run_sampler() command, which runs through the steps above until the input number of orbits has been accepted. Additionally, we can specify the number of orbits generated by prepare_samples each time the sequence is initiated with an argument called num_samples. Higher values for num_samples will output more accepted orbits, but may take longer to run since all initially prepared orbits will be run through the rejection step.

orbits = s.run_sampler(100, num_samples=1000)

Saving and Loading Results

Finally, we can save our generated orbits in a file that can be easily read for future use and analysis. Here we will walk through the steps of saving a set of orbits to a file in hdf5 format. The easiest way to do this is using orbitize.Results.save_results():

s.results.save_results('orbits.hdf5')

Now when you are ready to use your orbits data, it is easily accessible through the file we’ve created. One way to do this is to load the data into a new results object; in this way you can make use of the functions that we learn before, like plot_corner and plot_orbits. To do this, use the results module:

import orbitize.results
loaded_results = orbitize.results.Results() # create a blank results object to load the data
loaded_results.load_results('orbits.hdf5')

Alternatively, you can directly access the saved data using the h5py module:

import h5py
f = h5py.File('orbits.hdf5', 'r')
orbits = f['post']

print('orbits array dimensions: ', orbits.shape)
print('orbital elements for first orbit: ', orbits[0])

f.close()

orbits array dimensions:  (100, 8)
orbital elements for first orbit:  [42.56737306  0.18520168  2.50497349  1.46225095  3.29419715  0.60649612
 57.26589357  1.1478072 ]

And now we can easily work with the saved orbits that were generated by orbitize! Find out more about generating orbits in orbitize! with tutorials here.