How to write F2PY code that can call a Python function

Introduction

F2PY is a Fortran to Python interface generator contained within NumPy that allows Fortran code to be wrapped into a module that can be called from Python. This is particularly useful when you want to write fast numerical code, because Fortran code often runs considerably faster than the equivalent Python code.

If you are new to F2PY, I recommend this tutorial. The official documentation for F2PY can be found here.

One thing that I found difficult, when starting to use F2PY, was the comparative lack of documentation around how to write call-backs, in Fortran, to Python functions. Call-backs involving arrays were particularly badly documented. I have found very few working examples, and some outright misinformation, such as the commenter on this StackOverflow question who claims that it simply isn't possible to write a Fortran subroutine that takes arbitrary arrays as input and output and use it with F2PY. This is possible. In particular, it is possible to write a Fortran subroutine that calls an external (Python) function that takes an arbitrary array as input and output. A working example is included towards the end of this page.

A basic call-back

A call-back happens when we pass a function to another function. For example, in Python we could write:

def evaluate_at_0(func):
    return func(0)

This function takes another function func and returns the value of func at 0.

Suppose we wanted to write an evaluate_at_0 subroutine in Fortran, and have it be called from Python. The official documentation gives some examples of simple call-back subroutines in Fortran 70, but if you are writing your own Fortran code then there is a good chance you are using Fortran 90, so I will give an example that takes advantage of Fortran 90 formats.

subroutine evaluate_at_0(func, f_val)
    implicit none
! f2py intent(callback) func
    external :: func
    real :: func
    real, intent(out) :: f_val

    f_val = func(0)
end subroutine

If we save this as simple_callback.f90 then we can compile using F2PY by typing the following into a terminal:

python -m numpy.f2py -c simple_callback.f90 -m simple_callback

This creates a module that can be called in Python as follows:

import simple_callback

def func(x):
    return x * x

print(simple_callback.evaluate_at_0(func))

External functions with array input of fixed size

Building on our previous example, we might wish to write a Fortran subroutine that takes an external function with array input. For example, the following code can be used to evaluate func at x_val, which is an array of size 2:

subroutine vec_to_scalar(func, x_val, f_val)
    implicit none
! f2py intent(callback) func
    external :: func
    real :: func
    real, dimension(2), intent(in) :: x_val
    real, intent(out) :: f_val

    f_val = func(x_val)
end subroutine

We save this code as array_input.f90, and compile as before:

python -m numpy.f2py -c array_input.f90 -m array_input

Now we can call this code in Python:

import array_input
import numpy as np

def func(x):
    return x[0] + x[1]

x = np.array([1,5])
output = array_input.vec_to_scalar(func, x)
print(output)

We find that the output is 6.0, as expected.

External functions with array output of fixed size

Let us modify vec_to_scalar so that it will produce vector output. The most intuitive way to do so would be as follows:

subroutine this_has_a_bug(func, x_val, f_val)
  implicit none
! f2py intent(callback) func
  external :: func
  real :: func
  real, dimension(2), intent(in) :: x_val
  real, dimension(2), intent(out) :: f_val

  f_val = func(x_val)
end subroutine

As you might guess from the spoiler-y function name, this doesn't actually work. It will compile just fine, however, if we save it as array_input_and_output.f90 and type this in the command line:

python -m numpy.f2py -c array_input_and_output.f90 -m array_input_and_output

We can even call it from Python:

import array_input_and_output
import numpy as np

def func(x):
    return x

x = np.array([1,5])
output = array_input_and_output.this_has_a_bug(func, x)
print(output)

This code runs with no error messages! Unfortunately, the printed output is [1. 1.] rather than the expected [1. 5.]. Yikes. This is a serious problem. Note that it can also arise with functions that have scalar input and array output.

There is more than one way to invoke an external function in Fortran. Above, we have done it like this:

  f_val = func(x_val)

We could also have used the following syntax, however:

  call func(x_val, f_val)

When a function has array output, we need to use the latter way of invoking func. For technical reasons, this means we should no longer include the real :: func type declaration:

subroutine this_will_work(func, x_val, f_val)
  implicit none
! f2py intent(callback) func
  external :: func
  real, dimension(2), intent(in) :: x_val
  real, dimension(2), intent(out) :: f_val

  call func(x_val, f_val)
end subroutine

Compiling in the usual way, and replacing this_has_a_bug with this_will_work in the Python code above will give the desired output of [1. 5.].

Modifying the signature file

The this_will_work function defined above is very easy to break. In particular, suppose we wanted to modify f_val somehow, instead of just outputting it directly. For example, suppose we want to multiply the output by 2. We might do this by defining a new value f_temp to be our output from func, and then use this value to calculate our output f_val:

subroutine new_bug(func, x_val, f_val)
  implicit none
! f2py intent(callback) func
  external :: func
  real, dimension(2), intent(in) :: x_val
  real, dimension(2), intent(out) :: f_val
  real, dimension(2) :: f_temp

  call func(x_val, f_temp)
  f_val = 2 * f_temp

end subroutine

If you compile this, and call it in Python, in the way shown above, you will get ... actually, I cannot tell you what you will get. You'll get no error messages. You also probably won't get the same answer twice.

What has gone wrong?

A call-back that uses arrays of arbitrary size