他言語からJulia言語を呼び出すヒント(後編:Pythonから呼び出してみるよ)

本日は

他言語からJulia言語を呼び出すヒント(実行形式でC言語から呼び出すよ：中編) の続編です．

JuliaのコードをC言語を経由して共有ライブラリ形式までビルドし呼び出す準備ができました．

$ ls
Makefile external.jl  jl_application.dylib jl_application.c

jl_application.dylib を呼び出そうということで ctypes を用いて Python から呼び出しましょう．

注意

• PythonからJuliaを呼び出す場合は https://github.com/JuliaPy/pyjulia を使うのが賢い選択です．

Python から呼び出す

次のようなコードを用意しました

"""
usage

$ make
$ python call_shared_from.py
rotateZ(pi / 6, Float64[1, 0, 0]) = [0.866025, 0.5, 0.0]
Hello From Python World
[0.86602529 0.50000019 0.        ]
"""

import time
import math
import ctypes
from ctypes import cdll
import os
import numpy as np


class JL_App(object):

    def __init__(self, lib_path):
        self.lib = cdll.LoadLibrary(lib_path)
        self.lib.initialize()

    def __enter__(self):
        return self.lib

    def __exit__(self, type_, value, traceback):
        self.lib.finalize()
        return True


def main():
    # use context manager to control jl_init and jl_atexit
    with JL_App("jl_application.dylib") as jl_app:
        jl_app.using_Rotations()
        jl_app.rot_z.argtypes = [
            np.ctypeslib.ndpointer(flags="F"),
            ctypes.c_double,
            ctypes.c_double,
            ctypes.c_double,
            ctypes.c_double,
        ]
        theta = np.pi/6
        arr = np.empty(shape=(3), dtype=np.float64, order="F")
        jl_app.rot_z(arr, theta, 1., 0., 0.)
        print("Hello From Python World")
        print(arr)

if __name__ == '__main__':
    main()

解説

• 背伸びして with statement を用います．これでファイルのオープン/クローズと同様に Julia の init exit を自動的にすることができます.

• 呼び出すためにライブラリをロードします

from ctypes import cdll
cdll.LoadLibrary("/path/to/dynamic.dylib")

のようにして共有ライブラリを読み込みその中で定義されている関数 rot_z　を呼び出します．これは

int rot_z(double* arr, double theta, double x, double y, double z)

のように定義されているのでこれに合わせて argtypes を定義しておきます．

• See https://docs.python.org/3/library/ctypes.html#specifying-the-required-argument-types-function-prototypes

jl_app.rot_z.argtypes = [
    np.ctypeslib.ndpointer(flags="F"),
    ctypes.c_double,
    ctypes.c_double,
    ctypes.c_double,
    ctypes.c_double,
]

• Julia は Fortran と同様に column major なので flags="F" を指定しています．

  • それに応じて rot_z の第一引数に投げる arr の order オプションも "F" を指定しています.にしておきます.
  • 今回の記事の例ですと１次元なので大丈夫ですが，例えば，回転行列そのものを受け取るといった場合に
    arr = np.empty(shape=(3,3), dtype=np.float64) としておくと得られる結果は期待したそれの転置が得られてしまいます．

• jl_app.using_Rotations() の宣言も忘れずに(実装の内容は記事の中編を参考のこと)．

無事できたでしょうか？下記のような結果が得られるとOKです．お疲れ様です．

Hello From Python World
[0.86602529 0.50000019 0.        ]

呼び出しコスト重いわ・・・

Juliaのコードを呼び出して高速化できるでしょうか？
次のようなコードを書きます．

def rot_about_z(theta, vec):
    rotz = np.array([
        [math.cos(theta), -math.sin(theta), 0],
        [math.sin(theta),  math.cos(theta), 0],
        [0, 0, 1],
    ]).transpose()
    return vec.dot(rotz)

def bench():
    # use context manager to control jl_init and jl_atexit
    with JL_App("jl_application.dylib") as jl_app:
        jl_app.using_Rotations()
        jl_app.rot_z.argtypes = [
            np.ctypeslib.ndpointer(flags="F"),
            ctypes.c_double,
            ctypes.c_double,
            ctypes.c_double,
            ctypes.c_double,
        ]
        theta = np.pi/6
        import time
        stats = []
        arr = np.empty(shape=(3), dtype=np.float64, order="F")
        for i in range(1000):
            theta = 2*np.pi*np.random.random()
            vec = np.random.random(3)
            s = time.time()
            jl_app.rot_z(arr, theta, *vec)
            e = time.time()
            stats.append(1000*(e-s))
        print("Hello From Python World")
        print("mean elapsed", 1000*np.array(stats).mean(), "mu-sec")

    stats = []

    for i in range(1000):
        theta = 2*np.pi*np.random.random()
        vec = np.random.random(3)
        s = time.time()
        res = rot_about_z(theta, vec)
        e = time.time()
        stats.append(1000*(e-s))
    print("mean elapsed", 1000*np.array(stats).mean(), "mu-sec")

bench を呼び出した結果は次のようになります

# Julia
mean elapsed 151.30257606506348 mu-sec
# NumPy
mean elapsed 2.9592514038085938 mu-sec

うーーん．・・・pyjulia ではどうでしょう？

$ pip install julia # pyjulia をインストールする
$ ipython
Python 3.6.8 (default, Jul 12 2019, 20:53:32)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.6.1 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import julia
In [2]: import numpy as np
In [3]: j=julia.Julia()
In [4]: j.using("Rotations")
In [5]: j.eval("theta=pi/6; r=RotZ(theta)")
Out[5]:
[[0.8660254037844387, -0.49999999999999994, 0.0],
 [0.49999999999999994, 0.8660254037844387, 0.0],
 [0.0, 0.0, 1.0]]
In [6]: %timeit j.eval("r*rand(3)")
132 µs ± 34.6 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [7]: pyr=np.array(j.r).transpose()
In [8]: %timeit np.random.random(3).dot(pyr)
1.25 µs ± 4.11 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

なるほど・・・．やっぱり呼び出しのコストが実行時間に比べて大きいとかえって遅いんだなーってのがわかりました．各自で確かめてみてください．

まとめ

• まがりなりにも Julia のコードを共有ライブラリとして呼び出せるようにできました.
• 色々闇
• Pythonから呼び出すことを念頭におく場合は素直に pyjulia を使った方が良い.
• 簡易的でもいいからベンチマーク取っておこう

以上