Julia Advent Calendar 2017の23日目の記事です。
12/23が1日だけ空いていたので埋めてみます。

Juliaで最近流行りの機械学習というかDeep Learningをやりたいです。
というわけで、JuliaのDeep Learningフレームワークをいくつか触ったりして試行錯誤していたのですが、最近、PyCall.jl経由でPyTorchを使う、ってことに落ち着いてます。その紹介。

Julia製のDeep Learningフレームワークの紹介。

JuliaのDeep Learningフレームワークとしては、

[Mocha.jl](https://github.com/pluskid/Mocha.jl): Juliaでは古参のフレームワーク。現在も活発に開発されています。[Caffe](http://caffe.berkeleyvision.org/)に影響を受けたということで、ネットワーク定義の記述なんかがCaffeによく似ています。Define-and-Runです。
[MXNet.jl](https://github.com/dmlc/MXNet.jl): Apacheが開発している[MXNet](https://mxnet.apache.org/)のJuliaバインディングです。フレームワークの完成度は現状もっとも高いと思います。データサイエンティストが本気でデータ分析するならこれでしょうか。Define-and-Runです。
[Knet.jl](https://github.com/denizyuret/Knet.jl): Julia界隈では最近話題のフレームワークです。Define-by-Runです。GPUのカーネルをかなり頑張って書いているようで、他の有名なDLフレームワークなんかよりも[速い](https://github.com/ilkarman/DeepLearningFrameworks)場合もあるとか。
[Flux.jl](https://github.com/FluxML/Flux.jl): Julia界隈では、上のKnet.jlと並んで言及されることが多いフレームワークです。Define-by-Runです。現時点での完成度は全体としてはKnet.jlのほうが上だと思いますがが、特にRNN関係はFlux.jlのほうが充実しているようです。
[Merlin.jl](https://github.com/hshindo/Merlin.jl): 奈良先端科学技術大のチームが開発しているフレームワーク。Define-by-Runです。このAdvent Calendarでもいくつか記事を執筆されている[regonn](https://qiita.com/regonn)氏が、[別のAdvent Calendar](https://qiita.com/advent-calendar/2017/deep-learning-frameworks)で簡単な[解説](https://qiita.com/regonn/items/3d9a77d6be902f0f0ff6)を書いてくれています。現状、CPU専用でGPUは使えないようです。
[TensorFlow.jl](https://github.com/malmaud/TensorFlow.jl)、[Keras.jl](https://github.com/invenia/Keras.jl): 圧倒的な知名度を持つGoogleの機械学習フレームワーク[TensorFlow](https://www.tensorflow.org/)のJuliaバインディングです。ただ、TensorFlowは機械学習フレームワークであってDeep Learningフレームワークではないので、実際にDeep Learningするのは結構大変です。[Keras](https://keras.io/)が欲しい……と思ったら、公式パッケージには登録されていませんが [Keras.jl](https://github.com/invenia/Keras.jl)というのがあるようです（今、知りました）。ドキュメントを見ると、下で説明するように、まさに[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)を使ってPython版のKerasを呼んでいるようです。Define-and-Runです。

なんかがあります。
私はダイナミックなネットワークをいろいろ試してみたいと思っているので、Define-by-Runがいいわけですが、JuliaのDefine-by-Runフレームワークは、現時点では、正直、どれも完成度がいまいちなんですよね。
実装されている機能ブロック（CONVとかLSTMとか）が少ないので、チュートリアルから外れたちょっと複雑なことをやろうとすると、ほぼ自分で書かないといけなくなります。

やっぱり現時点ではDeep Learningフレームワークの充実度という点では、JuliaはPythonにかなわないです。
Define-by-Runの本家Chainerや、最近人気急上昇中のPyTorchを、Juliaでも使いたいです。

JuliaからPyCall.jl経由でPyTorchを使う

というわけで、人気急上昇中のPyTorchを、Juliaから使います。
Define-by-RunのDeep Learningフレームワークといえば日本製のChainerでしょ。なんで、PyTorchなの？っていう疑問は置いておきます。
PyCall.jlを使えば、そんなに大変ではありません。

インストール

まずは、PyCall.jl、PyTorch、torchvisionをインストールします。
Juliaとは別にインストールしたPython環境を使うこともできますが、JuliaのPkgディレクトリの中に専用のPython環境をインストールして使うほうが、いろいろと分かりやすいと思います。以下は、その前提でいきます。

PyCall.jlのインストール

Pkg.add("PyCall")

するだけです。Conda.jlによってPython環境がJuliaのPkgディレクトリ内にインストールされます。
デフォルトだと、Python 2.7がインストールされます。Python3を使いたい場合は、

ENV["CONDA_JL_VERSION"] = "3"
Pkg.add("PyCall")

と、Pkg.add("PyCall")の前に、環境変数CONDA_JL_VERSIONを設定しておきます。

PyTorch のインストール（conda経由）

PyTorch は、condaでのインストールが推奨されています。PyCall.jl環境であれば、pyimport_condaを使えばよいです。
GPUを使う場合には、あらかじめ、CUDAのツールキットをインストールしておきます。(cuDNNや、cuPy等のインストールは不要です）

具体的には、Linuxなら、

using PyCall
pyimport_conda("torch", "pytorch", "pytorch")

でOKです。

Windowsの場合は、

using PyCall
pyimport_conda("cuda", "cuda90", "peterjc123") # エラーは無視してOK
pyimport_conda("torch", "pytorch", "peterjc123")

で入ります。2行目のpyimport_conda("cuda", "cuda90", "peterjc123")の"cuda90"の部分は、CUDAのバージョンに合わせて変えてください。このコマンドを実行するとFailed to import required Python module cudaみたいなエラーが出ると思いますが、インストール自体はできているので無視して大丈夫です。

macでは、GPUを使いたい場合にはConda経由では駄目でソースからのコンパイルが必要とのことですが、私はよくわかっていないです。

torchvisionのインストール（pip経由）

torchvision のインストールはcondaではなくて、pipを使うのが一般的なようです。PyCall.jl環境でpipを使うには、pipをモジュールとしてimportして、main()を直接呼べばよいです。

using PyCall
@pyimport pip
pip.main(["install", "torchvision"])

PyYamlのインストールでエラーが出てしまった場合には、

using PyCall
pyimport_conda("yaml", "pyyaml")

とconda経由であらかじめPyYamlをインストールしておくのが手っ取り早いです。

使ってみる

ぐだぐだ説明するより、コードを見るのが早いでしょう。
PyTorchの最初のチュートリルである、Deep Learning with PyTorch: A 60 Minute Blitz
のNeural Networksの章をやってみます。
畳み込みネットを定義して、1ステップの学習をするという例題になります。

オリジナルのPythonのコードと、Juliaのコードを比較してみてください。

オリジナルのPythonコード

import torch
from torch.autograd import Variable
import torch.nn as nn
import torch.nn.functional as F

class Net(nn.Module):

    def __init__(self):
        super(Net, self).__init__()
        # 1 input image channel, 6 output channels, 5x5 square convolution
        # kernel
        self.conv1 = nn.Conv2d(1, 6, 5)
        self.conv2 = nn.Conv2d(6, 16, 5)
        # an affine operation: y = Wx + b
        self.fc1 = nn.Linear(16 * 5 * 5, 120)
        self.fc2 = nn.Linear(120, 84)
        self.fc3 = nn.Linear(84, 10)

    def forward(self, x):
        # Max pooling over a (2, 2) window
        x = F.max_pool2d(F.relu(self.conv1(x)), (2, 2))
        # If the size is a square you can only specify a single number
        x = F.max_pool2d(F.relu(self.conv2(x)), 2)
        x = x.view(-1, self.num_flat_features(x))
        x = F.relu(self.fc1(x))
        x = F.relu(self.fc2(x))
        x = self.fc3(x)
        return x

    def num_flat_features(self, x):
        size = x.size()[1:]  # all dimensions except the batch dimension
        num_features = 1
        for s in size:
            num_features *= s
        return num_features

net = Net()
print(net)

params = list(net.parameters())
print(len(params))
print(params[0].size())  # conv1's .weight

input = Variable(torch.randn(1, 1, 32, 32))
out = net(input)
print(out)

output = net(input)
target = Variable(torch.arange(1, 11))  # a dummy target, for example
criterion = nn.MSELoss()

loss = criterion(output, target)
print(loss)

net.zero_grad()     # zeroes the gradient buffers of all parameters

print('conv1.bias.grad before backward')
print(net.conv1.bias.grad)

loss.backward()

print('conv1.bias.grad after backward')
print(net.conv1.bias.grad)

import torch.optim as optim

# create your optimizer
optimizer = optim.SGD(net.parameters(), lr=0.01)

# in your training loop:
optimizer.zero_grad()   # zero the gradient buffers
output = net(input)
loss = criterion(output, target)
loss.backward()
optimizer.step()    # Does the update

対応するJuliaのコード

using PyCall
const torch = pyimport("torch")
const Variable = torch[:autograd][:Variable]
const nn = torch[:nn]
const F = torch[:nn][:functional]

@pydef type Net <: nn[:Module]

    __init__(self) = begin
        pybuiltin(:super)(Net, self)[:__init__]()
        self[:conv1] = nn[:Conv2d](1, 6, 5)
        self[:conv2] = nn[:Conv2d](6, 16, 5)
        self[:fc1] = nn[:Linear](16 * 5 * 5, 120)
        self[:fc2] = nn[:Linear](120, 84)
        self[:fc3] = nn[:Linear](84, 10)
    end

    forward(self, x) = begin
        # Max pooling over a (2, 2) window
        x =  F[:max_pool2d](F[:relu](self[:conv1](x)), (2, 2))
        # If the size is a square you can only specify a single number
        x = F[:max_pool2d](F[:relu](self[:conv2](x)), 2)
        x = x[:view](-1, self[:num_flat_features](x))
        x = F[:relu](self[:fc1](x))
        x = F[:relu](self[:fc2](x))
        x = self[:fc3](x)
    end

    num_flat_features(self, x) = begin
        size = x[:size]()[2:end]  # all dimensions except the batch dimension
        num_features = 1
        for s in size
            num_features *= s
        end
        return num_features
    end

end

net = Net()
println(net)

params = pybuiltin(:list)(net[:parameters]())
println(length(params))
println(params[1][:size]())  # conv1's .weight

input = Variable(torch[:randn](1, 1, 32, 32))
out = net(input)
println(out)

output = net(input)
target = Variable(torch[:arange](1, 11))  # a dummy target, for example
criterion = nn[:MSELoss]()

loss = criterion(output, target)
println(loss)

net[:zero_grad]()     # zeroes the gradient buffers of all parameters

println("conv1.bias.grad before backward")
println(net[:conv1][:bias][:grad])

loss[:backward]()

println("conv1.bias.grad after backward")
println(net[:conv1][:bias][:grad])

const optim = torch[:optim]

# create your optimizer
optimizer = optim[:SGD](net[:parameters](), lr=0.01)

# in your training loop:
optimizer[:zero_grad]()   # zero the gradient buffers
output = net(input)
loss = criterion(output, target)
loss[:backward]()
optimizer[:step]()

Pythonのa.bを、Juliaではa[:b]と書く、ということだけを気を付ければ、ほぼ1対1に対応している感じです。
ちなみに、上のコードではあえて元のPythonコードと1対1で対応するように書いているので、あまりJuliaっぽくない書き方になっている部分もあります。例えば、コード中のnum_flat_featuresメソッドは、

num_flat_features(self, x) = prod(x[:size]()[2:end])

と一行で書くこともできます。

PyTorchの`torch.FloatTensor`と、Julia の`Array`との相互変換

上のコードは、PyTorchのtorch.FloatTensorの機能を直接呼び出していて、JuliaのArrayを使っていないので、実質、Juliaで書いている意味がありません。

PyCall.jlには、numpy.arrayと、JuliaのArrayの間の相互変換を自動的に行ってくれる機能があるので、numpy.arrayを経由することで、PyTorchのtorch.FloatTensorと、Julia のArrayとの相互変換が可能です。

具体的には、PyTorchのtorch.FloatTensorの[:numpy]()メソッドでnumpy.arrayに変換すれば、PyCall.jlの自動変換機能によって、JuliaのArrayが得られます。例えば、

julia> output[:data][:cpu]()[:numpy]()
1×10 Array{Float32,2}:
 -0.0923741  -0.017475  0.0750015  …  0.0123179  -0.0106525  0.0577791

このコードでは、PyTorchのVaribleであるoutputから、output[:data]でtorch.FloatTensorを取り出して、さらに（GPU上にあるかもしれない）データを[:cpu]()でCPUに持ってきて、それから、[:numpy]()の呼び出しによって、numpy.array経由でJuliaのArrayを得ています。

逆に、JuliaのArrayをPyTorchに渡したい場合は、torch[:from_numpy]()を使います。例えば、上のコードでは、

input = Variable(torch[:randn](1, 1, 32, 32))

とPyTorchのtorch[:randn]を呼び出して、直接torch.FloatTensorを作っていますが、

julia> input = Variable(torch[:from_numpy](randn(Float32, (1, 1, 32, 32)))) 
PyObject Variable containing:
(0 ,0 ,.,.) = 
 -1.1789  0.2133  0.8533  ...   1.3377 -0.5419 -0.3846
  0.0239 -0.6832 -0.0486  ...  -0.6551  0.0996 -0.9503
 -0.6340  0.1009 -1.2199  ...   0.6817 -0.4266 -1.2444
           ...             ⋱             ...          
 -1.7544  0.5275  0.0504  ...  -1.3610  2.0705  0.9354
  0.8037 -1.4190 -2.2780  ...   1.3173 -0.3348  0.6090
 -1.2804  0.2406  0.7395  ...   0.3972 -0.1215  1.2111
[torch.FloatTensor of size 1x1x32x32]

のようにすれば、Juliaのrandn()関数で作ったJuliaのArrayから、torch.FloatTensor（やPyTorchのVariable）を得ることができます。

`a[:b]`が気持ち悪い

上のコードですが、a[:b]という書き方が気持ち悪いですね。（実は、個人的には、だいぶ見慣れてしまって、あまり違和感を感じなくなってきていますが）
Juliaでは（v0.6以前は）、.演算子をオーバーロードすることができせん。そのため、PyCall.jlでは、オーバロードが可能な[]演算子を使って、PyObject（PythonのオブジェクトのPyCall.jlによる表現)
が持つプロパティへのアクセスするようにしているわけです。

それでも、a[:b]じゃなくて、a.bと書きたいのであれば、Juliaの強力なマクロ機能を使うことが考えられます。これを行うのが、PySyntax.jlです。このパッケージでは、@pyというマクロを定義していて、a.bという記述を a[:b]に書き換えます。
例えば、上の、

output[:data][:cpu]()[:numpy]()

という気持ち悪い記述を、

@py output.data.cpu().numpy()

と書くことができます。
ただし、@pyマクロは、a.bという記述全てを機械的にa[:b]に書き換えてしまうので、例えば、Iterators.enumerate(iter)みたいな、Julia本来の.表記もIterators[:enumerate](iter)に書き換わってしまうので注意が必要です。
個人的には、自分でプログラムを書くのであれば、中途半端にPySyntax.jlを使うくらいなら、全てa[:b]で書いた方が良い気がしています。
ただ、github等にある最新論文のPyTorchモデル（当然、Pythonで書かれたもの）をちょっと試してみたい、なんかの時にはPythonの記述をほとんど変更することなくJuliaで実行できるPySyntax.jlは便利です。

`@pyimport`を使わない理由

PyCall.jlを使ったことがある人は、上のコードをみて、なぜ@pyimportマクロを使わないのかと思われるかもしれません。
個人的には、@pyimportは使わないで、pyimport()関数を使うほうがよいと思います。

理由その１

@pyimportは、Juliaのmoduleをその場（on-the-fly）で作ることで、a[:b]をa.bと記述することを可能にする仕組みです。
ただ、a.bと記述するのが可能なのは、@pyimportでインポートしたモジュールがもつメソッドだけです。Pythonのオブジェクト（PyObject）のメソッドはやっぱり、a[:b]で記述する必要があります。
その結果、a.bとa[:b]の２つの流儀の記述が入り混じって、よけいにワケワカラン状態になります。

理由その２

すぐ下に書いているように、Juliaの次期バージョン v0.7（≒v1.0）では、.演算子のオーバーロードが導入されます。
同時にPyCall.jlはv2.0になる予定ですが、そこで、@pyimportは削除される（少なくとも非推奨）になります。
近い将来、削除or非推奨になることが確実な機能を、新しいコードで使うのはやめるべきでしょう。

v0.7(≒v1.0)で導入される新機能 `.`演算子のオーバーロード

…などと、a[:b]という記法について書いてきたわけですが。。。

実は、タイムリーなことに、つい１週間ほど前に、Juliaの開発版masterに#24960がマージされました。
これによりコア言語の仕様として.演算子のオーバーロードが可能になります。
.演算子のオーバーロードの要望は、5年ほど前からissueに挙がっていたのですが、ついに実装されました。

今後、a.bという構文は、Base.getproperty(a, :b)の呼び出しとして解釈されることになります。したがって、Base.getproperty()を特定の型について定義することで、a.bという構文の振る舞いをプログラム上で変えることが可能になります。

PyCall.jlの開発者も、Julia v0.7にあわせて、.演算子のオーバーロードを全面的に使ってPyCallを書き直してPyCall v2.0をリリースすると表明しているので、近い将来、上のコードの醜悪なa[:b]をすべてa.bと記述することが可能になると思われます。

`.`演算子のオーバーロードがJulia界へ与える影響？

PyCall.jlを考えれば、.演算子のオーバーロードはうれしい限りです。
一方で、この新機能は、場合によっては、Juliaの生態系に大きな(悪)影響を与える可能性もあります。
実際、新機能の導入に否定的な意見（あるいは懸念）もかなりあったようです。

Julia言語の大きな特徴といえばmultiple dispatchが挙げられます。

いわゆるオブジェクト指向言語(Java、C++、Python、Ruby、…etc）では、メソッドは特定のオブジェクトに「所有」されています。
これは、通常、所有者のオブジェクトが、メソッドの第一引数として渡される、という形で実装されます。（Pythonであれば、メソッドの第一引数selfとして明示的に渡されます。JavaやC++ではthisという形で暗黙のうちに渡されます。）
一方、メソッドを呼び出すときには、object.method(arg1, arg2)のように、メソッドの所有者であるオブジェクト（＝メソッドの第一引数）は、引数リストではなくて、.の前に記述します。
つまり、.演算子というのは、メソッドの多数ある引数のうちで第一引数のみを特別扱いする仕組み（構文）なわけです。

これに対して、Juliaはmultiple dispatchな言語なので、Juliaの関数は特定のオブジェクトに「所有」されているわけではありません。したがって、Juliaの関数は全ての引数を平等に扱い、第一引数を特別扱いする仕組み（構文）はありません。
例えば、JuliaのVector、あるいはPythonのlistのようなコンテナcontainerに要素を追加したい場合、Pythonであれば、

container.append(1)

となるのに対して、Juliaでは、

push!(container, 1)

となります。
これが、Julian way です。

ところが、.演算子のオーバーロードは、オブジェクト指向の色が非常に強い要素をJuliaに持ち込むことになります。
世の中に、オブジェクト指向言語はものすごく多いわけで、メソッドがある特定のオブジェクトに「所有」されている（≒メソッドの第一引数を特別扱いする）という思想に慣れている人はとても多いでしょう。例えば、将来、Juliaがものすごく人気の言語になったとして、多くのデータサイエンティストがPythonからJuliaに移行してきたとしたら、Juliaでも、.演算子のオーバーロードを使用してmultiple dispatchではなくてオブジェクト指向なプログラムをバンバン書き出して、事実上、そっちが標準の書き方になってしまわないか…
実際に.演算子をオーバーロードするには、Base.getproperty、Base.setproperty! を新たに定義する必要があるわけでそれなりにJuliaの言語思想を理解している人でないと難しいでしょうけど、.演算子のオーバーロードを多用して見かけ上オブジェクト指向言語のように使えるパッケージが作られて、Pythonなどからの移行者を中心にそれを使うのが一般的になる可能性はあるのでは？

あるいは、極論ですけど、Julia初心者向けに、.juliarc.jlに

Base.getproperty(o, s::Symbol) = isdefined(o, s) ? getfield(o, s) : (args...) -> eval(s)(o, args...)

という「おまじない」を書くとよい、とすすめるサイトが現れたら？
（こうすることで、すべてのfunc(arg1, args...）という関数呼び出しを、arg1.func(args...)とオブジェクト指向的に記述することが可能になるはずです¹）

実際、GoとかLuaとかScalaとか、本来はオブジェクト指向言語ではないにも関わらず、構文糖衣みたいなものを使って、事実上、オブジェクト指向的に記述することが標準になっている言語もありますよね²。

というわけで、Julia界の醜いアヒルの子 .演算子のオーバーロードが白鳥となったとき、Julia界は果たして如何に？今後に注目！？

かなり長くなってしまった＆話が脱線しましたが、これで終わりです。

動くように修正しました。（v1.3で動作確認済） ↩
異論がある人はたくさんいると思います。そもそも私はこれら言語は全く素人です。 ↩