【解説】PythonからC言語の関数呼び出し・C言語からPythonの関数呼び出し

悩んでいる人

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))

それぞれの関数は、以下のような内容となります。

また、実行すると以下のような出力が得られます。

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関数」に振り分けると以下のようになります。

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ファイルを生成する必要があります。

ここでは、上記を踏まえたディレクトリ構成にしております。

具体的な全体構成は以下のようになっております。

./
|-- compile.sh
|-- newton-method-with-clib.py
|-- setup.py
`-- libs/
    |-- include/
    |   `-- wrapper_newton.h
    `-- src/
        |-- wrapper.c
        `-- newton_method.c

それぞれのファイルは、以下のような役割を持ちます。

以降では、以下に示すステップで解説をしていきます。

1. 必要なヘッダーの読み込み
2. 実装時の注意点
3. C言語から呼び出すPython関数の登録
4. C言語によるラッパー関数
5. newton_methodの前処理部分
6. newton_methodにおけるコールバック関数
7. newton_methodの後処理部分
8. モジュール定義部分（本体実装である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種類があります。

それぞれの違いは、以下のようになります。

また、一部のAPIには、参照を盗む（APIの内部でPy_DECREFが"勝手に"呼ばれるイメージ）ものも存在します。

ここでは、参照関係の理解のため、リストオブジェクトを題材に事例を取り上げていきたいと思います。

https://docs.python.org/ja/3/c-api/stable.html

リストオブジェクトによる事例

上記のリンクにアクセスすると、以下のAPIが存在することが確認できると思います。

それぞれ、以下のような関係になります。

また、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;
}

処理自体は、下記に示す通り、シンプルなものとなります。

1. 与えられた次元のnp.arrayオブジェクトを作成する。
2. C言語からの入力データをNumpyオブジェクトの該当インデックスにコピーする。
3. あらかじめ登録しておいたPython側の関数（py_objective_function）を引数とともに呼び出す。
4. 戻り値を確認後、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ステップで実現します。

1. モジュールに登録する関数の一覧を作成する。
2. 作成するモジュール名や内包する関数一覧との対応関係など、モジュールに関する情報を定義する。
3. モジュール呼び出し時の初期化処理を定義する。

モジュールに登録する関数一覧

C言語で定義した関数をモジュールに登録する際は、以下の形式の構造体を配列として定義します。

今回は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で呼び出す際のモジュールを定義します。

今回は、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 --inplace

Pythonからの呼び出し

一通り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などのライブラリを用いることでもっと手軽に処理を実装できます。

ガウスの消去法については、以下で解説しているので、興味がある方は参考にしてください。