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