Cythonで高速化#

Pythonで配列の要素に対して大量のループ計算を行う場合、Pythonの実行効率はC言語に比べて数百倍遅くなることがあります。まず、配列演算関数の例を通じて、Cythonがどのように配列の高速演算を実現するかを示します。

ベクトル集合の距離行列の計算#

一組のベクトルの中で、各ペアのベクトル間の距離を計算する関数を実装します。以下の配列Xの形状は(200, 3)で、200個の3次元空間中の点と見なすことができます:

import numpy as np

np.random.seed(42)
X = np.random.rand(200, 3)

各点間の距離を計算するために、NumPyのブロードキャスト機能を使用するか、SciPyのscipy.spatial.distance.pdist()を使用することができます。ただし、ここではPythonとCythonの性能を比較するために、三重ループを使用して各要素を個別に計算します:

def pairwise_dist_python(X):
    m, n = X.shape
    D = np.empty((m, m), dtype=np.float64)
    for i in range(m):
        for j in range(i, m):
            d = 0.0
            for k in range(n):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = D[j, i] = d**0.5
    return D

プログラムの結果が正しいことを確認するために、SciPyのpdist()を使用して同じ計算を行い、両者の計算速度を比較します。両者の計算速度は約300倍異なることがわかります。

from scipy.spatial.distance import pdist, squareform

%timeit squareform(pdist(X))
%timeit pairwise_dist_python(X)
np.allclose(squareform(pdist(X)), pairwise_dist_python(X))
204 μs ± 8.66 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
51 ms ± 834 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
True

%%cythonマジックコマンド#

Cythonプログラムはコンパイルが必要で、通常はsetup.pyプログラムを書く必要があります。これについては後で詳しく説明します。まず、Notebookの%%cythonマジックコマンドを使用してCythonプログラムを迅速にコンパイルして実行します。

%%cythonマジックコマンドはセルコマンドで、セル内のプログラム全体がCythonによって拡張モジュールにコンパイルされ、コンパイル後のモジュールからすべてのオブジェクトが自動的にロードされます。Cythonプログラムは独立したモジュールであるため、このモジュール内でnumpyライブラリを再ロードする必要があります。以下のプログラムはpairwise_dist_python()と完全に同じで、その計算速度もpairwise_dist_python()とほぼ同じで、わずかな改善しか見られません:

%load_ext helper.cython
%%cython
import numpy as np

def pairwise_dist_cython(X):
    m, n = X.shape
    D = np.empty((m, m), dtype=np.float64)
    for i in xrange(m):
        for j in xrange(i, m):
            d = 0.0
            for k in xrange(n):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = D[j, i] = d ** 0.5
    return D

以下でその計算速度をテストします:

%timeit pairwise_dist_cython(X)
np.allclose(pairwise_dist_cython(X), pairwise_dist_python(X))
40.3 ms ± 698 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
True

以下のプログラムでは、Cythonのすべての最適化手段を使用しています。cdefは変数の型を宣言するキーワードで、@cythonはCythonのコンパイルを指示するコマンドです。

%%cython
import numpy as np
import cython
from libc.math cimport sqrt

@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_dist_cython2(double[:, ::1] X):
    cdef int m, n, i, j, k
    cdef double tmp, d
    m, n = X.shape[0], X.shape[1]    
    cdef double[:, ::1] D = np.empty((m, m), dtype=np.float64)
    for i in range(m):
        for j in range(i, m):
            d = 0.0
            for k in range(n):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = D[j, i] = sqrt(d)
    return np.asarray(D)

Cythonでコンパイルされた関数は、SciPyのpdist()よりも少し速いです:

%timeit pairwise_dist_cython2(X)
np.allclose(pairwise_dist_cython2(X), pairwise_dist_python(X))
108 μs ± 1.61 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
True

Cythonで拡張モジュール作成#

NotebookでCythonプログラムの速度と正確性をテストした後、他のPythonプログラムから呼び出せるように拡張モジュールとしてコンパイルしたいと考えます。実際、%%cythonマジックコマンドは自動的にCythonコードを拡張モジュールにコンパイルし、そのモジュールからすべてのグローバルオブジェクトをロードします。以下では、sys.modulesを使用してpairwise_dist_cython2()関数が含まれるモジュールのファイルパスを見つけます。

import sys

sys.modules[pairwise_dist_cython2.__module__]
<module '_cython_magic_e995555757c69007b2c0463c54f193e7e1c94190' from 'C:\\Users\\ruoyu\\.ipython\\cython\\_cython_magic_e995555757c69007b2c0463c54f193e7e1c94190.cp313-win_amd64.pyd'>

Cythonプログラムを拡張モジュールとしてコンパイルするには、以下のようなsetup_fast_pdist.pyインストールスクリプトを書く必要があります。その中で、"fast_pdist"という名前の拡張モジュールを表すExtensionオブジェクトを定義します。これにはCythonソースプログラムfast_pdist.pyxが含まれており、このファイルの内容は前に定義したpairwise_dist_cython2()と同じです。NumPyの機能を使用しているため、numpy.get_include()を使用してNumPyのヘッダーファイルのパスを指定する必要があります。

%%writefile setup_fast_pdist.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np

ext_modules = [
    Extension("fast_pdist", ["fast_pdist.pyx"],
        include_dirs = [np.get_include()]),    
]

setup(
  name = 'a faster version of pdist',
  cmdclass = {'build_ext': build_ext},
  ext_modules = ext_modules
)
Overwriting setup_fast_pdist.py
%%writefile fast_pdist.pyx
# cython: language_level=3
import numpy as np
import cython
from libc.math cimport sqrt

@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_dist_cython2(double[:, ::1] X):
    cdef int m, n, i, j, k
    cdef double tmp, d
    m, n = X.shape[0], X.shape[1]    
    cdef double[:, ::1] D = np.empty((m, m), dtype=np.float64)
    for i in range(m):
        for j in range(i, m):
            d = 0.0
            for k in range(n):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = D[j, i] = sqrt(d)
    return np.asarray(D)
Overwriting fast_pdist.pyx

コマンドラインで以下のコマンドを実行すると、Cythonファイルがコンパイルされ、拡張モジュールfast_pdist.pydが生成されます。コマンドライン引数--inplaceは、生成された拡張モジュールがソースプログラムと同じパスに配置されることを意味します。

!python setup_fast_pdist.py build_ext --inplace
Compiling fast_pdist.pyx because it changed.
[1/1] Cythonizing fast_pdist.pyx

次に、この拡張モジュールをロードし、その中の関数を呼び出すことができます:

import fast_pdist

np.allclose(fast_pdist.pairwise_dist_cython2(X), pairwise_dist_python(X))
True

setup.pyを使用する以外に、cythonizeコマンドを使用することもできます。引数-iの意味は、前述の--inplaceと同じです。

Warning

Notebook環境で既にロードされている拡張モジュールを再コンパイルすると、ファイルが上書きできないエラーが発生することがあります。この場合、Notebookのカーネルプロセスを再起動して再コンパイルする必要があります。

!cythonize -i fast_pdist.pyx
error: could not delete 'fast_pdist.cp313-win_amd64.pyd': 拒\u7edd\u8bbf\u95ee。

Cythonの基本文法#

このセクションでは、Cythonの出力コードを確認しながら、基本的な文法について説明します。

C言語でのPythonオブジェクト型#

Cython言語で書かれたプログラムは、CythonコンパイラによってC言語プログラムにコンパイルされ、その後C言語コンパイラによってPythonから呼び出せる拡張モジュールにコンパイルされます。Cython言語では、Pythonのオブジェクト型とC言語の型を処理することができます。読者がCython言語をより深く理解するために、このセクションでは、C言語でPythonオブジェクトがどのように表現されているかを紹介します。

PythonはC言語で書かれており、すべてのオブジェクトはC言語の構造体で表現されます。どのようなオブジェクトの構造体でも、最初の2つのフィールドの意味は固定されています:

  • ob_refcnt:このオブジェクトの参照カウントで、参照カウントが0になると、この構造体が占めるメモリが解放されます。

  • ob_type:型オブジェクトへのポインタ。

PythonのC言語プログラムでは、PyObject構造体が定義されており、これには上記の2つのフィールドしかありません。他のオブジェクト型の構造体は、その後新しいフィールドを追加します。例えば、float型に対応する構造体は以下のようになります:

typedef struct {
    Py_ssize_t ob_refcnt;
    struct _typeobject *ob_type;
    double ob_fval;
} PyFloatObject;

64ビットシステムでは、Py_ssize_tとポインタは8バイトで表現され、double型の長さも8バイトです。したがって、Pythonのfloatオブジェクトは24バイトを占めます。以下では、sys.getsizeof()を使用して浮動小数点数オブジェクト1.0のバイト数を取得します:

import sys

sys.getsizeof(1.0)
24

すべてのPythonオブジェクトの最初の2つのフィールドの型は同じであるため、C言語ではPyObject *型のポインタを使用してPythonオブジェクトを表現します。Cythonでは、型が宣言されていないすべての変数はPythonオブジェクトを表すため、C言語にコンパイルされると、これらの変数はすべてPyObject *型のポインタになります。Pythonは、C言語で拡張モジュールを書くためのPython/C APIを提供しています。

以下では、%%cythonマジックコマンドに-a引数を追加して、Cythonプログラムとそのコンパイル後のC言語プログラムの関係を表示します。Notebookで以下のプログラムを実行すると、次のような折りたたみ可能なコードが出力されます。各行のCythonコードがコンパイルされたC言語コードの量を色で表示します。色が濃いほど、その行のコードの実行速度が遅くなる可能性があります。これらのC言語コードはすべてCythonによって自動生成されているため、読むのが難しいかもしれませんが、読者はこれらのコードを完全に理解する必要はなく、大まかな動作原理を理解するだけで十分です。

%%cython -a
a = 1.0
b = 2.0
c = a + b

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: a = 1.0
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_a, __pyx_float_1_0) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
/* … */
  __pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_4) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+2: b = 2.0
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_b, __pyx_float_2_0) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
+3: c = a + b
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_a); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_b); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_c, __pyx_t_4) < 0) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

上記のプログラムでは、変数abcはすべてPythonオブジェクトであるため、c = a + bは次のようなC言語コードにコンパイルされます。__Pyx_GOTREF__Pyx_DECREFは、ガベージコレクションが正常に動作することを保証するために、Pythonオブジェクトの参照カウンタを増減します。

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_a); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_b); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_c, __pyx_t_4) < 0) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

__pyx_n_s_aは文字列オブジェクト"a"です。__Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_a)は、グローバルディクショナリから"a"で表されるPythonオブジェクトを取得し、変数__pyx_t_1に代入します。PyNumber_Add()はPython/C API関数で、2つのPythonオブジェクトを数値として加算し、新しいPythonオブジェクト__pyx_t_3を返します。最後に、PyDict_SetItem()を呼び出して、__pyx_t_3をグローバルディクショナリ__pyx_dに追加し、対応するキーとして__pyx_n_s_cで表される文字列オブジェクト"c"を使用します。

したがって、単純な加算演算でも、グローバルディクショナリを2回検索して変数abで表されるオブジェクトを取得し、API関数を1回呼び出して加算演算を行い、結果をグローバルディクショナリに書き込む必要があります。

cdef を用いた変数型の宣言#

プログラムの実行速度を大幅に向上させるためには、Pythonプログラムで頻繁に使用される変数に対してcdefキーワードを使用して変数の型を宣言する必要があります。cdefキーワードを使用して変数の型を宣言すると、以下の2つの効果があります:

  • 変数のアクセス速度の向上:Pythonでは、グローバル変数やオブジェクトの属性とそれに対応する値の関係はディクショナリに保存されており、これらの変数や属性を読み書きするたびにディクショナリへのアクセスが必要です。cdefキーワードを使用して変数を宣言すると、その変数と値の対応関係はコンパイル時に決定されるため、ディクショナリの検索に必要な時間を節約できます。

  • 数値処理速度の向上:PythonではすべてのオブジェクトがC言語の構造体であり、それらを演算するためにはPythonが提供するAPI関数を呼び出す必要があります。値の型がわかっている場合、Cythonは可能な限りC言語が提供する演算機能を使用するため、計算速度が大幅に向上します。

以下の例では、cdefを使用してabcの3つの変数をdouble型として定義しています。生成されるC言語コードは以下のようになります:

static double a, b, c;
...
PyMODINIT_FUNC init(){
    ...
    a = 1.0;
    b = 2.0;
    c = a + b;
    ...
}

このコードは、3つのdouble型のグローバル変数を作成し、モジュールの初期化関数でこれらのグローバル変数に値を代入します。加算演算と変数の代入はC言語の操作であり、前節でAPI関数を呼び出すコードに比べて非常に簡潔です。

Warning

ここでcdefを使用して定義した3つのグローバル変数はC言語のグローバル変数であり、コンパイル後の拡張モジュールからPythonでこれらの値を取得することはできません。

%%cython -a
cdef double a = 1.0
cdef double b = 2.0
cdef double c = a + b

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: cdef double a = 1.0
  __pyx_v_54_cython_magic_a9eb7207565ddb950a621eff032ce32dee5975d5_a = 1.0;
/* … */
  __pyx_t_2 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+2: cdef double b = 2.0
  __pyx_v_54_cython_magic_a9eb7207565ddb950a621eff032ce32dee5975d5_b = 2.0;
+3: cdef double c = a + b
  __pyx_v_54_cython_magic_a9eb7207565ddb950a621eff032ce32dee5975d5_c = (__pyx_v_54_cython_magic_a9eb7207565ddb950a621eff032ce32dee5975d5_a + __pyx_v_54_cython_magic_a9eb7207565ddb950a621eff032ce32dee5975d5_b);

C言語の変数とPythonオブジェクトを演算する場合、C言語の変数をPythonオブジェクトに変換してからAPI関数を呼び出して演算を行います。以下の例では、sはC言語の変数で、aはPythonのfloatオブジェクトです。

saを直接加算すると、以下のようなC言語コードが実行されます。_v_sdouble型の変数で、_t_1から_t_3はPythonオブジェクトのポインタ、_s_aは文字列オブジェクト"a"です:

_t_1 = PyFloat_FromDouble(_v_s); //sをPythonオブジェクトに変換
__Pyx_GetModuleGlobalName(_t_2, _s_a); //グローバルディクショナリから"a"に対応するオブジェクトを取得
_t_3 = PyNumber_Add(_t_1, _t_2); //加算演算
_v_s = __pyx_PyFloat_AsDouble(_t_3); //加算演算の結果をdouble型に変換

<double>を使用してPythonオブジェクトをC言語のdouble型に強制変換するため、ここでの加算演算はC言語の加算演算子を使用します。以下のようなC言語コードが実行されます。_t_1はPythonオブジェクトのポインタ、_t_2double型です:

__Pyx_GetModuleGlobalName(_t_1, _s_a); 
_t_2 = __pyx_PyFloat_AsDouble(_t_1); 
_v_s = _v_s + _t_2;
%%cython -a
cdef double s = 0
a = 3.0
s = s + a #❶
s = s + <double>a #❷

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: cdef double s = 0
  __pyx_v_54_cython_magic_8c23bbe70ae22a6bfd9ef5ac40a39af6522c6c77_s = 0.0;
/* … */
  __pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_4) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+2: a = 3.0
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_a, __pyx_float_3_0) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
+3: s = s + a #❶
  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_54_cython_magic_8c23bbe70ae22a6bfd9ef5ac40a39af6522c6c77_s); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_a); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_5 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_5 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_54_cython_magic_8c23bbe70ae22a6bfd9ef5ac40a39af6522c6c77_s = __pyx_t_5;
+4: s = s + <double>a #❷
  __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_a); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_5 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_54_cython_magic_8c23bbe70ae22a6bfd9ef5ac40a39af6522c6c77_s = (__pyx_v_54_cython_magic_8c23bbe70ae22a6bfd9ef5ac40a39af6522c6c77_s + ((double)__pyx_t_5));

C言語のさまざまなデータ型に加えて、CythonはPythonの多くの組み込みデータ型、例えばリスト、辞書、タプルなどを認識できます。以下の例では、listを使用してclist変数をリストオブジェクトとして宣言し、int型のcindex変数をインデックスとしてリストの要素を取得します。比較のために、2つの動的変数pylistpyindexを使用して同じ操作を行います。

%%cython -a
cdef list clist = [1000, 2, 3]
cdef int cindex = 0
clist[cindex] #❶

pylist = [1000, 2, 3]
pyindex = 0
pylist[pyindex] #❷

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: cdef list clist = [1000, 2, 3]
  __pyx_t_2 = PyList_New(3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_INCREF(__pyx_int_1000);
  __Pyx_GIVEREF(__pyx_int_1000);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 0, __pyx_int_1000)) __PYX_ERR(0, 1, __pyx_L1_error);
  __Pyx_INCREF(__pyx_int_2);
  __Pyx_GIVEREF(__pyx_int_2);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 1, __pyx_int_2)) __PYX_ERR(0, 1, __pyx_L1_error);
  __Pyx_INCREF(__pyx_int_3);
  __Pyx_GIVEREF(__pyx_int_3);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 2, __pyx_int_3)) __PYX_ERR(0, 1, __pyx_L1_error);
  __Pyx_XGOTREF(__pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_clist);
  __Pyx_DECREF_SET(__pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_clist, ((PyObject*)__pyx_t_2));
  __Pyx_GIVEREF(__pyx_t_2);
  __pyx_t_2 = 0;
/* … */
  __pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_4) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+2: cdef int cindex = 0
  __pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_cindex = 0;
+3: clist[cindex] #❶
  if (unlikely(__pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_clist == Py_None)) {
    PyErr_SetString(PyExc_TypeError, "'NoneType' object is not subscriptable");
    __PYX_ERR(0, 3, __pyx_L1_error)
  }
  __pyx_t_2 = __Pyx_GetItemInt_List(__pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_clist, __pyx_v_54_cython_magic_30b42740991fe54adc6243f328268e0f7d234184_cindex, int, 1, __Pyx_PyInt_From_int, 1, 1, 1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
 4: 
+5: pylist = [1000, 2, 3]
  __pyx_t_2 = PyList_New(3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_INCREF(__pyx_int_1000);
  __Pyx_GIVEREF(__pyx_int_1000);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 0, __pyx_int_1000)) __PYX_ERR(0, 5, __pyx_L1_error);
  __Pyx_INCREF(__pyx_int_2);
  __Pyx_GIVEREF(__pyx_int_2);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 1, __pyx_int_2)) __PYX_ERR(0, 5, __pyx_L1_error);
  __Pyx_INCREF(__pyx_int_3);
  __Pyx_GIVEREF(__pyx_int_3);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_2, 2, __pyx_int_3)) __PYX_ERR(0, 5, __pyx_L1_error);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_pylist, __pyx_t_2) < 0) __PYX_ERR(0, 5, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+6: pyindex = 0
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_pyindex, __pyx_int_0) < 0) __PYX_ERR(0, 6, __pyx_L1_error)
+7: pylist[pyindex] #❷
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_pylist); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_pyindex); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 7, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

clistcindexの型が宣言されているため、clist[cindex]は以下のコードにコンパイルされ、リストオブジェクト内の特定のインデックスに対応する要素を取得するためにヘルパー関数__Pyx_GetItemInt_List()を直接呼び出します:

__Pyx_GetItemInt_List(clist, cindex, int, 1, __Pyx_PyInt_From_int, 1, 1, 1)

pylistpyindex変数には型宣言がないため、これらをPythonのオブジェクトとして扱います。pylist[pyindex]は以下のコードにコンパイルされます:

__Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_pylist);
__Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_pyindex);
__pyx_t_3 = __Pyx_PyObject_GetItem(__pyx_t_1, __pyx_t_2);

まず、__Pyx_GetModuleGlobalName()を呼び出してグローバル変数辞書からpylistpyindexに対応するオブジェクトを取得し、次に__Pyx_PyObject_GetItem()を呼び出してインデックスに対応する要素を取得します。

[]を使用したインデックスアクセス操作に加えて、Cythonは多くの内部型のメソッドと属性を認識できます。次のテーブルには、現在Cythonが認識できるPythonデータ型と関連する演算操作がリストされています。表内の既知の属性やメソッドに遭遇した場合、Cythonは効率的なコードを生成します。

データ型

認識可能な演算

bytesstrunicode

join()in

tuple

in

list

ininsert()reverse()append()extend()

dict

inget()has_key()keys()values()clear()copy()

set

inclear()add()pop()

slice

startstopstep

complex

cvalrealimag

前のプログラムにclist.append(pyindex)を追加して、append()がコンパイルされた後のプログラムを確認できます。未知の属性やメソッドに遭遇した場合、Pythonオブジェクトが提供する汎用APIインターフェースが呼び出されます。例えば、clist.count(0)などです。

def を用いた関数定義#

通常、プログラムの実行速度を向上させるために、Cythonでdefキーワードを使用して関数を定義し、Pythonからそれらを呼び出します。Pythonから呼び出す必要があるため、defで定義された関数の引数と戻り値はPythonオブジェクトです。しかし、Cythonでは関数の引数に型定義を追加することができ、Cythonはこれらの引数を型チェックし、対応するC言語の型に自動的に変換します。

以下のpy_square_add()関数の2つの引数はdouble型です。Pythonから呼び出す際には、2つのPythonのfloatオブジェクトが渡され、py_square_add()内部ではこれら2つのfloatオブジェクトがC言語のdouble型の変数に変換され、double型の結果を計算した後、それをPythonのfloatオブジェクトに変換して返します。

対応するC言語のコードは以下のようになります。ここで、xy_rはPythonのオブジェクトポインタ型であり、_v_x_v_yはC言語のdouble型変数です。API関数PyFloat_AsDouble()PyFloat_FromDouble()を使用して、floatオブジェクトとdouble型変数の間で変換を行います。

double _v_x = PyFloat_AsDouble(x);
double _v_y = PyFloat_AsDouble(y); 
PyObject * _r;
_r = PyFloat_FromDouble(((_v_x * _v_x) + (_v_y * _v_y)));
%%cython -a
def py_square_add(double x, double y):
    return x*x + y*y

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: def py_square_add(double x, double y):
/* Python wrapper */
static PyObject *__pyx_pw_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_1py_square_add(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static PyMethodDef __pyx_mdef_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_1py_square_add = {"py_square_add", (PyCFunction)(void*)(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_1py_square_add, __Pyx_METH_FASTCALL|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_1py_square_add(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  double __pyx_v_x;
  double __pyx_v_y;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("py_square_add (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_MACROS
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject **__pyx_pyargnames[] = {&__pyx_n_s_x,&__pyx_n_s_y,0};
  PyObject* values[2] = {0,0};
    if (__pyx_kwds) {
      Py_ssize_t kw_args;
      switch (__pyx_nargs) {
        case  2: values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = __Pyx_NumKwargs_FASTCALL(__pyx_kwds);
      switch (__pyx_nargs) {
        case  0:
        if (likely((values[0] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_x)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[0]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_y)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[1]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
        else {
          __Pyx_RaiseArgtupleInvalid("py_square_add", 1, 2, 2, 1); __PYX_ERR(0, 1, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        const Py_ssize_t kwd_pos_args = __pyx_nargs;
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values + 0, kwd_pos_args, "py_square_add") < 0)) __PYX_ERR(0, 1, __pyx_L3_error)
      }
    } else if (unlikely(__pyx_nargs != 2)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
      values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
    }
    __pyx_v_x = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
    __pyx_v_y = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_y == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("py_square_add", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 1, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_AddTraceback("_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1.py_square_add", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_py_square_add(__pyx_self, __pyx_v_x, __pyx_v_y);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_py_square_add(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x, double __pyx_v_y) {
  PyObject *__pyx_r = NULL;
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1.py_square_add", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(2, __pyx_n_s_x, __pyx_n_s_y); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_2 = __Pyx_CyFunction_New(&__pyx_mdef_54_cython_magic_6708ce801a2c4515da3cb93bf03e1b748ea510a1_1py_square_add, 0, __pyx_n_s_py_square_add, NULL, __pyx_n_s_cython_magic_6708ce801a2c4515da, __pyx_d, ((PyObject *)__pyx_codeobj__2)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_py_square_add, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+2:     return x*x + y*y
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = PyFloat_FromDouble(((__pyx_v_x * __pyx_v_x) + (__pyx_v_y * __pyx_v_y))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

Pythonの型(例えばlist)を使用して引数を宣言する場合、Cythonは型チェックを行い、認識可能な演算をより効率的なコードに変換します。以下の例では、alist引数をlist型として宣言しているため、len(alist)PyList_GET_SIZE(__pyx_v_alist)にコンパイルされ、alist[i]__Pyx_GetItemInt_List()関数の呼び出しにコンパイルされます:

%%cython -a
def sum_list(list alist): 
    cdef double s = 0
    cdef int i = 0
    for i in range(len(alist)):
        s += <double>alist[i]
    return s

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: def sum_list(list alist):
/* Python wrapper */
static PyObject *__pyx_pw_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_1sum_list(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static PyMethodDef __pyx_mdef_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_1sum_list = {"sum_list", (PyCFunction)(void*)(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_1sum_list, __Pyx_METH_FASTCALL|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_1sum_list(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_alist = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("sum_list (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_MACROS
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject **__pyx_pyargnames[] = {&__pyx_n_s_alist,0};
  PyObject* values[1] = {0};
    if (__pyx_kwds) {
      Py_ssize_t kw_args;
      switch (__pyx_nargs) {
        case  1: values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = __Pyx_NumKwargs_FASTCALL(__pyx_kwds);
      switch (__pyx_nargs) {
        case  0:
        if (likely((values[0] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_alist)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[0]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
        else goto __pyx_L5_argtuple_error;
      }
      if (unlikely(kw_args > 0)) {
        const Py_ssize_t kwd_pos_args = __pyx_nargs;
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values + 0, kwd_pos_args, "sum_list") < 0)) __PYX_ERR(0, 1, __pyx_L3_error)
      }
    } else if (unlikely(__pyx_nargs != 1)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
    }
    __pyx_v_alist = ((PyObject*)values[0]);
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("sum_list", 1, 1, 1, __pyx_nargs); __PYX_ERR(0, 1, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_AddTraceback("_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39.sum_list", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_alist), (&PyList_Type), 1, "alist", 1))) __PYX_ERR(0, 1, __pyx_L1_error)
  __pyx_r = __pyx_pf_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_sum_list(__pyx_self, __pyx_v_alist);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  goto __pyx_L0;
  __pyx_L1_error:;
  __pyx_r = NULL;
  __pyx_L0:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_sum_list(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_alist) {
  double __pyx_v_s;
  int __pyx_v_i;
  PyObject *__pyx_r = NULL;
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_AddTraceback("_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39.sum_list", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(3, __pyx_n_s_alist, __pyx_n_s_s, __pyx_n_s_i); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_2 = __Pyx_CyFunction_New(&__pyx_mdef_54_cython_magic_d4c804714e9354a89c70268af2fc7cd550bdae39_1sum_list, 0, __pyx_n_s_sum_list, NULL, __pyx_n_s_cython_magic_d4c804714e9354a89c, __pyx_d, ((PyObject *)__pyx_codeobj__2)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_sum_list, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+2:     cdef double s = 0
  __pyx_v_s = 0.0;
+3:     cdef int i = 0
  __pyx_v_i = 0;
+4:     for i in range(len(alist)):
  if (unlikely(__pyx_v_alist == Py_None)) {
    PyErr_SetString(PyExc_TypeError, "object of type 'NoneType' has no len()");
    __PYX_ERR(0, 4, __pyx_L1_error)
  }
  __pyx_t_1 = __Pyx_PyList_GET_SIZE(__pyx_v_alist); if (unlikely(__pyx_t_1 == ((Py_ssize_t)-1))) __PYX_ERR(0, 4, __pyx_L1_error)
  __pyx_t_2 = __pyx_t_1;
  for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) {
    __pyx_v_i = __pyx_t_3;
+5:         s += <double>alist[i]
    if (unlikely(__pyx_v_alist == Py_None)) {
      PyErr_SetString(PyExc_TypeError, "'NoneType' object is not subscriptable");
      __PYX_ERR(0, 5, __pyx_L1_error)
    }
    __pyx_t_4 = __Pyx_GetItemInt_List(__pyx_v_alist, __pyx_v_i, int, 1, __Pyx_PyInt_From_int, 1, 1, 1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 5, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_5 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_5 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 5, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_v_s = (__pyx_v_s + ((double)__pyx_t_5));
  }
+6:     return s
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_4 = PyFloat_FromDouble(__pyx_v_s); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 6, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_r = __pyx_t_4;
  __pyx_t_4 = 0;
  goto __pyx_L0;

cdef を用いたC関数定義#

def関数はPythonの呼び出しインターフェースを使用するため、Cythonプログラム内部でこれらの関数を呼び出す場合でも、引数と戻り値の型変換が必要です。大量のループでこのような関数を呼び出すと、大きなオーバーヘッドが発生します。cdefキーワードを使用して、Cythonプログラム内部でのみ呼び出せる関数を定義することができ、呼び出しのオーバーヘッドはC言語の関数と同じです。

以下のc_square_add()cdefを使用して定義され、引数と戻り値はdouble型です。a = c_square_add(1.0, 2.0)は以下のコードにコンパイルされます:

__pyx_v_a = __pyx_f_c_square_add(1.0, 2.0);

Tip

cdef関数の戻り値の型を宣言しない場合、その型はPythonオブジェクトになります。

%%cython -a
cdef double c_square_add(double x, double y):
    return x*x + y*y

cdef double a = c_square_add(1.0, 2.0)

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: cdef double c_square_add(double x, double y):
static double __pyx_f_54_cython_magic_3719867d67f48edc1093c858a4ef3d73694a22f6_c_square_add(double __pyx_v_x, double __pyx_v_y) {
  double __pyx_r;
/* … */
  /* function exit code */
  __pyx_L0:;
  return __pyx_r;
}
/* … */
  __pyx_t_3 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_3) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+2:     return x*x + y*y
  __pyx_r = ((__pyx_v_x * __pyx_v_x) + (__pyx_v_y * __pyx_v_y));
  goto __pyx_L0;
 3: 
+4: cdef double a = c_square_add(1.0, 2.0)
  __pyx_t_2 = __pyx_f_54_cython_magic_3719867d67f48edc1093c858a4ef3d73694a22f6_c_square_add(1.0, 2.0); if (unlikely(__pyx_t_2 == ((double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L1_error)
  __pyx_v_54_cython_magic_3719867d67f48edc1093c858a4ef3d73694a22f6_a = __pyx_t_2;

同じ関数をCython内で高速に呼び出し、かつPythonからも呼び出せるようにしたい場合、cpdefキーワードを使用できます。これにより、C言語関数とPythonから呼び出すためのラッパー関数の両方が生成されます。Cython内でcpdef関数を呼び出す場合、cdef関数よりも若干時間がかかりますが、def関数を呼び出すよりもはるかに高速です。

%%cython -a
cpdef double cp_square_add(double x, double y):
    return x*x + y*y

cp_square_add(1.0, 2.0)

Generated by Cython 3.0.12

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+1: cpdef double cp_square_add(double x, double y):
static PyObject *__pyx_pw_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static double __pyx_f_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_cp_square_add(double __pyx_v_x, double __pyx_v_y, CYTHON_UNUSED int __pyx_skip_dispatch) {
  double __pyx_r;
/* … */
  /* function exit code */
  __pyx_L0:;
  return __pyx_r;
}

/* Python wrapper */
static PyObject *__pyx_pw_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
static PyMethodDef __pyx_mdef_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add = {"cp_square_add", (PyCFunction)(void*)(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add, __Pyx_METH_FASTCALL|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  double __pyx_v_x;
  double __pyx_v_y;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cp_square_add (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_MACROS
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject **__pyx_pyargnames[] = {&__pyx_n_s_x,&__pyx_n_s_y,0};
  PyObject* values[2] = {0,0};
    if (__pyx_kwds) {
      Py_ssize_t kw_args;
      switch (__pyx_nargs) {
        case  2: values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = __Pyx_NumKwargs_FASTCALL(__pyx_kwds);
      switch (__pyx_nargs) {
        case  0:
        if (likely((values[0] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_x)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[0]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_y)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[1]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
        else {
          __Pyx_RaiseArgtupleInvalid("cp_square_add", 1, 2, 2, 1); __PYX_ERR(0, 1, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        const Py_ssize_t kwd_pos_args = __pyx_nargs;
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values + 0, kwd_pos_args, "cp_square_add") < 0)) __PYX_ERR(0, 1, __pyx_L3_error)
      }
    } else if (unlikely(__pyx_nargs != 2)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
      values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
    }
    __pyx_v_x = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
    __pyx_v_y = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_y == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L3_error)
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("cp_square_add", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 1, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_AddTraceback("_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c.cp_square_add", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_cp_square_add(__pyx_self, __pyx_v_x, __pyx_v_y);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_cp_square_add(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x, double __pyx_v_y) {
  PyObject *__pyx_r = NULL;
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = __pyx_f_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_cp_square_add(__pyx_v_x, __pyx_v_y, 0); if (unlikely(__pyx_t_1 == ((double)-1) && PyErr_Occurred())) __PYX_ERR(0, 1, __pyx_L1_error)
  __pyx_t_2 = PyFloat_FromDouble(__pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_r = __pyx_t_2;
  __pyx_t_2 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_AddTraceback("_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c.cp_square_add", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple_ = PyTuple_Pack(2, __pyx_n_s_x, __pyx_n_s_y); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);
/* … */
  __pyx_t_2 = __Pyx_CyFunction_New(&__pyx_mdef_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_1cp_square_add, 0, __pyx_n_s_cp_square_add, NULL, __pyx_n_s_cython_magic_e6d158981d6f99fd80, __pyx_d, ((PyObject *)__pyx_codeobj__2)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_cp_square_add, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* … */
  __pyx_t_2 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_2) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+2:     return x*x + y*y
  __pyx_r = ((__pyx_v_x * __pyx_v_x) + (__pyx_v_y * __pyx_v_y));
  goto __pyx_L0;
 3: 
+4: cp_square_add(1.0, 2.0)
  __pyx_t_3 = __pyx_f_54_cython_magic_e6d158981d6f99fd8029338235c32e9bcbcf298c_cp_square_add(1.0, 2.0, 0); if (unlikely(__pyx_t_3 == ((double)-1) && PyErr_Occurred())) __PYX_ERR(0, 4, __pyx_L1_error)

配列の効率的な処理#

科学計算プログラムでは、多くの配列操作が行われます。前節では、Cythonのメモリビュー(MemoryView)を使用して配列の要素に高速にアクセスする方法を簡単に紹介しました。本節では、その詳細な使い方を説明します。

Cythonのメモリビュー#

Cythonのメモリビューは以下の構文で宣言されます:

cdef 要素型[次元宣言] 変数名

次元宣言の形式は以下の通りで、:は軸を表し、::1は対応する軸上の要素が連続して格納されていることを示します:

  • [:]: 1次元配列

  • [::1]: 1次元連続配列

  • [:, :]:2次元配列

  • [:, ::1]:2次元配列、第1軸の要素が連続して格納

メモリビューはNumPy配列と同様に、shapebasestrideなどの情報を保持していますが、データストレージ自体は持たず、他のPythonオブジェクトやC言語の配列からデータストレージを取得します。

Cythonでは、メモリビューには2つの形式があり、使用法に応じてCythonが自動的に切り替えます。C言語レベルでは構造体ですが、Pythonオブジェクトとして使用する必要がある場合、Cythonは自動的にMemoryViewオブジェクトに変換します。以下に、メモリビューの使用例を示します。

以下のプログラムでは、❶bufはC言語のグローバル2次元配列で、❷viewはメモリビューであり、そのデータストレージをグローバル配列bufに初期化しています。

numpy.asarray()を使用して、MemoryViewオブジェクトをNumPy配列に変換し、配列のメソッドや関数を呼び出すことができます。❹また、MemoryViewオブジェクトは直接NumPyの関数に渡すこともできます。これらの関数内部では、asarray()のような関数が呼び出され、引数が配列に変換されます。

❺メモリビューをPython環境に返す場合、CythonはそれをMemoryViewオブジェクトに変換します。したがって、get_view()を使用して変換後のMemoryViewオブジェクトを取得し、その内容を確認することができます。

❻Cythonでメモリビューのshape属性を取得することは、メモリビュー構造体のshapeフィールドのデータを取得することに相当します。また、メモリビューの添字操作は、データストレージの対応するアドレスへの直接アクセスにコンパイルされるため、非常に効率的です。

%%cython
import numpy as np

cdef double buf[3][4]            #❶
cdef double[:, ::1] view = buf   #❷

def fill_value(double value):    
    np.asarray(view).fill(value) #❸

def sum_view():
    return np.sum(view)          #❹    
    
def get_view():
    return view                  #❺

def square_view():               #❻
    cdef int i, j
    for i in range(view.shape[0]):
        for j in range(view.shape[1]):
            view[i, j] *= view[i, j]

以下に、get_view()を使用してMemoryViewオブジェクトを取得し、その各種属性を確認します。これらはNumPy配列と全く同じです:

import helper.magics

view = get_view()
%C view.shape; view.strides; view.itemsize; view.nbytes
view.shape  view.strides  view.itemsize  view.nbytes
----------  ------------  -------------  -----------
(3, 4)      (32, 8)       8              96         

以下に、fill_view()square_view()sum_view()などの関数を使用して、C言語のグローバル配列を操作します。MemoryViewオブジェクトは添字演算もサポートしています:

fill_value(3.0)
print(sum_view())
square_view()
print(sum_view())
view[1, 2] = 10
print(sum_view())
36.0
108.0
109.0

以下に、MemoryViewオブジェクトをNumPy配列に変換します。これらはデータストレージを共有していることがわかります:

arr = np.asarray(view)
arr[1, 0] = 11
print(sum_view())
print(arr)
111.0
[[ 9.  9.  9.  9.]
 [11.  9. 10.  9.]
 [ 9.  9.  9.  9.]]

メモリビューのbase属性を使用して、実際のデータを保持するオブジェクトを取得できます。データはC言語のグローバル配列に保存されているため、Cythonはそのグローバル配列を表すarrayオブジェクトを作成し、その唯一のmemview属性を使用して新しいMemoryViewオブジェクトを取得できます。

print(view.base.__class__)
print(view.base.memview)
<class '_cython_magic_d99ab4ab18f2536ba0f0ee5152cd6ccca941a457.array'>
<MemoryView of 'array' object>

NumPy配列とメモリビュー#

以下に、メモリビューを使用してNumPy配列を操作する方法を示します:

double[:]を使用してxを1次元メモリビューとして宣言します。要素が連続して格納されていることを指定していないため、x[i]x.data + i * x.strides[0]にコンパイルされます。ここで、x.dataxのデータストレージの先頭アドレスで、単一バイトのポインタ型です。x.strides[0]は第0軸上の要素間のバイト間隔数です。

resの要素が連続して格納されていることを宣言しているため、res[i]((double *) res.data) + iにコンパイルされます。添字変数iと要素間のバイト間隔数の乗算が省略されるため、演算速度が向上します。

res.baseを使用してNumPy配列を返します。base属性はメモリビューのC言語構造体のフィールドではないため、CythonはまずそれをMemoryViewオブジェクトに変換し、Pythonのgetattr()関数を使用してbase属性を取得します。

関数内で配列の要素をループ処理する場合、boundscheckwraparoundの2つのコンパイルオプションを無効にすることができます。これにより、生成されるC言語コードは配列の添字の範囲チェックを行わず、負の添字もサポートしないため、配列要素へのアクセス速度が向上します。

%%cython
import numpy as np
import cython

@cython.boundscheck(False)
@cython.wraparound(False)
def square(double[:] x): #❶
    cdef int i
    cdef double[::1] res = np.empty_like(x)
    for i in range(x.shape[0]):
        res[i] = x[i] * x[i] #❷
    return res.base #❸

さらに、メモリビューはNumPy配列と同じスライス添字アクセス機能をサポートしています。メモリビューをスライスアクセスする場合、元のメモリビューとメモリを共有する新しい構造体が作成されるため、多少の演算オーバーヘッドがありますが、配列のスライス操作よりもはるかに高速です。

以下の例では、cpdefを使用して、Cython内で高速に呼び出せ、かつPythonからも呼び出せるnorm()関数を定義しています。これは1次元ベクトルをその場で正規化します。norm_axis()は2次元配列の指定された軸を正規化します。inplaceパラメータがTrueの場合、その場で正規化し、そうでない場合は新しい正規化された配列を返します。

❶メモリビューはスライス代入演算をサポートしており、この場合、一時的なメモリビュー構造体が作成され、2つのメモリビュー構造体間でデータがコピーされます。❷メモリビューをスライス添字で読み取る場合、新しいビューが作成され、このビューがnorm()関数に渡されます。

Tip

この例では、#cythonコメントを使用してboundscheckwraparoundコンパイルオプションを設定しています。

%%cython
#cython: boundscheck=False, wraparound=False
import numpy as np
import cython

cpdef norm(double[:] x):
    cdef double s
    cdef int i
    s = 0
    for i in range(x.shape[0]):
        s += x[i]*x[i]
    s = 1 / s**0.5
    for i in range(x.shape[0]):
        x[i] *= s 

def norm_axis(double[:, :] x, int axis=0, bint inplace=True):
    cdef int i
    cdef double[:, :] data
    if not inplace:
        data = np.empty_like(x)
        data[:] = x #❶
    else:
        data = x
        
    if axis == 1:
        for i in range(data.shape[0]):
            norm(data[i, :]) #❷
    elif axis == 0:
        for i in range(data.shape[1]):
            norm(data[:, i]) #❷ 
            
    return data.base

また、アドレス演算子を使用してメモリビューのデータのアドレスを取得し、それをC言語の関数にポインタとして渡すこともできます。以下の例では、❶cimportを使用してC言語の標準ヘッダーファイルstring.hからmemcpy()関数をインポートしています。この関数のシグネチャは以下の通りです:srcが指すアドレスを開始アドレスとする連続したnバイトを、dstが指すアドレスを開始アドレスとする領域にコピーします。

void *memcpy(void *dst, const void *src, size_t n);

❷アドレス演算子&を使用して、メモリビューの最初の要素を指すポインタを取得します。&dst[0]double *型のポインタを返します。memcpy()はバイト長を受け取るため、ここではsizeof(double)を使用して倍精度浮動小数点数のバイト数を計算し、メモリビューの第0軸の長さを掛けます。

ここで使用しているメモリビューは連続しているため、memcpy()を直接使用してメモリをコピーできます。連続していない場合、strides属性もC言語関数に渡す必要があり、これによりメモリビューの要素に正しくアクセスできます。

%%cython
from libc.string cimport memcpy #❶

def copy_memview(double[::1] src, double[::1] dst):
    memcpy(&dst[0], &src[0], sizeof(double)*dst.shape[0]) #❷
a = np.random.rand(10)
b = np.zeros_like(a)
copy_memview(a, b)
assert np.all(a == b)

また、型変換操作を使用してC言語のポインタをメモリビューに変換することもできます。例えば、addrがポインタの場合、<double[:10]>addrを使用して、長さ10の倍精度浮動小数点数のメモリビューに変換できます。長さは変数で指定することもできます。ただし、このような変換を行う際は、メモリの割り当てと解放に十分注意する必要があります。そうしないと、ダングリングポインタが発生し、プログラム全体がクラッシュする可能性があります。

PythonのオブジェクトとAPI#

Cythonは、Pythonの多くの組み込み型の一般的な操作を認識し、計算速度を向上させることができます。さらに、CythonでPython/C API関数を呼び出すことで、C言語でしか実現できない操作を実現することもできます。

listオブジェクトの操作#

Pythonでは、まず空のリストを作成し、append()メソッドを使用して要素を追加するか、[None]*nを使用してn個の要素を持つリストを作成し、インデックスを使用してリストの内容を設定することができます。Cythonでは、リストを操作するAPI関数を呼び出すことで、リストの作成速度を向上させることができます。

以下のコードは、API関数とappend()メソッドを使用してリストを作成する速度を比較しています。my_range()では、以下の3つのAPI関数を呼び出します:

  • object PyList_New(Py_ssize_t len): 指定されたサイズのリストを作成し、リストの要素をNULLに設定します。ここでのNULLはC言語の空アドレスであり、PythonのNoneオブジェクトではありません。つまり、このリストは作成されていますが、その内容は初期化されておらず、Pythonでは使用できません。

  • void PyList_SET_ITEM(object list, Py_ssize_t i, object o): リストlistの指定されたインデックスiの内容をoに設定します。これは実際にはC言語のマクロであり、要素がNULLのリストの要素を迅速に設定するために使用できます。

  • void Py_INCREF(object o): オブジェクトoの参照カウンタを1増やします。PyList_SET_ITEM()を呼び出してオブジェクトoをリストに追加する場合、リストはオブジェクトoの参照カウンタを増やす必要がありますが、PyList_SET_ITEM()は自動的に追加されたオブジェクトの参照カウンタを増やしません。そのため、Py_INCREF()を呼び出す必要があります。このような参照カウンタを増やさない関数は、Python/C APIのドキュメントで「参照を盗む」と注記されています。

%%cython --compile-args=-w
#cython: boundscheck=False, wraparound=False
from cpython.list cimport PyList_New, PyList_SET_ITEM #❶
from cpython.ref cimport Py_INCREF

def my_range(int n):
    cdef int i
    cdef object obj #❷
    cdef list result
    result = PyList_New(n)
    for i in range(n):
        obj = i
        PyList_SET_ITEM(result, i, obj)
        Py_INCREF(obj)
    return result

def my_range2(int n):
    cdef int i
    cdef list result
    result = []
    for i in range(n):
        result.append(i)
    return result

cimportを使用して、CythonのヘッダーファイルからPyList_NewPyList_SET_ITEMの2つの関数宣言を読み込みます。Cythonのヘッダーファイルは.pxd拡張子を持ち、C言語のヘッダーファイルと同様に、さまざまな関数と型の定義を含むことができます。

CythonのインストールディレクトリのIncludesディレクトリに、すべてのCythonヘッダーファイルがあります。numpyディレクトリにはNumPy関連の型とAPI関数の宣言が含まれ、libcディレクトリにはC言語標準ライブラリの関数の宣言が含まれ、cpythonディレクトリにはPython/C APIの関数宣言が含まれています。

obj変数は、C言語の整数変数iをPythonの整数オブジェクトに変換するために一時的に保存されます。

my_range()は直接目標サイズのリストを作成するため、段階的な拡張によるオーバーヘッドを省き、append()を使用する場合の1.5倍以上の速度を実現します。

%timeit list(range(100))
%timeit my_range(100)
%timeit my_range2(100)
857 ns ± 28.8 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
715 ns ± 4.23 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
1.46 μs ± 32 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

tupleオブジェクトの作成#

Pythonではtupleオブジェクトは不変ですが、CythonでAPI関数を呼び出してtupleオブジェクトを作成する際にその内容を設定することで、高速にtupleオブジェクトを作成することができます。以下のto_tuple_list()は、API関数を呼び出して2次元配列をタプルのリストに変換します。

リストと同様に、タプルにも対応する初期化関数があります:PyTuple_NewPyTuple_SET_ITEMで、その使用方法はリストと同じですので、ここでは繰り返しません。

%%cython --compile-args=-w
#cython: boundscheck=False, wraparound=False
from cpython.list cimport PyList_New, PyList_SET_ITEM
from cpython.tuple cimport PyTuple_New, PyTuple_SET_ITEM
from cpython.ref cimport Py_INCREF

def to_tuple_list(double[:, :] arr):
    cdef int m, n
    cdef int i, j
    cdef list result
    cdef tuple t
    cdef object obj
    
    m, n = arr.shape[0], arr.shape[1]
    result = PyList_New(m)
    for i in range(m):
        t = PyTuple_New(n)
        for j in range(n):
            obj = arr[i, j]
            PyTuple_SET_ITEM(t, j, obj)
            Py_INCREF(obj)
        PyList_SET_ITEM(result, i, t)
        Py_INCREF(t)
    return result

以下は、NumPy配列のtolist()メソッドとto_tuple_list()の速度を比較します:

import numpy as np

arr = np.random.randint(0, 10, (5, 2)).astype(np.double)
print(to_tuple_list(arr))

arr = np.random.rand(100, 5)
%timeit to_tuple_list(arr)
%timeit arr.tolist()
[(8.0, 9.0), (8.0, 8.0), (5.0, 7.0), (0.0, 9.0), (3.0, 0.0)]
14.4 μs ± 399 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
17 μs ± 116 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

array.arrayを動的配列として使用#

NumPyの章で紹介したように、Pythonの標準モジュールarrayarrayオブジェクトを1次元の動的配列として使用することができます。Cythonコードでも、array.arrayオブジェクトを動的配列として使用することができます。ただし、Cythonコードでarray.append()を呼び出して要素を追加する場合、Pythonの関数オブジェクトを呼び出す必要があるため、速度向上の効果はありません。Cythonのcpython/array.pxdヘッダーファイルには、arrayオブジェクトを迅速に操作するためのcdef関数が提供されています。

以下のin_circle()は、2次元座標配列pointsの中で、cxcyrで表される円の内部にあるすべての座標を収集します。ここで、(cx, cy)は円の中心座標、rは円の半径です。円の内部にある点の数が事前にわからないため、プログラムではarray.array動的配列を使用して条件を満たす点を逐次追加します。

%%cython -c-Ofast --compile-args=-w
#cython: boundscheck=False, wraparound=False
import numpy as np
from cpython cimport array

def in_circle(double[:, :] points, double cx, double cy, double r):
    cdef array.array[double] res = array.array("d") #❶
    cdef double r2 = r * r
    cdef double p[2] #❷
    cdef int i 
    for i in range(points.shape[0]):
        p[0] = points[i, 0]
        p[1] = points[i, 1]
        if (p[0] - cx)**2 + (p[1] - cy)**2 < r2:
            array.extend_buffer(res, <char*>p, 2) #❸
    return np.frombuffer(res, np.double).copy().reshape(-1, 2) #❹

cimportキーワードを使用してarrayのヘッダーファイルを読み込んだ後、その中のarray.array型を使用してCython変数resを定義し、新しいarray.arrayオブジェクトを作成してそれに割り当てます。"d"は要素の型が倍精度浮動小数点数であることを示します。❷pは2つの要素を持つC言語の配列で、現在処理中の点の座標を一時的に保存するために使用します。❸array.extend_buffer()を呼び出してpresに追加します。extend_buffer()の最初の引数はarray.arrayオブジェクトで、2番目の引数は追加するデータの先頭アドレスを指すchar*ポインタ、3番目の引数は追加する要素の数(バイト数ではない)です。❹最後に、numpy.frombuffer()を使用してresとメモリを共有するNumPy配列を作成し、Pythonがresオブジェクトをガベージコレクトできるように配列をコピーします。

以下は、in_circle()とNumPy関連メソッドの計算速度を比較します。ほとんどの点が円の外側にある場合、in_circle()の計算速度はより高速になります。

points = np.random.rand(10000, 2)
cx, cy, r = 0.3, 0.5, 0.05

%timeit points[(points[:, 0] - cx)**2 + (points[:, 1] - cy)**2 < r**2, :]
%timeit in_circle(points, cx, cy, r)
82.9 μs ± 1.4 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
25.5 μs ± 690 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)