f2py，Python函數，返回一個數組（向量值函數）

Question

在下面的Python中，我有五個函數包含在由func返回的數組中，這是我必須集成的。該代碼調用使用f2py產生的外部的Fortran模塊：

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例程是：

 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 

的問題是，Fortran的僅在func(x)集成第一功能。 參見印刷結果：

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] 

一種方法workarond即修改func(x)返回給定 位置的值的函數的陣列中：

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

，然後改變Fortran例程調用函數有兩個參數：

 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 

其中一期工程：

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] 

但這裏func被稱爲超過必要的5倍（在現實情況下func 具有高於300層的功能，所以它會被稱爲多300倍必要）。

• 有沒有人知道一個更好的解決方案，使Fortran識別出func(x)返回的所有數組？換句話說，將Fortran構建爲fx1 = f(x1)作爲一個數組，其中有5個元素對應於func(x)中的函數。

---

OBS：我使用f2py -c --compiler=mingw32 -m trapzdv trapzdv.f90

Answer 1

不幸的是，你不能返回從Python功能到Fortran的數組進行編譯。你需要一個子程序（意思是用call聲明調用），這是f2py不允許你做的事情。

在Fortran 90中，您可以創建返回數組的函數，但這不是f2py可以執行的功能，特別是因爲您的函數不是Fortran函數。

您唯一的選擇是使用循環解決方法，或重新設計python和Fortran之間的交互方式。

Answer 2

儘管這個答案沒有解決問題，但它是一種在Cython中做同樣的解決方法。這裏梯形法則和一個多項式積分器是爲矢量值函數實現的。下面我的代碼放在一個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 

以下測試：

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 

，並提供：

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] 

聚可以是任何順序的，這可能會給對於大多數情況更好的結果...