python.tex

\section{Python as a core computational language}
\label{sec:python}

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 1994\footnote{http://python-history.blogspot.com.es/2009/01/brief-timeline-of-python.html}. 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\cite{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\cite{Perez_2007} and Jupyter. Nowadays, it is one of the most used languages in fields like Astronomy\cite{2015arXiv150703989M} and small-to-medium Data Science, and heavily trusted for teaching undergraduate Computer Science in top universities\cite{guo2014python}.

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"\cite{Ray2014}. 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 \cite{cannon2005localized}.

% Add some examples here showing that Python is slow, comparing with numpy should be enough

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.

\subsection{Just-in-time compilation using numba}

As discussed earlier, while it is possible to use NumPy\cite{van_der_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 \verb|numba|, 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\cite{numba}.

% http://tex.stackexchange.com/a/86901/2488
\newsavebox\myv

\begin{lrbox}{\myv}\begin{minipage}{\textwidth}
\input{examples/annotations.tex}
\end{minipage}\end{lrbox}

% This is for code http://tex.stackexchange.com/a/34065/2488
\begin{figure*}
  \centering
    \resizebox{\textwidth}{!}{\usebox\myv}
  \caption{Fragment of numba annotations in poliastro.}
\end{figure*}

\subsection{Benchmarks against Fortran}

To test the suitability of Python and \verb|numba| 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 \cite{bate1971fundamentals} and \cite{vallado2001fundamentals} (from now on, referred to as BMW-Vallado algorithm) and the more recent algorithm by \cite{Izzo2014}, based on a Householder iteration scheme over a Lambert-invariant variable (see \cite{gooding1990procedure} 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\cite{Droettboom2016}\cite{Waskom2014}.

% May we add some math here?

Both algorithms implemented in Python and accelerated using \verb|numba| are present in \verb|poliastro|, 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. \verb|poliastro| documentation and source code are available online\footnote{https://github.com/poliastro/poliastro}\footnote{http://poliastro.readthedocs.org/en/v0.5.0/}.

\subsubsection{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 book\footnote{http://celestrak.com/software/vallado-sw.asp}, and the Fortran version was incorporated in \verb|poliastro| 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 \verb|poliastro| v0.3 the algorithm was translated to Python and accelerated with \verb|numba|.

\begin{lstlisting}[language=Python]
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
        break
    else:
        count += 1
        if tof_new <= tof:  # Bisection check
            psi_low = psi
        else:
            psi_up = psi
        psi = (psi_up + psi_low) / 2
\end{lstlisting}

% I should add a table and/or a figure here with these results

\subsubsection{Izzo algorithm}

For the preparation of this paper, we also implemented the Izzo algorithm.

\subsection{Gradual typing}

%We have already discussed the avantages of using \verb|numba| to accelerate numerical Python code. We have also seen that whether we manually specify the expected types of our functions or let the computer automatically infer them, in both cases there is some sort of type enforcement, which on the other hand is not part of the language itself and does not support its complete set of features. Python 3.5 introduced a new provisional module adding \textit{type hints} (also known as \textit{gradual typing}), which represent a similar idea. While this new feature has, to our knowledge, not been explored to help accelerating Python code yet, it can already be used by some Integrated Development Environments (IDEs) to supply more useful information to the developer and definitely has the potential to serve a similar purpose than \verb|numba|.

% Rewrite the above

\begin{lstlisting}[language=Python]
def greeting(name: str) -> str:
    return 'Hello ' + name
\end{lstlisting}