-
Notifications
You must be signed in to change notification settings - Fork 0
How to write F2PY code that can call a Python function
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 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 77, 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))
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.
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.].
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. We could want to multiply the output by 2, for example. 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, exactly as before, 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?
To answer this question, we will need to compile our code more carefully. We will use what the F2PY documentation refers to as the smart way. This involves first creating a signature file by typing the following into the command line:
python -m numpy.f2py array_input_and_output.f90 -m array_input_and_output -h array_input_and_output.pyf
This will produce a file array_input_and_output.pyf with the following contents:
! -*- f90 -*-
! Note: the context of this file is case sensitive.
python module new_bug__user__routines
interface new_bug_user_interface
subroutine func(x_val,f_temp) ! in :array_input_and_output:array_input_and_output.f90:new_bug:unknown_interface
real dimension(2),intent(in) :: x_val
real dimension(2) :: f_temp
end subroutine func
end interface new_bug_user_interface
end python module new_bug__user__routines
python module array_input_and_output ! in
interface ! in :array_input_and_output
subroutine new_bug(func,x_val,f_val) ! in :array_input_and_output:array_input_and_output.f90
use new_bug__user__routines
external func
real dimension(2),intent(in) :: x_val
real dimension(2),intent(out) :: f_val
end subroutine new_bug
end interface
end python module array_input_and_output
! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/
You may note that there are two parts to this signature file. The first, python module new_bug__user__routines, deals with the call-back. The second, python module array_input_and_output, deals with inputs and outputs of the function(s) defined in the file. Our problem lies in the first part of the file, the part that deals with the call-back. Zooming in a little:
subroutine func(x_val,f_temp) ! in :array_input_and_output:array_input_and_output.f90:new_bug:unknown_interface
real dimension(2),intent(in) :: x_val
real dimension(2) :: f_temp
end subroutine func
What does F2PY know about func? It knows that it has input x_val, and it knows that it has ... something ... f_temp. It does not know that f_temp is supposed to be output. That is what has caused this issue, and now that we have found the problem it is easily fixed. We simply modify the real dimension(2) :: f_temp line in the signature file so that it now says real dimension(2),intent(out) :: f_temp. Here is the new signature file, with that one change made:
! -*- f90 -*-
! Note: the context of this file is case sensitive.
python module new_bug__user__routines
interface new_bug_user_interface
subroutine func(x_val,f_temp) ! in :array_input_and_output:array_input_and_output.f90:new_bug:unknown_interface
real dimension(2),intent(in) :: x_val
real dimension(2), intent(out) :: f_temp
end subroutine func
end interface new_bug_user_interface
end python module new_bug__user__routines
python module array_input_and_output ! in
interface ! in :array_input_and_output
subroutine new_bug(func,x_val,f_val) ! in :array_input_and_output:array_input_and_output.f90
use new_bug__user__routines
external func
real dimension(2),intent(in) :: x_val
real dimension(2),intent(out) :: f_val
end subroutine new_bug
end interface
end python module array_input_and_output
! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/
We now compile F2PY via the signature file by typing the following into the command line:
python -m numpy.f2py -c array_input_and_output.pyf array_input_and_output.f90
Behold, our function new_bug has now been fixed!
In Fortran 90 we can have arrays of arbitrary size, declared by saying dimension(:). Sometimes F2PY can handle this and sometimes it cannot. The following subroutine was written by taking this_will_work, above, and converting dimension(2) to dimension(:):
subroutine sadly_this_will_not_work(func, x_val, f_val)
implicit none
! f2py intent(callback) func
external :: func
real, dimension(:), intent(in) :: x_val
real, dimension(:), intent(out) :: f_val
call func(x_val, f_val)
end subroutine
This would seem to be the obvious way to use arbitrary array sizes in this context, but, sadly, it fails.
Fortunately, this is not the only way to pass an array of arbitrary dimension to F2PY. The section of the F2PY documentation on array arguments gives the following example, written in Fortran 77:
C FILE: ARRAY.F
SUBROUTINE FOO(A,N,M)
C
C INCREMENT THE FIRST ROW AND DECREMENT THE FIRST COLUMN OF A
C
INTEGER N,M,I,J
REAL*8 A(N,M)
Cf2py intent(in,out,copy) a
Cf2py integer intent(hide),depend(a) :: n=shape(a,0), m=shape(a,1)
DO J=1,M
A(1,J) = A(1,J) + 1D0
ENDDO
DO I=1,N
A(I,1) = A(I,1) - 1D0
ENDDO
END
C END OF FILE ARRAY.F
Notably, the comments in this file declare to F2PY that the dimension integers N and M are intent(hide),depend(a) :: n=shape(a,0), m=shape(a,1). What does this mean? It means that this function can be called from Python using just one argument A, and F2PY will figure out for itself what N and M should be, based on the size of A.
We can use a similar technique to rescue our function above. Specifically, if we want the input and output to be the same (arbitrary) size, we can do this as follows:
subroutine arbitrary_size_array_callback(func, x_val, f_val, n)
implicit none
! f2py intent(callback) func
external :: func
real, dimension(n), intent(in) :: x_val
! f2py depend (n) :: x_val
real, dimension(n), intent(out) :: f_val
integer :: n
! f2py intent(hide) n
call func(x_val, f_val, n)
end subroutine
In this subroutine, we expect that f_val and x_val will both be of size n. The function will figure out for itself what n is based on the x_val that it is given when called from Python. It will then re-use that size with f_val. So, for example, if we save the above in our array_input_and_output.f90 file, we can call it from Python as follows:
import array_input_and_output
import numpy as np
def func(x):
return x
x = np.array([1,5])
output = array_input_and_output.arbitrary_size_array_callback(func, x)
print(output)
This will give the expected output of [1. 5.].