pyjuliaでPyTorch使ってDeep Learningする

昨年のアドベントカレンダーで、[PyCall.jlでPyTorch使ってDeep Learningする](https://qiita.com/yatra9/items/0a1a9a5ba19e9efe08c0)という投稿をしましたが、その続きです。（なんか、１年前と同じことやってますね。。）

昨年は、[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)を使ってJuliaから[PyTorch](https://pytorch.org/)を呼ぶという話でしたが、今年は、逆に、PyTorch（Python）から[pyjulia](https://github.com/JuliaPy/pyjulia)を使ってJuliaを呼ぶ、という話です。


# 愚痴
自分１人だけなのでであれば昨年書いたように、JuliaからPyTorchを呼ぶでいいんですが、いかんせん複数人でいろんなモデルを作って試すとなると、みんなPython使ってるわけで、やっぱり、現状ではPythonをメインに使わざるをえない部分があります。
実際、PyTorch便利ですよね。なんだかんだいって、[Flux.jl](https://github.com/FluxML/Flux.jl)や[Knet.jl](https://github.com/denizyuret/Knet.jl)などとは、フレームワークとしての完成度が違います。

なんですけど、PyTorch（というかPython）なんかで、深層学習なり数値計算するとなると、とにかくベクトル演算で一度に処理するように書かないといけない（For文は使ってはいけない）という強い縛りがあるわけです。なにしろ、PythonのFor文は死ぬほど遅いですからね。
もちろん、普通の、CNNとかLSTMとかなら、PyTorchが用意しているLayerを呼べばいいだけなんですが、いかんせん、ちょっと複雑なことをする独自Layerを作りたいとかなると、（PyTorchに用意されている関数では）どうやってもベクトル演算として書けないなんて処理もでてきます。

そういう場合、普通は、C++あるいはCUDAで[extention](https://pytorch.org/tutorials/advanced/cpp_extension.html)を作るということになるわけですが、C++ではなくてJuliaで書きたい。なにしろ、Juliaなら（[GPUArrays.jl](https://github.com/JuliaGPU/GPUArrays.jl)使えば）、CPUだけではなくて、GPUも簡単に使えるんですよ。

だいたい、そうでなくても、いったん Juliaに一度なれてしまうと、とにかく、Python（numpy、PyTorch）の全てをベクトル演算で書け（For文を使ったら負け）という縛りは、ほんとにキツイです。やってられないです。もしかして、私だけ？そんなことないですよね、みんなそう思いますよね。
人間の自然な思考はFor文であって、ベクトル演算ではないですよね。

愚痴がすぎました。

# pyjulia使って、PyTorchのTensorをjuliaでいじる
というわけで本題です。
[pyjulia](https://github.com/JuliaPy/pyjulia)を使うと、PyTorchのTensorをJuliaでいじれます。C/C++などでグルーコードを書く必要もありません。pythonのスクリプトの中にJuliaのスクリプトを文字列として埋め込むだけです。（もちろん、PythonとJuliaを別ファイルに分けることもできます）
PyTorchのカスタムレイヤーをJuliaで書く、それも、CPUにもGPUも同一のコードで対応できます。device agnosticというやつですね。

あるいは、カスタムレイヤーだけではなくて、結果の可視化なんかでも、ちょっと独自の可視化をやりたいとかなると、Pythonだとかなりつらかったりするんですが、Juliaの描画ライブラリ業界における期待の新星[Makie.jl](https://github.com/JuliaPlots/Makie.jl)なんかを使えば、大量のニューロンの状態を独自の方法で描画したい、とかも高速にできます。

# pyjuliaのインストール
pyjuliaの[公式ページ](https://github.com/JuliaPy/pyjulia)にしたがって、pyjuliaをインストールします。
注意点としては、公式ページにも書いてありますが、anaconda版のPythonではそのままでは動きません。
`python-jl`を使うか、あるいは、`pyenv`で`--enable-shared`オプションを使ってPythonをインストールする必要があります。

私は、anaconda版のPythonをインストールしたあと、`pyenv`に含まれる[python-build](https://github.com/pyenv/pyenv/tree/master/plugins/python-build)を使って、`--enable-shared`オプション版のPythonをanaconda版のpythonの上に上書きしてしまう、という強引な方法を使っています。
anaconda版のPythonとpython-buildするPythonのバージョンを合わせれば、とりあえずは、問題なく動きます。


# PyTorchからJuliaを呼ぶ
まずは、以下の内容のJuliaスクリプトを、カレントディレクトリに`PyTorchJulia.jl`という名前でつくります。
（本当は、真面目にPackageとして作るのがいいんでしょうけど、面倒なので）
これが、PyTorchのTensorをJuliaのArrayやCuArrayに変換するインタフェースのコードです。

```julia
module PyTorchJulia
export jl_array, jl_type
using GPUArrays
using CuArrays
using PyCall

const torch = pyimport("torch")
const is_cuarray_unsafe_wrap = :own in fieldnames(CuArray)
const py_str = pybuiltin("str")
const torch_types = Dict(
        "torch.float32" => Float32,
        "torch.float64" => Float64,
        "torch.float16" => Float16,
        "torch.uint8" => UInt8,
        "torch.int8" => Int8,
        "torch.int16" => Int16,
        "torch.int32" => Int32,
        "torch.int64" => Int64)

jl_type(dtype) = torch_types[pycall(py_str, String, dtype)]

function jl_array(x::PyObject)
    @assert x[:is_contiguous]()
    T = jl_type(x[:dtype])
    dims = reverse(x[:shape])
    device = x[:device]
    device_type = device[:type]
    if device_type == "cpu"
        ptr = Ptr{T}(x["data_ptr"]())
        return unsafe_wrap(Array, ptr, dims)
    elseif device_type == "cuda"
        if is_cuarray_unsafe_wrap
            ptr = Ptr{T}(x["data_ptr"]())
            return unsafe_wrap(CuArray, ptr, dims)
        else
            ctx = CuArrays.CuCurrentContext()
            @assert CuArrays.CUDAdrv.device(ctx) == CuArrays.CuDevice(device[:index])
            buf = CuArrays.Mem.Buffer(
                Ptr{Cvoid}(x["data_ptr"]()),
                x[:numel]() * sizeof(T),
                ctx)
            CuArrays.Mem.retain(buf)
            return CuArray{T, length(dims)}(buf, dims)
        end
    end
end

GPUArrays.backend(::Type{<:Array}) = GPUArrays.JLBackend()
end
```

で、これを使う、Python（PyTorch）のコードを書きます。

```Python
import torch
from julia import Main as julia

jl_script = '''
include("PyTorchJulia.jl")
using .PyTorchJulia

using GPUArrays
using PyCall
const torch = pyimport("torch")

Base.getproperty(o::PyObject, s::Symbol) = isdefined(o, s) ? getfield(o, s) : getindex(o, s)

f_impl!(x) = (x .*= 2)
f!(x) = f_impl!(jl_array(x))

g_impl!(out, x) = (@. out = x * 2)
function g(x)
    out = torch.empty_like(x)
    g_impl!(jl_array(out), jl_array(x))
    return out
end

function h_impl!(out, x, y)
    gpu_call(out, (out, x, convert(eltype(x), y))) do state, out, x, y
        idx = @linearidx(out)
        out[idx] = x[idx] + y
        return
    end
end
function h(x, y)
    out = torch.empty_like(x)
    h_impl!(jl_array(out), jl_array(x), y)
    return out
end
'''

julia.eval(jl_script)
x_org = torch.empty((2,3,4), dtype=torch.float32).uniform_()
x_cpu = x_org.clone()
julia.f_b(x_cpu)
assert torch.allclose(x_cpu, 2 * x_org)
x2_cpu = julia.g(x_cpu)
assert torch.allclose(x2_cpu, 2 * x_cpu)
x3_cpu = julia.h(x_cpu, 3.0)
assert torch.allclose(x3_cpu, x_cpu + 3)

if torch.cuda.is_available():
    x_gpu = x_org.cuda()
    julia.f_b(x_gpu)
    assert torch.allclose(x_gpu.cpu(), x_cpu)
    x2_gpu = julia.g(x_gpu)
    assert torch.allclose(x2_gpu.cpu(), x2_cpu)
    x3_gpu = julia.h(x_gpu, 3.0)
    assert torch.allclose(x3_gpu.cpu(), x3_cpu)
```

Pythonのスクリプト上で文字列として定義したJuliaの関数（`f!(x)`, `g(x)`, `h(x, y)`）を、
PyTorchのTensorである`x_cpu`を引数として
それぞれ、`julia.f_b(x_cpu)`，`x2_cpu = julia.g(x_cpu)`，`x3_cpu = julia.h(x_cpu, 3.0)` と呼び出しています。

あるいは、PyTorchのCUDA Tensorである、`x_gpu`を引数として、
それぞれ、`julia.f_b(x_gpu)`，`x2_gpu = julia.g(x_gpu)`，`x3_gpu = julia.h(x_gpu, 3.0)` と呼び出しています。

## 簡単な解説

- Pythonではビックリマーク`!`を識別名として使えないため、juliaの`f!`という名前の関数は、pyjuliaの仕様として、ビックリマーク`!`が`_b`に置換されて、Pythonからは`f_b`となります。

- `PyTorchJulia.jl`で定義している`jl_array()`という関数で、PyTorchのTensorをJuliaのArrayあるいはCuArrayに変換しています。
PyTorchのTensorのメモリをそのままJuliaのArrayやCuArrayとしてラップしているだけで、コピーすることなく高速に変換できます。

- ただし、PyTorchのTensorのメモリ配置がRow-mejorであるのに対して、JuliaのArrray（やCuArray）はColumn-majorです。jl_arrayではメモリのコピーをしないかわりに、PyTorchのTensorをJuliaのArray（やCuArray）に変換すると、次元の順番が逆になります。例えば、PyTorchの次元(2,3,4)のTensorを、jl_arrayでJuliaのArrayに変換すると、次元は(4,3,2)になります。
（Arrayを使う代わりに、例えば、[Strided.jl](https://github.com/Jutho/Strided.jl)を使えば、次元をひっくり返さないですみそうです。もし、Package化するとなれば、対応したいです）

- PyTorchのTensorがCPU上にある場合にはJuliaのArrayに、GPU上にある場合にはCuArrayに変換されます。上のコードでは、`f_impl!(x)`, `g_impl(x, y)`, `h_impl!(out, x, y)`はArrayかCuArrayかに関わらず同一コードで動作します。（device agnostic)
とくに、`h_impl!(out, x, y)`では、[GPUArrays.jl](https://github.com/JuliaGPU/GPUArrays.jl)の`gpu_call`を使って、CUDAのカーネルをPure Juliaで記述して呼んでいます。

- juliaのコードで、`f!(x)`と`f_impl!(x)`を別の関数に分けている（一つの関数にまとめていない）のは、[function barrier](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#kernel-functions-1)を使いたいためです。
`jl_array()`の戻り値は、PyTorchのTensorがCPUにあるかGPUにあるかによって、`Array`型か`CuArray`型になります。関数を分けることで、`f_impl!(x)`は、`x`の型に応じて特殊化されたコードがコンパイルされるので、型安定なコードになって、高速に計算できます。

- Juliaコード内で、新たなTensorを確保してPyTorch側に返したい場合には、`g(x)`の中で行っているように`out = torch.empty_like(x)`という形で、PyTorchのメソッドを使ってPyTorch側でTensorを確保する必要があります。

- 本当は、[昨年の記事](https://qiita.com/yatra9/items/0a1a9a5ba19e9efe08c0)で詳しく説明したように、juliaからPython（PyTorch）の関数を呼ぶ場合には、[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)の仕様により、`torch[:empty_like](x)`と書かなければなりません。これは、とても面倒なので、v1.0で導入された`Base.getproperty`をオーバーロードすることで、`out = torch.empty_like(x)`と呼べるようにしています。
なお、[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)公式でも、この記法を可能にする[PR](https://github.com/JuliaPy/PyCall.jl/pull/517)が提案されています。

- `PyTorchJulia.jl`の中で、`is_cuarray_unsafe_wrap`の値に応じて、`unsafe_wrap(CuArray, ptr, dims)`、あるいは、`CuArray{T, length(dims)}(buf, dims)`を呼んでいます。実は、この記事を書いている途中で、[CuArrays.jl](https://github.com/JuliaGPU/CuArrays.jl)に、unsafe_wrapを作ってほしいと[提案](https://github.com/JuliaGPU/CuArrays.jl/issues/203)したところ、[公式に対応していただきました](https://github.com/JuliaGPU/CuArrays.jl/pull/214)。すでに、masterにマージされていますが、まだ、Releaseにはなっていないので、旧版にも対応できるようにif文でわけています。

# JuliaからPyTorchを使う場合
Juliaをメインに使いたい場合、つまり、[昨年度の記事](https://qiita.com/yatra9/items/0a1a9a5ba19e9efe08c0)のように、JuliaからPyCall.jl経由でPyTorchを使う場合は、以下のJuliaスクリプトのようになります。
`Base.getproperty`をオーバーロードしたことで、PythonとJuliaのコードはほぼコピペで動くようになりました。

```julia

include("PyTorchJulia.jl")
using .PyTorchJulia

using GPUArrays
using PyCall
const torch = pyimport("torch")

Base.getproperty(o::PyObject, s::Symbol) = isdefined(o, s) ? getfield(o, s) : getindex(o, s)

f_impl!(x) = (x .*= 2)
f!(x) = f_impl!(jl_array(x))

g_impl!(out, x) = (@. out = x * 2)
function g(x)
    out = torch.empty_like(x)
    g_impl!(jl_array(out), jl_array(x))
    return out
end

function h_impl!(out, x, y)
    gpu_call(out, (out, x, convert(eltype(x), y))) do state, out, x, y
        idx = @linearidx(out)
        out[idx] = x[idx] + y
        return
    end
end
function h(x, y)
    out = torch.empty_like(x)
    h_impl!(jl_array(out), jl_array(x), y)
    return out
end

x_org = torch.empty((2,3,4), dtype=torch.float32).uniform_()
x_cpu = x_org.clone()
f!(x_cpu)
@assert torch.allclose(x_cpu, 2 * x_org)
x2_cpu = g(x_cpu)
@assert torch.allclose(x2_cpu, 2 * x_cpu)
x3_cpu = h(x_cpu, 3.0)
@assert torch.allclose(x3_cpu, x_cpu + 3)

if torch.cuda.is_available()
    x_gpu = x_org.cuda()
    f!(x_gpu)
    @assert torch.allclose(x_gpu.cpu(), x_cpu)
    x2_gpu = g(x_gpu)
    @assert torch.allclose(x2_gpu.cpu(), x2_cpu)
    x3_gpu = h(x_gpu, 3.0)
    @assert torch.allclose(x3_gpu.cpu(), x3_cpu)
end
```

## 今後の展望（野望）

JuliaからPyTorchを使う場合、PyCall.jlを使ってPython経由で使うということになるわけですが、せっかくJulia使っているんで、できればPythonのお世話になりたくないと思ったりします。
実は、PyTorchは、もうすぐReleaseされる予定のv1.0から、C++ (C++11以降)でフロントエンドを記述できるようになります。[公式ページ](https://pytorch.org/cppdocs/frontend.html)を見ていただくとわかりますが、C++のコードはPythonのコードとほぼ一対一に対応すると言ってよいです。
となると、このC++ Frontend APIを、最近Julia v1.0に向けて作業が進んでいる[Cxx.jl](https://github.com/Keno/Cxx.jl)を使ってラップすれば、Juliaから直接PyTorchを使うことができることになりそうです。できたら、やってみたいところです。