Nim でDeepLearningをしたい　（MNISTサンプルを触ります）

# 本日は

Nim で Deep Learning をしたい!
というモチベーションでArraymancerのにあるMNISTのサンプルを動かしてみましょう。

# 導入
最初は

```shell-session
$ nimble install arraymancer
```

でやっていたのですがMNISTのサンプルを読み込むロジックをインポートできなかったのでソースから直接インストールする方法をとります。

```shell-session
$ git clone https://github.com/mratsim/Arraymancer.git
$ cd Arraymancer
$ nimble install
```

# 学習セットの準備

学習データセットを導入します。

作業ディレクトリに `build` というディレクトリを作りその直下に下記のファイルを
http://yann.lecun.com/exdb/mnist/
から持ってきて配置します。

> train-images-idx3-ubyte.gz:  training set images (9912422 bytes) 
train-labels-idx1-ubyte.gz:  training set labels (28881 bytes) 
t10k-images-idx3-ubyte.gz:   test set images (1648877 bytes) 
t10k-labels-idx1-ubyte.gz:   test set labels (4542 bytes)

これらは予め解凍しておきましょう。

#コードを書く

コードは　

https://github.com/mratsim/Arraymancer/blob/master/examples/ex02_handwritten_digits_recognition.nim

から拝借します。

```nim:mnist.nim
import ../src/arraymancer, random

# This is an early minimum viable example of handwritten digits recognition.
# It uses convolutional neural networks to achieve high accuracy.
#
# Data files (MNIST) can be downloaded here http://yann.lecun.com/exdb/mnist/
# and must be decompressed in "./build/" (or change the path "build/..." below)
#

# Make the results reproducible by initializing a random seed
randomize(42)

let
  ctx = newContext Tensor[float32] # Autograd/neural network graph
  n = 32                           # Batch size

let
  # Training data is 60k 28x28 greyscale images from 0-255,
  # neural net prefers input rescaled to [0, 1] or [-1, 1]
  x_train = read_mnist_images("build/train-images.idx3-ubyte").astype(float32) / 255'f32

  # Change shape from [N, H, W] to [N, C, H, W], with C = 1 (unsqueeze). Convolution expect 4d tensors
  # And store in the context to track operations applied and build a NN graph
  X_train = ctx.variable x_train.unsqueeze(1)

  # Labels are uint8, we must convert them to int
  y_train = read_mnist_labels("build/train-labels.idx1-ubyte").astype(int)

  # Idem for testing data (10000 images)
  x_test = read_mnist_images("build/t10k-images.idx3-ubyte").astype(float32) / 255'f32
  X_test = ctx.variable x_test.unsqueeze(1)
  y_test = read_mnist_labels("build/t10k-labels.idx1-ubyte").astype(int)

# Configuration of the neural network
network ctx, DemoNet:
  layers:
    x:          Input([1, 28, 28])
    cv1:        Conv2D(x.out_shape, 20, 5, 5)
    mp1:        MaxPool2D(cv1.out_shape, (2,2), (0,0), (2,2))
    cv2:        Conv2D(mp1.out_shape, 50, 5, 5)
    mp2:        MaxPool2D(cv2.out_shape, (2,2), (0,0), (2,2))
    fl:         Flatten(mp2.out_shape)
    hidden:     Linear(fl.out_shape, 500)
    classifier: Linear(500, 10)
  forward x:
    x.cv1.relu.mp1.cv2.relu.mp2.fl.hidden.relu.classifier

let model = ctx.init(DemoNet)

# Stochastic Gradient Descent (API will change)
let optim = model.optimizerSGD(learning_rate = 0.01'f32)

# Learning loop
for epoch in 0 ..< 5:
  for batch_id in 0 ..< X_train.value.shape[0] div n: # some at the end may be missing, oh well ...
    # minibatch offset in the Tensor
    let offset = batch_id * n
    let x = X_train[offset ..< offset + n, _]
    let target = y_train[offset ..< offset + n]

    # Running through the network and computing loss
    let clf = model.forward(x)
    let loss = clf.sparse_softmax_cross_entropy(target)

    if batch_id mod 200 == 0:
      # Print status every 200 batches
      echo "Epoch is: " & $epoch
      echo "Batch id: " & $batch_id
      echo "Loss is:  " & $loss.value.data[0]

    # Compute the gradient (i.e. contribution of each parameter to the loss)
    loss.backprop()

    # Correct the weights now that we have the gradient information
    optim.update()

  # Validation (checking the accuracy/generalization of our model on unseen data)
  ctx.no_grad_mode:
    echo "\nEpoch #" & $epoch & " done. Testing accuracy"

    # To avoid using too much memory we will compute accuracy in 10 batches of 1000 images
    # instead of loading 10 000 images at once
    var score = 0.0
    var loss = 0.0
    for i in 0 ..< 10:
      let y_pred = model.forward(X_test[i*1000 ..< (i+1)*1000, _]).value.softmax.argmax(axis = 1).squeeze
      score += accuracy_score(y_test[i*1000 ..< (i+1)*1000], y_pred)

      loss += model.forward(X_test[i*1000 ..< (i+1)*1000, _]).sparse_softmax_cross_entropy(y_test[i*1000 ..< (i+1)*1000]).value.data[0]
    score /= 10
    loss /= 10
    echo "Accuracy: " & $(score * 100) & "%"
    echo "Loss:     " & $loss
    echo "\n"
```

前半はデータセットのロードとネットワークの定義です。後半は学習(Validationも含む)のループを回しています。
Trainerを使わないChainerを知っているとそこまで抵抗なくソースを読むことができますね。

# 動かしましょう

これをビルドして実行します。

```shell-session
$ nim c -d:release mnist.nim
$ ./mnist
Elapsed   0.01471120000000001[sec]
Epoch is: 0
Batch id: 0
Loss is:  194.3992156982422
Elapsed   7.5248887[sec]
Epoch is: 0
Batch id: 200
Loss is:  2.858644962310791
Elapsed   14.9354357[sec]
Epoch is: 0
Batch id: 400
Loss is:  1.602403879165649
...
```

学習は進んでいるように見えます。
ちなみに動作が非常に遅い・CPUが仕事をしていない時はBLASとかが入っていないのかもしれません。

http://verifiedby.me/adiary/091 

の記事をそのままつかって導入しました。

```shell-session
$ sudo apt install libatlas-base-dev
$ sudo apt install libatlas-doc
$ sudo apt install libopenblas-base
$ sudo apt install libopenblas-dev
```

#わかっていないこと
## GPU で実行したい

`Arraymancer/nim.cfg` をながめるとおそらく

```
nim c -d:release -d:cuda mnist.nim
```

でいけるのではないでしょうか。私の環境ではソースのコンパイルで失敗します。

## CPU のバックエンドをMKLにしたい

たぶん
```
nim c -d:release -d:mkl mnist.nim
```
でいけるはず・・・。~~MKLはインストールしたつもりですが `-d:mkl` で動作速度が変わらないので効いていないのかも・・・。~~

~~わかっていないことが多いですが、動かせたので良いことにします。~~




# おまけ

Chainer 版はこちら。

```python3:pymnist.py
import chainer
import chainer.functions as F
import chainer.links as L
from chainer import training, datasets, iterators, optimizers
from chainer.training import extensions
import numpy as np

DEVICE = -1
BATCH_SIZE = 32


class MNISTCONV(chainer.Chain):

    def __init__(self):
        super(MNISTCONV, self).__init__()
        with self.init_scope():
            ksize = 5
            self.cv1 = L.Convolution2D(1, 20, ksize=ksize)
            linear_size1 = (28 - (ksize - 1)) // 2
            self.cv2 = L.Convolution2D(20, 50, ksize=ksize)
            linear_size2 = (linear_size1 - (ksize - 1)) // 2
            self.hidden = L.Linear(linear_size2 * linear_size2 * 50, 500)
            self.classifier = L.Linear(500, 10)

    def __call__(self, x):
        h = self.cv1(x)
        h = F.relu(h)
        h = F.max_pooling_2d(h, 2)
        h = self.cv2(h)
        h = F.relu(h)
        h = F.max_pooling_2d(h, 2)
        h = self.hidden(h)
        h = F.relu(h)
        h = self.classifier(h)
        return h


def train():
    model = L.Classifier(MNISTCONV())
    if DEVICE >= 0:
        chainer.cuda.get_device_from_id(DEVICE).use()
        chainer.cuda.check_cuda_available()
        model.to_gpu()

    optimizer = optimizers.SGD(lr=0.01)
    optimizer.setup(model)

    train, test = chainer.datasets.get_mnist(ndim=3)
    train_iter = chainer.iterators.SerialIterator(
        train, BATCH_SIZE, shuffle=True)
    test_iter = chainer.iterators.SerialIterator(
        test, BATCH_SIZE, repeat=False, shuffle=False)

    updater = training.updaters.StandardUpdater(
        train_iter, optimizer, device=DEVICE)
    trainer = training.Trainer(updater, (5, 'epoch'), out='result')

    log_interval = (1, 'epoch')
    trainer.extend(extensions.Evaluator(test_iter, model,
                                        device=DEVICE), trigger=log_interval)
    trainer.extend(extensions.LogReport(trigger=log_interval))
    trainer.extend(extensions.PrintReport(
        ['iteration', 'main/loss', 'validation/main/loss',
         'main/accuracy', 'validation/main/accuracy', 'elapsed_time']), trigger=log_interval)

    trainer.extend(extensions.ProgressBar())

    with chainer.using_config('train', True):
        trainer.run()

    chainer.serializers.save_npz('result/mnistconv.npz', model)
```

CPUの実行速度で比べるとNumpyのバックエンドをOpenBLASにしているときはNimのほうが若干速いです。
(MKLをバックエンドにしてiDepp を導入するととても速くなりChainerが勝ちます。)

# 追記

`-d:mkl` が仕事しないので `Arraymancer/nim.cfg` にあるように

```shell-session

$ nim c -d:release \
 --define:"openmp" \
 --define:"blas=mkl_intel_lp64" \
 --define:"lapack=mkl_intel_lp64" \
 --clibdir:"/opt/intel/mkl/lib/intel64" \
 --passL:"/opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a" \
 --passL:"-lmkl_core" \
 --passL:"-lmkl_gnu_thread" \
 --passL:"-lgomp" \
 --dynlibOverride:"mkl_intel_lp64" \
 mnist.nim
```

と愚直にコンパイルオプションを書いたら実行の速度が改善されました。ChainerでのMKL+iDeepの組み合わせといい勝負をしています。悶々が解決できてよかった・・・。