Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

# お題はこちら

## Rotations.jl について

$$\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)
|__/                   |

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


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 を使うメリット

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 もあります．

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

おまたてしました．最低限の 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\")")
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


$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/Pythonの記事を書いています．
Why not register and get more from Qiita?
1. We will deliver articles that match you
By following users and tags, you can catch up information on technical fields that you are interested in as a whole
2. you can read useful information later efficiently
By "stocking" the articles you like, you can search right away