{'fraction_non_alphanumeric': 0.033829907599472, 'fraction_numerical': 0.011238921365264944, 'mean_word_length': 4.862480654432899, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
orbitize!: A Comprehensive Orbit-fitting Software Package for the High-contrast Imaging Community

orbitize! is an open-source, object-oriented software package for fitting the orbits of directly-imaged objects. It packages the Orbits for the Impatient (OFTI) algorithm and a parallel-tempered Markov Chain Monte Carlo (MCMC) algorithm into a consistent and intuitive Python API. orbitize! makes it easy to run standard astrometric orbit fits; in less than 10 lines of code, users can read in data, perform one fit using OFTI and another using MCMC, and make two publication-ready figures. Extensive pedagogical tutorials, intended to be navigable by both orbit-fitting novices and seasoned experts, are available on our documentation page. We have designed the orbitize! API to be flexible and easy to use/modify for unique cases. orbitize! was designed by members of the exoplanet imaging community to be a central repository for algorithms, techniques, and know-how developed by this community. We intend for it to continue to expand and change as the field progresses and new techniques are developed, and call for community involvement in this process. Complete and up-to-date documentation is available at orbitize.info.


INTRODUCTION
By repeatedly imaging exoplanets, we can directly observe them moving along their orbits. The physics behind orbits is well-established and straightforward to computationally model, and can reveal much about the properties and dynamical histories of planetary systems. For example, misalignment of the orbital plane with the stellar spin axis can indicate disturbances early in the lifetime of a system before planets formed (Bate et al. 2010;Maire et al. 2019). In systems with circumstellar dust, a planet can warp a disk when the two are mutually inclined (Lagrange et al. 2010;Dawson et al. 2011) or scatter comets into the inner planetary system (Thébault & Beust 2001;Zieba et al. 2019). Planet-planet scattering and resonant migration in a gas disk can excite observable eccentricities (Yu & Tremaine 2001;Scharf & Menou 2009), potentially implying the presence of unseen planets in planetary systems and constraining formation scenarios. Improved astrometric precision (e.g., Gravity Collaboration et al. 2019) could soon lead to the detection of unseen planets based on perturbations in the orbits we observe. In the future, accurate orbit models of exo-Earths from future space imaging missions will be critical to properly assess their climates and habitability (Williams & Pollard 2002).
While the physics of orbits is straightforward, orbitfitting is challenging, especially for the current generation of directly imaged planets. The wide angular separations needed to detect these objects translate to decades-long or longer orbital periods, which means that the parameter space of possible orbits is often large. Orbit-fitting packages such as ExoFast (Eastman et al. 2013) and PyAstrOFit (Wertz et al. 2017) explore these parameter spaces using a Markov-Chain Monte Carlo (MCMC) approach, but such algorithms often converge slowly in this regime. Several specialized techniques have been developed in the past few years to address some of these difficulties (e.g. Blunt et al. 2017, O'Neil et al. 2019, Brandt et al. 2018), but it is often left as an exercise to the reader to implement, debug, and combine these techniques.
In this paper, we present orbitize!, an open-source orbit-fitting software package inspired by radvel (Fulton et al. 2018), designed to meet the needs of the highcontrast imaging community. orbitize! is designed to consolidate the algorithms, techniques, and knowhow of the high-contrast orbit-fitting community in one place. It is fast and robust, but also clearly written, well-documented, and easy to use. We have designed orbitize! to be flexible and easily modifiable so it can grow with the field of high-contrast orbit-fitting. We believe orbitize! to be a code base with a comparatively low barrier to understanding and contribution. We seek to remove obstacles to becoming an expert in directimaging orbit-fitting, enabling the field to advance more rapidly.
This paper is organized as follows: in Section 2, we review Bayesian orbit-fitting for directly imaged astrometry and discuss the various algorithms used to perform this procedure. In Section 3, we outline the design of orbitize! and give examples of how to use the code. In Section 4, we discuss community involvement guidelines and provide a list of items on the orbitize! to-do list, and we conclude in Section 5.

Defining the Orbit
Orbits in orbitize! are parametrized using the angle conventions from Green (1985): semi-major axis (a), eccentricity (e), inclination angle (i), argument of periastron of the secondary's orbit (ω), longitude of ascending node (Ω), epoch of periastron passage (τ ), parallax, and total mass. Note that we express epoch of periastron passage (τ ) as a fraction of the orbital period past a specified reference date t ref (default January 1, 2020): where t 0 is the time of periastron and P is the orbital period. We chose to fit in τ rather than t 0 because the prior bounds for τ are straightforward (between 0 and 1) no matter the orbital period. A screen-capture from an interactive module intended to help readers visualize Screen-capture of an animated, interactive Python module intended to help users visualize orbits and the orbitize! coordinate system. Top: RA vs decl. The orbit's line of nodes is shown as a dashed line, and the primary's location as a pink star. The dark purple arc shows the portion of the orbit in front of the sky plane, and the lighter purple arc shows the portion that is behind the sky plane. Middle: stellar radial velocity versus time. Bottom: interactive sliders that users can drag to set the values of orbital elements in the interactive version. This figure is available as a video in the online version of this article. In addition, readers are encouraged to download and interact with this visualization at this link.
these parameters is shown in Figure 1, and an interactive version is available online 1 . In the orbitize! coordinate system, motion along the positive z-direction causes a redshift. The positive xdirection is in the positive R.A. direction offset from the primary star, and the positive y-direction is in the positive decl. direction. The orbital elements are defined as usual within this reference frame, with i=0 • corresponding to a face-on orbit. We caution, however, that users must be careful comparing the outputs of RV orbitsolving packages like radvel with those of orbitize!, since radvel fits the orbit of the star, and orbitize! fits the orbit of the planet. In practice, this just means adding π to any argument of periastron (ω) values returned by an RV code. See Figure 1 of Fulton et al. (2018) for a visualization of this difference. a Expressed as fraction of orbital period past Jan 1, 2020

Bayesian Orbit Fitting
For high-contrast imaging, Bayesian orbit-fitting is the process of converting time series measurements of a planet's location relative to its primary to a posterior over its orbital parameters. The inputs to this process are times at which measurements were taken and relative position measurements, most often expressed as a planet-star separation (ρ) and position angle (θ), or as offsets in R.A. (α) and decl. (δ).
To compute the posterior over orbital parameters, we use Bayes' theorem: where p(m|D) is the posterior, L(m|D) is the likelihood (the probability of the orbit model given the data), and p(m) is the prior (the probability of the orbit model itself). By default, orbitize! uses a Gaussian likelihood: where N is the number of observation epochs with observations measured in terms of ρ and θ, M is the number of observation epochs with observations measured in terms of α and δ, t i is the epoch of the ith observation, ρ m (t i ) is the separation predicted by the model orbit at t i , ρ o (t i ) is the observed separation at t i , and σ ρo(ti) is the observational uncertainty on the observed separation. The inverse tangent function in Equation 4 accounts for angle wrapping near 0 and 2π to ensure that the difference between the model and observed value is are determined by solving Kepler's equation, discussed further in Section 3.1.3. Total mass and parallax are either included as free parameters and assigned priors motivated by observations or held fixed in this analysis. In orbitize!, users can select either option.

Algorithms for Orbit-fitting
Prior to more modern Bayesian techniques, wellconstrained visual orbits were typically fit with leastsquare minimization techniques (e.g. Binnendijk 1960, Sozzetti et al. 1998). While such methods are computationally efficient, and they effectively find the maximum likelihood orbit with estimates of the 1D errors and 2D covariances, they do not produce full posteriors. This is particularly problematic for poorly constrained orbits, when the posteriors are distinctly non-Gaussian.
A more computationally intensive method is Least-Square Monte Carlo (LSMC), which generates input orbits by drawing randomly from a multi-dimensional parameter space of orbital parameters. Each of the input orbits is run through an iterative χ 2 minimizer until a stopping criterion is reached (e.g. Mugrauer et al. 2012). Several recent publications (Ginski et al. 2013;Vigan et al. 2016;Maire et al. 2019) have demonstrated that that while LSMC is capable of finding families of plausible orbits, its output is often significantly different from the Bayesian posterior. Another Monte Carlo based method was used in Wagner et al. (2016) and Wagner et al. (2018) to fit short orbital arcs. This method used a multi-step grid-search algorithm to find the maximum likelihood orbit, and Gaussian point estimates to approximate confidence intervals.
Both LSMC and grid search algorithms have difficulty fitting systems with data covering only short orbital arcs, in part because both prioritize finding the set of orbital parameters corresponding to the maximum likelihood solution. In a well-constrained orbit where the posterior is a six-dimensional Gaussian, either method should find the best-fitting parameters and recover the marginalized uncertainty in each parameter. When the posteriors are decidedly non-Gaussian, however (which is much more often the case for orbit fitting in direct imaging), both techniques fail to adequately derive the shape of the posterior. Additionally, neither of these methods appears to offer a gain in computational efficiency compared to Orbits for the Impatient (OFTI; see below) or MCMC.
In the following two subsections, we describe OFTI and MCMC, the two backend algorithms available in orbitize!. These algorithms represent efficient Bayesian methods for producing plausible sets of orbital parameters that represent the full multi-dimensional posterior of orbital parameters. These posteriors are robust probability density functions with confidence intervals that are an accurate reflection of the constraints imposed by the data.

Orbits for the Impatient (OFTI)
The OFTI algorithm is described in detail in Blunt et al. (2017), but we review it briefly here. Trial orbits are drawn randomly from priors, and are "scaledand-rotated" to match the data. Scaling-and-rotating involves modifying the semimajor axis and position angle of nodes of the trial orbit to match the most constrained astrometric data point within the observational uncertainties, cutting down the large parameter space of possible orbits. For each scaled-and-rotated trial orbit, Kepler's equation is solved (See Section 3.1.3) and a likelihood is computed. Finally, each orbit is either accepted or rejected by comparing the likelihood probability to a uniform random number.
The OFTI algorithm is most efficient when the orbital posteriors are similar to the priors, or in other words when the parameter space is relatively unconstrained. This occurs most often for short orbital arcs, when the data span only a small fraction of the total orbital period.
In the case of OFTI, each individual orbit considered is uncorrelated with the rest, and so the only stopping criterion is the number of samples desired, which varies by application. To plot plausible orbit tracks, ∼100 sets of orbital parameters are sufficient. ∼1000 may be needed for accurate medians and 1σ confidence intervals on each parameter. Plots of 1D marginalized posteriors are wellsampled with ∼ 10 4 sets, and 2D parameters with ∼ 10 6 .

Affine-Invariant & Parallel Tempered MCMC
Markov-chain Monte Carlo (MCMC) algorithms are commonly used to sample the posterior of planetary orbits. orbitize! makes use of two such algorithms: the Affine-invariant sampler from emcee (Foreman-Mackey et al. 2013), and the parallel-tempered MCMC (PTM-CMC) sampler from ptemcee (Vousden et al. 2016). Given the complex covariances and often multi-modal posteriors of the orbits of directly-imaged planets, the Affine-invariant sampler alone generally fails to fully sample the posterior without fine-tuned starting locations for the walkers. We offer the use of ptemcee, which runs multiple Affine-invariant samplers with different likelihood weights, as an alternative to overcome this difficulty.
MCMC algorithms generally have similar run times for orbits with partial phase coverage (Blunt et al. 2017). However, convergence time is cut down significantly if the orbital elements are well-constrained, such that the posterior is close to a multivariate Gaussian.
Unlike OFTI, the MCMC algorithms coded in orbitize! require an initial period for the walkers to fully converge before they sample the posterior in an unbiased fashion. Convergence is assessed using a combination of the walkers' autocorrelation time as recommended by Foreman-Mackey et al. (2013) and visual inspection of the walker positions over time to determine if they are fully exploring parameter space. Unconverged posteriors typically appear "lumpy," as different chains are still in different regions of parameter space, compared to generally smooth (though not necessarily Gaussian) posteriors for converged chains. See Figure 4 of Blunt et al. (2017) for an illustration. The multimodal posteriors of ω and Ω can also illustrate convergence. Without a measurement of radial position or velocity, values 180 • apart in both ω and Ω are degenerate. Thus, the chains are more likely to have converged when the one-dimensional marginalized posteriors on ω and Ω show symmetric equal peaks.

CODE DESIGN
3.1. API orbitize! comprises several modules, each a separate code file. This section refers to orbitize! version 1.6.0, but the most up-to-date documentation can be found online 2 . The primary modules are:

Input Data
The orbitize.read input module is designed to read astrometric measurements as input to orbitize! in any file format supported by astropy.io.ascii.read(), including csv, fixed width, cds, and LaTeX formats.
The main method in this module is orbitize.read input.read file().
This module also contains the method orbitize.read input.write orbitize input() which takes a table of measurements formatted in the orbitize! format and writes it as an ASCII file in any file format supported by astropy.io.ascii.read().

Priors
Priors in orbitize! are represented as subclasses of the abstract super class orbitize.priors.Prior. Each prior class can draw random values from its distribution as well as compute the probability of drawing particular values from that prior distribution. Other parts of the code use these two methods to interface with the priors. Currently, orbitize! implements uniform, Gaussian, log-uniform (Jeffreys), sine, and linear priors. The default priors listed in Table 1 are used for orbit fits, but users can easily replace priors by instantiating their own priors and replacing them in the orbitize.system.System instance. 3

Kepler's Equation Solver
The orbitize.kepler module converts orbital elements to position and radial velocity measurements by generating an orbital ellipse, placing the companion at the appropriate phase of the orbit, and rotating the orbit based on viewing geometry. The instantaneous position of a companion along an orbital ellipse is defined by the true anomaly ν, the angle between the location of the companion and periastron, as measured from the focus of the orbital ellipse. ν cannot be analytically calculated from a given epoch and set of orbital elements. It is instead calculated from the eccentric anomaly E, the angle between the projection of the companion's location onto a circle of radius a that intersects the true orbit at periastron, which is calculated in turn from the mean anomaly M , the fraction of the period that has elapsed since periastron passage. See Seager (2010) for an excellent derivation of the relevant equations from first principles and visualizations of these quantities.
The conversion between M and E, typically called "Kepler's equation," is given by: This equation cannot be solved analytically. In orbitize!, we solve Kepler's equation using one of two different methods depending on the eccentricity of the orbit. For e < 0.95, we use Newton's method to estimate E with a default tolerance of |E n − E n−1 | < 10 −9 . We use E 0 = M as a starting point and typically achieve the required tolerance in three or four iterations. If a solution to E is not found within a default of fifty iterations, the procedure is restarted with E 0 = π. If a solution is still not found, we instead use the numerical approximation described below to estimate E.
For highly-eccentric orbits (e ∼ 1) the number of iterations required to reach the required tolerance increases dramatically, especially as M approaches 0 or 3 http://orbitize.info/en/latest/tutorials/Modifying Priors.html 2π, significantly reducing the speed of this solver (Fig. 2,  left panel). This can be especially problematic for the parallel-tempered MCMC sampler described in Section 3.1.6, where the high-temperature walkers will explore the full range of eccentricities allowed by the prior distribution. To mitigate this slowdown, we instead use the numerical approximation for Kepler's equation described by Mikkola (1987) for orbits with e ≥ 0.95. Mikkola's algorithm invokes a cubic approximation to Kepler's equation and a numerically intensive (but single iteration) high-precision correction formula.
We verified the accuracy of Mikkola's method by comparing the value of E calculated for a range of (M, e) combinations using Newton's method (with a tolerance of 10 −18 ) and using Mikkola's method. We find a maximum absolute difference between the eccentric anomaly computed via both methods of 1.6 × 10 −15 (Fig. 2, middle panel), well below the nominal tolerance on the eccentric anomaly solver of |E n − E n−1 | < 10 −9 . While Mikkola's method is strictly more accurate than applying Newton's method with a tolerance of 10 −9 , it is much more computationally expensive, especially for low-eccentricity orbits where it is a factor of ten slower (Fig. 2, right panel). For the vast majority of applications for orbitize!the default tolerance of 10 −9 corresponds to a position angle error of 3.6 × 10 −7 deg for a face-on circular orbit, well below the precision of current astrometric measurements.
These algorithms are implemented in C for maximum computational efficiency. A Python version of this module is also included in orbitize!, and the package will revert to the Python version if the C-solver fails to compile on a user's machine.

Planetary System
The main component of the orbitize.system module is the orbitize.system.System class, which stores observational data, priors, and posteriors for a star-planet system.
To initialize an instance of orbitize.system.System, a user inputs a data table, created using the orbitize.read input module, a total system mass, and a parallax. The init () method automatically initializes a list of orbitize.priors.Prior objects in a standard order, detailed in the online documentation. orbitize.system.System also has a compute model() method that takes in an array of potential orbital parameters and computes model predictions to be compared against the data. This is the method that calls the Keplerian orbit-solver (Section 3.1.3). Note that this abstraction makes the sampler agnostic of physics; a user could use any model computation code as a drop-in replacement. Finally, this class has an attribute results, a list of orbitize.results.Results Figure 2. Performance summary of our two Keplerian solvers. Left: Number of iterations by the Newton solver required to achieve a tolerance of 1e-9 as a function of eccentricity and mean anomaly. The Mikkola solver is used for orbits with e ≥ 0.95 to mitigate this slowdown at high eccentricities. Middle: absolute difference between the eccentric anomaly values computed by the Mikkola and Newton solvers. A tolerance of 1e-18 was used in this experiment. The absolute difference between the two solvers is negligible compared with current astrometric precision. Right: time needed for both solvers to achieve a tolerance of 1e-18. The Mikkola solver achieves significant performance gains over the Newton solver for high eccentricities. objects (Section 3.1.5). In version 1.6.0, orbitize! is limited to single-planet systems. Users can create orbitize.system.System objects for each planet in a multi-planet system, but orbit-fits using this framework currently do not take into account the dynamical effects of other planets on the orbit fit.

Results
The orbitize.results module contains code for the orbitize.results.Results class. Each instance of this class represents a posterior (whether calculated using OFTI or MCMC). The orbitize.sampler.Sampler objects add their output to an instance of this class. This module also controls saving, loading, and visualizing results.
The run sampler() method of the orbitize.sampler.OFTI class generates a posterior of permissible orbits. This method iteratively runs a series of methods that define the process of generating and statistically vetting potential orbits until a desired number of orbits are accepted. Accepted orbits are then added to an instance of the orbitize.Results class. This process is arbitrarily parallelizible, and users can easily set the number of CPU cores available for fits.
orbitize.sampler.MCMC uses MCMC algorithms to generate a Markov Chain representing the posterior. Setting the attribute num temps to 1 invokes the Affine-Invariant sampler, and setting it to a number greater than 1 invokes the parallel tempered sampler. Users also choose the number of walkers, the number of threads, and whether they would like to fix certain parameters so that they are not sampled by the MCMC algorithm. A convenient API for inspecting chains to assess convergence is available, and is explained in the online tutorials and documentation.

Driver
Finally, the orbitize.driver module automates the creation of the data table and orbitize.system.System objects.
It is a convenient short-cut for standard orbit-fits, and allows new users to begin using the code relatively quickly. While using this object makes it simple to run orbitize! in a standard way, users are encouraged to learn about the underlying API in order to learn how to, e.g. modify priors, customize plots, and set parameters specific to OFTI or MCMC.

Example Usage
To illustrate the ease with which a user can run an orbit fit with orbitize!, we provide a minimal code example below. One of the strengths of orbitize! is its customizability, however, and this example is merely intended to show how easy it is to run a "standard" orbit fit, not to give a sense of the full scale of the code's capabilities. We encourage users to peruse the online tutorials 4 for more in-depth examples of how to use and modify orbitize!.
In the example below, an orbitize.driver.Driver object is initialized, and the OFTI algorithm runs until 10,000 orbits have been accepted. To use MCMC instead of the OFTI algorithm, the user just needs to switch out the keyword in line 6. Inline comments are provided to aid understanding.

Community Involvement Guidelines
In order to realize our goals of consolidating direct imaging orbit-fitting best practices in one place and continually adapting orbitize! to serve the direct imaging community, we require substantial community involvement. We strongly encourage anyone interested in using orbitize! to contribute code. We maintain a document with specific contributor guidelines 5 on our GitHub page.
In addition to directly contributing code, we encourage community members to request features, report bugs, and provide feedback through raising issues on GitHub. This is the most efficient way to reach the entire development team.

Future Work
orbitize! has a long and exciting list of planned updates. A version 2.0 release, intended to implement many of these, is planned for 2020. Key upgrades in version 2.0 will be enabling dynamical mass measurements by jointly fitting RV and stellar astrometry datasets, as well as fitting multi-planet systems. To summarize, the major upcoming features are: In this paper, we have presented orbitize!, an opensource Python package for fitting the orbits of directlyimaged planets. orbitize! uses OFTI and MCMC, two efficient, Bayesian methods for computing posteriors. We aim to encourage community contributions, remove barriers to becoming an expert in orbit-fitting, and provide an open-source development environment in order to meet the orbit-fitting needs of the high-contrast exoplanet imaging community.
The authors thank those at academic and telescope facilities whose labor maintains spaces for scientific inquiry, particularly those whose communities are excluded from the academic system. This package was born and developed at the winter 2018 and 2019 AAS Hack Days, and the authors wish to thank the organizers of these events. The authors would also like to thank Dillon Dong, Jasmine Garani, Melisa Tallis, and Daniel Yahalomi for their time and initial work on orbitize!, and Junellie Gonzalez, Kelly Kosmo O'Neil, Ryan Rubenzahl, and Jean-Baptiste Ruffio for participating in our hackathons and for their anticipated future contributions to orbitize! Finally, we thank the anonymous individual who unknowingly named orbitize! at the 2018 AAS Hack Day.
Yes, you need to use the exclamation point!