Pythonでプログラムを実装したが、処理速度の観点でボトルネックになっている部分があるため、C言語を用いて高速化したい。
また、C言語からPythonで定義した関数も呼び出したい。
これらを実現する方法について教えて欲しい。
こんなお悩みを解決します。
今回は、求根アルゴリズムの1つとして有名なニュートン法を題材に、PythonからC言語の関数を呼び出す方法・C言語からPythonの関数を呼び出す方法について解説します。
Python、C言語の両方から関数を呼び出す方法について解説するので、興味がある方はぜひ最後までご覧ください。
全体構成
下記の2パターンについて解説する上で、全体構成を説明したいと思います。
- PythonからC言語の関数を呼び出すケース
- C言語からPythonの関数を呼び出すケース
PythonからC言語の関数を呼び出すケース
Pythonで実装された関数のうち、C言語に置き換えたい部分を以下のように表現します。
上記の「C言語に置き換えたいPythonのコード」を置き換えて呼び出す場合、呼び出し関係は以下のようになります。
C言語からPythonの関数を呼び出すケース
同様に、C言語から呼び出したいPythonの関数を以下のように表現します。
上記の「C言語から呼び出したいPythonの関数」をC言語から呼び出す場合、呼び出し関係は以下のようになります。
今回は、非線形連立方程式をNewton法で解くことになるため、C言語からPythonの関数を呼び出す必要があります。
したがって、以降では「C言語からPythonの関数を呼び出すケース」を対象に解説を進めます。
対象とする非線形連立方程式
ここでは、以下の非線形連立方程式を解くことを考えます。
$$
\begin{eqnarray}
\left\{
\begin{array}{lcc}
x + 2y + 1 & = & 0\\
x^{2} + 2y^{2} - 3 & = & 0\\
\end{array}
\right.
\end{eqnarray}
$$
また、この解析解は、以下のようになります。
\((x, y) = (1, -1), \left(-\dfrac{5}{3}, \dfrac{1}{3} \right)\)
Pythonによる実装
Pythonによる実装例は、以下のようになります。
今回は、numpy等のライブラリを用いて実装しています。
また、Github上の以下のファイルと対応します。
https://github.com/yuruto-free/python-c-api-newton/blob/v0.1.0/original/newton_method.py
import numpy as np
def objective_function(vec):
"""
Parameters
----------
vec : numpy.ndarray
vector of x
vec[i]: x_{i}
Returns
-------
fv : numpy.ndarray
vector of f(vec)
"""
x, y = vec
fv = np.array([
x + 2 * y + 1.0,
x * x + 2 * y * y - 3.0,
])
return fv
def calc_Jacobian_matrix(vec, delta=1e-3):
"""
Parameters
----------
vec : numpy.ndarray
vector of x
vec[i]: x_{i}
delta : double
fixed difference value
Returns
-------
Jf : numpy.ndarray
Jacobian matrix
Jf[i, j] = \dfrac{\partial f_{i}}{\partial x_{j}}
"""
ndim = vec.size
Jf = np.empty((ndim, ndim), dtype=np.float64)
for idx in np.arange(ndim):
val = vec[idx]
# calculate f(x + delta)
vec[idx] = val + delta
right_f= objective_function(vec)
# calculate f(x - delta)
vec[idx] = val - delta
left_f = objective_function(vec)
# calculate 0.5 * (f(x + delta) - f(x - delta)) / delta
diffs = 0.5 * (right_f - left_f) / delta
Jf[:, idx] = diffs
vec[idx] = val
return Jf
def newton_method(vec, max_iter, tol=1e-10, delta=1e-3):
"""
Parameters
----------
vec : numpy.ndarray
vector of x
vec[i]: x_{i}
max_iter : int
maximum iteration
tol : double
relative tolerance
Returns
-------
hat : numpy.ndarray
vector of \hat{x}
hat[i]: \hat{x}_{i}
is_convergent : boolean
convergent status
True: convergent
False: divergence
"""
hat = np.copy(vec)
is_convergent = False
for _ in np.arange(max_iter):
# convergence test
if np.linalg.norm(hat) < tol:
is_convergent = True
break
# Step1: calculate Jacobian matrix
Jf = calc_Jacobian_matrix(hat, delta=delta)
# Step2: calculate function value
fv = objective_function(hat)
# Step3: calculate delta_v
delta_v = np.linalg.solve(Jf, -fv)
# Step4: update vec
hat += delta_v
return [hat, is_convergent]
if __name__ == '__main__':
max_iter = 1000
tol = 1e-10
# ===============
# solve pattern 1
# ===============
print('Solve pattern 1')
vec = np.array([2.0, 1.0])
exact_vec = np.array([1.0, -1.0])
print('[init value] x: {:.5f}, y: {:.5f}'.format(*vec))
hat, _ = newton_method(vec, max_iter, tol=tol)
err = np.linalg.norm(hat - exact_vec)
print('[estimated] x: {:.5f}, y: {:.5f} ({:.5e})'.format(*hat, err))
print('')
# ===============
# solve pattern 2
# ===============
print('Solve pattern 2')
vec = np.array([-1.0, 1.0])
exact_vec = np.array([-5.0/3.0, 1.0/3.0])
print('[init value] x: {:.5f}, y: {:.5f}'.format(*vec))
hat, _ = newton_method(vec, max_iter, tol=tol)
err = np.linalg.norm(hat - exact_vec)
print('[estimated] x: {:.5f}, y: {:.5f} ({:.5e})'.format(*hat, err))それぞれの関数は、以下のような内容となります。
| 関数名 | 内容 |
|---|---|
objective_function | 目的関数 今回の解く対象となる非線形連立方程式が定義されている。 |
calc_Jacobian_matrix | ヤコビ行列の計算用関数 |
newton_method | ニュートン法により近似解を求める関数 |
また、実行すると以下のような出力が得られます。
Solve pattern 1
[init value] x: 2.00000, y: 1.00000
[estimated] x: 1.00000, y: -1.00000 (0.00000e+00)
Solve pattern 2
[init value] x: -1.00000, y: 1.00000
[estimated] x: -1.66667, y: 0.33333 (2.28878e-16)C言語に置き換える関数・C言語から呼び出すPython関数
定義したPython関数を「C言語に置き換える関数」、「C言語から呼び出すPython関数」に振り分けると以下のようになります。
| 関数名 | 呼び出し関係 |
|---|---|
objective_function | C言語から呼び出すPython関数 |
calc_Jacobian_matrix | C言語に置き換える関数 |
newton_method | C言語に置き換える関数 |
C言語による実装
ここから、C言語による関数の置き換え方について説明していきます。
また、Github上に一通り格納しているので、参考にしてください。
https://github.com/yuruto-free/python-c-api-newton/tree/v0.1.0/c-api
全体構成
C言語による実装結果をPythonから呼び出すために、C言語による実装結果を共有ライブラリ（.so）に変換後、Pythonからロードする、という方法を取ります。
このため、以下に示すように、wrapper_newtonlib.soファイルを生成する必要があります。
wrapper_newtonlib.soファイル作成時の対象範囲ここでは、上記を踏まえたディレクトリ構成にしております。
具体的な全体構成は以下のようになっております。
./
|-- compile.sh
|-- newton-method-with-clib.py
|-- setup.py
`-- libs/
|-- include/
| `-- wrapper_newton.h
`-- src/
|-- wrapper.c
`-- newton_method.c
それぞれのファイルは、以下のような役割を持ちます。
| ディレクトリ名 | ファイル名 | 役割 |
|---|---|---|
| . | compile.sh | C言語ファイルをコンパイルするためのshell script |
| newton-method-with-clib.py | C言語での実装結果を用いて非線形連立方程式をNewton法で解くスクリプト | |
| setup.py | C言語ファイルをコンパイルするためのPython script | |
| libs/include | wrapper_newton.h | C言語でNewton法を実装する際のAPI定義 |
| libs/src | wrapper.c | 前処理、関数呼び出し、後処理部分を実装したファイル |
| newton_method.c | Newton法本体 |
以降では、以下に示すステップで解説をしていきます。
- 必要なヘッダーの読み込み
- 実装時の注意点
- C言語から呼び出すPython関数の登録
- C言語によるラッパー関数
newton_methodの前処理部分newton_methodにおけるコールバック関数newton_methodの後処理部分- モジュール定義部分（本体実装である
newton_method.cの解説は省略）
必要なヘッダーの読み込み
PythonのオブジェクトをC言語で処理するためには、Pythonのデータ型を扱う変数やC言語のデータ型に変換するための関数が必要になります。
これらは、Pythonのライブラリとして提供されており、以下の形式で利用できるようになります。
#define PY_SSIZE_T_CLEAN
#include <Python.h> // Python用のヘッダー
#include <numpy/arrayobject.h> // Numpy用のヘッダー実装時の注意点
C言語では、プログラム実装時に用いる変数を定義する必要があるため、実装者が変数を管理します。
一方、Pythonは変数の管理はPython側で良きに処理しています。
具体的には、Pythonスクリプトにおいて、どこからも参照されなくなったオブジェクトは自動的に割り当てが解除されます。
言い換えると、オブジェクトの参照状態が適切に管理されていない場合、解除忘れのオブジェクトがメモリ上に残り続けるまたは、解除済みのオブジェクトにアクセスすることになります。
これは、メモリリークやメモリアクセス違反に繋がります。
C言語でPythonオブジェクトを管理するために、Py_INCREFやPyDECREFといった「参照カウント」を増減させるマクロが用意されています。
これらのマクロを適切に呼び出しつつ、Python/C APIを利用することになります。
上記の実装をする上で、Python/C APIが返却する変数のパターンを知っておく必要があります。
以降では、これらのパターンについて解説します。
Python/C APIが返却する変数のパターン
Python/C APIが返却する変数のパターンとして、「New reference」と「Borrowed reference」の2種類があります。
それぞれの違いは、以下のようになります。
| 対象 | 内容 |
|---|---|
| New reference | 返却するPythonのオブジェクトに所有権を付与した状態で返す（所有参照）。 不要になったら、 Py_DECREFで参照カウントを減らす必要がある。 |
| Borrowed reference | 返却するPythonオブジェクトに所有権を付与しない状態で返す（借用参照）。 参照カウントを減らしてはいけない。 |
また、一部のAPIには、参照を盗む（APIの内部でPy_DECREFが"勝手に"呼ばれるイメージ）ものも存在します。
ここでは、参照関係の理解のため、リストオブジェクトを題材に事例を取り上げていきたいと思います。
https://docs.python.org/ja/3/c-api/stable.html
リストオブジェクトによる事例
上記のリンクにアクセスすると、以下のAPIが存在することが確認できると思います。
それぞれ、以下のような関係になります。
| API名 | 参照関係 | 使用時の注意点 |
|---|---|---|
PyList_New | New reference | 所有権が付与されたオブジェクトを返却するため、参照カウントの管理が必要になる。 |
PyList_GetItem | Borrowed reference | 所有権のないオブジェクトを返却する。参照カウントを減らしてはならない。 |
PyList_SetItem | - | 所有権を盗み取るAPIのため、API呼出し後もオブジェクトを利用する場合、呼び出す前にPy_INCREFで参照カウントを増やしておく必要がある。 |
また、Pythonで下記のようなコードをC言語に置き換える場合も注意が必要となります。
src_data = [3, 1, 4, 1, 6]
dst_data = [6, 2, 8, 2, 2]
item = src_data[2] # 4
dst_data[2] = item # 8 -> 4上記の内容を素直にC言語で実装すると以下のようになりますが、オブジェクトを又貸ししている状態となるため、実装時には注意が必要です。
// 不適切な実装例
PyObject *item = PyList_GetItem(src_data, 2); // 所有権の無いオブジェクトを返却する
PyList_SetItem(dst_data, 2, item); // 所有権を盗み取る
// このまま、以下を実行した場合、該当オブジェクトがないためエラーとなる
// PyList_GetItem(src_data, 2);以降では、今回の実装内容について解説していきます。
C言語から呼び出すPython関数の登録
C言語からPythonの関数を呼び出すためには、Pythonの関数を登録するC言語の関数が必要になります。
ここでは、グローバル変数に該当する関数を登録し、Newton法を実行する際に参照する方針とします。
実装結果は、以下のようになります。
//! define objective function in python code
static PyObject *py_objective_function = NULL;
/**
* @fn static PyObject *update_function(PyObject *self, PyObject *args)
* @brief set objective function used by Newton's method
* @param[in] self python object
* @param[in] args arguments from python script
* @return Py_None
*/
static PyObject *update_function(PyObject *self, PyObject *args) {
PyObject *tmp;
if (!PyArg_ParseTuple(args, "O:set_callback", &tmp)) {
PyErr_SetString(PyExc_TypeError, "Invalid argument");
goto EXIT_SET_OBJECTIVE_FUNCTION;
}
if (!PyCallable_Check(tmp)) {
PyErr_SetString(PyExc_TypeError, "parameter must be callable");
goto EXIT_SET_OBJECTIVE_FUNCTION;
}
Py_XINCREF(tmp); //! add a reference to new callback
Py_XDECREF(py_objective_function); //! dispose of previous callback
py_objective_function = tmp;
EXIT_SET_OBJECTIVE_FUNCTION:
Py_RETURN_NONE;
}また、利用方法としては、以下のようになります。
import wrapper_newtonlib as newtonlib
def objective_function(vec):
pass
newtonlib.set_objective_function(objective_function)上記を例に引数の処理方法、参照カウントの使い方を解説していきます。
引数の処理方法
Pythonでキーワード引数を利用しない場合、C言語側での宣言方法は以下のようになります。キーワード引数を用いる例は別途、解説します。
static PyObject *called_function_from_python(PyObject *self, PyObject *args)ここで、selfはモジュールオブジェクト、argsが関数の引数となります。
今回の例の場合、selfにはset_objective_functionなどが登録されています。ほとんど使うことはないかと思います。
引数の処理には、PyArg_ParseTuple関数を利用します。
下記にある通り、フォーマットを指定してデータを読み込みます。
https://docs.python.org/ja/3/c-api/arg.html
今回、Pythonで定義した関数（=オブジェクト）を引数に取るため、下記のような処理となります。
// Pythonオブジェクトとしてそのまま読み込む
if (!PyArg_ParseTuple(args, "O:set_callback", &tmp)) {
PyErr_SetString(PyExc_TypeError, "Invalid argument");
goto EXIT_SET_OBJECTIVE_FUNCTION;
}
// 渡されたオブジェクトが呼び出し可能かを確認
if (!PyCallable_Check(tmp)) {
PyErr_SetString(PyExc_TypeError, "parameter must be callable");
goto EXIT_SET_OBJECTIVE_FUNCTION;
}参照カウントの増減
変数取得後は、参照カウントを増減させ、引数で受け取った関数を登録します。
ここで、Py_XINCREF、Py_XDECREFは、NULLを許容するバージョンのPy_INCREF、Py_DECREFとなります。
Py_XINCREF(tmp); //! 参照カウントを増加させる
Py_XDECREF(py_objective_function); //! 参照カウントを減少させる
py_objective_function = tmp;最後に、Py_RETURN_NONEマクロを用いて関数の戻り値を返却します。
Py_RETURN_NONEは、戻り値としてNoneを返すマクロであり、以下と同じ働きをします。
PyObject *ret = Py_None;
Py_INCREF(ret);
return ret;C言語によるラッパー関数
以降では、下記のPython関数をC言語に置き換え、Pythonから呼び出すための実装方法について解説します。
newton_method(vec, max_iter, tol=1e-10, delta=1e-3)C言語による実装を以下に示します。
#define NEWTONLIB_DELTA (1e-3)
#define NEWTONLIB_TOL (1e-10)
// 中略
/**
* @fn static PyObject *newtonlib(PyObject *self, PyObject *args)
* @brief wrapper function of newton_method
* @param[in] self python object
* @param[in] args arguments from python script
* @param[in] keywds keyword arguments from python script
* @return py_output
*/
static PyObject *newtonlib(PyObject *self, PyObject *args, PyObject *keywords) {
PyObject *py_vec, *py_output;
PyArrayObject *np_arr;
int32_t idx;
int32_t func_val;
int32_t ndim;
double *vec, *hat;
static char *kwlist[] = {"vec", "max_iter", "tol", "delta", NULL};
struct NewtonLib_ArgParam_t params;
struct NewtonLib_Output_t output;
//! initialize
py_vec = NULL;
py_output = PyList_New(0);
vec = NULL;
hat = NULL;
params.tol = (double)NEWTONLIB_TOL;
params.delta = (double)NEWTONLIB_DELTA;
//! Step1: check arguments
if (!PyArg_ParseTupleAndKeywords(args, keywords, "Oi|dd", kwlist, &py_vec, ¶ms.max_iter, ¶ms.tol, ¶ms.delta)) {
PyErr_SetString(PyExc_TypeError, "Invalid argument");
goto EXIT_NEWTONLIB;
}
//! Step2: check object type
if (!PyArray_Check(py_vec)) {
PyErr_SetString(PyExc_TypeError, "'vec' needs to numpy array");
goto EXIT_NEWTONLIB;
}
//! Step3: get vector size
np_arr = (PyArrayObject *)py_vec;
ndim = (int32_t)PyArray_SIZE(np_arr);
//! Step4: malloc
vec = (double *)malloc(sizeof(double) * ndim);
hat = (double *)malloc(sizeof(double) * ndim);
if ((NULL == vec) || (NULL == hat)) {
PyErr_SetString(PyExc_MemoryError, "lack of memory");
goto EXIT_NEWTONLIB;
}
//! Step5: copy data
for (idx = 0; idx < ndim; idx++) {
vec[idx] = *((double *)PyArray_GETPTR1(np_arr, idx));
}
//! Step6: call newton_method function
params.ndim = ndim;
params.vec = vec;
params.callback = alternative_callback_function;
output.hat = hat;
func_val = NewtonLib_newton_method((const struct NewtonLib_ArgParam_t *)¶ms, &output);
if ((int32_t)NEWTONLIB_SUCCESS != func_val) {
switch ((int)func_val) {
case NEWTONLIB_JACOBIAN:
PyErr_SetString(PyExc_RuntimeError, "[Newton's method]Failed to calculate Jacobian matrix.");
break;
case NEWTONLIB_OBJFUNCVAL:
PyErr_SetString(PyExc_RuntimeError, "[Newton's method]Failed to calculate objective function value.");
break;
case NEWTONLIB_SOLVE_SLE:
PyErr_SetString(PyExc_RuntimeError, "[Newton's method]Failed to solve simultaneous linear equations.");
break;
default:
break;
}
goto EXIT_NEWTONLIB;
}
//! Step7: store result
store_result(&py_output, ndim, hat, output.is_convergent);
EXIT_NEWTONLIB:
//! finalize
if (NULL != vec) {
free(vec);
vec = NULL;
}
if (NULL != hat) {
free(hat);
hat = NULL;
}
return py_output;
}上記のStep1～Step5までが前処理部分、Step6がコールバック関数の処理部分、Step7が後処理部分となります。
以降で、該当箇所の解説をしていきます。
newton_methodの前処理部分
Step1～Step5に対する解説を行っていきます。
Step1・Step2：引数処理と引数チェック
ここでは、キーワード引数を用いる場合の引数処理について解説します。
該当する処理のみを記載すると以下のようになります。
#define NEWTONLIB_DELTA (1e-3)
#define NEWTONLIB_TOL (1e-10)
// 中略
static PyObject *newtonlib(PyObject *self, PyObject *args, PyObject *keywords) {
PyObject *py_vec;
static char *kwlist[] = {"vec", "max_iter", "tol", "delta", NULL};
struct NewtonLib_ArgParam_t params;
//! initialize
py_vec = NULL;
params.tol = (double)NEWTONLIB_TOL;
params.delta = (double)NEWTONLIB_DELTA;
//! Step1: check arguments
if (!PyArg_ParseTupleAndKeywords(args, keywords, "Oi|dd", kwlist, &py_vec, ¶ms.max_iter, ¶ms.tol, ¶ms.delta)) {
PyErr_SetString(PyExc_TypeError, "Invalid argument");
goto EXIT_NEWTONLIB;
}
//! Step2: check object type
if (!PyArray_Check(py_vec)) {
PyErr_SetString(PyExc_TypeError, "'vec' needs to numpy array");
goto EXIT_NEWTONLIB;
}
// 以下省略
}キーワード引数を用いる場合、C言語側（呼び出し先）では受け取る際のキーワード一覧を用意しておく必要があります。
今回の場合、以下が該当します。
static char *kwlist[] = {"vec", "max_iter", "tol", "delta", NULL};上記の設定により、Python側（呼び出し元）では以下のような振る舞いを持つ関数として利用できます。
import numpy as np
from wrapper_newtonlib import newton_method
vec = np.array([1.0, -0.5])
max_iter = 100
tol = 1e-11
delta = 1e-4
# [Pattern1]error: invalid argument
newton_method()
# [Pattern2]error: invalid argument
newton_method(vec)
# [Pattern3]tol: 1e-10, delta: 1e-3
newton_method(vec, max_iter)
# [Pattern4]tol: 1e-11, delta: 1e-4
newton_method(vec, max_iter, tol=tol, delta=delta)
# [Pattern5]vec: np.array([-1.0, 1.0]), max_iter: 300, tol: 1e-10, delta: 1e-3
newton_method(vec=np.array([-1.0, 1.0]), max_iter=300)また、今回のケースでは、vecはNumpyオブジェクトのため、期待通りのオブジェクトになっているか、あわせて確認しています。
PyArray_Check(py_vec)Numpy/C APIは、以下の公式サイトでまとめられています。
https://numpy.org/doc/stable/reference/c-api/index.html
Step3～Step5：領域確保・データの読み込み
引数処理後は、Numpyオブジェクトからデータを取得し、C言語で確保したメモリ領域にコピーします。
データ領域を無駄に消費しますが、処理自体はC言語だけに閉じるため高速になります。
//! Step3: get vector size
np_arr = (PyArrayObject *)py_vec;
ndim = (int32_t)PyArray_SIZE(np_arr);
//! Step4: malloc
vec = (double *)malloc(sizeof(double) * ndim);
hat = (double *)malloc(sizeof(double) * ndim);
if ((NULL == vec) || (NULL == hat)) {
PyErr_SetString(PyExc_MemoryError, "lack of memory");
goto EXIT_NEWTONLIB;
}
//! Step5: copy data
for (idx = 0; idx < ndim; idx++) {
vec[idx] = *((double *)PyArray_GETPTR1(np_arr, idx));
}また、Numpyオブジェクトは、1次元配列を仮定しているため、PyArray_SIZEやPyArray_GETPTR1でNumpyオブジェクトから情報を取得しています。
newton_methodにおけるコールバック関数
C言語のデータをPythonで扱う為には、データ型の変換が必要になります。
今回は、Python側で定義された関数を呼び出す際に、データ型の変換を行う代替コールバック関数（alternative_callback_function）を定義しています。
alternative_callback_functionの実装は以下のようになります。
/**
* @fn static int32_t alternative_callback_function(int32_t ndim, const double *input, double *output);
* @brief call python function from c library
* @param[in] ndim number of dimensions
* @param[in] input input vector
* @param[out] output outpu vector
* @return NEWTONLIB_SUCCESS success
* NEWTONLIB_FAILED failed
*/
static int32_t alternative_callback_function(int32_t ndim, const double *input, double *output) {
int32_t ret = (int32_t)NEWTONLIB_FAILED;
int32_t idx;
PyObject *numpy_list, *py_data, *py_out;
PyArrayObject *np_arr = NULL;
npy_intp dims[1] = {0};
if (PyCallable_Check(py_objective_function)) {
//! set argument
dims[0] = ndim;
numpy_list = PyArray_SimpleNew(1, &dims[0], NPY_DOUBLE);
np_arr = (PyArrayObject *)numpy_list;
//! store input to numpy array
for (idx = 0; idx < ndim; idx++) {
py_data = PyFloat_FromDouble(input[idx]);
PyArray_SETITEM(np_arr, PyArray_GETPTR1(np_arr, idx), py_data);
Py_DECREF(py_data);
}
//! call function
py_out = PyObject_CallFunctionObjArgs(py_objective_function, np_arr, NULL);
Py_DECREF(numpy_list);
//! check result
if (py_out && PyArray_Check(py_out)) {
np_arr = (PyArrayObject *)py_out;
//! store output from numpy array
for (idx = 0; idx < ndim; idx++) {
output[idx] = *((double *)PyArray_GETPTR1(np_arr, idx));
}
Py_DECREF(py_out);
ret = (int32_t)NEWTONLIB_SUCCESS;
}
}
return ret;
}処理自体は、下記に示す通り、シンプルなものとなります。
- 与えられた次元の
np.arrayオブジェクトを作成する。 - C言語からの入力データをNumpyオブジェクトの該当インデックスにコピーする。
- あらかじめ登録しておいたPython側の関数（
py_objective_function）を引数とともに呼び出す。 - 戻り値を確認後、Numpyオブジェクトの内容をC言語の出力データの該当インデックスにコピーする。
C言語形式の入力データ
=> Numpy形式の入力データ
=> 関数呼び出し
=> Numpy形式の出力データ
=> C言語形式の出力データnewton_methodの後処理部分
推定値が得られたら、呼び出し元であるPythonにデータを返却できるように後処理を行います。
返却するデータを作成する部分は、以下のようになります。
/**
* @fn static void store_result(PyObject **py_output, int32_t ndim, double *hat, int32_t is_convergent)
* @brief convert output data of c format to python object
* @param[out] py_output output python object
* @param[in] ndim dimension
* @param[in] hat estimated vector of x
* @param[in] is_convergent convergent status
* @return none
*/
static void store_result(PyObject **py_output, int32_t ndim, double *hat, int32_t is_convergent) {
PyObject *numpy_list, *py_data, *py_convergent;
PyArrayObject *np_arr;
int32_t idx;
npy_intp dims[1] = {0};
dims[0] = ndim;
numpy_list = PyArray_SimpleNew(1, &dims[0], NPY_DOUBLE);
np_arr = (PyArrayObject *)numpy_list;
for (idx = 0; idx < ndim; idx++) {
py_data = PyFloat_FromDouble(hat[idx]);
PyArray_SETITEM(np_arr, PyArray_GETPTR1(np_arr, idx), py_data);
Py_DECREF(py_data);
}
py_convergent = ((int32_t)NEWTONLIB_CONVERGENT == is_convergent) ? Py_True : Py_False;
Py_INCREF(py_convergent);
PyList_Append(*py_output, numpy_list);
PyList_Append(*py_output, py_convergent);
}ここでも、参照カウントの考慮が必要になるため、適切なカウント処理を挟みます。
また、PythonのTrue、FalseをC言語でも扱えますが、これらも参照カウントの管理が必要になるため、次のような実装により参照カウントを増やしています。
py_convergent = ((int32_t)NEWTONLIB_CONVERGENT == is_convergent) ? Py_True : Py_False;
Py_INCREF(py_convergent);モジュール定義部分
最後に、PythonからC言語の関数を呼び出すための"つなぎ"の部分について解説します。
モジュール定義時は、以下の3ステップで実現します。
- モジュールに登録する関数の一覧を作成する。
- 作成するモジュール名や内包する関数一覧との対応関係など、モジュールに関する情報を定義する。
- モジュール呼び出し時の初期化処理を定義する。
モジュールに登録する関数一覧
C言語で定義した関数をモジュールに登録する際は、以下の形式の構造体を配列として定義します。
| 構造体メンバ | C言語での型定義 | 内容 |
|---|---|---|
ml_name | const char * | モジュールから呼び出す際のメソッド名 |
ml_meth | PyCFunction | C言語で定義した関数へのポインタ |
ml_flags | int | 関数の呼び出し方を示すフラグビット |
ml_doc | const char * | C言語で定義した関数に対するdocstringの内容を示すポインタ |
今回は2つの関数があるため、それぞれ構造体を定義します。
定義結果は以下のようになります。
//! register c functions to use python script
static PyMethodDef wrapper_methods[] = {
{
/**
* @defgroup PublicSetterAPI
* @brief setter module name of objective function
* @{
*/
"set_objective_function",
/** @} */ // end of PublicSetterAPI
//! target c function
update_function,
//! argument type
METH_VARARGS,
//! output message of __doc__ function
"Set objective function used by Newton's method",
},
{
/**
* @defgroup PublicWrapperAPI
* @brief newton method module name
* @{
*/
"newton_method",
/** @} */ // end of PublicWrapperAPI
//! target c function
(PyCFunction)newtonlib,
//! argument type
METH_VARARGS | METH_KEYWORDS,
//! output message of __doc__ function
"Calculate approximate solutions using Newton's method",
},
{NULL, NULL, 0, NULL}, // 番兵
};モジュールに関する情報の定義
以下の形式の構造体を用いて、Pythonで呼び出す際のモジュールを定義します。
| 構造体メンバ | C言語での型定義 | 内容 |
|---|---|---|
m_base | PyModuleDef_Base | 常にPyModuleDef_HEAD_INITで初期化する |
m_name | const char * | モジュール名 |
m_doc | const char * | このモジュールに対するdocstringの内容を示すポインタ |
m_size | Py_ssize_t | モジュールごと個別のメモリ領域を持たせる場合に非負の値を指定 今回は、コールバック関数をグローバル変数で管理するため、 -1を指定 |
m_methods | PyMethodDef * | 「モジュールに登録する関数一覧」で定義した変数を指定 |
今回は、Newton法のラッパー関数として定義するため、「wrapper_newtonlib」というモジュール名にします。
具体的な定義結果は以下のようになります。
//! defile python module
static struct PyModuleDef wrapper_newtonlib = {
PyModuleDef_HEAD_INIT,
"wrapper_newtonlib",
NULL,
-1,
wrapper_methods,
};モジュール呼び出し時の初期化処理
最後に、モジュール呼び出し時の処理を記述します。上記のモジュール定義をPyModule_Create関数の引数に渡し、呼び出せばOKです。
また、今回はNumpyを利用するため、import_array関数も実行しています。
/**
* @fn PyMODINIT_FUNC PyInit_wrapper_newtonlib(void)
* @brief PyInit function
*/
PyMODINIT_FUNC PyInit_wrapper_newtonlib(void) {
import_array();
return PyModule_Create(&wrapper_newtonlib);
}ライブラリの生成
C言語で実装したものをPythonから呼び出せるように共有ライブラリを生成する準備をします。
共有ライブラリは、setuptoolsを用いると簡単に生成できるため、setup.pyを作成後、gccによるコンパイルにより共有ライブラリを生成します。
今回のケースにおける、setup.pyの実装例を以下に示します。
from setuptools import setup, Extension
import numpy as np
setup(
name='newton_method',
version='0.1.0',
author='yuruto',
ext_modules=[
Extension(
name='wrapper_newtonlib',
sources=[
'libs/src/wrapper.c',
'libs/src/newton_method.c',
],
include_dirs=['./libs/include', np.get_include()],
extra_compile_args=['-O3'],
),
],
zip_safe=False,
)また、コンパイル時は下記のコマンドを実行します。
# compile.shの内容
python setup.py build_ext --inplacePythonからの呼び出し
一通りC言語での実装が完了したため、Pythonから共有ライブラリを読み込み、モジュールを利用する方法について説明します。
全体像は以下のようになります。
import sys
sys.path.append('.')
sys.path.append('..')
import numpy as np
import wrapper_newtonlib as newtonlib
def objective_function(vec):
"""
Parameters
----------
vec : numpy.ndarray
vector of x
vec[i]: x_{i}
Returns
-------
fv : numpy.ndarray
vector of f(vec)
"""
x, y = vec
fv = np.array([
x + 2 * y + 1.0,
x * x + 2 * y * y - 3.0,
])
return fv
if __name__ == '__main__':
max_iter = 1000
tol = 1e-10
# set objective function
newtonlib.set_objective_function(objective_function)
# ===============
# solve pattern 1
# ===============
print('Solve pattern 1')
vec = np.array([2.0, 1.0])
exact_vec = np.array([1.0, -1.0])
print('[init value] x: {:.5f}, y: {:.5f}'.format(*vec))
hat, _ = newtonlib.newton_method(vec, max_iter, tol=tol)
err = np.linalg.norm(hat - exact_vec)
print('[estimated] x: {:.5f}, y: {:.5f} ({:.5e})'.format(*hat, err))
print('')
# ===============
# solve pattern 2
# ===============
print('Solve pattern 2')
vec = np.array([-1.0, 1.0])
exact_vec = np.array([-5.0/3.0, 1.0/3.0])
print('[init value] x: {:.5f}, y: {:.5f}'.format(*vec))
hat, _ = newtonlib.newton_method(vec, max_iter, tol=tol)
err = np.linalg.norm(hat - exact_vec)
print('[estimated] x: {:.5f}, y: {:.5f} ({:.5e})'.format(*hat, err))このうち、今回の対応事項に関係する部分は以下のようになります。
# 作成したモジュールの読み込み
import wrapper_newtonlib as newtonlib
# 目的関数の定義
def objective_function(vec):
# 内容は省略
pass
# =====
# 中略
# =====
# 目的関数の設定
newtonlib.set_objective_function(objective_function)
# Newton法の実行
hat, _ = newtonlib.newton_method(vec, max_iter, tol=tol)
# 以下、省略上記のスクリプトを実行すると以下のような結果が得られます。
Solve pattern 1
[init value] x: 2.00000, y: 1.00000
[estimated] x: 1.00000, y: -1.00000 (0.00000e+00)
Solve pattern 2
[init value] x: -1.00000, y: 1.00000
[estimated] x: -1.66667, y: 0.33333 (2.28878e-16)Pythonのみで実装した結果と同じ結果が得られていることが分かります。
まとめ
今回は、Python/C APIを用いて、PythonからC言語の関数を呼び出す方法・C言語からPythonの関数を呼び出す方法について解説しました。
参照カウントの管理や所有権の所在等、注意する部分はありますが、C言語で実装することにより、処理速度の改善が期待できるので、処理時間を気にしている方は一度試してみてはどうでしょうか。
また、今回はガウスの消去法を自前で実装しましたが、BLASやLAPACKなどのライブラリを用いることでもっと手軽に処理を実装できます。
ガウスの消去法については、以下で解説しているので、興味がある方は参考にしてください。
-
-
【解説＆実装】ガウスの消去法
続きを見る