データサイエンティスト(?)の青田です。これは Nikkei Advent Calendar 2021 の25日目の記事です。

はじめに

本記事ではCythonを用いるときに発生しがちなつまずきポイントについて実用的な具体例を示す。つまずきポイントをほぼコピペで乗り越えられることを意識して執筆した。

すでに数多くのブログや公式ドキュメントがある中でこれを書いたモチベーションがある。 Cythonはコンパイルしようとするだけでも4種類のやり方が存在し、型宣言の仕方は3種類の作法があり、numpyとの連携方法は2種類存在する。

このように同じことをやろうとしたときの選択肢の多さが混乱を招いているように感じた。いろんな流派が存在するものの、ここでは自分の方法を目的別に示す。これにより、利用者の選択の時間を削減し、Cythonを道具として使いやすくなるだろう。

本記事は以下のトピックについての例を示す。

• pyxファイルのコンパイル方法
  • 高速化のおまじないを添えて
• 高速化されているかの確認方法
• CythonでC++のstdの使用例
  • ここではstd::vectorを用いてランレングス符号化を実装する
• Cythonとnumpyの連携
  • cimport numpy を用いる方法
  • Typed Memoryviews を用いる方法
• Cythonを含んだmoduleのinstall例
  • 自作packageをpip installできるようにする

なお、他の目的についてはすでにある文献に任せる。

• jupyter上でCythonを書く
  • https://qiita.com/pashango2/items/45cb85390193d97523ca#jupyteripython%E3%81%A7%E5%A7%8B%E3%82%81%E3%82%8Bcython
• 宣言の仕方やクラスの書き方
  • https://qiita.com/pashango2/items/fed0a1c06c7077639edc
• Cythonの並列化
  • https://www.isus.jp/products/python-distribution/thread-parallelism-in-cython/
  • 個人的には多くの場合において、Cythonで定義される関数は1スレッドでの処理にしてmultiprocessingで並列処理を書くほうが良いと思う

また前提知識は以下になる。

• コマンドラインでlsやcd等の基本的な操作ができること
• Pythonの実行環境や各種ライブラリのinstallを自力でできること
• Pythonの基本的文法や型の概念を理解していること (tutorialの6章まで終わっていれば本記事では十分です)

本記事執筆時に用いたPythonのversionは3.9.6であり、Cythonのversionは0.29.24である。

pyxファイルのコンパイル方法

まず多くの人がつまずくポイントがコンパイル方法だと思う。筆者のおすすめの方法はcythonizeコマンドを用いた方法である。

具体例

細かい理屈は抜きにまずはとにかくやってみよう。次の例にしたがって、引数を受け取り何もせずに返す関数をCythonで実装してみよう。

まず以下の内容のファイルを作成する。名前はとりあえず、cythonized_functions.pyxとでもしておこう。

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

cpdef identity_map(x):
    return x

そして次にこれをコンパイルしてみる。cythonized_functions.pyxのディレクトリで次のコマンドを打つと直ちにコンパイルされる。

cythonize -3 -a -i cythonized_functions.pyx

おそらく、同一ディレクトリの中にcythonized_functions.cpython-39-x86_64-linux-gnu.soという隠しファイルが存在するだろう。環境によってファイル名は若干異なるもののこれがコンパイルされた関数の本体である。

先程定義した関数はimportできるようになる。同一リポジトリ内だったら以下のように呼び出せる。

from cythonized_functions import identity_map  # import
identity_map('pen' + 'pineapple' + 'apple' + 'pen')  # 呼び出し
# -> 'penpineappleapplepen' #このような結果が得られるはずです。

もちろん、これだけではCythonの恩恵を得ることはできないが、どうやってコンパイルしてPythonから呼び出すかの流れは把握できただろう。

解説

解説パートでは「話が細かい！」と思われる方向けに読み飛ばしても後に支障は出ないように心がけて執筆した。

まずは、pyxファイルの冒頭に記述したおまじないから解説する。

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

1行目は、Cythonが変換する言語の指定をしている。Cythonはpyxで定義された関数を一度CかC++に変換してからコンパイルしている。デフォルトはCであり、ここではわざわざC++を指定している。なぜこんなことをするかというと、C++で用意されている便利なデータ構造を用いたいときがあるからだ。詳しくは、後の"CythonでC++のstdの使用例"で見ていく。

2行目は、高速化のためのコンパイル最適化オプションである。

3行目が一番高速化に重要な行である。language_level=3はpyxファイルがPython3系の文法に基づいたファイルであることを示している。boundscheck=Falseは範囲外参照チェックを無効にする。配列外の要素を指定したときにもerrorにならないので、用いるのは諸刃の剣だが高速化にはかなり効く。wraparound=FalseはPythonでやりがちな負のindex指定を無効にする。

4行目は割り算の挙動を指定している。これもC側の処理に合わせることで、Python側の諸々のチェックを無効にし高速化が達成されると理解すればOKだろう。割り算を含む関数では特に重宝するだろう。

次に関数を定義するときのcpdefだが、これはCython独特の定義の仕方である。他にcdefもdefもあるが、とりあえず、cpdefで定義しておけば問題はなさそう。defで定義すると高速化が効かず、cdefで定義するとPython内から呼び出せないという違いがある。これより細かい話は公式ドキュメントに譲ることにする。

そしてコンパイル時の

cythonize -3 -a -i cythonized_functions.pyx

のオプションについては、前からPython3系、高速化状況のhtml生成、その場でコンパイル を意味している。 特に"高速化状況のhtml生成"については次の章で見ていくファイルとなる。

高速化されているかの確認方法

「Cythonを使っても早くならなかった」そういう声をしばしば聞く。ほとんどの場合、Cythonをうまく書けていないだけであり、原因はPython Objectを取り扱うときの型チェック等のオーバーヘッドであることが多い。しかし、なまじ動いてしまうだけにどこが悪いのがわからず途方にくれてしまう人が多いだろう。Cythonで定義した関数のどこがコンパイルされて、どこがPythonのままなのかを判断するのは経験がないと難しい。

実は、先のコンパイルコマンドで生成されるhtmlファイルを見れば、改善点がわかるのである。ここでは、いくつか関数を実装してみて、htmlファイルにどう出力されるか観察してみる。

具体例

まずはcythonized_functions.pyxとして以下の内容のファイルを作成して、コンパイルコマンドcythonize -3 -a -i cythonized_functions.pyxでコンパイルしてみよう。

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

cpdef identity_map(x):
    return x

cpdef int identity_map_typed(int x):
    return x

cdef int identity_map_typed_cdef(int x):
    return x

cpdef int sum_to(int x):
    num_list = list(range(x + 1))
    ret = sum(num_list)
    return ret

cpdef int sum_to_for(int x):
    cdef int i
    cdef int ret = 0
    for i in range(x + 1):
        ret += i
    return ret

すると以下のようなcythonized_functions.htmlが生成されるだろう。

この図は、Pythonが呼び出されている行ほど黄色くハイライトされている。つまり黄色い行が少ない関数ほど高速に処理される。

たとえば3番目の関数cdef int identity_map_typed_cdef(int x)は一切黄色い行がなく、この関数のすべてがPythonを経由することなく実行されることを示す。ただこれはcdefで定義しているため、Pythonからは呼び出せないことに注意だ。

基本的には最後のcpdef int sum_to_for(int x)の様に、関数の1行目以外は黄色くハイライトされないのが理想だ。関数の1行目は、Pythonとデータを受け渡しするためのコードが入るのでcpdefを使っている限り黄色くハイライトされるだろう。

解説

ここではもっと細かく、先の例の関数5つについて解説していく。

まずは始めの3つ、cpdef identity_map(x), cpdef int identity_map_typed(int x), cdef int identity_map_typed_cdef(int x)を比較する。 これらはどれも受け取った引数をそのまま返す関数だが、型付けが異なる。 最初のcpdef identity_map(x)は型付けを全くしていないため、受け取るときも返すときもPythonによる型チェックが入り関数呼び出しのオーバーヘッドが大きい。

次に残りの2つ、cpdef int sum_to(int x), cpdef int sum_to_for(int x)を比較する。これらは0からxまでの総和を取る関数である。 前者はPythonの作法でlistやsumなどを用いて実装していて、後者は愚直にfor文で実装している。Cythonにおいては後者の様に愚直に実装するほうが良い。 前者ではrange(x + 1)だけで済むところをlist(range(x + 1))と冗長に書いている。この例ではlistを用いなくても良いだろうが、一般的にはlistなどに(連番でない)数列等を格納して処理することが多いと思い、例としてlistで囲った。重要なポイントはlistなどを用いるとPython Objectを扱うことになりPythonによる諸々のチェックが混入する点である。またsumというのもPythonの組み込み関数であるため遅くなる。そのため、後者の様に愚直にfor文で回したほうが早いのである。もちろん、各変数には適切に型をつけることが必要である。

速度を向上させるためとはいえ、listやdict等の便利なデータ構造を用いることができないのは非常に厄介な制約だと感じるだろう。 これに関してはC++のvectorやunordered_map等を用いることにより解決可能である。詳しくは次章で見ていこう。

CythonでC++のstdの使用例

前章で述べたようにCythonでlistやdictといったようなPython Objectを扱うとパフォーマンスが下がる。 ここではC++のstdを用いることで、パフォーマンスを落とすことなくPythonの様に柔軟に実装できる具体例を示す。

具体例

ランレングス符号化をPythonとCythonを用いて実装してみる。

ランレングス符号化とは例えば以下のような入出力に対応する処理である。日本語だと連長圧縮といったほうが理解しやすいかもしれない。 競技プログラミングでよく使う関数から例を取ってきたので一般的な出力でないかもしれないが、この形式でお付き合い願いたい。

# 入力
[1, 1, 0, 0, 0, 2, 1, 1, 1, 1, 1] #例えば長さ10の配列

# 出力
([1, 0, 2, 1], # 連続する要素を1つの要素に
 [2, 3, 1, 5], # 各要素の連続する回数 #もともと長さ10の配列なので、この配列も合計は10になる
 [0, 2, 5, 6]) # 各要素の入力配列における開始index

これをPythonで実装すると以下のようになる。

def run_length_encoding(s: iter):
    '''
    ランレングス符号化(連長圧縮)を行う
    s (iterable object) ... 圧縮したい数列

    return
    ----------
    s_composed, s_num, s_idx
     それぞれ、圧縮後の配列、連続する個数、その要素が始まるidx
    '''
    s_composed = []
    s_num = []
    s_idx = [0]
    pre = s[0]
    cnt = 1
    for i, ss in enumerate(s[1:], start=1):
        if pre == ss:
            cnt += 1
        else:
            s_num.append(cnt)
            s_composed.append(pre)
            s_idx.append(i)
            cnt = 1
            pre = ss
    s_num.append(cnt)
    s_composed.append(pre)
    return s_composed, s_num, s_idx

一方Cythonで実装すると以下のようになる。ただしlistの代わりにstd::vectorを用いる。

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

from libcpp.vector cimport vector #ここでcppのvectorを呼び出している。

# 以下は型のエイリアス的なもの。
ctypedef long long LL
ctypedef vector[LL] vec

cpdef vector[vec] run_length_encoding_cython(LL[:] s): #次章で解説するが、numpy配列前提の引数。また型が大事。
    '''docstringは省略'''
    cdef:
        vec s_composed, s_num, s_idx
        LL pre, cnt, i, ss
        vector[vec] ret

    s_idx.push_back(0) #これがlist.appendに相当する操作
    pre = s[0]
    cnt = 1

    for i in range(1, s.shape[0]): 
        ss = s[i]
        if pre == ss:
            cnt += 1
        else:
            s_num.push_back(cnt)
            s_composed.push_back(pre)
            s_idx.push_back(i)
            cnt = 1
            pre = ss
    s_num.push_back(cnt)
    s_composed.push_back(pre)
    # 以下は戻り値用の操作 #Pythonの様に return s_composed, s_num, s_idx でも動くがここでは型付けに準拠する
    ret.push_back(s_composed)
    ret.push_back(s_num)
    ret.push_back(s_idx)
    return ret

PythonとCythonで定義した関数の実行時間の違いを測ってみよう。

import numpy as np
rand = np.random.randint(0, 10, size=10**7) #乱数列を生成

%timeit run_length_encoding(rand)
# -> 1.96 s ± 108 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit run_length_encoding_cython(rand)
# -> 871 ms ± 11.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

というわけで、Cythonに書き直すことで、実行時間が半分以下になった。今回はCythonで定義した関数を1回しか呼び出さなかったが、関数を複数回呼び出すときはもっと高速化される。Pythonは関数を呼び出すときのオーバーヘッドも大きいためである。場合によっては数十~百倍高速化することも珍しくない。

解説

冒頭のおまじないは上記で見たとおりである。

# distutils: language=c++

次に、

from libcpp.vector cimport vector #ここでcppのvectorを呼び出している。

の行でC++のstd::vectorを使えるようにimportしている。cimportというのはCython独特のimportである。

そして以下の部分でlistに相当する変数を宣言している。

cdef:
    vec s_composed, s_num, s_idx
    LL pre, cnt, i, ss
    vector[vec] ret

s_composed, s_num, s_idxはvector[long long]という型になるのだが、いちいちこの長い型を書いているのはめんどくさいので、冒頭に

ctypedef long long LL
ctypedef vector[LL] vec

を書いて、vecと書くだけで型宣言できるようにしている。

空のvecに要素を追加するには.push_backメソッドを呼び出す。ここらへんはC++の作法になる。Pythonみたいに.appendとしてもエラーになるので気をつけよう。

s_idx.push_back(0)

そして最後の戻り値は

ret.push_back(s_composed)
ret.push_back(s_num)
ret.push_back(s_idx)
return ret

というふうにretにまとめている。型宣言(vector[vec])からわかるようにこれはvectorの入れ子である。 これはC++が複数の戻り値を取れないという制約があるためこうしている(できないことはないのですが...)。vectorはPythonの世界に戻るときにlistに自動変換されるので特に我々が書く際に意識する必要はない。

Cythonとnumpyの連携

Cythonを使う場合は基本的にnumpyと連携することが前提となる。とくに引数が配列の場合はそうだ。 一方でnumpyとの連携の仕方は大きく2種類存在している。

1. cimport numpyを用いた方法
2. Typed Memoryviewsを用いた方法

どちらも一長一短があるが、推奨されている方法は後者である。ここで前者を紹介したのは前者のほうが楽な場合も多いからだ。また前者と後者の方法は共存できる。

ここでは、まず理解をなしに、1,2のそれぞれについて具体例を示す。細かいことを気にしなければ2の方法をとりあえず真似することをおすすめする。

具体例

配列sを受け取り合計する関数を1,2の方法で実装してみる。

まずは方法1から

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# distutils: include_dirs= ["/opt/conda/lib/python3.9/site-packages/numpy/core/include"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

cimport numpy as cnp
import numpy as np

ctypedef long long LL

cpdef LL sum_cython_numpy(cnp.ndarray[cnp.int64_t, ndim=1] s):
    cdef LL ret
    ret = s.sum()
    return ret

上記の例は3行目でnumpyにpathを通しているがここは人によって違うのでコンパイルエラーが出るかもしれない。もしそうなった場合は、次章の"Cythonを含んだmoduleのinstall例"を参考にしてほしい。大事なのは関数の引数の型をcnp.ndarray[cnp.int64_t, ndim=1]と指定している点である。 型を指定されたsはnp.ndarrayであり、numpyのメソッドを使って短く書けるのがこの方法のメリットである。

次に方法2の例を示す。

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

ctypedef long long LL

cpdef LL sum_typed_memoryviews(LL[:] s):
    cdef LL ret = 0
    cdef LL l = s.shape[0]
    cdef LL i
    for i in range(l):
        ret += s[i]
    return ret

大事な点は引数の型にLL[:]を指定している点である。Cythonはnumpy配列のメモリを読み取り、自動的にTyped Memoryviewsとして扱う。これはnp.ndarrayと同じメモリを参照しているが、numpyのメソッドを用いることはできない点に注意だ。

方法1も2も同じ実行結果である。速度的には後者のほうが早いこともあって後者の方法が推奨されている。 しかしコードの長さからわかるように、cimport numpyを用いたほうがnumpyのメソッドを呼び出せるため記述が簡単になり、バグを埋め込みにくいという利点もある。よって筆者は1,2の方法を使い分けたり、場合によっては同時に使ったりすることが重要だと考える。

解説

方法1をする際に、Cythonでnumpyをcimportする必要がある。そのためにまずpathを通す必要がある。

# distutils: include_dirs= ["/opt/conda/lib/python3.9/site-packages/numpy/core/include"]

の行が、それに相当する。これはお使いのnumpyでnumpy.get_include()を実行することで確認できる。

cimport numpy as cnp

をしたあとは、numpyの型を使えるようになる。たとえば、上記の例では

cpdef LL sum_cython_numpy(cnp.ndarray[cnp.int64_t, ndim=1] s):

という行で、引数がcnp.ndarray[cnp.int64_t, ndim=1]という型である情報を与えている。

型の定義は以下のファイルで確認できる。

https://github.com/cython/cython/blob/0.28.x/Cython/Includes/numpy/__init__.pxd#L730

次に方法2であるが、引数にLL[:]という型を付けている。繰り返しになるがCythonはnumpy配列のメモリを読み取り、自動的Typed Memoryviewsとして扱う。これはnp.ndarrayと同じメモリを参照しているが、numpyのメソッドを用いることはできない点に注意だ。

一応速度比較を行ってみよう。

from cythonized_functions import sum_cython_numpy, sum_typed_memoryviews
import numpy as np
rand = np.random.randint(0, 100, size=100)

%timeit sum_cython_numpy(rand)
# -> 1.38 µs ± 21.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%timeit sum_typed_memoryviews(rand)
# -> 275 ns ± 2.16 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Typed Memoryviewsを用いたほうが5倍以上早いことが読み取れる。この理由は方法2の関数呼び出し等のオーバーヘッドの少なさからくる。逆にrandの要素数が大きければ両者の速度は近づいていく。

Cythonを含んだmoduleのinstall例

最後にCythonを含んだmoduleをインストールできるようにしよう。 以下の具体例を参考にして書くことによって, pip install -e .で自作moduleがインストールできるようになる。またPyPIに公開するときも参考になるだろう。

具体例

ひとまず以下の構造にファイルを配置してみよう。

.
├── setup.py
├── MANIFEST.in
└── cython_module_example
    ├── __init__.py
    └── cythonized_functions.pyx

setup.pyは以下

import setuptools

try:
    import numpy as np
except ImportError:
    from setuptools import dist
    dist.Distribution().fetch_build_eggs(['numpy==1.21.2']) #ご使用のversionに
    import numpy as np

try:
    from Cython.Build import cythonize
except ImportError:
    from setuptools import dist
    dist.Distribution().fetch_build_eggs(['Cython==0.29.24']) #ご使用のversionに
    from Cython.Build import cythonize

setuptools.setup(
    name="cython_module_example", #ご自由に名前を付けてください
    version="0.0.1",
    author="Your name",
    author_email="hogehoge@example.com",
    url="hogehoge",
    description="自分で作ったやつ",
    packages=['cython_module_example'], # __init__.pyがおいてあるディレクトリ名を相対的に指定します。
    classifiers=[
        "Programming Language :: Python :: 3.9.6", #ここも環境に合わせて書いてください
        "Operating System :: OS Independent",
    ],
    install_requires=['numpy==1.21.2',
                      'Cython==0.29.24'],
    setup_requires=['setuptools>=18.0',
                    'Cython>=0.29.24',
                    'numpy>=1.21.2'],
    ext_modules=cythonize(["cython_module_example/cythonized_functions.pyx"]), #ここが味噌です
    include_dirs=np.get_include(), #ここでcythonにnumpyのpathを通します
    include_package_data=True #これがあると配布時にpyxのファイルを含めることができます。下記のMANIFEST.inを定義してください。
)

MANIFEST.inは以下

include cython_module_example/*.pyx

cythonized_functions.pyxは以下

# distutils: language=c++
# distutils: extra_compile_args = ["-O3"]
# cython: language_level=3, boundscheck=False, wraparound=False
# cython: cdivision=True

ctypedef long long LL

cpdef LL sum_typed_memoryviews(LL[:] s):
    cdef LL ret = 0
    cdef LL l = s.shape[0]
    cdef LL i
    for i in range(l):
        ret += s[i]
    return ret

__init__.pyは空で問題ない。

ここまで準備できたところで、setup.pyのあるディレクトリに戻って、以下のコマンドを実行しよう。

pip install -e .

すると、Pythonのインタラクティブシェルで以下のimportができるようになるだろう。

from cython_module_example.cythonized_functions import sum_typed_memoryviews

これで、どこでも自作moduleが使えるようになった。

解説

すべてについて解説するのは長くなってしまうので、ここではCythonに関連する部分を解説する。setup.pyに書き方は他のブログや公式ドキュメントに譲ることにする。

まずは最初に気になるのはこの部分だろう。

try:
    import numpy as np
except ImportError:
    from setuptools import dist
    dist.Distribution().fetch_build_eggs(['numpy==1.21.2']) #ご使用のversionに
    import numpy as np

try:
    from Cython.Build import cythonize
except ImportError:
    from setuptools import dist
    dist.Distribution().fetch_build_eggs(['Cython==0.29.24']) #ご使用のversionに
    from Cython.Build import cythonize

Cythonを用いたmoduleのinstallにはnumpyやCythonが必要であり、それらが事前にないとエラーがでる。そのため、"なかったらインストールする"という処理を最初に加えている。

次に、installするときにコンパイルすることを命令する部分。

    ext_modules=cythonize(["cython_module_example/cythonized_functions.pyx"]), #ここが味噌です
    include_dirs=np.get_include(), #ここでcythonにnumpyのpathを通します
    include_package_data=True #これがあると配布時にpyxのファイルを含めることができます。下記のMANIFEST.inを定義してください。

ここはおまじないである。

最後に、必須ではないが、MANIFEST.inも大事である。このファイルを定義しておくことで、配布する(python setup.py sdist)ときにpyxファイルも一緒に配布する。 pyxファイル以外のファイルを含めるときにもMANIFEST.inに追加しよう。例えば小規模なサンプルデータや、パラメーターなどを記述したjson, yamlファイルを含む場合に便利である。

さいごに

長い技術ブログになってしまったが、使いたいところからコピペではじめてだんだん理解していくのが良さそうである。 自分もそうだったし、今も知らないことばかりである。 最後に繰り返しになるが、このブログでは以下のようなことについてコピペで動くサンプルと、簡単な解説を提供した。

• pyxファイルのコンパイル方法
  • 高速化のおまじないを添えて
• 高速化されているかの確認方法
• CythonでC++のstdの使用例
  • ここではstd::vectorを用いてランレングス符号化を実装する
• Cythonとnumpyの連携
  • cimport numpy を用いる方法
  • Typed Memoryviews を用いる方法
• Cythonを含んだmoduleのinstall例
  • 自作packageをpip installできるようにする

以上のことは自分が実務で高速化するときも用いている事柄である。最近リリースした紙面ビューアーのハイライト機能では、画像処理の高速化をする際に用いた(関連記事 既存OCRを超える精度の文字領域の検出と、それに基づくiOS版「紙面ビューアーアプリ」のUXの向上 )。 この記事によってPythonを高速化したい人の一助になれたなら幸いである。

以上 Nikkei Advent Calendar 2021 の25日目の記事でした。他にも多種多様な記事が公開されているので、年末年始の情報収集にいかがでしょうか？