8

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?

petervanya
  • 306
  • 4
  • 11

2 Answers2

9

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.

DavidW
  • 29,336
  • 6
  • 55
  • 86
  • Is there any documentation for option 3? Did you test it? The `allocatable,dimension(n_b+1)` is clearly against Fortran rules but perhaps f2py allows it for its callback interfaces. However, `callable` has still implicit interface inside `create_hist`. – Vladimir F Героям слава Jan 10 '16 at 17:53
  • I have tested option 3 and it does work. It's a horrendous hack of my own invention so there's no documentation. I think you're right though - I should delete the "allocatable" in the specification (I'll edit the answer to do that). – DavidW Jan 10 '16 at 22:16
  • "However, callable has still implicit interface inside create_hist" - yes, that's almost certainly true. I don't really know Fortran well enough to fix that, but I'd gratefully accept a suggestion (or edit) that does. – DavidW Jan 10 '16 at 22:17
  • Right - to follow this up the proper way to specify the callback interface looks to be to use interface blocks, as described here http://www.fortran90.org/src/best-practices.html#callbacks (which I'm sure Vladimir F already knows, but might be useful to someone else). Unfortunately f2py does seem to be smart enough to read it, so you still have to create and modify the pyf file – DavidW Jan 10 '16 at 23:05
  • Yes, without the `allocatable` in the interface it should be better. – Vladimir F Героям слава Jan 11 '16 at 12:14
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.

  • Would a large enough assumed-shape argument work, also? And should there also be a "returned size" argument so the caller knows how much of the array is valid (depends on use case, I'm sure)? [Asking in ignorance of things more than anything.] – francescalus Jan 03 '16 at 19:32
  • Yes, assumed-shape would work too. You are also right about the returned size. Otherwise it would have to be inferred from the values of `bins` (they could be initialized to a negative number, or something. – Vladimir F Героям слава Jan 03 '16 at 21:24
  • @VladimirF do you know of some fortran library/subroutine that can build 2d histograms. something similar to hist2d from numpy i.e. with user defined bins. – Alexander Cska Aug 22 '20 at 20:07