Fortranで記述して、f2pyを使うのもよいのですが、Cで記述した関数が手元にあるもので、そちらを使いたいと思ったしだいです。
Python APIを使ったり、SWIGやCythonを利用したりと、いろいろ方法はある様子ですが、ctypesモジュールを利用して、呼び出す方法を試してみました。
コンパイラは、Visual Studioではなくて、MinGWのgccを使います。関数をコンパイルするだけなら、コンソールからちょっとコマンドを打つだけなので、こちらのほうが楽そうです。
簡単なものから順番に試していきたいと思います。
まずは、一つの整数を渡して二倍して返す関数です。
/* doublenumber.c */
#include
int doublenumber(int n)
{
return n * 2;
}
コンパイルします
gcc -Wall -fPIC -c doublenumber.c
windowsでは、-fPICは無視されるとの警告がでました。
どうやらデフォルトでPICオプションが適用される様子です。
オブジェクトファイルdoublenumber.oができました。
共有ライブラリを作成します。
gcc -shared -o doublenumber.so doublenumber.o
共有ライブラリdoublenumber.soファイルができました。
pythonから利用してみようと思います。
import ctypes
mydll = ctypes.CDLL('doublenumber.so')
mydll.doublenumber(5)
10
とちゃんと利用できています。
実数を引数として与えると、ctypes.ArgumentErrorが発生します。
さて、次に実数型を利用してみましょう。
C関数の引数をdouble vにして、同様に二倍した値を返すものにします。
/* doubledouble.c */
#include
double doublenumber(double v)
{
return v * 2.0;
}
コンパイルして、共有ライブラリにします。
gcc -Wall -c doubledouble.c
gcc -shared -o doubledouble.so doubledouble.o
実数を渡してもエラーが出ます。
型をc言語用に変換して、渡します。
mydll.doublenumber(ctypes.c_double(5.3))
1075131187
あれ、エラーは出ませんが、返り値がおかしいですね。
どうやらデフォルトで、関数はc_intを返すと仮定されているようです。
関数オブジェクトのrestype属性を設定するとよいようです。
mydll.doublenumber.restype = ctypes.c_double
mydll.doublenumber(ctypes.c_double(5.3))
10.6
ちゃんと値が得られました。
さて、引数でc_double指定をするのも大変です。
argtypesを設定してみます。
mydll.doublenumber.argtypes = [ctypes.c_double]
mydll.doublenumber(5.3)
10.6
簡単に5.3を渡すことができました。
整数5を渡しても、自動的に変換されるようです。
次は、文字列を渡す関数を試してみます。
/* charcount.c */
#include
int charcount(char *s)
{
int c;
c = 0;
while (s[c] != '\0')
c++;
return c;
}
そういえばstrlen()なんて関数もありましたね。
コンパイルして共有ライブラリを作成します。
gcc -Wall -c charcount.c
gcc -shared -o charcount.so charcount.o
テストするPythonスクリプトは以下の通り
C言語の文字列は、文字配列のポインタになります。
文字型のポインタは、ctypes.c_char_pで指定できます。
import ctypes
mydll = ctypes.CDLL("charcount.so")
mydll.charcount.argtypes = [ctypes.c_char_p]
mydll.charcount.restype = ctypes.c_int
print mydll.charcount("Hello")
実行すると文字列の長さの 5 がプリントされました。
次に、文字列を渡して、それを変更して返すような関数を作ってみます。
/* editstring.c */
#include
#include
char *editstring(char *s)
{
char buff[1024];
sprintf(buff, "*** %s ***", s);
strcpy(s, buff);
return s;
}
バッファ内で文字列を編集して、それを渡された文字列領域にコピーして、
そのポインタを返す関数です。
渡した文字列領域が編集されることになります。
Pythonの文字列は変更できないので、編集のための文字列バッファを用意するらしいです。
Pythonスクリプトは以下の通りです。
import ctypes
mydll = ctypes.CDLL("editstring.so")
mydll.editstring.argtypes = [ctypes.c_char_p]
mydll.editstring.restype = ctypes.c_char_p
s = ctypes.create_string_buffer("Hello", 1024)
result = mydll.editstring(s)
print '[%s]' % result
実行すると、
[*** Hello ***]
と表示され、編集されたのが確認できます。
返り値の型をc_char_pにして受け取ったresultは、
そのまま文字列としてプリントされますが、sはそうはいきません。
print sとすると c_char_arrayのオブジェクトであると表示されます。
その中の文字列をプリントするときは、
print '%s' % s.value
とするとよいようです。
でも、普通にbuffのポインタを返したらよいかな。
/* newstring.c */
#include
char *newstring(char *s)
{
char buff[1024];
sprintf(buff, "*** %s ***", s);
return buff;
}
こちら、コンパイルすると「関数が局所変数のアドレスを返します」
と警告が出ます。それはそうですよね。
この関数を抜けると、このメモリ領域は、解放されてしまいそうですが、、、。
Pythonスクリプトは以下の通り。
import ctypes
mydll = ctypes.CDLL("newstring.so")
mydll.newstring.argtypes = [ctypes.c_char_p]
mydll.newstring.restype = ctypes.c_char_p
string = mydll.newstring("Hello")
print string
ちゃんと*** Hello ***と出力されました。
このような書き方のほうが書きやすいですね。
直接渡した文字列をいじる必要がなければ、create_string_buffer()
は使わなくてもよいのでしょうか。
type(string)すると、
となります。
でもちょっとメモリが心配。
Python側で文字列バッファを用意したほうがいいかな。
次は、複数の実数値が入った配列をPythonから渡して、その和を返すようなC関数を試してみます。
/* sumvalues.c */
#include
double sumvalues(double *v, int n)
{
int i;
double sum;
sum = 0.0;
for (i=0; i {
sum = sum + v[i];
printf("%f\n", v[i]);
}
return sum;
}
配列の長さはCの関数から分かりませんので、その長さも渡すようにしています。
コンパイルして、共有ライブラリを作ります。
gcc -c sumvalues.c
gcc -shared -o sumvalues.so sumvalues.o
Pythonスクリプトは以下の通り
import ctypes
mydll = ctypes.CDLL("sumvalues.so")
mydll.sumvalues.restype = ctypes.c_double
n = 10
nvalues = ctypes.c_double * n
vlist = nvalues(*[1.0+i*1.0 for i in range(n)])
print vlist, len(vlist)
v = mydll.sumvalues(vlist, ctypes.c_int(n))
print v
共有ライブラリをを呼び出すところと、返り値を指定するのは同じです。
配列を使うときは、まず長さ分の型を作ります。
nvaluesという型で、長さがnのctypes.cdoubleの型を作ります。
その型でvlistというn個の実数値が配列に入ったものを作成します。
len()を用いてvlistの長さを調べることができます。
vlistは、c_double_Array_10オブジェクトになっています。
vlistは、ctypesの型になっていますので、あらためてargtypesで設定する必要もありません。長さだけ、ctypes.c_int()の型にして渡します。
1から10までの足し算をすることができました。
実際に利用するときは、Pythonの関数でラッピングするとよいかと思いました。
そのように書いたPythonスクリプトは以下の通りです。
import ctypes
mydll = ctypes.CDLL("sumvalues.so")
mydll.sumvalues.restype = ctypes.c_double
def sumvalues(vlist):
n = len(vlist)
nvalues = ctypes.c_double * n
arg1 = nvalues(*vlist)
v = mydll.sumvalues(arg1, ctypes.c_int(n))
return v
if __name__ == '__main__':
vlist = [1.0+i*1.0 for i in range(10)]
v = sumvalues(vlist)
print v
次は、複数の値を返してもらう関数を書いてみます。
関数は、基本的に一つの値しか返しません。
ですので、引数に値を入れてもらうための変数（へのポインタ）を渡すことにします。
ソースファイル
/* summean.c */
#include
int summean(double *v, int n, double *sum, double *mean)
{
int i;
double s;
s = 0.0;
for (i=0; i {
s = s + v[i];
}
*sum = s;
*mean = s/n;
return 0;
}
コンパイルして共有ライブラリを作成します。
gcc -c summean.c
gcc -shared -o summean.so summean.o
Pythonスクリプトは以下の通り
import ctypes
mydll = ctypes.CDLL("summean.so")
def summean(vs):
n = len(vs)
nvalues = ctypes.c_double * n
vlist = nvalues(*vs)
s = ctypes.c_double(0.0)
m = ctypes.c_double(0.0)
mydll.summean(vlist, ctypes.c_int(n), ctypes.pointer(s), ctypes.pointer(m))
return (s.value, m.value)
if __name__ == '__main__':
vs = [1.0+i*1.0 for i in range(10)]
s, m = summean(vs)
print s, m
前と同様にPythonスクリプトでsummean関数としてラッピングしました。
結果を保存する変数sおよびmをc_double(0.0)で用意します。
関数summeanには、そのポインタを渡す必要があるため、ctypes.pointer()関数を利用しています。
Pythonの関数の返り値として、s.valueおよびm.valueをタプルにして返します。
つぎは、いよいよ構造体です。まずは、シンプルな構造体を渡して、その内容を参照した計算結果が返ってくる関数を作ってみたいと思います。
ここでは、四角の左、右、下、上の座標を格納する構造体です。
ソースファイル
/* squarearea.c *\
#include
typedef struct {
double left;
double right;
double bottom;
double top;
} SQUARE;
double squarearea(SQUARE sq)
{
double area;
double width;
double height;
width = sq.right - sq.left;
height = sq.top - sq.bottom;
area = width * height;
return area;
}
コンパイルして、共有ライブラリを作成します。
gcc -c squarearea.c
gcc -shared -o squarearea.so squarearea.o
次に、Pythonスクリプトです。
Pythonスクリプト側で構造体を取り扱うためには、ctypes.Structureを継承したクラスを作ります。
そのクラスに、_fields_属性を設定することで、Cの構造体と一致させるということです。
_fields_属性は、フィールド名 と フィールド型 を持つ 2要素タプル のリストです。
ctypes.Structureで定義されているコンストラクタが、このフィールドを参照して、初期化するようです。
Pythonスクリプトは、以下のようになりました。
import ctypes
mydll = ctypes.CDLL("squarearea.so")
class SQUARE(ctypes.Structure):
_fields_ = [('left', ctypes.c_double),
('right', ctypes.c_double),
('bottom', ctypes.c_double),
('top', ctypes.c_double)]
if __name__ == '__main__':
sq = SQUARE(0.3, 0.8, 0.4, 3.4)
mydll.squarearea.restype = ctypes.c_double
v = mydll.squarearea(sq)
print v
幅0.5、高さ3ですので、1.5が返ってきます。
返り値の型を設定するのを忘れないようにしましょう。
C言語で返り値として構造体を使うと、複数の値を返すことができて便利ですよね。
先ほどの関数で、高さを1/2したものを返すようにしてみましょう。
ソースファイル
/* squarehalf.c */
#include
typedef struct {
double left;
double right;
double bottom;
double top;
} SQUARE;
SQUARE squarehalf(SQUARE sq)
{
double area;
double width;
double height;
SQUARE hsq;
hsq = sq;
height = sq.top - sq.bottom;
hsq.top = hsq.bottom + height/2.0;
return hsq;
}
hsqにsqを入れて、そのtopの値を、半分の高さになるように設定しなおしています。
コンパイルして、共有ライブラリを作ります。
gcc -c squarehalf.c
gcc -shared -o squarehalf.so squarehalf.o
Pythonスクリプト
import ctypes
mydll = ctypes.CDLL("squarehalf.so")
class SQUARE(ctypes.Structure):
_fields_ = [('left', ctypes.c_double),
('right', ctypes.c_double),
('bottom', ctypes.c_double),
('top', ctypes.c_double)]
if __name__ == '__main__':
sq = SQUARE(0.3, 0.8, 0.4, 3.4)
v = SQUARE()
mydll.squarehalf.restype = SQUARE
v = mydll.squarehalf(sq)
print sq.left, sq.right, sq.bottom, sq.top
print v.left, v.right, v.bottom, v.top
返り値をSQUAREにすることと、返り値を入れる変数vをSQUARE()クラスのインスタンスにしただけです。
0.3 0.8 0.4 3.4
0.3 0.8 0.4 1.9
と出力され、半分の高さの四角を表すデータになりました。
構造体の中に、別の構造体がある場合はどう書けるのでしょうか。
先ほど面積を求めた四角を、左下、右上の座標データPOINTで表してみます。
/* squarearea2.c */
#include
typedef struct {
double x;
double y;
} POINT;
typedef struct {
POINT leftbottom;
POINT righttop;
} SQUARE;
double squarearea(SQUARE sq)
{
double area;
double width;
double height;
width = sq.righttop.x - sq.leftbottom.x;
height = sq.righttop.y - sq.leftbottom.y;
area = width * height;
return area;
}
コンパイルして共有ライブラリにします。
gcc -c squarearea2.c
gcc -shared -o squarearea2.so squarearea2.o
Pythonスクリプトは以下のようになりました。
import ctypes
mydll = ctypes.CDLL("squarearea2.so")
class POINT(ctypes.Structure):
_fields_ = [('x', ctypes.c_double),
('y', ctypes.c_double)]
class SQUARE(ctypes.Structure):
_fields_ = [('leftbottom', POINT),
('righttop', POINT)]
if __name__ == '__main__':
sq = SQUARE(POINT(0.3, 0.4), POINT(0.8, 3.4))
v = SQUARE()
mydll.squarearea.restype = ctypes.c_double
v = mydll.squarearea(sq)
print v
それぞれ対応するクラスを書いただけで、意外なところはないですね。
ただ、クラス定義する順番で、POINTをSQUAREの後で行うと、
SQUAREで利用しているPOINTが定義できないというエラーが起こります。
構造体の定義では、その中に自分自身をもつものもありますね。
再帰的な構造は、中に同じ自分自身の型へのポインタを入れたりします。
不完全型と呼ばれているようです。
簡単のため、
(v1, (v2, (v3, (v4, (v5, NULL)))))
のようなリストに似た形のデータを考えてみます。
ソースファイル
/* showlist.c */
#include
struct CELL {
int value;
struct CELL *next;
};
int showlist(struct CELL *c, int flg)
{
printf("(%d, ", c->value);
if ((struct CELL *)NULL != c->next)
showlist(c->next, 0);
else
printf("NULL");
if (1 == flg)
printf(")\n");
else
printf(")");
return 0;
}
CELLという型の定義の中に、CELLへのポインタが含まれています。
前の例のように、typedef struct { ... } CELL; のような書き方をすると、コンパイル時にエラーがでます。使われるより前の方に、CELLというワードが入るような記述をする必要があるということです。
showlistは、再帰的に呼び出される関数で、渡されたCELLの値を出力したのちに、続くCELLの描画を行うため、showlistを呼び出します。nextに、NULLポインタが入っていれば、そこで終了です。
flgは、最初に呼び出すときに1をセットすると、それ以降とを区別することができます。
コンパイルして、共有ライブラリにします。
gcc -c showlist.c
gcc -shared -o showlist.so sholist.o
このような不完全な構造体定義の場合、Pythonスクリプトは、次のようにするようです。
import ctypes
mydll = ctypes.CDLL("showlist.so")
class CELL(ctypes.Structure):
pass
CELL._fields_ = [('value', ctypes.c_int),
('next', ctypes.POINTER(CELL))]
if __name__ == '__main__':
c = CELL(0, None)
for i in range(5):
c = CELL((i+1), ctypes.pointer(c))
mydll.showlist(ctypes.pointer(c), ctypes.c_int(1))
クラスCELLの定義では、なにも行わずに、あとから改めて_fields_属性を設定するようです。
ctypes.POINTER(CELL)で、CELLへのポインタ型という指定ができます。
実際に利用するときに、CELL型のcというもののポインタは、ctypes.pointer(c)で得ます。紛らわしいですね。
Pythonスクリプトでcの値を設定していき、showlistにcへのポインタと整数1を渡すと、
(5, (4, (3, (2, (1, (0, NULL))))))
と表示されます。この表示は、C関数内で行われたものです。
このようにctypesを使って共有ライブラリを呼び出す方法だと、CのソースファイルにPythonで利用するという配慮をせずに、Pythonから呼び出すことができました。
Pythonスクリプト側で、ちょっと使いやすいように包んであげれば呼び出しやすくなるかと思います。
f2pyを使ってFortranプログラムで計算するのもいいのですが、すでにCの関数を作っているようなら、この方法は便利ですね。
ということで、自作のガウス分布をフィットするルーチンをPythonから呼び出してみました。速度の比較のために、Pythonスクリプトでも書いてみました。numpyを利用せずに、Cのソースをほぼ1対1に直したものです。
#--- C func ---
#1.02713 0.0241351
#4.93666 0.0580069
#-0.938185 0.023901
#2.49769 0.0360503
# time_c = 0.000507835 s
#--- Python func ---
#1.02713 0.0241351
#4.93666 0.0580069
#-0.938185 0.023901
#2.49769 0.0360503
# time_py = 0.00944325 s
#-------------------
# time_py/time_c = 18.5951
#-------------------
200点のデータをフィットしました。やっぱりCで記述されたものは早いですね。
Pythonスクリプトの方は、numpyの機能を利用して実装すると、どこまで速くなるのでしょうか。
