Spicing up your Python with a little bit of Fortran - part 1

22 Aug 2018 on Python, Fortran, Hpc

How do you integrate your beautiful Python code with a little bit rough but also incredibly speedy Fortran routines?

Whether you like it or not, for a long time Fortran has been dominating language in scientific computing community and there are still many maintained projects with large Fortran codebases. There are also new projects developed using Fortran from the ground up. Why is it so? Is it a power of habit?

While it is certainly true that some people code in Fortran just because they are used to it, one cannot deny that code generated by most Fortran compilers cannot be outperformed (in terms of speed) by most other languages. Combine it with native support for multidimensional arrays and you can see why scientist and engineers still use this ancient programming language.

Personally, I find Fortran to be a terrible language to code in. Its array indexing starts from 1 (what kind of monster came up with this?), it has some peculiar syntactic constructs (why do you need to call a subroutine but can invoke a function like in any other language?) and so on and so forth. Writing any nontrivial Fortran program, especially one with some user interface and/or IO handling is undoubtedly a daunting task. It would be cool if we could leverage Fortran’s speed in some higher level language - and fortunately we can.

This post is meant as a first part of a series on combining Python with Fortran. I found that resources on the subject are somewhat scattered out there on the Internet, making it hard to find a single guide or series of guides covering everything from beginner level and gradually through more advanced topics. Hopefully the series I’m planning to write could serve as one.

In today’s part we’ll start with the following tasks:

• wrapping a Fortran subroutine using f2py,
• choosing a compiler used by f2py and configuring its flags.

The above should be enough to get you started with some simple routines.

Prerequisites

For the examples showed below I used Python 3.6 but everything should work just fine for 2.7. Besides Python you will also need to have the following installed:

• Python packages: numpy and imageio
• some Fortran compiler

What we would like to achieve

For demonstration purposes we will implement a simple program for generating images of Mandelbrot set. I’ve choosen this particular problem since it is well known and easy to understand, yet requires handling of two-dimensional arrays.

The idea is to code the computational-heavy part, i.e. computing each point’s escape iteration, in Fortran and then write some script in Python to invoke our routine.

If you don’t remember a thing about Mandelbrot here is a quick reminder: consider point $p = p_x + i p_y$ in a complex plane and define a sequence: \[z_n = \begin{cases} p & n = 0\\ z_{n-1}^2 + p & n > 0 \end{cases}\]

The point $p$ lies in the Mandelbrot set if and only if $|z_n| < 2$ for every $n$. For a point $p$ not belonging to the Mandelbrot set, the first $n$ for which the above inequality is violated can be used for computing nice color of this point. This is another whole topic though, we will just use this value of $n$ as a number on a scale of grey.

If preceeding description does not make sense to you, you still can make sense of the forecoming sections - just treat a Fortran code as a black box. You will miss some fun though.

Step 1: writing our Fortran routines

The below is the implementation of the required functions and subroutines:

• function color: given point’s coordinates and maximum iterations compute index of iteration after which it beceomes evident that it is not in the Mandelbrot set. After reaching maximum number of iterations we simply return maxiter.
• subroutine compute_colors: basically run color for every point in the grid and store result in output array. Note that the grid is defined by its upper left and bottom right corners as well as number of pixels to compute - which is the same as shape of the output array.

! File mandelbrot.f90
module mandelbrot
  implicit none
contains

  integer(8) function color(p_x, p_y, maxiter)
    real(8), intent(in) :: p_x, p_y
    integer(8), intent(in) :: maxiter
    real(8) :: x_tmp, y_tmp, x, y

    color = 0
    x = p_x
    y = p_y

    do while (x * x + y * y < 4 .and. color < maxiter)
       x_tmp = x * x - y * y + p_x
       y_tmp = 2 * x * y + p_y
       x = x_tmp
       y = y_tmp
       color = color + 1
    end do
  end function color

  subroutine compute_colors( &
       top, left, bot, right, & ! coordinates of corners
       out_arr, cols,rows, & ! grid size
       maxiter) ! maximum number of iteratoins
    real(8), intent(in) :: top, left, bot, right
    integer(8), intent(in) :: rows, cols
    integer(8), dimension(rows, cols), intent(out) :: out_arr
    integer(8), intent(in) :: maxiter

    integer(8) :: i, j
    real(8) :: x, y

    do i = 1, cols
       ! Compute abscissa of point corresponding to
       ! (i, j element in grid) - basically a convex
       ! combination.
       x =  left * (cols-i)/(cols-1.0) + (i-1)/(cols-1.0) * right
       do j = 1, rows
          ! As above but for the ordinate
          y = (rows-j)/(rows-1.0) * top + (j-1)/(rows-1.0) * bot
          out_arr(j, i) = color(x, y, maxiter)
       end do
    end do
  end subroutine compute_colors
end module mandelbrot

Nothing special here. Now let’s move on to wrapping our routines with f2py.

Step 2: wrapping Fortran routines with f2py

The f2py script (and acoompanying Python modules) are part of numpy. The script basically creates Python extension that calls your Fortran routines. There are three ways of using f2py, described in its documentation. Today we will only use the simple scenario, leaving the more advanced options for another time.

Basic usage of the f2py goes as follows:

f2py -m mandelbrot -c mandelbrot.f90

This should create file mandelbrot.cpython-36m-x86_64-linux-gnu.so in current directory (exact name may differ depending on your OS and Python version). Breaking down the above command:

• the -m argument controls how the created Python module will be named. It is the name that you will import to use your Fortran functions.
• the -c tells f2py to compile given files.

The newly created module should be importable. Exceute the following Python code to examine its contents:

import mandelbrot

help(mandelbrot)

The beginning of output should be similar to the following:

NAME
    mandelbrot

DESCRIPTION
    This module 'mandelbrot' is auto-generated with f2py (version:2).
    Functions:
    Fortran 90/95 modules:
      mandelbrot --- color(),compute_colors().

So our mandelbrot module has a mandelbrot Fortran module that contains two functions. If you didn’t get why we get a module inside module keep in mind that we could pass multiple Fortran files to f2py, each of them containing separate Fortran module. So the top level one is named following -m argument and its children correspond to Fortran modules you just compiled.

Now let’s look what are the signatures of functions that f2py generated. Displaying their __doc__ attributes should give something similar to the output below.

color = color(p_x,p_y,maxiter)

Wrapper for ``color``.

Parameters
----------
p_x : input float
p_y : input float
maxiter : input long

Returns
-------
color : long

Wrapper for ``compute_colors``.

Parameters
----------
top : input float
left : input float
bot : input float
right : input float
cols : input long
rows : input long
maxiter : input long

Returns
-------
out_arr : rank-2 array('q') with bounds (rows,cols)

There are no surprises for the color function but look what happend with compute_colors. The f2py detected that we have a single output argument and wrapped our subroutine as a function that returns an array. Behind the scenes an array will be created for us at runtime and passed to Fortran routine and then returned back wrapped nicely as numpy.ndarray. Creating arrays to be filled by subroutines is OK, but having them delivered to you is so much better!

Step 3: write a Python script that uses it

We’ll go with a simple script that will just pass through command line arguments to the compute_colors.

from argparse import ArgumentParser
import imageio
import numpy as np
from mandelbrot import mandelbrot

def main():
    parser = ArgumentParser()
    parser.add_argument('--width', type=int, default=1920,
                        help='width of the image')
    parser.add_argument('--height', type=int, default=1080,
                        help='height of the image')
    parser.add_argument('-t', type=float, default=1,
                        help='ordinate of top-left corner')
    parser.add_argument('-l', type=float, default=-2,
                        help='abscissa of top-left corner')
    parser.add_argument('-b', type=float, default=-1,
                        help='ordinate of bottom-right corner')
    parser.add_argument('-r', type=float, default=1,
                        help='abscissa of bottom-right corner')
    parser.add_argument('-i', type=float, default=100,
                        help='maximum number of iterations')
    parser.add_argument('output', type=str, default=None,
                        help='path to the output file')

    args = parser.parse_args()
    arr = mandelbrot.compute_colors(args.t, args.l, args.b, args.r,
                                    args.width, args.height, args.i)
    arr = arr / np.max(arr)
    imageio.imwrite(args.output, arr)

if __name__ == '__main__':
    main()

Executing the script like this should produce nice image of Mandelbrot set in shades of gray:

python generate.py --width 1920 --height 1080 mandelbrot.jpg

Example Mandelbrot image generated by our script

Step 4 (optional): additional seasoning

So, suppose that for whatever reason you are planning to generate a lot of big ass images of Mandelbrot set. You are really concerned about performance so you are considering tweaking Fortran compiler flags and maybe throwing some parallelization into the mix. Or maybe something is not working as expected and you wish to include some debug information so you may solve the problem. To accomplish any of those tasks you need to specify compiler flags used by f2py. You can do so in one of two, mostly equivalent, ways:

• by specifying --f90flags command line arguments
• by specifying F90FLAGS environmental variable

So for instance, assuming that you are using gfortran as your Fortran compiler, you might specify -Ofast and -finline-functions optimization flags in one of the following ways:

f2py --f90flags="-Ofast -finline-functions" -c mandelbrot.f90 -m mandelbrot

F90FLAGS="-Ofast -finline-functions" f2py  -c mandelbrot.f90 -m mandelbrot.o

But wait, “assuming that you are using gfortran”?. How do I tell f2py to use specific Fortran compiler? Again, you have basically two options:

• speficy --fcompiler command line argument
• specify F90 environmental variable

This time however things are a bit more tricky. The --fcompiler argument accepts some predefined values representing compilers that f2py can work with. For instance:

• “pg”: Portland Group Fortran compiler
• “intelem”: ifort
• “gnu95”: gfortran

On the other hand, the F90 env variables has to be set to an executable used to compile the code. So for compiling our extension using PG fortran compiler we could do one of the following:

f2py --fcompiler=pg -c mandelbrot.f90 -m mandelbrot

F90=pgfortran f2py -c mandelbrot.f90 -m mandelbrot

Why isn’t using the --fcompiler as intuitive as using F90 environmental variable? Most probably everything boils down to the fact that f2py can work only with predefiend compilers it knows how to handle - you may see that f2py itself adds some compiler flags which differ depending on compiler you use. Still, the names they chose for the possible values are unintuitive at least.

Conclusions

Hopefully after reading this post you can wrap you Fortran routines into nice Python modules and alter copmiler flags used by f2py to suite your needs. There are still some interesting questions that we shall answer someday else, for instance:

• if we would like to include additional libraries in the linking process, how do we do that?
• what if it was not my intention to convert compute_colors to function returning an array and I’d rather pass the target array myself?
• how do I package Python project with Fortran extensions?

We will look into the above topics in the next post in the series.