PyCall.jlでPyTorch使ってDeep Learningする

[Julia Advent Calendar 2017](https://qiita.com/advent-calendar/2017/julialang)の23日目の記事です。
12/23が1日だけ空いていたので埋めてみます。


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

# Julia製のDeep Learningフレームワークの紹介。
JuliaのDeep Learningフレームワークとしては、
<dl>
  <dt>[Mocha.jl](https://github.com/pluskid/Mocha.jl)</dt>
  <dd>Juliaでは古参のフレームワーク。現在も活発に開発されています。[Caffe](http://caffe.berkeleyvision.org/)に影響を受けたということで、ネットワーク定義の記述なんかがCaffeによく似ています。Define-and-Runです。</dd>
  <dt>[MXNet.jl](https://github.com/dmlc/MXNet.jl)</dt>
  <dd>Apacheが開発している[MXNet](https://mxnet.apache.org/)のJuliaバインディングです。フレームワークの完成度は現状もっとも高いと思います。データサイエンティストが本気でデータ分析するならこれでしょうか。Define-and-Runです。</dd>
  <dt>[Knet.jl](https://github.com/denizyuret/Knet.jl)</dt>
  <dd>Julia界隈では最近話題のフレームワークです。Define-by-Runです。GPUのカーネルをかなり頑張って書いているようで、他の有名なDLフレームワークなんかよりも[速い](https://github.com/ilkarman/DeepLearningFrameworks)場合もあるとか。</dd>
  <dt>[Flux.jl](https://github.com/FluxML/Flux.jl)</dt>
  <dd>Julia界隈では、上のKnet.jlと並んで言及されることが多いフレームワークです。Define-by-Runです。現時点での完成度は全体としてはKnet.jlのほうが上だと思いますがが、特にRNN関係はFlux.jlのほうが充実しているようです。</dd>
  <dt>[Merlin.jl](https://github.com/hshindo/Merlin.jl)</dt>
  <dd>奈良先端科学技術大のチームが開発しているフレームワーク。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は使えないようです。</dd>
  <dt>[TensorFlow.jl](https://github.com/malmaud/TensorFlow.jl)、[Keras.jl](https://github.com/invenia/Keras.jl)</dt>
  <dd>圧倒的な知名度を持つ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です。</dd>
</dl>

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

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

# Juliaから[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)経由で[PyTorch](http://pytorch.org/)を使う
というわけで、人気急上昇中の[PyTorch](http://pytorch.org/)を、Juliaから使います。
Define-by-RunのDeep Learningフレームワークといえば日本製の[Chainer](https://chainer.org/)でしょ。なんで、[PyTorch](http://pytorch.org/)なの？っていう疑問は置いておきます。
[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)を使えば、そんなに大変ではありません。

## インストール
まずは、[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)、[PyTorch](http://pytorch.org/)、[torchvision](http://pytorch.org/docs/master/torchvision/)をインストールします。
Juliaとは別にインストールしたPython環境を使うこともできますが、JuliaのPkgディレクトリの中に専用のPython環境をインストールして使うほうが、いろいろと分かりやすいと思います。以下は、その前提でいきます。

### [PyCall.jl](https://github.com/JuliaPy/PyCall.jl)のインストール
```julia
Pkg.add("PyCall")
```
するだけです。[Conda.jl](https://github.com/JuliaPy/Conda.jl)によってPython環境がJuliaのPkgディレクトリ内にインストールされます。
デフォルトだと、Python 2.7がインストールされます。Python3を使いたい場合は、

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

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

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

具体的には、Linuxなら、

```julia
using PyCall
pyimport_conda("torch", "pytorch", "pytorch")
```
でOKです。

Windowsの場合は、

```julia
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](http://pytorch.org/docs/master/torchvision/)のインストール（pip経由）
torchvision のインストールはcondaではなくて、pipを使うのが一般的なようです。[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)環境でpipを使うには、pipをモジュールとしてimportして、main()を直接呼べばよいです。

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

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

```julia
using PyCall
pyimport_conda("yaml", "pyyaml")
```
とconda経由であらかじめPyYamlをインストールしておくのが手っ取り早いです。

## 使ってみる
ぐだぐだ説明するより、コードを見るのが早いでしょう。
PyTorchの最初のチュートリルである、[Deep Learning with PyTorch: A 60 Minute Blitz](http://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html)
の[Neural Networks](http://pytorch.org/tutorials/beginner/blitz/neural_networks_tutorial.html)の章をやってみます。
畳み込みネットを定義して、1ステップの学習をするという例題になります。

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

### オリジナルのPythonコード
```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のコード
```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`メソッドは、

```julia
num_flat_features(self, x) = prod(x[:size]()[2:end])
```
と一行で書くこともできます。

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

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

[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)には、`numpy.array`と、Juliaの`Array`の間の相互変換を自動的に行ってくれる機能があるので、`numpy.array`を経由することで、PyTorchの`torch.FloatTensor`と、Julia の`Array`との相互変換が可能です。

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

```jlcon
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]()`を使います。例えば、上のコードでは、

```julia
input = Variable(torch[:randn](1, 1, 32, 32))
```
とPyTorchの`torch[:randn]`を呼び出して、直接`torch.FloatTensor`を作っていますが、

```jlcon
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](https://github.com/JuliaPy/PyCall.jl)では、オーバロードが可能な`[]`演算子を使って、`PyObject`（Pythonのオブジェクトの[PyCall.jl](https://github.com/JuliaPy/PyCall.jl)による表現)
が持つプロパティへのアクセスするようにしているわけです。

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

```julia
output[:data][:cpu]()[:numpy]()
```
という気持ち悪い記述を、

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

##`@pyimport`を使わない理由
[PyCall.jl](https://github.com/JuliaPy/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](https://github.com/JuliaPy/PyCall.jl)はv2.0になる予定ですが、そこで、`@pyimport`は削除される（少なくとも非推奨）になります。
近い将来、削除or非推奨になることが確実な機能を、新しいコードで使うのはやめるべきでしょう。

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

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

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

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

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

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

Julia言語の大きな特徴といえば[multiple dispatch](http://nbviewer.jupyter.org/gist/StefanKarpinski/b8fe9dbb36c1427b9f22)が挙げられます。

いわゆるオブジェクト指向言語(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であれば、

```python
container.append(1)
```
となるのに対して、Juliaでは、

```julia
push!(container, 1)
```
となります。
これが、**Julian way** です。

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

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

```julia
Base.getproperty(o, s::Symbol) = isdefined(o, s) ? getfield(o, s) : (args...) -> eval(s)(o, args...)
```
という「おまじない」を書くとよい、とすすめるサイトが現れたら？
（こうすることで、すべての`func(arg1, args...）`という関数呼び出しを、`arg1.func(args...)`とオブジェクト指向的に記述することが可能になるはずです[^1]）

[^1]: 動くように修正しました。（v1.3で動作確認済）

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

[^2]: 異論がある人はたくさんいると思います。そもそも私はこれら言語は全く素人です。

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

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