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

# 本日は

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

- [Embedding Julia](https://docs.julialang.org/en/v1/manual/embedding/#Embedding-Julia-1)
- [https://github.com/JuliaLang/julia/tree/master/test/embedding](https://github.com/JuliaLang/julia/tree/master/test/embedding)

## やりたいこと

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



# お題はこちら

実応用を考えると，自分で書いたコードや外部パッケージを呼び出したいといふような要請は当然としてあるのでそれを想定したシナリオを考えます．ここでは回転行列の実装をサポートする [`Rotation.jl`](https://github.com/FugroRoames/Rotations.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 でも動きます）

```julia
               _
   _       _ _(_)_     |  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
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
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` を用います

- See https://github.com/JuliaCI/BenchmarkTools.jl
- See https://github.com/JuliaArrays/StaticArrays.jl

```julia
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` と適当につけてみます．

```julia: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` という名前で記述しておきます．

```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
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: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` することで共有ライブラリができます．

```console
$ 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
```

# 呼び出しは？

[後編](https://qiita.com/SatoshiTerasaki/items/0c4a4241387131a3dcdd) に続きます．