Main Data History
Show Index Toggle 0 comments
  •  Quick Edit
  • poliastro: An Astrodynamics library written in Python with Fortran performance


    Python is a fast-growing language both for astronomic applications and for educational purposes, but it is often criticized for its suboptimal performance and lack of type enforcement. In this paper we present poliastro, a pure Python library for Astrodynamics that overcomes these obstacles and serves as a proof of concept of Python strengths and its suitability to model complex systems and implement fast algorithms.

    poliastro features core Astrodynamics algorithms (such as resolution of the Kepler and Lambert problems) written in pure Python and compiled using numba, a modern just-in-time Python-to-LLVM compiler. As a result, preliminary benchmarks suggest a performance increase close to the reference Fortran implementation, with negligible impact in the legibility of the Python source. We analyze the effects of these tools, along with the introduction of new ahead-of-time compilers for numerical Python and optional type declarations, in the interpreted and dynamic nature of the language.

    poliastro relies on well-tested, community-backed libraries for low level astronomical tasks, such as astropy and jplephem. We comment the positive outcomes of the new open development strategies and the permissive, commercial-friendly licenses omnipresent in the scientific Python ecosystem.

    While recent approaches involve writing Python programs which are translated on the fly to lower level code, traditional Python libraries for scientific computing have succeeded because they leverage computing power to compiled languages. We briefly present tools to build wrappers to Fortran, C/C++, MATLAB and Java, which can be also useful for validation and verification, reusability of legacy code and other purposes.



    The Python programming language has seen broad recognition in the last decade among scientists and academics, being one of the most popular languages in astronomy(Robitaille et al., 2013) and for educational purposes(Guo, 2014). It is highly trusted in corporative environments as well as a tool for scripting, automating tasks and creating high level APIs and wrappers.

    However, the current situation of scientific codes is still delicate: most scientists and engineers that write numerical software (which often features a strong algorithmic component and has tight performance requirements) usually do not have any formal training on computer programming, let alone software engineering best practices(Wilson et al., 2014). In fact, in the Aerospace industry there is a sad track record of software failures(Albee et al., 2000; Lions et al., 1996) that could have been avoided by following better software and systems engineering practices.

    When selecting a certain programming language for a specific problem, we as engineers have the obligation to consider as much information as possible and make an informed decision based on technical grounds. For example, if defect density were to be selected as the single figure to rank the contenders, well-established languages for space applications such as FORTRAN or C would perform worse than functional languages such as Haskell or Erlang(Ray et al., 2014). Another common misconception is to assume that each language features certain properties, while languages are abstract specifications and language implementations are the concrete systems we can measure. Other metrics that could be taken into account are readability and programmer productivity, specially considering that “programmers write roughly the same number of lines of code per unit time regardless of the language they use”(Wilson et al., 2014).

    In this paper we claim that the Python programming language, with the aid of both young projects and solid, well tested libraries, can be an optimal solution for the prototyping stage of the development and a fair complement to traditional alternatives in the production stage, in terms of performance, availability, maturity and maintainability. As a demonstrator we selected basic problems in Astrodynamics and compared the performance of existing FORTRAN or C++ implementations with our new Python implementations, comparing them in terms of code complexity and performance in section \ref{sec:python}. Our Python code is available as part of the f2py package, an open source Python library for Astrodynamics and Orbital Mechanics focused on interplanetary applications and released under the MIT license(Cano Rodríguez 2016). As a complement, in section \ref{sec:interface} we present an overview of the techniques that can be used to use code written in Fortran, C/C++, Java and MATLAB from Python, and their advantages and tradeoffs.

    To conclude, we remark that JCC and many other software packages would not be possible without the vast number of open source projects that lay the foundations for present and future work. Other authors have highlighted the potential of open source in the aerospace industry in terms of software reusability and collaboration between academia and private companies(Ziemer et al., 2012). In section \ref{sec:development} we comment the practical outcomes of this philosophy, the challenges that need to be solved and its potential ramifications for the Aerospace industry.

    Python as a core computational language


    The Python programming language was started by Guido van Rossum in 1989 as a successor to the ABC language, and v1.0 was released in 19941. It is therefore not new, and in fact it predates the Java language, first released in 1996. On the other hand, Python first uses for scientific purposes appeared as early as 1995, with the creation of a special interest group on numerical arrays(Millman 2011). However, in recent times the ecosystem has greatly improved, with the application of Open Development principles (see \ref{sec:development} for further discussion), the increasing interest and involvement of private companies and the generous funds given to projects like IPython(Perez 2007) and Jupyter. Nowadays, it is one of the most used languages in fields like Astronomy(Momcheva 2015) and small-to-medium Data Science, and heavily trusted for teaching undergraduate Computer Science in top universities(Guo 2014).

    One of the most important differences between Python and compiled languages like Fortran or C is its dynamic typing nature. The variety of type systems across has traditionally been a major source of debate among programmers, and in fact some studies suggest that there is “a small but significant relationship between language class and defects”(Ray 2014). Languages featuring dynamic typing, as it is the case with Python, are often easier to write and read but more difficult to debug, as there are no guarantees about the types of the arguments. While there is an increasing interest in developing type inference systems (see for instance the Julia and Scala languages), these are extremely difficult to set up for languages like Python (Cannon 2005).

    In the following sections we analyze some alternatives that have been put in place to overcome these limitations without greatly affecting the philosophy of the language.

    Just-in-time compilation using numba

    As discussed earlier, while it is possible to use NumPy(Walt 2011) to vectorize certain kinds of numerical operations, there might be other cases where this may not be feasible and where the dynamic nature of Python leads to a performance penalty, specially when the algorithm involves several levels of nested looping. To overcome these limitations we used poliastro , an open source library which can infer types for array-oriented and math-heavy Python code and generate optimized machine instructions using the LLVM compiler infrastructure(Team 2016).

    Benchmarks against Fortran

    To test the suitability of Python and for writing expensive mathematical algorithms in terms of performance and legibility we compare two algorithms for solving Lambert’s problem: a simple bisection iteration over an universal variable as presented in (Bate 1971) and (Vallado 2001) (from now on, referred to as BMW-Vallado algorithm) and the more recent algorithm by (Izzo 2014), based on a Householder iteration scheme over a Lambert-invariant variable (see (Gooding 1990) for the definition). These two algorithms are very different in nature: the former favors a simple approach that is robust and simple to implement, while the latter employs clever analytical transformations and a higher order root finding method to converge in very few iterations, hence using fewer function evaluations. Also, the former works only for single revolution solutions, whereas the latter can also find solutions corresponding to multiple revolutions.

    The data was plotted using matplotlib and seaborn(Giuca; 2016)(Allan; 2014).

    Both algorithms implemented in Python and accelerated using are present in , a MIT-licensed open source library dedicated to problems focused on interplanetary Astrodynamics problems, such as orbit propagation, solution of Lambert’s problem, conversion between position and velocity vectors and classical orbital elements and orbit plotting. documentation and source code are available online23.

    Bate-Mueller-White-Vallado algorithm

    Several implementations of the BMW-Vallado algorithm are freely available on the Internet as the companion software of Vallado’s book4, and the Fortran version was incorporated in v0.2 and distributed under the same terms with explicit permission of the author (see \ref{sec:interface} for discussion on calling Fortran functions from Python). In v0.3 the algorithm was translated to Python and accelerated with .

    while count < numiter:
        y = norm_r0 + norm_r + A * (psi * c3(psi) - 1) / c2(psi)**.5
        # ...
        xi = np.sqrt(y / c2(psi))
        tof_new = (xi**3 * c3(psi) + A * np.sqrt(y)) / np.sqrt(k)
        if np.abs((tof_new - tof) / tof) < rtol:  # Convergence check
            count += 1
            if tof_new <= tof:  # Bisection check
                psi_low = psi
                psi_up = psi
            psi = (psi_up + psi_low) / 2