Advanced Plotting

by Henry Ngo (2018)

The results.py module contains several plotting functions to visualize the results of your orbitize orbit fit. Basic uses of these functions are covered in the OFTI and MCMC tutorials. Here, we will examine these plotting functions more deeply. This tutorial will be updated as more features are added to orbitize.

1. Test orbit generation with OFTI

In order to have sample data for this tutorial, we will use OFTI to generate some orbits for a published dataset on the GJ 504 system. The following code block is from the OFTI Tutorial , with 10000 orbits generated. Please see that tutorial for details.

Note: If you have already run this tutorial and saved the computed orbits, you may skip to Section 3 and load up your previously computed orbits instead of running this block below.

import orbitize.driver

myDriver = orbitize.driver.Driver(
    "{}/GJ504.csv".format(orbitize.DATADIR),  # relative or absolute 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,  # parallax of system [mas]
    mass_err=0.08,  # mass error [M_sun]
    plx_err=0.26,  # parallax error [mas]
)
s = myDriver.sampler
orbits = s.run_sampler(1000)

1000/1000 orbits found

2. Accessing a Results object with computed orbits

After computing your orbits from either OFTI or MCMC, they are accessible as a Results object for further analysis and plotting. This object is an attribute of s, the sampler object defined above.

myResults = (
    s.results
)  # array of MxN array of orbital parameters (M orbits with N parameters per orbit)

It is also useful to save this Results object to a file if we want to load up the same data later without re-computing the orbits.

myResults.save_results("plotting_tutorial_GJ504_results.hdf5")

For more information on the Results object, see below.

myResults?

Type:        Results
String form: <orbitize.results.Results object at 0x7fcd20f3c2e0>
File:        ~/Documents/GitHub/orbitize/orbitize/results.py
Docstring:
A class to store accepted orbital configurations from the sampler

Args:
    system (orbitize.system.System): System object used to do the fit.
    sampler_name (string): name of sampler class that generated these results
        (default: None).
    post (np.array of float): MxN array of orbital parameters
        (posterior output from orbit-fitting process), where M is the
        number of orbits generated, and N is the number of varying orbital
        parameters in the fit (default: None).
    lnlike (np.array of float): M array of log-likelihoods corresponding to
        the orbits described in ``post`` (default: None).
    version_number (str): version of orbitize that produced these results.
    data (astropy.table.Table): output from ``orbitize.read_input.read_file()``
    curr_pos (np.array of float): for MCMC only. A multi-D array of the
        current walker positions that is used for restarting a MCMC sampler.

Written: Henry Ngo, Sarah Blunt, 2018

API Update: Sarah Blunt, 2021

Note that you can also add more computed orbits to a results object with myResults.add_samples():

myResults.add_samples?

Signature: myResults.add_samples(orbital_params, lnlikes, curr_pos=None)
Docstring:
Add accepted orbits, their likelihoods, and the orbitize version number
to the results

Args:
    orbital_params (np.array): add sets of orbital params (could be multiple)
        to results
    lnlike (np.array): add corresponding lnlike values to results
    curr_pos (np.array of float): for MCMC only. A multi-D array of the
        current walker positions

Written: Henry Ngo, 2018

API Update: Sarah Blunt, 2021
File:      ~/Documents/GitHub/orbitize/orbitize/results.py
Type:      method

3. (Optional) Load up saved results object

If you are skipping the generation of all orbits because you would rather load from a file that saved the Results object generated above, then execute this block to load it up. Otherwise, skip this block (however, nothing bad will happen if you run it even if you generated orbits above).

import orbitize.results

if "myResults" in locals():
    del myResults  # delete existing Results object
myResults = orbitize.results.Results()  # create empty Results object
myResults.load_results("plotting_tutorial_GJ504_results.hdf5")  # load from file

4. Using our Results object to make plots

In this tutorial, we’ll work with two plotting functions: plot_corner() makes a corner plot and plot_orbits() displays some or all of the computed orbits. Both plotting functions return matplotlib.pyplot.figure objects, which can be displayed, further manipulated with matplotlib.pyplot functions, and saved.

%matplotlib inline
import matplotlib.pyplot as plt

4.1 Corner plots

This function is a wrapper for corner.py and creates a display of the 2-D covariances between each pair of parameters as well as histograms for each parameter. These plots are known as “corner plots”, “pairs plots”, and “scatterplot matrices”, as well as other names.

corner_figure = myResults.plot_corner()

Choosing which parameters to plot

Sometimes, the full plot with all parameters is not what we want. Let’s use the param_list keyword argument to plot only semimajor axis, eccentricity and inclination. Here are the possible string labels for this fit that you can enter for param_list and the corresponding orbit fit parameter:

Label	Parameter name
sma1	semimajor axis
ecc1	eccentricity
inc1	inclination
aop1	argument of periastron
pan1	position angle of nodes (aka longitude of ascending node)
tau1	epoch of periastron passage (expressed as a fraction of orbital period past a specified offset)
mtot	system mass
plx	system parallax

Note: for labels with numbers, the number corresponds to the companion (sma1 is the first object’s semimajor axis, sma2 would be the second object, etc)

corner_figure_aei = myResults.plot_corner(param_list=["sma1", "ecc1", "inc1"])

Limiting which samples to display

By picking out the panels we show, the plot can be easier to read. But in this case, we see that the plotted x-range on semimajor axis does show the posterior very well. This function will pass on all corner.corner() keywords as well. For example, we can use corner’s range keyword argument to limit the panels to only display 95% of the samples to avoid showing the long tails in the distribution.

corner_figure_aei_95 = myResults.plot_corner(
    param_list=["sma1", "ecc1", "inc1"], range=(0.95, 0.95, 0.95)
)

For other keywords you can pass to corner, see the corner.py API.

Making single variable histogram plots

One use of the param_list keyword is to just plot the histogram for the distribution of one single orbit fit parameter. We can do this by just providing one single parameter.

histogram_figure_sma1 = myResults.plot_corner(param_list=["sma1"], range=(0.95,))

Axis label customization

The axis labels seen on the above plots are the default labels that orbitize passes to corner.py to make these plots. We can override these defaults by passing our own labels through the labels keyword parameter as per the corner.py API.

Note that the length of the list of labels should match the number of parameters plotted.

# Corner plot with alternate labels
corner_figure_aei_95_labels = myResults.plot_corner(
    param_list=["sma1", "ecc1", "inc1"],
    range=(0.95, 0.95, 0.95),
    labels=("SMA (AU)", "eccen.", "inc. (deg)"),
)

Overplotting best fit (“truth”) values

One feature of corner.py is to overplot the contours and histograms with a so-called “truth” value, which we can use for many purposes. For example, if we are sampling from a known distribution, we can use it to plot the true value to compare with our samples. Or, we can use it to mark the location of the mean or median of the distribution (the histogram and contours make it easy to pick out the most likely value or peaks of the distribution but maybe not the mean or median). Here is an example of overplotting the median on top of the distribution for the full corner plot.

import numpy as np

median_values = np.median(myResults.post, axis=0)  # Compute median of each parameter
range_values = (
    np.ones_like(median_values) * 0.95
)  # Plot only 95% range for each parameter
corner_figure_median_95 = myResults.plot_corner(
    range=range_values, truths=median_values
)

Overall, we find that some of the parameters have converged well but others are not well constrained by our fit. As mentioned above, the output of the plot_corner() methods are matplotlib Figure objects. So, if we wanted to save this figure, we might use the savefig() method.

corner_figure_median_95.savefig("plotting_tutorial_corner_figure_example.png")

4.2 Orbit Plot

The plot_orbits method in the Results module allows us to visualize the orbits sampled by orbitize!. The default call to plot_orbits will draw 100 orbits randomly chosen out of the total orbits sampled (set by parameter num_orbits_to_plot). In addition, to draw each of these orbits, by default, we will sample each orbit at 100 evenly spaced points in time throughout the orbit’s orbital period (set by parameter num_epochs_to_plot). These points are then connected by coloured line segments corresponding to the date where the object would be at that point in the orbit. By default, orbits are plotted starting in the year 2000 (set by parameter start_mjd) and are drawn for one complete orbital period. We usually choose to begin plotting orbits at the first data epoch, using this keyword as illustrated below.

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)
)
orbit_figure.savefig("example_orbit_figure.png", dpi=250)

<Figure size 1008x432 with 0 Axes>

Customizing orbit plot appearence

In the above figure, we see that the x and y axes are RA and Dec offsets from the primary star in milliarcseconds. By default the axes aspect ratio is square, but we can turn that off. This is normally not recommended but for some situations, it may be desirable.

(Note: Each call to plot_orbits selects a different random subset of orbits to plot, so the orbits shown in the figures here are not exactly the same each time).

orbit_figure_non_square = myResults.plot_orbits(start_mjd=epochs[0], square_plot=False)

<Figure size 1008x432 with 0 Axes>

The colorbar shows how the line segment colors correspond to the date, beginning at the first data epoch. We can also turn off the colourbar.

orbit_figure_no_colorbar = myResults.plot_orbits(
    start_mjd=epochs[0], show_colorbar=False
)

<Figure size 1008x432 with 0 Axes>

We can also modify the color and ending-epoch of the separation/position angle panels as follows:

orbit_figure_no_colorbar = myResults.plot_orbits(
    start_mjd=epochs[0], sep_pa_color="lightblue", sep_pa_end_year=2100.0
)

<Figure size 1008x432 with 0 Axes>

Choosing how orbits are plotted

Plotting one hundred orbits may be too dense. We can set the number of orbits displayed to any other value.

orbit_figure_plot10 = myResults.plot_orbits(start_mjd=epochs[0], num_orbits_to_plot=10)

<Figure size 1008x432 with 0 Axes>

We can also adjust how well we sample each of the orbits. By default, 100 evenly-spaced points are computed per orbit, beginning at start_mjd and ending one orbital period later. Decreasing the sampling could lead to faster plot generation but if it is too low, it might not correctly sample the orbit, as shown below.

orbit_figure_epochs10 = myResults.plot_orbits(
    start_mjd=epochs[0], num_epochs_to_plot=10, num_orbits_to_plot=10
)

<Figure size 1008x432 with 0 Axes>

In this example, there is only one companion in orbit around the primary. When there are more than one, plot_orbits will plot the orbits of the first companion by default and we would use the object_to_plot argument to choose a different object (where 1 is the first companion).

4.3 Residual plots

It is often useful to look at the residuals to the orbit fits, to more carefully understand how well the model and data agree. This can be used to spot systematics in the data, or fitting issues in getting the model to fit the data. There are additional plotting features in the orbitize.plot module, including plot_residuals(). Currently, plot_residuals() can be used to look at the residuals to the relative astrometry data in separation and PA space. Submit a new issue if you want to see other residual plots.

In the example before, we modify the start_mjd and sep_pa_end_year keywords to zoom in on the segment with data.

import orbitize.plot

orbitize.plot.plot_residuals(myResults, start_mjd=epochs[0] - 365, sep_pa_end_year=2014)

As you can see, the plot shows several models randomly drawn from the posterior. So what are the residuals with respect to? The residuals are with respect to the median sep-vs-time and PA-vs-time trajectories from the random draws. This is true for both the data points and for the models. The residuals of the model with respect to the median model is another way to visualize the spread in the model predictions. If you only plotted one model (by setting the keyword num_orbits_to_plot=1), the data residuals would by with respect to that model, and the model residuals would be 0 throughout.

4.4 Working with matplotlib Figure objects

The idea of the Results plotting functions is to create some basic plots and to return the matplotlib Figure object so that we can do whatever else we may want to customize the figure. We should consult the matplotlib API for all the details. Here, we will outline a few common tasks.

Let’s increase the font sizes for all of the text (maybe for a poster or oral presentation) using matplotlib.pyplot. This (and other edits to the rcParams defaults) should be done before creating any figure objects.

plt.rcParams.update({"font.size": 16})

Now, we will start with creating a figure with only 5 orbits plotted, for simplicity, with the name orb_fig. This Figure object has two axes, one for the orbit plot and one for the colorbar. We can use the .axes property to get a list of axes. Here, we’ve named the two axes ax_orb and ax_cbar. With these three objects (orb_fig and ax_orb, and ax_cbar) we can modify all aspects of our figure.

orb_fig = myResults.plot_orbits(start_mjd=epochs[0], num_orbits_to_plot=5)
ax_orb, ax_sep, ax_pa, ax_cbar = orb_fig.axes

<Figure size 1008x432 with 0 Axes>

First, let’s try to add a figure title. We have two options. We can use the Figure’s suptitle method to add a title that spans the entire figure (including the colorbar).

orb_fig.suptitle("5 orbits from GJ-504 fits")  # Adds a title spanning the figure
orb_fig

Alternatively, we can just add the title over the Ra/Dec axes instead.

orb_fig.suptitle("")  # Clears previous title
ax_orb.set_title("5 orbits from GJ-504 fits")  # sets title over Axes object only
orb_fig

We can also change the label of the axes, now using matplotlib.Axes methods.

ax_orb.set_xlabel("$\Delta$(Right Ascension, mas)")
ax_orb.set_ylabel("$\Delta$(Declination, mas)")
orb_fig

If we want to modify the colorbar axis, we need to access the ax_cbar object

ax_orb.set_xlabel(‘\(\Delta\)RA [mas]’) ax_orb.set_ylabel(‘\(\Delta\)Dec [mas]’) # go back to what we had before ax_cbar.set_title(‘Year’) # Put a title on the colorbar orb_fig

We may want to add to the plot itself. Here’s an exmaple of putting an star-shaped point at the location of our primary star.

ax_orb.plot(0, 0, marker="*", color="black", markersize=10)
orb_fig

And finally, we can save the figure objects.

orb_fig.savefig("plotting_tutorial_plot_orbit_example.png")