本日は

他言語からJulia言語を呼び出すヒント(前編 : 実行形式でC言語から呼び出すよ) の続きを書きます．前半では，C言語の main のなかでJuliaのコードを実行する方法を公式ドキュメントのExampleを基にして書きました．

やりたいこと

今回は手元のMacで ほにゃらら.dylib (Linux だと ほにゃらら.so)　の形にして Julia を呼び出してみたいと思います．

お題はこちら

実応用を考えると，自分で書いたコードや外部パッケージを呼び出したいといふような要請は当然としてあるのでそれを想定したシナリオを考えます．ここでは回転行列の実装をサポートする Rotation.jl に依存したコードを呼び出すぞーな目標で進んでいきます．

`Rotations.jl` について

例えば３次元空間における $z$ 軸の周りでの回転行列は RotZ で得ることができます．

$$
\operatorname{RotZ}(\theta)=
\begin{bmatrix}
\cos\theta & -\sin\theta & 0 \\
\sin\theta & \cos\theta & 0 \\
0 & 0 & 1
\end{bmatrix}
$$

Rotations.jl は公式パッケージなのでパッケージモードで add すれば入手できます．上記の式をJuliaのREPL上で実現してみます．（Julia の開発版で行なっていますが安定版の Julia1.1.1 でも動きます）

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.3.0-DEV.547 (2019-07-15)
 _/ |\__'_|_|_|\__'_|  |  Commit ba039209aa (5 days old master)
|__/                   |

(v1.3) pkg> add Rotations
julia> using Rotations
julia> θ=π/6
0.5235987755982988
julia> r=RotZ(θ);
julia> r
3×3 RotZ{Float64} with indices SOneTo(3)×SOneTo(3)(0.523599):
 0.866025  -0.5       0.0
 0.5        0.866025  0.0
 0.0        0.0       1.0

回転行列が定義されたので適当なベクトルに作用させてみましょう

julia> φ = 5π/6
2.6179938779914944
julia> vec=[cos(φ),sin(φ),0]
3-element Array{Float64,1}:
 -0.8660254037844387
  0.49999999999999994
  0.0
julia> r*vec
3-element StaticArrays.SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 -1.0
  0.0
  0.0
julia> # 実はSymPy の要素も作用させられる
julia> using SymPy
julia> @vars x y z
(x, y, z)
julia> r*[x,y,z]
3-element StaticArrays.SArray{Tuple{3},Sym,1,3} with indices SOneTo(3):
 0.866025403784439*x - 0.5*y
 0.5*x + 0.866025403784439*y
                       1.0*z

余談ですが $\sin$ と $\cos$ の値を両方取得する方法として sincos という Julia の関数が提供されています．

julia> sincos(θ)
(0.49999999999999994, 0.8660254037844387)
help?> sincos
search: sincos MissingException

  sincos(x)

  Simultaneously compute the sine and cosine of x, where the x is in radians.

  ────────────────────────────────────────────────────────────────────────────────────────────────

  sincos(A::AbstractMatrix)

  Compute the matrix sine and cosine of a square matrix A.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> S, C = sincos(fill(1.0, (2,2)));

  julia> S
  2×2 Array{Float64,2}:
   0.454649  0.454649
   0.454649  0.454649

  julia> C
  2×2 Array{Float64,2}:
    0.291927  -0.708073
   -0.708073   0.291927

Rotation.jl を使うメリット

回転行列ぐらいなら自力で書くこともできますが，StaticArrays.jl で提供されている型を利用しますのでパフォーマンスがよくなります．パフォーマンス計測ツールとして BenchmarkTools.jl を用います

julia> using Rotations
julia> using BenchmarkTools
julia> θ=π/6
0.5235987755982988
julia> rotz=RotZ(θ)
3×3 RotZ{Float64} with indices SOneTo(3)×SOneTo(3)(0.523599):
 0.866025  -0.5       0.0
 0.5        0.866025  0.0
 0.0        0.0       1.0
julia> @benchmark rotz*rand(3)
BenchmarkTools.Trial:
  memory estimate:  144 bytes
  allocs estimate:  2
  --------------
  minimum time:     67.636 ns (0.00% GC)
  median time:      68.834 ns (0.00% GC)
  mean time:        77.553 ns (9.67% GC)
  maximum time:     27.189 μs (99.71% GC)
  --------------
  samples:          10000
  evals/sample:     978
julia> # 自前で書いた場合
julia> myrotz=[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1]
3×3 Array{Float64,2}:
 0.866025  -0.5       0.0
 0.5        0.866025  0.0
 0.0        0.0       1.0
julia> @benchmark myrotz*rand(3)
BenchmarkTools.Trial:
  memory estimate:  224 bytes
  allocs estimate:  2
  --------------
  minimum time:     104.570 ns (0.00% GC)
  median time:      106.462 ns (0.00% GC)
  mean time:        115.148 ns (5.48% GC)
  maximum time:     28.184 μs (99.58% GC)
  --------------
  samples:          10000
  evals/sample:     937

なお，一般的に $Ax+b$ みたいな変換を提供する CoordinateTransformations.jl もあります．

See https://github.com/FugroRoames/CoordinateTransformations.jl

お題はわかったのでサンプルはよ

おまたてしました．最低限の Rotations.jl のチュートリアルをしましたのでここで使われた RotZ を呼び出すコードを書いてみます．ファイル名を external.jl と適当につけてみます．

external.jl

import Rotations

function rotateZ(θ::Float64,v::Array{Float64})
    r = Rotations.RotZ(θ)
    return convert(Array{Float64},r*v)
end

return 式で convert(Array{Float64},r*v) がありますがこれは大人の事情です.

C の方のコードも書くぞいわ

~~「今日も１日がんばるぞい」の「ぞい」をフランス語読みするとゾワになります Galois 理論とかもそうですね．~~

external.jl と同じ階層に下記のような共有ライブラリのコードを書きます．ファイル名は何でもいいですがここでは jl_application.c という名前で記述しておきます．

jl_application.c

#include <julia.h>
#include <stdio.h>

void initialize() {
    jl_init();
}

void finalize() {
    jl_atexit_hook(0);
}

//実装と名前が噛み合ってませんが，皆さんのセンスに合わせて適宜修正してください
void using_Rotations() {
    /*
    Technically, you can include your julia script.
    Note that, you can't use jl_eval_string("include(\"external.jl\")")
    use "Base.include(Main, \"external.jl\")" instead
    See https://github.com/JuliaLang/julia/issues/28825
    */
    jl_eval_string("Base.include(Main, \"external.jl\")");
}

int rot_z(double* arr, double theta, double x, double y, double z) {
    char eq[256];
    sprintf(eq, "rotateZ(%lf,Float64[%lf,%lf,%lf])", theta, x, y, z);
    jl_value_t *rotz = jl_eval_string(eq);
    double* ret = (double*)jl_array_data(rotz);
    for (int i = 0; i < 3; ++i)
    {
        arr[i] = ret[i];
    }
    return 0;
}

解説

rot_z という関数を呼び出すことを目標にします.それを呼び出す前に jl_init でJulia本体を起動しないといけません.
rot_z の実装が external.jl に依存しているのでそれもロードしておく必要があります．JuliaのREPL上だと

julia> include("external.jl")

ですませられますが Base.include(Main, "external.jl") を evaluate させないといけないようですね．詳しくはこちらのIssue https://github.com/JuliaLang/julia/issues/28825 をどうぞ

jl_call1 で自前で定義した関数を呼ぶとなぜかセグフォしてしまう（辛い）ので jl_eval_string に突っ込むために sprintf で文字列を作成し rotateZ(theta, Float64[x,y,z]) のような式を実行させるようにしています（闇）．あとは arr に値を渡します.
external.jl で定義した rotateZ の戻り値を Array{Float64} にコンバートしないと (*double) へキャストできないせいか（呼び出しの時点で）セグフォがおきます．先ほど述べたとおりRotations.RotZ(theta) は特殊な型を持つのでそれが原因かなと思います．jl_value_t って結構やわなのね(´・ω・｀)．
ちょっと闇闇してますが良さげな回避策があるとJuliaコミュニティがとても喜びます．

ビルドしたいです

コードが書けたのでコンパイルしましょう. Makefile も書いておきましょう．

Makefile

JL_SHARE = $(shell julia -e 'print(joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia"))')
CFLAGS   += $(shell $(JL_SHARE)/julia-config.jl --cflags)
CXXFLAGS += $(shell $(JL_SHARE)/julia-config.jl --cflags)
LDFLAGS  += $(shell $(JL_SHARE)/julia-config.jl --ldflags)
LDLIBS   += $(shell $(JL_SHARE)/julia-config.jl --ldlibs)

TARGET_NAME=jl_application
C_PROGRAM_FILE=jl_application.c

all: run

shared_example: $(C_PROGRAM_FILE)
    gcc $< -o $(TARGET_NAME).dylib -O2 -shared $(CFLAGS) $(LDFLAGS) $(LDLIBS)

run: shared_example
    python call_shared_from.py # あとで出てくる

clean:
    rm -f ./$(TARGET_NAME).dylib

下記のように make することで共有ライブラリができます．

$ ls
Makefile jl_application.c external.jl
$ make shared_example
gcc jl_application.c -o jl_application.dylib -O2 -shared -std=gnu99 -I'/Users/terasaki/julia_dev/include/julia' -DJULIA_ENABLE_THREADING=1 -fPIC -L'/Users/terasaki/julia_dev/lib' -Wl,-rpath,'/Users/terasaki/julia_dev/lib' -ljulia
# jl_application.dylib が生成される
$ Makefile             external.jl          jl_application.dylib jl_application.c

呼び出しは？

後編に続きます．

他言語からJulia言語を呼び出すヒント(中編：共有ライブラリ呼び出すよ)

本日は