CFFIで外部関数を呼び出す

CFFIで外部関数を呼び出す#

Pythonの標準ライブラリには、外部関数を呼び出すためのctypesモジュールが提供されています。しかし、その使用方法が煩雑であるため、本書ではCFFI拡張ライブラリを使用して外部関数を呼び出す方法を紹介します。CFFIはC言語外部関数インターフェース（C Foreign Function Interface）の略で、C言語と同じ構文を使用してPythonとCの間の呼び出しインターフェースを宣言します。CFFIはダイナミックリンクライブラリ(DLL)を直接ロードしてその中の関数を呼び出すことができ、またC言語コンパイラを使用してインターフェースを拡張モジュールとしてコンパイルすることもできます。

import numpy as np
import sys
from cffi import FFI

2つの呼び出しモード#

CFFIは、外部関数を呼び出すための2つのモードを提供しています：ABIモードとAPIモードです。ABIモードでは、CFFIはctypesと同じ機能を実現できます。まず、外部関数の呼び出しインターフェースを宣言し、次にDLLをロードし、最後にC言語のデータを作成して外部関数を呼び出します。APIモードでは、CFFIはC言語コンパイラを使用してC言語コードとその呼び出しインターフェースを一緒にコンパイルし、Pythonの拡張モジュールとして生成します。

ABIモード#

ABIはアプリケーションバイナリインターフェース（Application Binary Interface）の略で、ユーザープログラムがソースコードなしでDLL内の関数を使用できることを意味します。以下の例では、normal.cにsnd_pdfとgnd_pdf()が定義されており、それぞれ標準正規分布と一般正規分布の確率密度関数を計算します。

%%writefile normal.c
#include <math.h>

double snd_pdf(double x)
{
    return exp(-0.5 * x*x) / sqrt(2*M_PI);
}

double gnd_pdf(double x, double mu, double sigma)
{
    return snd_pdf((x - mu) / sigma) / sigma;
}

Overwriting normal.c

ヘッダーファイルnormal.hには、これらの関数のシグネチャが含まれています：

%%writefile normal.h
double snd_pdf(double x);
double gnd_pdf(double x, double mu, double sigma);

Overwriting normal.h

gccコンパイラの-sharedオプションを使用して、ソースプログラムをDLLとしてコンパイルできます。これにより、ソースプログラム内のすべての非静的関数がエクスポートされ、ユーザーは生成されたDLLを使用してこれらのエクスポートされた関数を呼び出すことができます。

!gcc -Ofast -shared -o normal.dll normal.c

CFFIを使用してDLL内の外部関数を呼び出す基本的な手順は以下の通りです：

FFIオブジェクトffiを作成します。
ffi.cdef()を呼び出して、各種の型と関数を宣言します。
ffi.dlopen()を呼び出してDLLをロードし、DlLibraryオブジェクトlibを取得します。
libの各属性を通じてDLL内の関数を呼び出します。

以下では、CFFIを使用してDLLをロードし、その中のgnd_pdf()関数を呼び出します。❶まず、FFIのインスタンスffiを作成します。❷ffi.cdef()を使用してDLL内の関数呼び出しインターフェースを定義します。関数の定義方法はC言語の関数宣言と同じです。❸ffi.dlopen()を呼び出してDLLをロードします。これにより、FFILibraryのインスタンスlibが返されます。❹lib.gnd_pdf()を通じてDLL内のgnd_pdf()関数を呼び出します。CFFIはPythonの数値オブジェクトを関数パラメータに対応する数値型に変換し、関数の戻り値をPythonのオブジェクトに変換します。❺最後に、ffi.dlclose()を呼び出してロードされたDLLを解放します。解放された後でないと、再度コンパイルコマンドを実行して既存のDLLファイルを上書きすることはできません。

from cffi import FFI

defines = """
double snd_pdf(double x);
double gnd_pdf(double x, double mu, double sigma);
"""

ffi = FFI()  # ❶
ffi.cdef(defines)  # ❷
lib = ffi.dlopen("normal.dll")  # ❸
print(lib.gnd_pdf(1, 4.5, 10))  # ❹
ffi.dlclose(lib)  # ❺

0.03752403469169379

APIモード#

CFFIのAPIモードは、C言語コンパイラを使用してC言語ファイルと関数呼び出しインターフェースをPythonの拡張モジュールとしてコンパイルします。関数呼び出しインターフェースを定義した後、ffi.set_source()を使用して拡張モジュールのコンパイルに必要な情報を設定し、ffi.compile()を呼び出してコンパイルを実行します。具体的に実行されるコンパイルコマンドを確認したい場合は、verboseパラメータをTrueに設定できます。CFFIはdistutilsモジュールを使用して拡張ライブラリをコンパイルします。C言語コンパイラの設定方法については、前章の関連内容を参照してください。

set_source()の最初のパラメータは出力される拡張モジュール名で、2番目のパラメータは拡張モジュールのソースプログラムに埋め込まれるC言語コードです。extra_objectsパラメータは、一緒にコンパイルする必要がある他のファイルです。extra_objects以外にも、以下のようなよく使用されるパラメータがあります。これらはすべて文字列のリストです：

include_dirs: ヘッダーファイルの検索パス
library_dirs: 静的リンクライブラリの検索パス
extra_compile_args: 追加のコンパイル引数
extra_link_args: 追加のリンク引数

完全なパラメータの説明は、%PYTHON%\Lib\distutils\extension.pyにあります。

ffi = FFI()
ffi.cdef(defines)
ffi.set_source("normal_math", '#include "normal.h"', extra_objects=["normal.c"])
ffi.compile();

CFFIでコンパイルされた拡張モジュールには、ffiとlibの2つのオブジェクトのみが含まれます。これらは、前述のDLLをロードした際に作成されたffiとlibと同じように使用できます。

from normal_math import ffi, lib

lib.gnd_pdf(1, 4.5, 10)

0.03752403469169379

また、埋め込まれたC言語コード内で直接C言語関数を定義することもできます。これにより、extra_objectsを使用して追加のソースファイルを指定する必要がなくなります。例えば、以下のプログラムもnormal_math拡張モジュールを正しくコンパイルします。

ffi.set_source("normal_math", """
#include <math.h>

double snd_pdf(double x)
{
    return exp(-0.5 * x*x) / sqrt(2*M_PI);
}
...
""")

コンパイラを呼び出して拡張ライブラリをコンパイルする前に、CFFIは自動的に拡張ライブラリのソースファイルnormal_math.cを生成します。このファイルには、snd_pdf()の呼び出しインターフェース関数が含まれています。その内容は以下の通りです：

static PyObject *
_cffi_f_snd_pdf(PyObject *self, PyObject *arg0)
{
  double x0;
  double result;
  PyObject *pyresult;

  x0 = (double)_cffi_to_c_double(arg0);
  if (x0 == (double)-1 && PyErr_Occurred())
    return NULL;

  Py_BEGIN_ALLOW_THREADS
  _cffi_restore_errno();
  { result = snd_pdf(x0); }
  _cffi_save_errno();
  Py_END_ALLOW_THREADS

  (void)self; /* unused */
  pyresult = _cffi_from_c_double(result);
  return pyresult;
}

拡張モジュール内のlibを使用してsnd_pdf()関数を呼び出すと、実際にはインターフェース関数_cffi_f_snd_pdf()が呼び出されます。そのパラメータarg0はPythonから渡された数値オブジェクトで、_cffi_to_c_double()を使用してdouble型の数値x0に変換された後、snd_pdf(x0)を呼び出して計算結果resultを取得します。その型はdoubleです。最後に、_cffi_from_c_double()を呼び出してdouble型の数値をPythonの数値オブジェクトに変換します。

CFFIは、cdef()で定義された各関数呼び出しインターフェースに対してインターフェース関数を作成します。これらの関数は、PythonオブジェクトをC言語のデータ型に変換し、対応するC言語関数を呼び出し、関数の戻り値をPythonオブジェクトに変換します。実際の計算を行う関数の実行時間が短い場合、インターフェース関数の処理時間が占める割合が増加します。この場合、計算をベクトル化することで実行速度を向上させることができます。

cffiマジックコマンド#

Notebook環境でcffiを使用してC言語関数を呼び出すのを容易にするために、本書では%%cffiマジックコマンドを提供しています。まず、%load_ext helper.cffiを使用してこのマジックコマンドをロードし、次に%%cffiを使用して拡張モジュールとしてコンパイルする必要があるC言語プログラムを記述します。これにより、cdef()に必要な呼び出しインターフェースが自動的に生成され、重複しない拡張モジュールとしてコンパイルされます。ffiとlibのオブジェクトを指定されたグローバル変数にラップして保存します。デフォルトのグローバル変数名はcですが、%%cffiの-nパラメータを使用してグローバル変数名を指定できます。例えば、%cffi -n tmpとします。

%load_ext helper.cffi

%%cffi
#include <math.h>

double snd_pdf(double x)
{
    return exp(-0.5 * x*x) / sqrt(2*M_PI);
}

double gnd_pdf(double x, double mu, double sigma)
{
    return snd_pdf((x - mu) / sigma) / sigma;
}

上記のセルのプログラムを実行すると、グローバル変数cはコンパイル結果のffiとlibをラップしたFlattenAttrオブジェクトになります。そのffiとlib属性を使用して、以前と同じようにC言語で定義された関数を呼び出すことができます：

c.lib.snd_pdf(1.0)

0.24197072451914337

属性名がffiとlibでない場合、FlattenAttrオブジェクトは自動的にffiとlibから対応する属性名を探します。したがって、以下のように外部関数を呼び出すこともできます：

c.gnd_pdf(1.0, 2.0, 0.3)

0.005140929987637018

C言語のデータ型#

外部関数を呼び出す際、CFFIは自動的にPythonオブジェクトを対応するC言語のデータ型に変換します。例えば、floatオブジェクトはdouble型に変換されます。より複雑なデータ型を使用する場合は、ffi.new()を使用してC言語データをラップするオブジェクトを作成する必要があります。ffi.new()はCDataオブジェクトを返し、作成されたC言語データを解放する責任があります。ガベージコレクション時に、free()を呼び出して管理されているC言語データを解放します。

配列とポインタ#

C言語の配列はffi.new("type[]", init)を使用して作成します。ここで、typeは配列の要素型で、initは配列の長さまたは初期値を表すシーケンスです。配列要素へのポインタはffi.addressof(arr, offset)を使用して作成します。配列とポインタはどちらも添字演算をサポートし、ステップ1のスライス添字をサポートしますが、負の添字はサポートしません。配列をPythonのリストに変換するにはlist()を使用し、ポインタをリストに変換するにはffi.unpack()を使用します。2番目のパラメータは要素の数です。

buf1 = ffi.new("double[]", 5)
buf1[0:5] = range(5)
buf2 = ffi.new("double[]", list(range(5)))

pbuf1 = ffi.addressof(buf1, 2)
pbuf1[0] = 3.14
print(list(buf1))
print(ffi.unpack(pbuf1, 3))

[0.0, 1.0, 3.14, 3.0, 4.0]
[3.14, 3.0, 4.0]

C言語の配列とメモリを共有するNumPy配列を作成できます。まず、ffi.buffer()を使用してバッファインターフェースをサポートするオブジェクトを返し、次にnp.frombuffer()を使用してバッファインターフェースオブジェクトをNumPy配列に変換します。Pythonのバッファインターフェースはデータのメモリアドレスと長さのみを定義し、データ型は定義しないため、dtypeパラメータを使用して配列要素の型を指定する必要があります。作成されたNumPy配列はC言語の配列とメモリを共有するため、arr_buf1[0]とbuf1[0]はメモリ内で同じアドレスにあります。

arr_buf1 = np.frombuffer(ffi.buffer(buf1), dtype=np.double)
arr_buf1[0] = 100
print(list(buf1))

[100.0, 1.0, 3.14, 3.0, 4.0]

ffi.from_buffer()を使用して、NumPy配列の先頭アドレスを指すポインタオブジェクトを取得できます。ポインタの型は最初のパラメータによって決定されます。型がdouble []の場合、返されるオブジェクトは配列の長さを保持し、list()を使用してリストに変換できます。型がdouble *の場合、ffi.unpack()を使用して長さを指定する必要があります。

arr = np.array([1, 2, 3], dtype=np.int16)
parr1 = ffi.from_buffer("short []", arr)
parr2 = ffi.from_buffer("short *", arr)
print(list(parr1))
print(ffi.unpack(parr2, 3))

[1, 2, 3]
[1, 2, 3]

型を指定しない場合、返されるオブジェクトの型はchar[]です。以下は、arrのデータストレージ領域内のすべてのバイトを表示します。

print(list(ffi.from_buffer(arr)))

[b'\x01', b'\x00', b'\x02', b'\x00', b'\x03', b'\x00']

以下の例では、CFFIを使用してMKLライブラリの配列演算関数を呼び出します。MKLはIntel Math Kernel Libraryの略称で、高度に最適化され、マルチスレッド処理された数学ルーチンを提供し、性能が非常に重要な科学、工学、金融などの分野のアプリケーションに向いています。vmdCos()には4つのパラメータがあり、それぞれ配列の長さ、入力配列、出力配列、計算モードです。以下ではmkl_rt.dllをロードし、vmdCos()の関数呼び出しインターフェースを宣言します。

文字列とバイト列#

C言語には文字列型がなく、通常はchar *型のポインタを使用して文字列のデータを指し、文字列はバイト0x00で終了します。char *はPythonのバイト列型bytesに対応します。以下ではC言語標準ライブラリのstrlen()関数を呼び出します。そのパラメータはPythonのbytesオブジェクトです。

ffi = FFI()
ffi.cdef("size_t strlen(const char *);")
lib = ffi.dlopen("msvcrt.dll")
lib.strlen(b"123")

wchar_t *型のポインタはPythonの文字列型strに対応します：

ffi.cdef("size_t wcslen(const wchar_t *);")
lib.wcslen("abc")

Pythonのバイト列と文字列オブジェクトはどちらも不変であるため、それらを内容を変更する関数に渡すことはできません。例えばC言語の文字列連結関数strcat()とwcscat()は、その最初のパラメータが指すメモリを変更します。この場合、ffi.new()を使用してchar []型またはwchar_t []型の配列オブジェクトを作成できます。Pythonで配列の内容を取得する必要がある場合は、ffi.string()を使用できます。これは配列の要素型に応じてバイト列オブジェクトまたは文字列オブジェクトを返します。

ffi.cdef(
    """
char * strcat(char *, char *);
wchar_t * wcscat(wchar_t *, wchar_t *);
"""
)
buf_char = ffi.new("char []", 1024)
lib.strcat(buf_char, b"123")
lib.strcat(buf_char, b"abcd")

buf_wchar = ffi.new("wchar_t []", 1024)
lib.wcscat(buf_wchar, "123")
lib.wcscat(buf_wchar, "abcd")

print(ffi.string(buf_char), ffi.string(buf_wchar))

b'123abcd' 123abcd

構造体#

C言語では、構造体を使用して関連するデータのグループを保存し、オブジェクト指向プログラミングのカプセル化機能に似た機能を実現することがよくあります。以下のiir.hとiir.cは、2次IIRフィルタを実装し、構造体IIR2を使用してフィルタの係数と状態を保存します。プログラムでは単精度浮動小数点数floatを使用して演算を行います。iir2_init()はフィルタの状態をリセットし、iir2_step()は単一のサンプル値をフィルタリングし、フィルタリング後の結果を返します。iir2_run()は配列内の複数のサンプル値をフィルタリングし、結果を出力配列に書き込みます。

%%writefile iir.h
typedef struct _IIR2
{
    float b0, b1, b2;
    float a1, a2;
    float z0, z1;
} IIR2;

void iir2_init(IIR2 * self);
float iir2_step(IIR2 * self, float x);
void iir2_run(IIR2 * self, float * x, float * y, size_t n);

Overwriting iir.h

%%writefile iir.c
#include <stddef.h>
#include "iir.h"

void iir2_init(IIR2 * self)
{
    self->z0 = 0;
    self->z1 = 0;
}

float iir2_step(IIR2 * self, float x)
{
    float y;
    y = self->b0 * x + self->z0;
    self->z0 = self->b1 * x + self->z1 - self->a1 * y;
    self->z1 = self->b2 * x - self->a2 * y;
    return y;
}

void iir2_run(IIR2 * self, float * x, float * y, size_t n)
{
    size_t i;
    for(i=0;i<n;i++){
        y[i] = iir2_step(self, x[i]);
    }
}

Overwriting iir.c

以下ではAPIモードを使用してコードを拡張モジュールiir_filterにコンパイルします。ここではiir.hファイルの内容を直接使用して関数呼び出しインターフェースを定義します。ffi.cdef()はほとんどのC言語構文をサポートしているため、ヘッダーファイルを使用して関数呼び出しインターフェースを定義することで、コードの重複を避け、エラーの可能性を減らすことができます。ただし、C言語のマクロのサポートは限られているため、ヘッダーファイルに#include、#ifdefなどのマクロコマンドが含まれている場合、そのヘッダーファイルを直接使用してインターフェースを定義することはできません。

ffi = FFI()
with open("iir.h") as f:
    ffi.cdef(f.read())
ffi.set_source(
    "iir_filter",
    '#include "iir.h"',
    extra_objects=["iir.c"],
    extra_compile_args=["-Ofast"],
)
ffi.compile();

使用しやすくするために、以下ではIIR2クラスを使用してIIR2構造体および関連する関数をラップします。❶初期化メソッド__init__()では、IIR2 *型を使用してCDataオブジェクトを作成し、辞書を使用して構造体の各フィールドを初期化します。また、リストを使用してフィールドの順序に従って構造体を初期化することもできます。

iir2_run()がパラメータxとyの要素にアクセスする際のインデックスは連続的に変化するため、これらのパラメータに対応するNumPy配列がC言語連続メモリであることを保証する必要があります。❷numpy.ascontiguousarray()を呼び出すことで、C言語連続メモリを保証し、dtypeパラメータを使用して要素型を指定できます。run()メソッドを呼び出す際のパラメータxがC言語連続メモリ型でない場合、または要素型が単精度浮動小数点数でない場合、すべてのデータがコピーされます。

from iir_filter import ffi, lib


class IIR2:
    def __init__(self, b, a):
        self.iir = ffi.new(
            "IIR2 *", dict(b0=b[0], b1=b[1], b2=b[2], a1=a[1], a2=a[2])
        )  # ❶
        lib.iir2_init(self.iir)

    def run(self, x):
        x = np.ascontiguousarray(x, dtype=np.float32)  # ❷
        y = np.zeros_like(x)
        lib.iir2_run(
            self.iir,
            ffi.from_buffer("float []", x),
            ffi.from_buffer("float []", y),
            len(x),
        )
        return y

以下ではscipy.signal.butter()を使用して2次のバターワースローパスフィルタを設計し、フィルタの係数を使用してIIR2オブジェクトを作成します。signal.lfilter()の出力との最大誤差を比較します。結果から、単精度浮動小数点数の演算が十分に正確であることがわかります。フィルタの入力はnumpy.random.normal()によって作成された正規分布の乱数で、その要素型は倍精度浮動小数点数です。ascontiguousarray()によって単精度浮動小数点数配列に変換された後、lib.iir2_run()に渡されてフィルタリング演算が行われます。

from scipy import signal

b, a = signal.butter(2, 0.1)
iir = IIR2(b, a)
x = np.random.normal(size=200)
y = iir.run(x)
y2 = signal.lfilter(b, a, x)
np.abs(y - y2).max()

np.float64(4.368061093384945e-07)

Pythonでは、構造体のフィールドは属性アクセスを使用してアクセスできます。例えば：

st_iir = iir.iir
print(st_iir.z0, st_iir.z1)

-0.32050031423568726 0.14441882073879242

構造体を辞書に変換する必要がある場合、dir()を使用して属性リストを取得し、getattr()を使用して各属性の値を取得できます。以下では辞書内包表記を使用してst_iirを辞書に変換します。

{name: getattr(st_iir, name) for name in dir(st_iir)}

{'a1': -1.5610181093215942,
 'a2': 0.6413515210151672,
 'b0': 0.020083365961909294,
 'b1': 0.04016673192381859,
 'b2': 0.020083365961909294,
 'z0': -0.32050031423568726,
 'z1': 0.14441882073879242}

CFFIのメモリ管理#

Pythonのガベージコレクション機構により、Pythonプログラムを書く際にメモリの割り当てと解放について注意する必要はありません。しかし、C言語で書かれた外部関数を呼び出す際には、メモリの割り当てと解放に注意する必要があります。CFFIを使用する場合、そのメモリ管理メカニズムを理解することで、メモリリークや不正なポインタアクセスなどのエラーを回避できます。

Pythonでは、すべてのC言語データはCDataオブジェクトを使用して管理されます。このオブジェクトは内部にC言語データのメモリを指すポインタを持っています。これを使用してメモリを管理する場合、以下の3つの状況に分けられます：

ffi.new()を呼び出してC言語データのメモリを割り当て、そのメモリを管理するCDataOwnオブジェクトを返します。これがガベージコレクションされると、C言語データも同時に解放されます。したがって、外部関数がまだC言語データを使用している間は、このCDataOwnオブジェクトがガベージコレクションされないようにする必要があります。
外部関数を呼び出してメモリを割り当て、Pythonに返す場合、割り当てられたメモリのアドレスを保持するCDataオブジェクトが作成されます。しかし、CDataオブジェクトがガベージコレクションされても、malloc()によって割り当てられたメモリは解放されません。CDataオブジェクトがガベージコレクションされる前に、メモリを解放する外部関数を呼び出す必要があります。
ffi.from_buffer()を使用して、メモリビュー機能をサポートするPythonオブジェクトからメモリポインタを取得する場合、CDataOwnGCオブジェクトが得られます。これはターゲットオブジェクトの参照カウントを増やし、ガベージコレクションされないようにします。

以下ではffi.new()を使用して整数配列を作成し、返されるオブジェクトの型がCDataOwnであることを確認します。これはC言語データの解放権を持っていることを示し、これがガベージコレクションされると、対応するC言語データのメモリも同時に解放されます。

arr = ffi.new("int[]", 100)
type(arr)

_cffi_backend.__CDataOwn

C言語でffi.new()を使用して作成された配列を使用する場合、対応するCDataOwnオブジェクトがメモリ回収されないように注意する必要があります。以下の例では、Points構造体は2つのポインタを使用して、2つの1次元配列を指し、それぞれ点のX軸とY軸の座標値を保存します。

ffi = FFI()
ffi.cdef(
    """
typedef struct _Points{
    double * x;
    double * y;
    int n;
} Points;
"""
)

以下ではPoints構造体を作成し、そのxとyフィールドをffi.new()を使用して作成された2つの配列に設定します。❶xフィールドに代入する前後にxの参照カウントを表示します。出力から、代入文がこのオブジェクトの参照カウントを増やしていないことがわかります。❷作成された配列はPython変数で参照されていないため、フィールドの代入文が実行された後、この配列のメモリはガベージコレクションされます。points構造体のyフィールドは不正なポインタになり、それが指すメモリは後続のプログラムで使用される可能性があります。

n = 10
points = ffi.new("Points *", {"n": n})
x = ffi.new("double []", n)
print("refcount before assign:", sys.getrefcount(x))
points.x = x  # ❶
print("refcount after assign:", sys.getrefcount(x))
points.y = ffi.new("double []", n)  # ❷ 注意：このように使用しないでください

refcount before assign: 2
refcount after assign: 2

以下ではC言語関数malloc()を呼び出してメモリを割り当て、それが返されると、CFFIは割り当てられたメモリのアドレスを保持するCDataオブジェクトを作成します。このオブジェクトがガベージコレクションされても、malloc()によって割り当てられたメモリは解放されません。free()を呼び出す必要があります。

ffi = FFI()
ffi.cdef(
    """
void * malloc(size_t size);
void free(void * ptr);
"""
)
lib = ffi.dlopen("msvcrt.dll")

buf = lib.malloc(100)
print(type(buf))
lib.free(buf)

<class '_cffi_backend._CDataBase'>

メモリ解放を忘れないようにするために、ffi.gc(cdata, destructor)を使用できます。これはCDataGCPオブジェクトを返し、このオブジェクトがガベージコレクションされると自動的にdestructorが呼び出され、外部関数によって割り当てられたメモリが解放されます。

buf = ffi.gc(lib.malloc(100), lib.free)
type(buf)

_cffi_backend.__CDataGCP

以下のプログラムは、Pythonオブジェクトからメモリポインタを取得する場合の例を示しています。ffi.from_buffer()を呼び出す前、abytesの参照カウントは3です。呼び出した後、参照カウントは1増加し、bufを削除した後、参照カウントは3に戻ります。

import sys

abytes = b"abcdef"
print("refcount before call from_buffer():", sys.getrefcount(abytes))
buf = ffi.from_buffer(abytes)
print(type(buf))
print("refcount after call from_buffer():", sys.getrefcount(abytes))
del buf
print("refcount after release buf:", sys.getrefcount(abytes))

refcount before call from_buffer(): 3
<class '_cffi_backend.__CDataFromBuf'>
refcount after call from_buffer(): 4
refcount after release buf: 3

Pythonのオブジェクト構造#

Pythonのオブジェクトは、C言語レベルではすべて構造体として表現されています。CFFIを使用することで、Pythonオブジェクトとメモリを共有する構造体を作成し、その中身を読み書きすることが可能です。本節では、Pythonのfloat、tuple、およびlistオブジェクトのメモリにCFFIを使ってアクセスし、その仕組みを紹介します。

`float`オブジェクト#

浮動小数点数（float）オブジェクトの構造体は以下のようになります。

ob_refcnt はガベージコレクション用の参照カウンタであり、ob_type はこのオブジェクトの型へのポインタです。ob_fval フィールドには実際の浮動小数点値が格納されます。

typedef struct {
    Py_ssize_t ob_refcnt;   // 参照カウント
    PyTypeObject *ob_type;  // オブジェクトの型
    double ob_fval;         // float値（Cのdouble）
} PyFloatObject;

次に、PyFloatObject 構造体を定義し、float 型の変数 x に対して id() で取得したアドレスを ffi.cast() を使って PyFloatObject * 型のポインタにキャストします。簡便のため、ob_type フィールドはポインタ型ではなく、ポインタサイズの整数型として宣言しています。

ffi = FFI()
ffi.cdef(
    """
typedef struct {
ssize_t ob_refcnt;
ssize_t ob_type;
double ob_fval;
} PyFloatObject;
"""
)

x = 3.14
x_st = ffi.cast("PyFloatObject *", id(x))

次に、sys.getrefcount() と ob_refcnt を使って参照カウントを取得してみます。getrefcount() を呼び出すと、その引数として渡されたオブジェクトの参照カウントが一時的に 1 増加します。そのため、x_st.ob_refcnt の値は sys.getrefcount(x) よりも 1 少なくなります。

sys.getrefcount(x), x_st.ob_refcnt

(2, 1)

次に、ob_type フィールドの値が float 型のアドレスと一致することを確認します。これにより、ob_type が float 型（の型オブジェクト）へのポインタであることが分かります。

x_st.ob_type == id(float)

True

実際の数値の値は、ob_fval フィールドから取得できます。

x_st.ob_fval

3.14

Python では float 型は変更不可能（イミュータブル）な型ですが、C言語の構造体を経由することで、その値を強制的に書き換えることができます。以下のコードでは、まず x の参照である y を作成し、x_st を通じて ob_fval の値を変更します。その結果、x と y の両方の値が書き換えられます。

Warning

これはあくまでデモ目的のコードです。実際のコードではこのようにオブジェクトの内部を直接書き換えることは推奨されません。

y = x
x_st.ob_fval = 6.28
print(x, y)

6.28 6.28

`tuple`オブジェクト#

tupleオブジェクトのC言語における構造体は以下のとおりです。ob_refcntとob_typeについては前述のとおりです。ob_sizeはこのtupleオブジェクトの要素数を表し、ob_itemは要素のアドレスを格納する配列です。

一見するとob_itemの長さは1で、1つのアドレスしか格納できないように見えます。しかし実際には、tupleオブジェクトのメモリを確保する際に、その長さに応じた十分なメモリが確保されます。C言語では配列の長さチェックが行われないため、十分なメモリさえ確保されていれば、たとえばob_item[3]のようにアクセスしても問題ありません。

typedef struct {
    Py_ssize_t ob_refcnt;
    PyTypeObject *ob_type;
    Py_ssize_t ob_size; 
    PyObject *ob_item[1];
} PyTupleObject;

次に、(1, 2, 3)というtupleオブジェクトのサイズを確認してみましょう。ヘッダー部分の3つの8バイト（ob_refcnt、ob_type、ob_size）と、要素のポインタ3つ分の8バイトずつを合わせて、合計で48バイトになります。

a = (1, 2, 3)
a.__sizeof__()

CFFIでは、次のようにtupleの構造体を宣言します。CFFIでは配列のインデックス範囲がチェックされるため、ob_item[0]のように宣言しておくと、そのチェックが無効になります。次に、ob_typeフィールドを使って、そのオブジェクトがtupleであるかどうかを確認します。

ffi = FFI()
ffi.cdef(
    """
typedef struct {
    ssize_t ob_refcnt;
    ssize_t ob_type;
    ssize_t ob_size;
    ssize_t ob_item[0];
} PyTupleObject;
"""
)

a_st = ffi.cast("PyTupleObject *", id(a))
a_st.ob_type == id(tuple)

True

次にob_sizeフィールドでこのtupleオブジェクトのサイズを取得します。

a_st.ob_size

次に、ob_itemフィールドに格納されているアドレスを確認します。

for i in range(len(a)):
    print(i, a_st.ob_item[i] == id(a[i]))

True
True
True

tupleはPythonでは変更不可能（イミュータブル）なオブジェクトですが、CFFIを経由すれば、その要素を変更することが可能です。

次のコードでは、まずxの参照カウントを出力し、その後xのアドレスをob_item[1]に書き込みます。このような直接メモリを操作する方法では、xの参照カウントは更新されないため、書き込んだ後もカウントは変化しません。しかし、tupleの1番目の要素がxに置き換わっていることが確認できます。

ただし、この例では参照カウントの更新が行われていないため、xがガベージコレクション（GC）された後、tupleの1番目の要素は解放済みメモリ（いわゆる“ごみ”）を指すことになります。

x = "abcd"
print(sys.getrefcount(x))
a_st.ob_item[1] = id(x)
print(a)
print(sys.getrefcount(x))
a

3
(1, 'abcd', 3)
3

(1, 'abcd', 3)

`list`オブジェクト#

上のtupleは変更不可能なオブジェクトであるため、メモリを確保する際にob_itemフィールドに十分なサイズを割り当てておけば問題ありません。しかし、listの場合は作成後に要素の追加や削除が可能なため、構造は次のようになります。

ob_sizeフィールドはリストの現在の長さを保持します。ob_itemフィールドには、動的に確保された配列のアドレスが格納されます。そして、allocatedフィールドには、この動的配列に確保されたメモリのサイズ（要素数）が保持されます。

typedef struct {
    Py_ssize_t ob_refcnt;
    PyTypeObject *ob_type;
    Py_ssize_t ob_size; 
    PyObject **ob_item;
    Py_ssize_t allocated;
} PyListObject;

次に、空のリストのサイズを確認してみましょう。PyListObject構造体には5つのフィールドがあり、その合計サイズは40バイトです。

x = []
x.__sizeof__()

次に、要素が3つあるリストのサイズを確認します。__sizeof__()の結果が72バイトで、空リストのサイズ40バイトとの差は32バイトになります。これは、4つ分のポインタを格納できるob_item配列が確保されていることを示しています。

x = [1, 2, 3]
x.__sizeof__()

次のコードでは、CFFIを使ってPyListObject構造体を宣言し、リストxをその構造体にキャストしてx_stとして扱います。

ffi = FFI()
ffi.cdef(
    """
typedef struct {
    ssize_t ob_refcnt;
    ssize_t ob_type;
    ssize_t ob_size; 
    ssize_t *ob_item;
    ssize_t allocated;
} PyListObject;
"""
)

x_st = ffi.cast("PyListObject *", id(x))

ob_sizeとallocatedフィールドの値を確認します。確保された配列は4つの要素を格納できるサイズになっていますが、実際には3つの要素が格納されていることがわかります。

x_st.ob_size, x_st.allocated

(3, 4)

次に、ob_itemフィールドにはリストの要素のアドレスが格納されていることを確認します。

for i in range(len(x)):
    print(i, x_st.ob_item[i] == id(x[i]))

True
True
True

ob_itemにはまだ1つの空きスペースがあるため、1つの要素を追加しても、メモリの再分配は発生しません。

print(x_st.ob_item)
x.append("abc")
print(x_st.ob_item)

<cdata 'ssize_t *' 0x00000170410CB450>
<cdata 'ssize_t *' 0x00000170410CB450>

さらに1つの要素を追加すると、現在のob_item配列がいっぱいになり、メモリが再分配されて配列のアドレスが変わることが確認できます。また、今回は8要素を格納できる新しい配列が確保されました。このように、配列が満杯になるたびにメモリの再分配が行われ、リストにどんどん要素を追加することができます。再分配は、現状の要素数に対して一定の伸び率でob_itemの長さを決定するため、再分配の回数はそれほど多くなく、append()操作はO(1)の計算量として扱うことができます。

print(x_st.ob_item)
x.append("xyz")
print(x_st.ob_item)
print(x_st.ob_size, x_st.allocated)

<cdata 'ssize_t *' 0x00000170410CB450>
<cdata 'ssize_t *' 0x0000017040ECC1B0>
5 8

9軸センサー融合アルゴリズム#

9軸センサーとは、加速度センサー、角速度センサー、および磁気センサーのことで、スマートフォンやVRデバイスで最も一般的なセンサーです。加速度センサーはデバイスの3軸に沿った加速度を検出し、角速度センサーはデバイスの3軸周りの回転速度を検出し、磁気センサーは3軸方向の磁場強度を検出します。センサーの融合アルゴリズムを使用すると、これらのセンサーのデータに基づいて、デバイスの地球座標系における回転姿勢を計算することができます。

ABIモード#

以下のコマンドを実行して動的リンクライブラリファイルfusion.dllにコンパイルします。

!gcc -shared Fusion/Fusion/*.c -I. -o fusion.dll

以下は、C言語でFusionを使用する例です。まず、FusionAhrsInitialise()を呼び出して構造体を初期化します。次に、FusionAhrsSetSettings()を呼び出して各種設定を行います。ループ内で各センサー構造体のデータを設定し、FusionAhrsUpdate()を呼び出して融合アルゴリズムの状態を更新します。その計算結果は、四元数を表すquaternionフィールドに保存されます。四元数は3次元空間の回転を表すために使用され、回転軸と回転角を簡単に提供することができます。

#include "../../Fusion/Fusion.h"
#include <stdbool.h>
#include <stdio.h>
#include <time.h>

#define SAMPLE_RATE (100) // replace this with actual sample rate
#define SAMPLE_PERIOD (0.01f) // replace this with actual sample period

int main() {
    FusionAhrs ahrs;
    FusionAhrsInitialise(&ahrs);

    // Set AHRS algorithm settings
    const FusionAhrsSettings settings = {
            .convention = FusionConventionNwu,
            .gain = 0.5f,
            .gyroscopeRange = 2000.0f, /* replace this with actual gyroscope range in degrees/s */
            .accelerationRejection = 10.0f,
            .magneticRejection = 10.0f,
            .recoveryTriggerPeriod = 5 * SAMPLE_RATE, /* 5 seconds */
    };
    FusionAhrsSetSettings(&ahrs, &settings);

    // This loop should repeat each time new gyroscope data is available
    while (true) {
        const clock_t timestamp = clock(); // replace this with actual gyroscope timestamp
        FusionVector gyroscope = {0.0f, 0.0f, 0.0f}; // replace this with actual gyroscope data in degrees/s
        FusionVector accelerometer = {0.0f, 0.0f, 1.0f}; // replace this with actual accelerometer data in g
        FusionVector magnetometer = {1.0f, 0.0f, 0.0f}; // replace this with actual magnetometer data in arbitrary units

        // Update gyroscope AHRS algorithm
        FusionAhrsUpdate(&ahrs, gyroscope, accelerometer, magnetometer, SAMPLE_PERIOD);
    }
} 

C言語のヘッダーファイルでは、FusionVectorは共用体であり、arrayまたはaxisフィールドを通じて3軸のデータにアクセスできます。

typedef union {
    float array[3];

    struct {
        float x;
        float y;
        float z;
    } axis;
} FusionVector;

Pythonではx, y, zフィールドのみを使用してデータにアクセスするため、以下のように簡略化できます。

typedef struct {
    float x;
    float y;
    float z;
} FusionVector;

以下は、ヘッダーファイルでの各データ型および関数の定義です。FusionAhrsUpdate()の3つのFusionVector型のセンサーパラメータは値渡しパラメータであることに注意してください。

fusion_types = """
typedef struct {
    float x;
    float y;
    float z;
} FusionVector;

typedef struct {
    float w;
    float x;
    float y;
    float z;
} FusionQuaternion;

typedef enum {
    FusionConventionNwu, /* North-West-Up */
    FusionConventionEnu, /* East-North-Up */
    FusionConventionNed, /* North-East-Down */
} FusionConvention;

typedef struct {
    FusionConvention convention;
    float gain;
    float gyroscopeRange;
    float accelerationRejection;
    float magneticRejection;
    unsigned int recoveryTriggerPeriod;
} FusionAhrsSettings;

typedef struct {
    FusionAhrsSettings settings;
    FusionQuaternion quaternion;
    FusionVector accelerometer;
    bool initialising;
    float rampedGain;
    float rampedGainStep;
    bool angularRateRecovery;
    FusionVector halfAccelerometerFeedback;
    FusionVector halfMagnetometerFeedback;
    bool accelerometerIgnored;
    int accelerationRecoveryTrigger;
    int accelerationRecoveryTimeout;
    bool magnetometerIgnored;
    int magneticRecoveryTrigger;
    int magneticRecoveryTimeout;
} FusionAhrs;
"""

fusion_functions = """
void FusionAhrsInitialise(FusionAhrs *const ahrs);
void FusionAhrsSetSettings(FusionAhrs *const ahrs, const FusionAhrsSettings *const settings);
void FusionAhrsUpdate(
    FusionAhrs *const ahrs, 
    const FusionVector gyroscope, 
    const FusionVector accelerometer, 
    const FusionVector magnetometer, 
    const float deltaTime);
"""

以下でFFIオブジェクトを作成し、各種構造体と関数シグネチャを宣言し、動的リンクライブラリfusion.dllをロードします。

import cffi

ffi = cffi.FFI()
ffi.cdef(fusion_types)
ffi.cdef(fusion_functions)
lib = ffi.dlopen("fusion.dll")

使いやすくするために、以下のFusionクラスで外部関数をオブジェクトとしてラッピングします。❶四元数構造体fusion.quaternionの値を読み取りやすくするために、構造体をオブジェクトの属性qとして保存します。この属性はfusion.quaternionのメモリアドレスと同じです。❷FusionAhrsUpdate()の3つのセンサーパラメータはすべて値渡しで渡されるため、まず添字演算を使用してポインタが指す構造体を取得します。❸dataパラメータのセンサーの各値を対応する構造体のフィールドに書き込んだ後、FusionAhrsUpdate()を呼び出して計算を完了します。❹最後に、四元数構造体の各フィールドを含むタプルオブジェクトを返します。

class Fusion:
    def __init__(self, period, gain):
        self.period = period
        self.fusion = ffi.new("FusionAhrs *")
        self.settings = ffi.new("FusionAhrsSettings *")
        self.settings.gain = gain
        lib.FusionAhrsInitialise(self.fusion)
        lib.FusionAhrsSetSettings(self.fusion, self.settings)
        self.gyroscope = ffi.new("FusionVector *")
        self.accelerometer = ffi.new("FusionVector *")
        self.magnetometer = ffi.new("FusionVector *")

        self.q = self.fusion.quaternion  # ❶
        self.gyr = self.gyroscope[0]  # ❷
        self.acc = self.accelerometer[0]
        self.mag = self.magnetometer[0]

    def update(self, data):
        self.gyr.x = data["gx"]
        self.gyr.y = data["gy"]
        self.gyr.z = data["gz"]
        self.acc.x = data["ax"]
        self.acc.y = data["ay"]
        self.acc.z = data["az"]
        self.mag.x = data["mx"]
        self.mag.y = data["my"]
        self.mag.z = data["mz"]
        lib.FusionAhrsUpdate(
            self.fusion, self.gyr, self.acc, self.mag, self.period  # ❸
        )
        q = self.q
        return (q.w, q.x, q.y, q.z)  # ❹

以下でCSVファイルからセンサーデータを読み込み、update()メソッドをループしてセンサーの姿勢を計算します。

import polars as pl
import numpy as np

df = pl.read_csv("data/imu_sample_data.csv").select(pl.col("*").name.map(str.strip))

quaternions = []
fusion = Fusion(period=0.01, gain=0.5)
for row in df.iter_rows(named=True):
    quaternions.append(fusion.update(row))
quaternions = np.array(quaternions)
quaternions

array([[ 1.00007558e+00, -3.36902449e-04, -1.18124088e-04,
        -9.39859636e-03],
       [ 1.00001717e+00, -5.45008166e-04, -2.44000985e-04,
        -1.44591173e-02],
       [ 9.99948204e-01, -8.69491778e-04, -3.46390007e-04,
        -1.86079293e-02],
       ...,
       [ 9.98859704e-01,  4.99216840e-02, -4.33470076e-03,
        -3.36576160e-03],
       [ 9.98873293e-01,  4.96541448e-02, -4.23847046e-03,
        -3.45249730e-03],
       [ 9.98875380e-01,  4.96137887e-02, -4.14189976e-03,
        -3.53873521e-03]])

APIモード#

Pythonで外部関数をループ呼び出しする場合、PythonオブジェクトをC言語のデータに変換する必要があり、この呼び出しインターフェースは時間がかかります。計算速度を向上させるために、C言語で配列内の各要素に対してFusionAhrsUpdate()をループ呼び出しする関数FusionAhrsUpdateArray()を記述することができます。以下のプログラムを実行して、APIモードでFusionAhrsUpdateArray()とFusionライブラリのC言語コードを一緒にコンパイルし、拡張モジュール_fusionを作成します。

❶PythonではFusionAhrs構造体を作成する必要がありますが、そのフィールドにアクセスする必要はないため、フィールドの宣言を省略できます。C言語コンパイラは、ヘッダーファイルの宣言に基づいてフィールド情報を補完します。❷C言語では、共用体FusionVectorのarrayフィールドを使用してデータにアクセスできるため、プログラムの記述が容易になります。

from pathlib import Path

ffi = cffi.FFI()
ffi.cdef(
    """
typedef enum {
    FusionConventionNwu, /* North-West-Up */
    FusionConventionEnu, /* East-North-Up */
    FusionConventionNed, /* North-East-Down */
} FusionConvention;

typedef struct {
    FusionConvention convention;
    float gain;
    float gyroscopeRange;
    float accelerationRejection;
    float magneticRejection;
    unsigned int recoveryTriggerPeriod;
} FusionAhrsSettings;

typedef struct {
...;  //{1}
} FusionAhrs;

void FusionAhrsInitialise(FusionAhrs * fusionAhrs);
void FusionAhrsSetSettings(FusionAhrs *const ahrs, const FusionAhrsSettings *const settings);
void FusionAhrsUpdateArray(FusionAhrs * fusionAhrs, 
        float samplePeriod, float * in_data, float * out_data, int length);
"""
)

ffi.set_source(
    "_fusion",
    """
#include "FusionAhrs.h"

void FusionAhrsUpdateArray(FusionAhrs * fusionAhrs, 
        float samplePeriod, 
        float * in_data, float * out_data, int length) {
    FusionVector gyroscope;
    FusionVector accelerometer;
    FusionVector magnetometer;
    int i=0, j=0, n, k;
    for(n=0;n<length;n++){
        for(k=0;k<3;k++) gyroscope.array[k] = in_data[i++]; //{2}
        for(k=0;k<3;k++) accelerometer.array[k] = in_data[i++];
        for(k=0;k<3;k++) magnetometer.array[k] = in_data[i++];
        
        FusionAhrsUpdate(fusionAhrs, 
                         gyroscope, accelerometer, magnetometer, 
                         samplePeriod);

        for(k=0;k<4;k++) out_data[j++] = fusionAhrs->quaternion.array[k];
    }
}
""",
    include_dirs=[str(Path("Fusion/Fusion").absolute())],
    extra_objects=["Fusion/Fusion/*.c"],
)
ffi.compile();

以下でFusionクラスを使用して、上記の外部関数をオブジェクトとしてラッピングします。FusionAhrsUpdateArray()は、入力配列がC言語の連続メモリであり、要素が単精度浮動小数点数である必要があるため、呼び出す前にnp.ascontiguousarray()を使用して入力配列を指定された形式に変換します。

from _fusion import ffi, lib


class Fusion2:
    def __init__(self, period, gain):
        self.period = period
        self.fusion = ffi.new("FusionAhrs *")
        self.settings = ffi.new("FusionAhrsSettings *")
        self.settings.gain = gain
        lib.FusionAhrsInitialise(self.fusion)
        lib.FusionAhrsSetSettings(self.fusion, self.settings)

    def update(self, data):
        in_data = np.ascontiguousarray(data, dtype=np.float32)
        out_data = np.zeros((len(in_data), 4), dtype=np.float32)
        lib.FusionAhrsUpdateArray(
            self.fusion,
            self.period,
            ffi.from_buffer("float *", in_data),
            ffi.from_buffer("float *", out_data),
            len(in_data),
        )
        return out_data

CSVファイルからセンサーデータを配列sensor_dataに変換し、この配列をupdate()メソッドに渡します。C言語のループで9軸センサーの姿勢を計算し、その結果をABIモードで計算した結果と比較します。

sensor_data = df.select("gx", "gy", "gz", "ax", "ay", "az", "mx", "my", "mz").to_numpy()

f = Fusion2(period=0.01, gain=0.5)
quaternion2 = f.update(sensor_data)
np.allclose(quaternions, quaternion2)

True

データの可視化#

最後に、quaternion ライブラリで得られた姿勢を表すクォータニオン配列を回転行列に変換し、これらの回転行列を用いて箱を回転させ、plotly で 3D の動画を作成します。

import plotly.io as pio

pio.renderers.default = "plotly_mimetype+notebook_connected"

import pyvista as pv
from helper.plotly import polydata_to_lines
from plotly import graph_objects as go
import quaternion as Q

cube = pv.Cube(z_length=0.3).extract_all_edges()
scatter = polydata_to_lines(cube)
points = np.vstack([scatter.x, scatter.y, scatter.z])

q = Q.from_float_array(quaternions)
rotated_points = Q.as_rotation_matrix(q) @ points

frames = []

for i, (x, y, z) in enumerate(rotated_points[::5]):
    frames.append(
        go.Frame(
            data=[
                go.Scatter3d(
                    x=x, y=y, z=z, mode="lines", line=dict(color="blue", width=2)
                )
            ],
            name=str(i),
        )
    )

fig = go.Figure(
    data=frames[0].data,
    layout=go.Layout(
        title="Fusion Result",
        height=500,
        width=500,
        scene=dict(xaxis_title="X", yaxis_title="Y", zaxis_title="Z"),
        updatemenus=[
            dict(
                type="buttons",
                showactive=False,
                buttons=[
                    dict(
                        label="Play",
                        method="animate",
                        args=[None, {"frame": {"duration": 30, "redraw": True}}],
                    ),
                    dict(
                        label="Pause",
                        method="animate",
                        args=[[None], {"mode": "immediate", "frame": {"duration": 0}}],
                    ),
                ],
            )
        ],
    ),
    frames=frames,
)

fig.layout.yaxis = dict(scaleanchor="x", scaleratio=1)
fig.layout.scene.aspectratio = dict(x=1, y=1, z=1)
fig.layout.scene.aspectmode = "manual"
rng = 1.5
fig.layout.scene.xaxis.range = -rng, rng
fig.layout.scene.yaxis.range = -rng, rng
fig.layout.scene.zaxis.range = -rng, rng
fig.show()

CFFIで外部関数を呼び出す

Contents

CFFIで外部関数を呼び出す#

2つの呼び出しモード#

ABIモード#

APIモード#

cffiマジックコマンド#

C言語のデータ型#

配列とポインタ#

文字列とバイト列#

構造体#

CFFIのメモリ管理#

Pythonのオブジェクト構造#

floatオブジェクト#

tupleオブジェクト#

listオブジェクト#

9軸センサー融合アルゴリズム#

ABIモード#

APIモード#

データの可視化#

`float`オブジェクト#

`tuple`オブジェクト#

`list`オブジェクト#