f2py error with allocatable arrays

Question

I have a Fortran subroutine that I would like to use in Python.

subroutine create_hist(a, n, dr, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(out), allocatable :: hist(:)
    real(8), intent(out), allocatable :: bins(:)

    n_b = n_bins(a, n, dr)  ! a function calculating the number of bins
    allocate(bins(n_b+1))
    allocate(hist(n_b))
    hist = 0

    !... Create the histogram hist by putting elems of a into bins
end subroutine

It is a simple programme to get an array of numbers a and create a histogram based on given bin size dr. First, it gets the number of bins using function n_bins and then accordingly allocates the space for the arrays bins and hist.

While gfortran compiles this code fine, f2py raises an error:

/tmp/tmpY5badz/src.linux-x86_64-2.6/mat_ops-f2pywrappers2.f90:32.37:
      call create_hist_gnu(a, n, dr, bins, hist)
Error: Actual argument for 'bins' must be ALLOCATABLE at (1)

As I understand it, f2py does not tolerate allocating space for arrays at runtime (no idea why, as this seems to be a natural scientific need).

Could anyone please suggest a way to allocate Fortran arrays at run-time so that f2py is happy with that?

Answer 1

There's a few valid options/work-rounds that you can use.

1. Probably the simplest would be to arrange for python to call the function n_bins, and then use the result from that to call (a slightly modified version of) the create_hist function with non-allocatable output arrays of the correct sizes.

i.e. in Python:

n_b = n_bins(a,dr) # don't need to pass n - inferred from a
bins,hist = create_hist(a,dr,n_b) # bins and hist are auto-allocated

with the Fortran interface to create_hist now defined as

subroutine create_hist(a, n, dr, n_b, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(in) :: n_b
    integer, intent(out) :: hist(n_b)
    real(8), intent(out) :: bins(n_b+1)

    ! code follows

! you also need to specify the function n_bins...

---

That only works in cases where you can call n_bins cheaply and from outside create_hist. I know of 2 methods for simulating allocatable arrays for cases where that doesn't apply (i.e. the code to work out the array sizes is expensive and can't easily be separated).

2. The first is to use module level allocatable arrays (described in the documentation here). This is essentially an allocatable global variable - you call your function, it saves the data into the global variable and then you access that from Python. The disadvantage is that it isn't thread-safe (i.e. it's bad if you call create_hist simultaneously in parallel)

module something
    real(8), allocatable :: bins(:)
    integer, allocatable :: hist(:)
contains
    subroutine create_hist(a,n,dr)
        integer, intent(in) :: n
        real(8), intent(in) :: a(n)
        real(8), intent(in) :: dr
        integer :: n_b

        n_b = n_bins(a,n,dr)

        allocate(bins(n_b+1))
        allocate(hist(n_b))
        ! code follows
    end subroutine
end module

The Python call then looks like

something.create_hist(a,n,dr)
bins = something.bins # or possible bins = copy.copy(something.bins)
hist = something.hist # or possible hist = copy.copy(something.hist)

3. The other way which I quite like is to only allocate the arrays within the function (i.e. do not pass them in/out as parameters). However, what you do pass is a Python callback function that is called at the end and saves the arrays. This is thread-safe (I believe).

The fortran code then looks something like

subroutine create_hist(a,n,dr,callback)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    external callable ! note: better to specify the type with an interface block (see http://www.fortran90.org/src/best-practices.html#callbacks)
    integer :: n_b
    real(8), allocatable :: bins(:)
    integer, allocatable :: hist(:)

    n_b = n_bins(a,n,dr)

    allocate(bins(n_b+1))
    allocate(hist(n_b))
    ! code follows
    call callable(bins,hist,n_b)
end subroutine

unfortunately it's then slightly more involved. You need to create a signature file with the command f2py -m fortran_module -h fortran_module.pyf my_fortran_file.f90 (this is to build a Python module called fortran_module - change the names as appropriate), and then modify the relevant lines of it to clarify the dimensions of the callback function:

python module create_hist__user__routines 
    interface create_hist_user_interface 
        subroutine callable(bins,hist,n_b) ! in :f:my_fortran_file.f90:stuff:create_hist:unknown_interface
            real(kind=8), dimension(n_b+1) :: bins
            integer, dimension(n_b) :: hist
            integer :: n_b
        end subroutine callable
    end interface create_hist_user_interface
end python module create_hist__user__routines

compile that with f2py -c fortran_module.pyf my_fortran_file.f90

and then a short Python wrapper (to give you an easy interface) that looks like

def create_hist(a,dr):
    class SaveArraysCallable(object):
        def __call__(self,bins,hist):
            self.bins = bins.copy()
            self.hist = hist.copy()

    f = SaveArrayCallable()
    fortran_module.create_hist(a,dr,f)
    return f.bins, f.hist

---

Summary

For many cases option 1 is probably the best. Otherwise option 2 or option 3. I prefer option 3 because there aren't multithreading pitfalls to avoid (but really this is a corner-case that you'll likely never see). Option 2 is easier to implement.

Answer 2

As far as I know f2py does not support dummy arguments (sometimes known as function parameters) which have the attribute ALLOCATABLE. An alternative would be to make it a large enough explicit shape array. You would then use just a certain part of the array.

subroutine create_hist(a, n, dr, n_b_max, bins, hist)
    integer, intent(in) :: n
    real(8), intent(in) :: a(n)
    real(8), intent(in) :: dr
    integer, intent(in) :: n_b_max
    integer, intent(out) :: hist(n_b_max)
    real(8), intent(out) :: bins(n_b_max+1)

By the way, using real(8) is not a portable way how to specify 64-bit precision and will fail with certain compilers. Even good old double precision is better as it will always compile to something and often to 64 bit.