天天看點

如何在Fortran中調用Python

Python是機器學習領域不斷增長的通用語言。擁有一些非常棒的工具包,比如scikit-learn,tensorflow和pytorch。氣候模式通常是使用Fortran實作的。那麼我們應該将基于Python的機器學習遷移到Fortran模型中嗎?資料科學領域可能會利用HTTP API(比如Flask)封裝機器學習方法,但是HTTP在緊密耦合的系統(比如氣候模式)中效率太低。是以,可以選擇直接從Fortran中調用Python,直接通過RAM傳遞氣候模式的狀态,而不是通過高延遲的通信層,比如HTTP。

有很多方法可以實作通過Python調用Fortran,但是從Fortran調用Python的方法卻很少。從Fortran調用Python,可以看作是将Python代碼嵌入到Fortran,但是Python的設計并不是像嵌入式語言Lua。可以通過以下三種方法實作從Fortran調用Python:

•Python的C語言API。這是最常用的方式,但需要實作大量的C封裝代碼。•基于Cython。Cython用于從Python中調用C語言,但也可以實作從C調用Python。•基于CFFI。CFFI提供了非常友善的方法可以嵌入Python代碼。

無論選擇哪種方法,使用者都需要将可執行Fortran檔案連結到系統Python庫,比如通過添加-lpython3.6到Fortran模式的Makefile檔案。下面通過Hello World示例示範如何通過Fortran調用Python。Fortran代碼儲存在test.f90檔案,如下:

! test.f90
program call_python
  use, intrinsic :: iso_c_binding
  implicit none
  interface
     subroutine hello_world() bind (c)
     end subroutine hello_world
  end interface

  call hello_world()

end program call_python
           

複制

首先導入Fortran 2003内部定義的和C語言類型互通的子產品iso_c_binding。但使用CFFI時,我們不需要寫任何C代碼,CFFI會生成C類型的打包接口。下一行則定義了一個C函數hello_world接口,這可以在C語言中實作,但是這裡我們使用Python和CFFI。最後,調用hello_world。

為了使用hello_world,我們需要建構CFFI标注,并儲存在builder.py中,此代碼用于建立可以連結Fortran程式的動态庫:

import cffi
ffibuilder = cffi.FFI()

header = """
extern void hello_world(void);
"""

module = """
from my_plugin import ffi
import numpy as np

@ffi.def_extern()
def hello_world():
    print("Hello World!")
"""

with open("plugin.h", "w") as f:
    f.write(header)

ffibuilder.embedding_api(header)
ffibuilder.set_source("my_plugin", r'''
    #include "plugin.h"
''')

ffibuilder.embedding_init_code(module)
ffibuilder.compile(target="libplugin.dylib", verbose=True)           

複制

‍首先,我們導入cffi包,并且聲明了外部函數接口(FFI)對象。這看起來似乎比較奇怪,這隻是CFFI實作這種目的的方式。下一步,header字元串中包含了需要調用的函數接口的定義。module字元串中包含了真正需要執行的Python程式。裝飾器@ffi.def_extern用于标記hello_world函數。my_plugin用于擷取ffi對象。def_extern裝飾器用于處理C類型,指針等。然後,ffibuilder.embedding_api(header)定義了API,embedding_init_code定義了Python代碼。

看起來比較奇怪的是在字元串中定義Python代碼,但CFFI需要以這種方式将Python代碼建構為共享庫對象。ffibuilder.set_source來設定源代碼資訊(?)。

然後執行以下語句建立共享庫libplugin.dylib:

python builder.py           

複制

然後使用下列指令編譯Fortran程式:

gfortran -o test -L./ -lplugin test.f90           

複制

以上是在Mac OSX上建立的共享庫,如果在Linux上,共享庫應該以.so結尾。如果一切沒有問題,那麼就可以執行檔案了:

./test
hello world           

複制

以上示範了如何使用CFFI從Fortran中調用Python程式,而不需要寫任何C程式。

FAQ

必須将所有Python代碼寫入header字元串嗎

不需要這樣。你可以直接在不同的Python子產品中定義Python代碼(比如my_module.py),然後在module字元串的開頭導入即可。比如,builder.py可以改為如下形式:

...

module = """
from my_plugin import ffi
import my_module

@ffi.def_extern()
def hello_world():
    my_module.some_function()

"""
...           

複制

這将在Python中使用可導入的形式使用Python程式。在添加到Fortran中之前,你也可以通過python -c "import my_module"測試一下。如果失敗了,你可能需要将包含my_module子產品的路徑添加到Python的sys.path變量中。

如何傳遞Fortran數組給Python

stack overflow page回答了此問題。下面是一個示例,将代碼定義在一個子產品檔案中,比如my_module.py:

# my_module.py

# Create the dictionary mapping ctypes to np dtypes.
ctype2dtype = {}

# Integer types
for prefix in ('int', 'uint'):
    for log_bytes in range(4):
        ctype = '%s%d_t' % (prefix, 8 * (2**log_bytes))
        dtype = '%s%d' % (prefix[0], 2**log_bytes)
        # print( ctype )
        # print( dtype )
        ctype2dtype[ctype] = np.dtype(dtype)

# Floating point types
ctype2dtype['float'] = np.dtype('f4')
ctype2dtype['double'] = np.dtype('f8')


def asarray(ffi, ptr, shape, **kwargs):
    length = np.prod(shape)
    # Get the canonical C type of the elements of ptr as a string.
    T = ffi.getctype(ffi.typeof(ptr).item)
    # print( T )
    # print( ffi.sizeof( T ) )

    if T not in ctype2dtype:
        raise RuntimeError("Cannot create an array for element type: %s" % T)

    a = np.frombuffer(ffi.buffer(ptr, length * ffi.sizeof(T)), ctype2dtype[T])\
          .reshape(shape, **kwargs)
    return a           

複制

asarray函數使用CFFI的ffi對象轉換指針ptr為給定形狀的numpy數組。可以使用如下形式在builder.py中的module字元串中調用:

module = """
import my_module

@ffi.def_extern()
def add_one(a_ptr)
    a = my_module.asarray(a)
    a[:] += 1
"""           

複制

add_one也可以定義在my_module.py中。最後,我們需要定義與函數相關的頭檔案資訊,并且添加到builder.py的header字元串中:

header  = """
extern void add_one (double *);
"""           

複制

最後,在Fortran中以如下形式調用:

program call_python
  use, intrinsic :: iso_c_binding
  implicit none
  interface
    subroutine add_one(x_c, n) bind (c)
        use iso_c_binding
        integer(c_int) :: n
        real(c_double) :: x_c(n)
    end subroutine add_one
  end interface
  
  real(c_double) :: x(10)

  print *, x
  call add_one(x, size(x))
  print *, x

end program call_python           

複制

這一部分,我們介紹了如何在Fortran中嵌入Python代碼塊,以及如何傳遞數組給Fortran或從Fortran傳遞數組給Python。然後,有些方面還是不太友善。

必須要在三個不同的區域定義python函數簽名嗎

任何要傳遞給Fortran的Python函數,都必須要要在三個區域進行定義。

•首先,必須在header.h中進行C頭檔案聲明•然後,執行函數必須要在builder.py的module字元串中,或一個外部子產品中•最後,Fortran代碼中必須包含定義子程式的interface塊(接口塊)

這對于改變Python函數來說就顯得有些麻煩。比如,我們寫了一個Python函數:

def compute_precipitation(T, Q):
    ...           

複制

必須要在所有區域進行聲明。如果我們想添加一個垂直渦度W作為輸入參數,我們必須要修改builder.py以及調用Fortran的程式。顯而易見,對于大的工程來說,這就變得極為麻煩。

對于一般通信而言,采用了一小部分fortran/python代碼封裝器。主要依賴于一個外部python子產品:

# module.py
import imp

STATE = {}


def get(key, c_ptr, shape, ffi):
    """Copy the numpy array stored in STATE[key] to a pointer"""

    # wrap pointer in numpy array
    fortran_arr = asarray(ffi, c_ptr, shape)
    
    # update the numpy array in place
    fortran_arr[:] = STATE[key]
    

def set(key, c_ptr, shape, ffi):
    """Call python"""
    STATE[key] = asarray(ffi, c_ptr, shape).copy()


def call_function(module_name, function_name):

    # import the python module
    import importlib
    
    mod = importlib.import_module(module_name)
    
    # the function we want to call
    fun = getattr(mod, function_name)
    
    # call the function
    # this function can edit STATE inplace
    fun(STATE)           

複制

全局變量STATE是一個包含了函數需要的所有資料的Python字典。get和set函數的功能主要就是将Fortran數組傳遞給STATA或者從STATE中取出Fortran數組。如果這些函數使用了Fortran/CFFI封裝器,那麼可以使用如下方式從Fortran中調用Python函數cumulus.compute_precipitation(state_dict):

call set("Q", q)
call set("T", temperature)
call set("ahother_arg", 1)

call call_function("cumulus", "compute_precipitation")

call get("prec", prec)           

複制

如果需要傳遞更多的資料給compute_precipitation,那麼需要添加更多的call set指令,這可能會改變Python函數的執行。我們就不需要改變builder.py中的任何代碼。

結論

上面描述了如何傳遞Fortran資料給Python函數,然後再擷取計算輸出。為了解決頻繁更改接口的問題,我們将fortran資料放到了Python子產品的字典中。通過調用給定的名稱來擷取資料,并且将計算結果也存儲到相同的字段中,然後,Fortran代碼通過索引字典中正确的關鍵詞來擷取結果。Cython中使用了類似的架構,但CFFI更為友善。最重要的是,從C語言中調用Cython需要導入Python.h頭檔案,還要運作Py_initialize和init_my_cython_module函數。然而,CFFI會在背景完成這些操作。

這篇文章隻是起到一個簡單的訓示性作用,有很多問題都沒有讨論,比如如何傳遞Fortran字元給Python。更多的代碼資訊,見Github。

感興趣的也可以看一下Forpy[2]這個包。

References

[1]

: https://www.noahbrenowitz.com/post/calling-fortran-from-python/

[2]: https://github.com/ylikx/forpy#using-forpy-with-anaconda