F2py, a Python function that returns an array (vector-valued function)

Question

In the following Python, I have five functions contained in the array returned by func , which I need to integrate. The code calls the Fortran plug-in generated with f2py :

 import numpy as np from numpy import cos, sin , exp from trapzdv import trapzdv def func(x): return np.array([x**2, x**3, cos(x), sin(x), exp(x)]) if __name__ == '__main__': xs = np.linspace(0.,20.,100) ans = trapzdv(func,xs,5) print 'from Fortran:', ans print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.), exp(20.)]) 

Fortran routine:

  subroutine trapzdv(f,xs,nf,nxs,result) integer :: I double precision :: x1,x2 integer, intent(in) :: nf, nxs double precision, dimension(nf) :: fx1,fx2 double precision, intent(in), dimension(nxs) :: xs double precision, intent(out), dimension(nf) :: result external :: f result = 0.0 do I = 2,nxs x1 = xs(I-1) x2 = xs(I) fx1 = f(x1) fx2 = f(x2) result = result + (fx1+fx2)*(x2-x1)/2 enddo return end 

The problem is that Fortran only integrates the first function in func(x) . See print result:

 from Fortran: [ 2666.80270721 2666.80270721 2666.80270721 2666.80270721 2666.80270721] exact: [ 2.66666667e+03 4.00000000e+04 9.12945251e-01 -4.08082062e-01 4.85165195e+08] 

One way that func(x) should change to return the value of a given position in an array of functions:

 def func(x,i): return np.array([x**2, x**3, cos(x), sin(x), exp(x)])[i-1] 

And then modify the Fortran procedure to call the function with two parameters:

  subroutine trapzdv(f,xs,nf,nxs,result) integer :: I double precision :: x1,x2,fx1,fx2 integer, intent(in) :: nf, nxs double precision, intent(in), dimension(nxs) :: xs double precision, intent(out), dimension(nf) :: result external :: f result = 0.0 do I = 2,nxs x1 = xs(I-1) x2 = xs(I) do J = 1,nf fx1 = f(x1,J) fx2 = f(x2,J) result(J) = result(J) + (fx1+fx2)*(x2-x1)/2 enddo enddo return end 

What works:

 from Fortran: [ 2.66680271e+03 4.00040812e+04 9.09838195e-01 5.89903440e-01 4.86814128e+08] exact: [ 2.66666667e+03 4.00000000e+04 9.12945251e-01 -4.08082062e-01 4.85165195e+08] 

But here func is called 5 times larger than necessary (in the real case, func has more than 300 functions, so it will be called 300 times bigger than necessary).

• Does anyone know a better solution to make Fortran recognize the entire array returned by func(x) ? In other words, make Fortran build fx1 = f(x1) as an array with 5 elements corresponding to functions in func(x) .

---

OBS: I compile with f2py -c --compiler=mingw32 -m trapzdv trapzdv.f90

Answer 1

Unfortunately, you cannot return an array from a python function in Fortran. To do this, you need a subroutine (which means that it is called using the call statement), and this is what f2py does not allow you to do.

In Fortran 90 you can create functions that return arrays, but again this is not what f2py can do, especially since your function is not a Fortran function.

Your only option is to use a workaround workaround or redesign of how you want to interact with python and Fortran.

Answer 2

Although this answer does not solve the question, this workaround does the same in Cython. Here, a trapezoidal rule and a polynomial integrator are implemented for a vector function. In the code below, I entered integratev.pyx :

 import numpy as np from numpy.linalg import inv cimport numpy as np FLOAT = np.float32 ctypedef np.float_t FLOAT_t def trapzv(f, np.ndarray xs, int nf): cdef int nxs = xs.shape[0] cdef np.ndarray ans = np.zeros(nf, dtype=FLOAT) cdef double x1, x2 for i in range(1,nxs): x1 = xs[i-1] x2 = xs[i] ans += (f(x2)+f(x1))*(x2-x1)/2. return ans def poly(f, np.ndarray xs, int nf, int order=2): cdef int nxs = xs.shape[0] cdef np.ndarray ans = np.zeros(nf, dtype=FLOAT) cdef np.ndarray xis = np.zeros(order+1, dtype=FLOAT) cdef np.ndarray ais if nxs % (order+1) != 0: raise ValueError("poly: The size of xs must be a multiple of 'order+1'") for i in range(order,nxs,order): xis = xs[i-order:i+1] X = np.concatenate([(xis**i)[:,None] for i in range(order+1)], axis=1) ais = np.dot( inv(X), f(xis).transpose() ) for k in range(1,order+2): ans += ais[k-1,:]/k * (xis[-1]**k - xis[0]**k) return ans 

The following test was used:

 import numpy as np from numpy import cos, sin , exp import pyximport; pyximport.install() import integratev from subprocess import Popen def func(x): return np.array([x**2, x**3, cos(x), sin(x), exp(x)]) if __name__ == '__main__': xs = np.linspace(0.,20.,33) print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.)+1, exp(20.)-1]) ans = integratev.trapzv(func,xs,5) print 'trapzv:', ans ans = integratev.poly(func,xs,5,2) print 'poly:', ans 

Donation:

 exact: [ 2.66666667e+03 4.00000000e+04 9.12945251e-01 5.91917938e-01 4.85165194e+08] trapzv: [ 2.66796875e+03 4.00390625e+04 8.83031547e-01 5.72522998e-01 5.00856448e+08] poly: [ 2.66666675e+03 4.00000000e+04 9.13748980e-01 5.92435718e-01 4.85562144e+08] 

Poly can be of any order, which is likely to give the best results for most cases ...