2017/11/26

Coulomb GAN

25個のガウス分布で生成した2次元の点の分布の学習が普通のGANでは困難だったので、Coulomb GAN [1]を使ってみました。

学習したい分布

Geometric GAN [2]の実験に使われている混合ガウス分布を使います。 25個のガウス分布を縦横それぞれに-21, -10.5, 0, 10.5, 21の位置の組合せで配置します。分散は0.1です。 下図のような分布です。

25000点サンプリングして描画しています。 点が集まっている部分の周りの楕円は、25個の混合ガウスモデル(GMM)でフィットさせることで求めています。

この分布に従ったサンプルを生成する生成器\(G\)を学習します。

普通のGAN

まずは、普通のGANで試します。生成器\(G\)のネットワークはKerasのコードで示すと、

def generator_model_elu():
    model = Sequential()
    model.add(Dense(128, input_dim=4))
    model.add(ELU())
    model.add(Dense(128))
    model.add(ELU())
    model.add(Dense(128))
    model.add(ELU())
    model.add(Dense(2))
    return model
です。4つの乱数を入力すると、2次元空間での点の座標を出力します。識別器\(D\)は
def discriminator_model_elu():
    model = Sequential()
    model.add(Dense(128, input_dim=2))
    model.add(ELU())
    model.add(Dense(128))
    model.add(ELU())
    model.add(Dense(128))
    model.add(ELU())
    model.add(Dense(1))
    return model
です。2次元空間での点の座標を入力として、真の分布から生成された点なら1、\(G\)が生成した点なら0を出力します。

どちらのオプティマイザにもAdamを使います。また、Lossはbinary_crossentropyで計算します。学習させた結果が以下です。 失敗集です。網羅的に試しているわけではありません。

ParameterValueImages
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
512
0.5
0.0001
0.001
0.5
0.0001
0.001
129エポックまで実行しましたが、同じ状況で変わりません。下の図は13エポック目の図です。
同パラメータで再実行してみましたが、下の図のように形状は変わるものの進みません。
さらに、同パラメータで再実行してみましたが、下の図のように形状は変わるもののやはり進みません。
さらに、同パラメータで再実行してみましたが、下の図のように形状は変わるものの進みません。
さらに、同パラメータで再実行してみましたが、下の図のように形状は変わるものの進みません。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
128
0.5
0.000025
0.001
0.5
0.000025
0.001
広がったものの、形状に変化がみられないので中断。中断したエポックでの結果が以下です。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
128
0.5
0.001
0.0005
0.5
0.001
0.0005
下の図は5エポック目のものですが、296エポックまで実行しても状況は変わりません。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
128
0.5
0.000025
0.001
0.5
0.000025
0.001
変化しないので、中断。中断したエポックでの結果が以下です。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
512
0.5
0.00001
0.01
0.5
0.00001
0.01
これも変化しないので、中断。中断したエポックでの結果が以下です。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
512
0.5
0.00001
0.0001
0.5
0.00001
0.0001
20エポック目くらいまでは分布が広がるものの、
最終的には、中心に集まってしまいます。残念。

全く同じパラメータで再実行してみましたが、おかしな所に点が集まり学習が進まなくなりました。

batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
512
0.5
0.000001
0.00001
0.5
0.000001
0.00001
棒がくるくる回っているのですが、時間がかかりそうなので中断。中断したエポックでの結果が以下です。
batch_size
disc \(\beta_1\)
disc decay
disc lr
gen \(\beta_1\)
gen decay
gen lr
512
0.1
0.00001
0.00001
0.5
0.00001
0.0002
途中までは比較的順調に進んだ結果が以下です。途中までは順調でしたが、途中から\(G\)が踊り始めます。クリックするとAPNGによるアニメーションが表示されます(32MB)。999エポックまで実行しました。
discは識別器のAdamのパラメータ、genは生成器のAdamのパラメータを表します。

この他にも色々な値で試しましたが、残念ながら再現したい25個のガウス分布っぽいものに行き着くことができませんでした。

Coulomb GAN

というわけで、Coulomb GAN [1]です。 各点を電荷に見立てることで、間違った分布を生成器が学習することを防ぎます。 真の分布から得られた点を正の電荷(+)、生成器から得られた点を負の電荷(-)として、 負の電荷に対するポテンシャルを考えると、下図のようになります。

電荷(サンプル点)の位置 マイナス電荷に対するポテンシャル

このポテンシャルをそのままLossとして生成器を学習します。 そうすると、ポテンシャルの低いほうへ負の電荷が移動する方向に生成器が表す分布が変形していきます。 言い換えると、負の電荷、すなわち、生成器が生成した偽の点\(\bm{a}^{\rm fake}\)同士は離れる方向に、かつ、真の分布から得られた真の点\(\bm{a}^{\rm real}\)には近づくように学習することになります。識別器\(D\)はCoulomb GANにおいては偽の点と真の点を識別するように学習するわけではなく、代わりにポテンシャルの分布を学習していきます。したがって、\(D\)のLossは \[ \mathcal{L}_D \propto \sum_i \left( D(\bm{a}_i) - \Phi(\bm{a}_i)\right)^2 \] となります。\(\Phi(\bm{a}_i)\)は点\(i\)の位置\(\bm{a}_i\)でのポテンシャルです。\(\bm{a}_i\)は偽の点と真の点の両方の位置を表しますが、計算するポテンシャルは偽の点としてのポテンシャルです。ポテンシャルの分布を直接学習できるわけではありませんが、サンプル点におけるポテンシャルの値を模擬できるように\(D\)を学習することで任意の点のポテンシャルの値を予測できるようになります。このようにポテンシャルの形状に近い形を学習した\(D\)を使うと、任意の点の勾配を計算できるようになります。

\(D\)を使うと、生成器\(G\)のLossを \[ \mathcal{L}_G \propto \sum_i D(G(z_i)) \] と計算できます。ポテンシャルが低くなる方向に偽の点が生成されやすくなるように\(G\)が学習されます。

真の点が作り出す偽の点に対するポテンシャルは次の式で計算します。 \[ \phi(r) = \frac{-1}{(r^2+\varepsilon^2)^{(d-2)/2}} \] 偽の点が作り出す偽の点に対するポテンシャルは符号が逆になるだけです。ここで、\(r\)は、ある点\(a_i\)からのユークリッド距離、\(d\)は次元数で、今回は3のみを使います。\(\varepsilon\)はソフトニングのためのパラメータです(Plummer radiusとも言う)。グラフの形状は下図のようになります。

Gnuplot Produced by GNUPLOT 5.0 patchlevel 7 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 -30 -20 -10 0 10 20 30 Potential x ε=0.25 ε=0.25 ε=0.5 ε=0.5 ε=1 ε=1 ε=2 ε=2 ε=4 ε=4

バッチサイズを\(N_B\)、真の点を\(\bm{a}^{\rm real}_i\)、偽の点を\(\bm{a}^{\rm fake}_i\)とすると、あるミニバッチで計算される、位置\(\bm{x}\)でのポテンシャルは \[ \Phi(\bm{x}) = \sum_i^{N_B} \phi(|\bm{a}_i^{\rm real}-\bm{x}|) - \sum_i^{N_B} \phi(|\bm{a}_i^{\rm fake}-\bm{x}|)\] となります。

[1]では、\(\varepsilon\)を指数関数的に減少させているようなので、それに倣って減少させるようにしました。減少する速さを想像しやすいように、半減期(Half life)の式で計算することとしました。単位はイテレーションです。半減期が5000なら、5000イテレーション後の\(\varepsilon\)が半分になります。

また、上記の\(\Phi(\bm{x})\)の定義では、ポテンシャルを発生させている点自身のポテンシャル(\(r=0\)のときの値)を計算してしまいます。その結果、\(\varepsilon\)が小さくなるにつれ、各点でのポテンシャルの絶対値が \[ \phi(t) = \frac{-1}{((\varepsilon_0 \times 2^{\frac{-t}{t_{\rm hl}}})^2)^{(d-2)/2}} \] となります。d=3とすれば、 \[ \phi(t) = \frac{-1}{\varepsilon_0} \times 2^{\frac{t}{t_{\rm hl}}} \] です。半減期\(t_{\rm hl}\)=1000、\(\varepsilon\)の初期値\(\varepsilon_0\)=1、1エポックが200イテレーション、100エポック処理すると仮定すると、\(\phi(20000)\sim -1\times 10^6\)となります。\(D\)がポテンシャルの形状を全く追随できていないとすると、MSEは\(\sim 1\times 10^{12}\)になります。floatの上限は\(3.402823\times 10^{+38}\)なので、まだ余裕がありますが、最終的には、\(D\)のLossはinfやNaNになります。

これを避けるため、以下の実験では、自分自身のポテンシャルを計算から除くようにしました。つまり、 \[ \Phi(\bm{x}) = \sum_{i,a_i^{\rm real} \neq x}^{N_B} \phi(|\bm{a}_i^{\rm real}-\bm{x}|) - \sum_{i,a_i^{\rm fake} \neq x}^{N_B} \phi(|\bm{a}_i^{\rm fake}-\bm{x}|) \] で計算します。

実験 (バッチサイズ512)

それでは、Coulomb GANで25個のガウス分布を学習してみます。

各種ハイパーパラメータの値は下表の通りです。ニューラルネットワークは普通のGANのときと同じものを使います。\(D\)の出力は、\(D\)に入力された位置のポテンシャルで、それは1次元の値なので使いまわせます。

Parameter nameValue
Samples per gaussian4000 (Total 100k)
Batch size512
\(d\)3
Adam lr for \(D\) and \(G\)0.001
Adam \(\beta_1\) for \(D\) and \(G\)0.5
Adam decay for \(D\) and \(G\)0.0001
\(\varepsilon\)3.0
Half life of \(\varepsilon\)5000 iteration
[1]ではバッチサイズを128にして実験してうまく学習できているようですが、今回の実験では学習に失敗することが多かったため、512に増やしています。ポテンシャルの計算をするときに使う真の点と偽の点の数はバッチサイズと等しいため、バッチサイズが小さいとポテンシャルの計算が雑になり、うまく\(D\)を学習できないのでしょう。ガウス分布が25個あるので、128だと真の点が平均で~5個しかなく、運悪く、ある分布に対する点が0個になってしまうことが起きていそうです。

500エポック後の結果は以下の通りです。まず、500エポック目の\(G\)で25000点を生成したときの分布が下の図です。

さらに、上の図にガウシアン25個のGMM(Gaussian Mixture Model)をフィットさせた結果が下の図です。全体にフィットしているガウス分布が1個あるため、背景が灰色になっています。ちょうど(-10,-10)の位置に集まっている点にフィットできていません。見た目ではまあまあ元の分布に近くなっています。
次の図は、500エポック目の\(D\)が推定するポテンシャル分布を表しています。真の点と偽の点がほぼ同じ位置に分布しているため、互いのポテンシャルが相殺し、大きな値は見られません。

一方、170エポック目では下の図のようにまだ(0, 0)と(10, 0)の位置に点が集まっていません。

このときの\(D\)が推定するポテンシャル分布(下図)を見てみると、その位置のポテンシャルの値が小さく(赤く)なっていることが分かります。真の点しかないので、ポテンシャルがその位置で小さくなります。他に小さくなっている部分をみてみると、(-20, -10)と(0, 20), (20, 0)があります。よーく上の図の分布を見てみると、その位置には点があるものの他の部分と違い点が集中していません。

このように、偽の点がない部分のポテンシャルが小さくなるので、ポテンシャルの小さい位置に偽の点が生成されるように\(G\)が学習されていきます。

下のAPNGによるアニメーションは、\(G\)の分布の発展の様子です。クリックすると再生されます(約34MB)。

少しずつ、抜けている部分に点が集まっていく様子が見て取れます。

実験 (バッチサイズ128)

バッチサイズ128、\(\varepsilon\)の半減期20000、Adam decay \(2.5\times 10^{-5}\)で試してみたときの、500エポック目の\(G\)が生成した点が下の図になります。イテレーションの回数が4倍に増えるので、半減期とAdam decayを調整していますが、それ以外は同じです。

中央(0, 0)の分布を再現できていないことが分かります。\(D\)が作るポテンシャルの分布は下の図のようになっていて、
見て分かるとおり、中央のポテンシャルは十分低いものの、その周辺にポテンシャルの壁がそびえ立っているため、中に点が入っていけない状況に陥っています。このケースでは70エポック目あたりから中央の周囲にポテンシャルの壁が出来始めていました。

というわけで、

今回の実験では、何かの違いにより[1]に書かれているハイパーパラメータではCoulomb GANによる学習を進めることができず、色々ハイパーパラメータを探すことになりましたが、見つかってしまえば何度実行しても学習できていました。ガウス分布の学習に関しては普通のGANより学習しやすいことは確かなようです。普通のGANよりバリエーションが豊かであることを確認してみたいのですが、まだまだ先は長そうです。

コード

今回の実験で使ったコードは以下の通りです。 Coulomb GANの著者らのTensorFlowベースのコードと[1]を参考にしつつ実装しています。

importしているモジュールはここに置いています。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
# -*- coding: utf-8 -*-
# Coulomb GAN
import sys, os, math
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../keras-examples/src')
import numpy as np
from sklearn import mixture
from util.history import ExperimentHistory
from plot.gmm_plot import (plot_points_and_gmm, plot_potential,
                           plot_real_or_fake, plot_generator_distribution)
from gan.models_2dpoints import gan_model_relu, gan_model_elu
from gan.coulomb import CoulombPotentials
from gan.gaussian_mixture.datagen import gen_2D_samples_from_5x5gm, RandomSampler
from keras.models import Sequential
from keras.optimizers import Adam
import keras.backend as K

GENERATED_IMAGE_PATH = 'coulomb_images/'
SAMPLES_PER_GAUSSIAN = 4000
BATCH_SIZE = 128*4
NUM_EPOCH = 1000

def raw_loss(y_true, y_pred):
    return K.mean(y_pred, axis=-1)

class DiscriminatorLabelNormal:
    def __call__(self, points_real, points_fake):
        assert len(points_real) == BATCH_SIZE
        assert len(points_fake) == BATCH_SIZE
        return [1]*len(points_real) + [0]*len(points_fake)

    def get_eps(self):
        return 0

class DiscriminatorLableCoulomb:
    def __init__(self, eh):
        self.total_num_batches = 0
        self.cp = CoulombPotentials(eh.plummer_kernel_dim, eh.plummer_kernel_eps,
                                    eh.plummer_kernel_ignore_self_potential)
        self.eps_half_life = eh.plummer_kernel_eps_half_life

    def __call__(self, points_real, points_fake):
        self.cp.eps = eh.plummer_kernel_eps * math.pow(2, -self.total_num_batches/self.eps_half_life)
        self.total_num_batches += 1
        potential_real, potential_fake = self.cp(points_real, points_fake)
        return np.concatenate((potential_real, potential_fake))

    def get_eps(self):
        return self.cp.eps

def train(eh):
    # Generate training data and draw points
    X_train = gen_2D_samples_from_5x5gm(SAMPLES_PER_GAUSSIAN)
    gmm = mixture.GaussianMixture(n_components=25, covariance_type='full').fit(X_train)
    plot_points_and_gmm(X_train, gmm.predict(X_train), gmm.means_, gmm.covariances_,
                        'Gaussian Mixture', GENERATED_IMAGE_PATH+"gmm.png")

    # Random sampler for each batch and for plotting points of G
    rs = RandomSampler(4, "normal" if eh.random_with_normal_dist else "uniform")

    # Make an object for making correct labels of D
    dl = DiscriminatorLableCoulomb(eh) if eh.coulomb_gan else DiscriminatorLabelNormal()

    # Make generator G and discriminator D
    generator, discriminator = gan_model_elu()
    d_opt = Adam(lr=eh.disc_Adam_lr, beta_1=eh.disc_Adam_beta_1, decay=eh.disc_Adam_decay)
    g_opt = Adam(lr=eh.gen_Adam_lr,  beta_1=eh.gen_Adam_beta_1,  decay=eh.gen_Adam_decay)
    discriminator.compile(loss='mse' if eh.coulomb_gan else 'binary_crossentropy', optimizer=d_opt)
    discriminator.trainable = False
    cgan = Sequential([generator, discriminator]) # G+D with fixed weights of D
    cgan.compile(loss=raw_loss if eh.coulomb_gan else 'binary_crossentropy', optimizer=g_opt)

    num_batches = int(X_train.shape[0] / BATCH_SIZE)
    print('Number of batches:', num_batches)
    eh.write(GENERATED_IMAGE_PATH+"history.log", {"Generator":generator, "Discriminator":discriminator},
             {"Generator opt":g_opt, "Discriminator opt":d_opt}, None)
    for epoch in range(NUM_EPOCH):
        if eh.X_train_is_shuffled:
            np.random.shuffle(X_train)
        for index in range(num_batches):
            noise = np.array(rs(BATCH_SIZE))
            points_real = X_train[index*BATCH_SIZE:(index+1)*BATCH_SIZE]
            points_fake = generator.predict(noise, verbose=0)

            # Update discriminator
            X = np.concatenate((points_real, points_fake))
            Y = dl(points_real, points_fake)
            d_loss = discriminator.train_on_batch(X, Y)

            # Update generator
            noise = np.array(rs(BATCH_SIZE))
            g_loss = cgan.train_on_batch(noise, [1]*BATCH_SIZE) # labels are ignored if Coulomb GAN
            print("epoch: %d, batch: %d, g_loss: %e, d_loss: %e, eps: %e" %
                  (epoch, index, g_loss, d_loss, dl.get_eps()))

        plot_generator_distribution(generator, epoch, GENERATED_IMAGE_PATH, rs)
        if eh.coulomb_gan:
            plot_potential(discriminator, "Epoch={0:d}".format(epoch),
                           GENERATED_IMAGE_PATH+"p{0:03d}".format(epoch))
        else:
            plot_real_or_fake(discriminator, "Epoch={0:d}".format(epoch),
                              GENERATED_IMAGE_PATH+"p{0:03d}".format(epoch))

if __name__ == '__main__':
    eh = ExperimentHistory()
    eh.batch_size = BATCH_SIZE
    eh.samples_per_gaussian = SAMPLES_PER_GAUSSIAN
    eh.random_with_normal_dist = False
    eh.X_train_is_shuffled = True
    eh.plummer_kernel_dim = 3.0
    eh.plummer_kernel_eps = 3.0
    eh.plummer_kernel_eps_half_life = 5000.0
    eh.plummer_kernel_ignore_self_potential = True
    eh.disc_Adam_decay = 1e-4
    eh.disc_Adam_lr = 1e-3
    eh.disc_Adam_beta_1 = 0.5
    eh.gen_Adam_decay = eh.disc_Adam_decay
    eh.gen_Adam_lr = eh.disc_Adam_lr
    eh.gen_Adam_beta_1 = eh.disc_Adam_beta_1
    eh.coulomb_gan = True
    if not eh.coulomb_gan:
        GENERATED_IMAGE_PATH = 'generated_images/'
    if not os.path.exists(GENERATED_IMAGE_PATH):
        os.mkdir(GENERATED_IMAGE_PATH)
    train(eh)
    print("finish")

2017/11/23

N次元の複数点間の距離を一括で計算 (その2)

N次元の複数点間の距離を一括で計算の計算は冗長でした。 np.tileを使う必要がありません。コードは以下のようになります。

# Make 3 points in 2D (embedded into 1x3 matrix)
p1 = np.array([[[1,2],[3,4],[5,6]]])

# Make 5 points in 2D (embedded into 5x1 matrix)
p2 = np.array([[[0,1],[1,2],[2,3],[1,1],[2,2]]]).reshape(5,1,-1)

print(np.sqrt(np.sum((p1 - p2)**2, axis=2)))
結果は、
[[ 1.41421356  4.24264069  7.07106781]
 [ 0.          2.82842712  5.65685425]
 [ 1.41421356  1.41421356  4.24264069]
 [ 1.          3.60555128  6.40312424]
 [ 1.          2.23606798  5.        ]]
となります。

2017/11/19

N次元の複数点間の距離を一括で計算

目的

Python numpyを使って、N次元空間内のK個の位置とL個の位置の間の距離を一括で計算します。 2重ループは使いません。

計算方法

最初にnumpyをimportします。

import numpy as np
N次元の位置を要素としてK個含む行ベクトルを作ります。ここでは、N=2, K=3としました。
# Make 3 points in 2D (embedded into 1x3 matrix)
p1 = np.array([[[1,2],[3,4],[5,6]]])
p1.shapeをprintすると、
p1.shape = (1, 3, 2)
となります。数式っぽく書くと、 \[ p_1 = \begin{pmatrix} (1, 2) && (3, 4) && (5, 6) \end{pmatrix} \] です。

同様に、N次元の位置を要素としてL個含む列ベクトルを作ります。ここでは、N=2, L=3としました。

# Make 5 points in 2D (embedded into 5x1 matrix)
p2 = np.array([[[0,1],[1,2],[2,3],[1,1],[2,2]]]).reshape(5, 1, -1)
p2.shapeをprintすると、
p2.shape = (5, 1, 2)
となります。数式っぽく書くと、 \[ p_1 = \begin{pmatrix} (0,1) \\ (1,2) \\ (2,3) \\ (1,1) \\ (2,2) \end{pmatrix} \] です。

行ベクトルp1を列方向に5個コピーして、5x3の行列m1を作ります。 同様に、列ベクトルp2を行方向に3個コピーして、5x3の行列m2を作ります。

# Make 5x3 matrix by copying elements
m1 = np.tile(p1, (5, 1, 1))
m2 = np.tile(p2, (1, 3, 1))
これを数式っぽく書くと、 \[ m_1 = \begin{pmatrix} (1, 2) && (3, 4) && (5, 6) \\ (1, 2) && (3, 4) && (5, 6) \\ (1, 2) && (3, 4) && (5, 6) \\ (1, 2) && (3, 4) && (5, 6) \\ (1, 2) && (3, 4) && (5, 6) \end{pmatrix},\ \ m_2 = \begin{pmatrix} (0,1) && (0,1) && (0,1) \\ (1,2) && (1,2) && (1,2) \\ (2,3) && (2,3) && (2,3) \\ (1,1) && (1,1) && (1,1) \\ (2,2) && (2,2) && (2,2) \end{pmatrix} \] となります。

p2に含まれる位置を始端、p1に含まれる位置を終端とするベクトルを全組合せについて一括で計算します。

# Make vectors from p2 to p1
mv = m1 - m2
計算結果は、 \[ M_v = \begin{pmatrix} (1,1)&&(3,3)&&(5,5)\\ (0,0)&&(2,2)&&(4,4)\\ (-1,-1)&&(1,1)&&(3,3)\\ (0,1)&&(2,3)&&(4,5)\\ (-1,0)&&(1,2)&&(3,4)\end{pmatrix} \] となります。

各要素について二乗距離を計算します。

# Calculate square distances
ms = np.sum(np.square(mv), axis=2)
計算結果は、
ms = [[ 2 18 50]
 [ 0  8 32]
 [ 2  2 18]
 [ 1 13 41]
 [ 1  5 25]]
数式で書くと、 \[ M_s = \begin{pmatrix} 2&&18&&50\\ 0&&8&&32\\ 2&&2&&18\\ 1&&13&&41\\ 1&&5&&25\end{pmatrix} \] となります。

あとは各要素の平方根を求めれば、各位置間の距離が得られます。

# Calculate distances
md = np.sqrt(ms)
計算結果は、
md = [[ 1.41421356  4.24264069  7.07106781]
 [ 0.          2.82842712  5.65685425]
 [ 1.41421356  1.41421356  4.24264069]
 [ 1.          3.60555128  6.40312424]
 [ 1.          2.23606798  5.        ]]
となります。数式で書くと、 \[ M_d = \begin{pmatrix} 1.41421356 && 4.24264069 && 7.07106781 \\ 0 && 2.82842712 && 5.65685425 \\ 1.41421356 && 1.41421356 && 4.24264069 \\ 1 && 3.60555128 && 6.40312424 \\ 1 && 2.23606798 && 5 \end{pmatrix} \] です。

コードまとめ

上記のコードをまとめておきます。

import numpy as np

# Make 3 points in 2D (embedded into 1x3 matrix)
p1 = np.array([[[1,2],[3,4],[5,6]]])

# Make 5 points in 2D (embedded into 5x1 matrix)
p2 = np.array([[[0,1],[1,2],[2,3],[1,1],[2,2]]]).reshape(5, 1, -1)

# Make 5x3 matrix by copying elements
m1 = np.tile(p1, (5, 1, 1))
m2 = np.tile(p2, (1, 3, 1))

# Make vectors from p2 to p1
mv = m1 - m2

# Calculate square distances
ms = np.sum(np.square(mv), axis=2)

# Calculate distances
md = np.sqrt(ms)

print("p1.shape =", p1.shape)
print("p2.shape =", p2.shape)
print("ms =", ms)
print("md =", md)

2017/11/16

Boston Dynamics の Meet Spot の紹介動画

四足歩行+アーム(首と顔っぽいが)のMeet Spotと名づけられたロボットの紹介動画。
https://www.youtube.com/watch?v=AO4In7d6X-c

動画自体は2017年8月14日にYoutubeに公開されている。

TEDでCEOのMarc Raibert氏がデモと解説していて、
段差を移動したり、操作なしで飲み物を取ってきたり、
2足でジャンプしたりできる。

動きが、なんというかガッションガッションと動く昔ながらのロボットではなく、
スムーズ。動物の動きにどんどん近づいているようだ。

2017/10/29

CRISPR (その1)

クリスピーじゃなくてCRISPR-Cas9(クリスパーキャスナイン)という遺伝子編集技術をつかった様々な結果をこれまで目にしていたのですが、なにか良く分からないが遺伝子を編集できるんだな、という以上の知識はありませんでした。偶然、書店で平積みされていた「CRISPR 究極の遺伝子編集技術の発見(ISBN978-4-16-390738-3)」というタイトルのハードカバーの書籍を見つけ、CRISPR-Cas9を見つけた主要な科学者が著者だったので、購入。翻訳本なので、1年くらい前に元の書籍は発売されていた模様。

全体が物語風味で構成されていますが、専門家ではない人にも技術の内容が分かりやすい構成になっています。第3章までが基礎技術部分の解説で、CRISPR-Cas9がDNAを目的の位置で正確に切断する原理が説明されています。第4章は、CRISPR-Cas9の切断能力をベースにした遺伝子編集技術の解説、それ以降は、すでに実施された様々な応用例と、倫理的な問題についての議論が書かれています。

以下はこの書籍を読んでいる最中にとった専門外の人間のメモなので、そのまま信じないように。ググって見つけたものを真に受けてはいけません。

遺伝子関連の言葉

まず、
ゲノム∋遺伝子∋塩基対
という関係があります。また、
塩基={A, G, C, T}≒ヌクレオチド
塩基対={A=T, C=G}
ヌクレオチドの列=遺伝子=タンパク質=アミノ酸の列
ゲノムの大きさ=塩基対の長さ
とのこと。塩基に色々くっつくとヌクレオチドになるらしい。
遺伝子が翻訳されてタンパク質になるので、1対1に対応している。

遺伝子のサイズ

遺伝子の大体のサイズが記載されていたので、表にしてみました。

名前 遺伝子の数 塩基対の数
ウイルス 数個 数千個
細菌 約4,000個 数百万個
ハエ 1万4千個 数億個
ヒト・ネズミ・蛙 21,000個(*1) 32億個 (*2)
サラマンダー (*1)の10倍 (*2)の10倍
一部の植物 (*1)の100倍 (*2)の100倍
ミトコンドリア ? 16,000個
(*2)このうちの8%(2億5千万個)が人類の祖先がレトロウイルスに感染した名残

おおよそ、ひとつの遺伝子は数千個~1万数千個の塩基対で構成されているようです。アミノ酸が何個の塩基対で作られるかは、この書籍には記載がありませんでした。本題からはずれるのでなくても不思議ではありませんが。


さて、ここからCRISPRに行き着くまでが遠いのですが、遺伝子治療と呼ばれていた遺伝子を挿入する技術とCRISPR-Cas9を使った編集技術は全く異なり、少なくとも細菌の遺伝子を書き換えるときの正確性は極めて高い技術だということです。ただし、人間やサルやマウスの場合は、CRISPR-Cas9を個々の細胞にまで届けないといけないので、色々課題が残っているようです。

MNISTでのmixup法の効果

数日前にarXivに投稿されていたmixup法(https://arxiv.org/abs/1710.09412)をMNISTで試してみました。 手法自体は簡単で、ニューラルネットワークの学習時に使う入力と出力を2つずつランダムに選び、それを比率\(\lambda\)で混ぜるだけです。式で書くと、ランダムに選んだ入力と出力を\((x_i,y_i), (x_j,y_j)\)のように2つ選ぶとき、学習時に渡す入力\(\tilde{x}\)と出力\(\tilde{y}\)はそれぞれ、 \[ \tilde{x} = \lambda x_i + (1-\lambda) x_j \] \[ \tilde{y} = \lambda y_i + (1-\lambda) y_j \] となります。\(\lambda\)はサンプルごとにベータ分布からランダムに決めます。

今回はMNISTで試すので、数字が半透明で重なります。ベータ分布の\(\alpha=\beta=2\)で32個生成した結果が下図です。

\(\alpha=\beta\)の下で、\(\alpha\)を色々変えて学習した結果が下表です。学習は30エポックのみです。

\(\alpha\) Train acc Validation acc
n/a 0.999933 0.991600
0.1 0.999767 0.991800
0.2 0.999483 0.990600
0.3 0.999283 0.990300
0.4 0.998983 0.991600
0.5 0.998633 0.989500
1 0.997483 0.989700
2 0.996383 0.990000

1行目は通常の方法で混ぜずに学習させた場合の結果です。Train accの計算には混ぜていない学習用データを用いています。\(\alpha\)を増やすと半々で混ざるサンプルの確率が増え、Train accが徐々に下がっていきます。一方、Validation accは下がっているようにも見えますが、\(\alpha=0.1\)の場合は若干上がっています。差が小さすぎるので誤差かもしれません。

論文に書かれているような効果は今回のMNIST実験では見られませんでしたが、こういう混ぜ方をしても学習ができることは確認できました。

以下、今回の実験で使ったコードです。122行目のeh.alphaの値を変えることで\(\alpha\)を変更します。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# -*- coding: utf-8 -*-
import sys,os,math
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../keras-examples/src')
import numpy as np
import keras
from keras.models import Sequential, Model, Input, model_from_json
from keras.layers import Dense, Activation, Reshape, Flatten, Dropout
from keras.layers.normalization import BatchNormalization
from keras.layers.convolutional import Conv2D
from keras.layers.advanced_activations import LeakyReLU
from keras.utils import to_categorical
from util.history import ExperimentHistory
from keras.datasets import mnist
from PIL import Image, ImageDraw, ImageFont

TRAIN_EPOCH = 30
GENERATED_IMAGE_PATH = 'generated_images/'
FONT_DIR = "C:/Windows/Fonts"

# Image size is the same as MNIST
IMAGE_WIDTH = 28
IMAGE_HEIGHT = 28

def make_model():
    model = Sequential()
    model.add(Conv2D(64, (5, 5), strides=(2, 2), padding='same',
                     input_shape=(IMAGE_WIDTH, IMAGE_HEIGHT, 1), 
                     data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    model.add(Conv2D(128, (5, 5), strides=(2, 2), data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    model.add(Flatten())
    model.add(Dense(256))
    model.add(LeakyReLU(0.2))
    model.add(Dropout(0.5))
    model.add(Dense(10))
    model.add(Activation('softmax'))
    print(model.summary())
    return model

def combine_images(generated_images, margin=20):
    total = len(generated_images)
    cols = int(math.sqrt(total))
    rows = math.ceil(float(total)/cols)
    width, height = generated_images[0].shape[0:2]
    print(generated_images[0].shape)
    combined_image = np.ones(((height+margin)*rows, width*cols), dtype=generated_images[0].dtype)*-0.5
    for index, image in enumerate(generated_images):
        i = int(index/cols)
        j = index % cols
        combined_image[(height+margin)*i:(height+margin)*i+height, width*j:width*(j+1)] = image[:, :, 0]
    return combined_image, cols, rows, width, height, margin

# Data generator
def data_gen(x, y, batch_size, alpha, num_classes):
    # Generate data eternally
    assert len(x) == len(y)
    epoch = 0
    index = 0
    font = ImageFont.truetype(FONT_DIR+"/{}".format("Tahoma.ttf"), 10) # For drawing numbers
    while True:
        index += 1
        data = []
        labels = []
        for i in range(batch_size):
            i1 = np.random.randint(0, len(x)) 
            i2 = np.random.randint(0, len(x)) 
            lam = np.random.beta(alpha, alpha)
            # Mixup two samples
            xm = lam * x[i1] + (1.0 - lam) * x[i2]
            ym = lam * y[i1] + (1.0 - lam) * y[i2]
            data.append(xm)
            labels.append(ym)

        ## Draw mixed-up images
        # image, cols, rows, width, height, margin = combine_images(data)
        # image = image*127.5 + 127.5
        # if not os.path.exists(GENERATED_IMAGE_PATH):
        #     os.mkdir(GENERATED_IMAGE_PATH)
        # img = Image.fromarray(image.astype(np.uint8))
        # draw = ImageDraw.Draw(img)
        # for i in range(len(data)):
        #     yi = int(i/cols)
        #     xi = i % cols
        # img.save(GENERATED_IMAGE_PATH+"%04d_%04d.png" % (epoch, index))

        yield np.array(data), np.array(labels).reshape((batch_size, num_classes))

def run(eh):
    # Prepare MNIST data
    (x_train, y_train), (x_val, y_val) = mnist.load_data()
    x_val = (x_val.astype(np.float32) - 127.5)/127.5
    x_val = x_val.reshape(x_val.shape[0], x_val.shape[1], x_val.shape[2], 1)
    y_val = keras.utils.to_categorical(y_val, num_classes=10)

    x_train = (x_train.astype(np.float32) - 127.5)/127.5
    x_train = x_train.reshape(x_train.shape[0], x_train.shape[1], x_train.shape[2], 1)
    y_train = keras.utils.to_categorical(y_train, num_classes=10)

    # Make a model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model = make_model()
    model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy'])

    # Train the model
    dg = data_gen(x_train, y_train, eh.batch_size, eh.alpha, 10)
    hist = model.fit_generator(dg, steps_per_epoch=len(x_train)/eh.batch_size, epochs=TRAIN_EPOCH,
                               validation_data=(x_val, y_val)).history

    # Evaluate the model by train and validation data
    score_train = model.evaluate(x_train, y_train, batch_size=eh.batch_size)
    score_val = model.evaluate(x_val, y_val, batch_size=eh.batch_size)

    comments = "#Final model loss_train={0:10.6f} acc_train={1:10.6f} loss_val={2:10.6f} acc_val={3:10.6f}"\
               .format(score_train[0], score_train[1], score_val[0], score_val[1])
    eh.write("history.log", model, opt, hist, comments)

if __name__ == '__main__':
    eh = ExperimentHistory()
    eh.batch_size = 32
    eh.alpha = 2
    run(eh)

2017/10/15

機械学習における著作権

投稿時点(2017/10/15)の「日本の法律」では、機械学習に放り込む目的で著作物を使う場合、著作権が制限される条件(著作権法47条の7)があるため、商用・非商用にかかわらず自由に利用することができるようだ。

ITmediaの記事
http://www.itmedia.co.jp/news/articles/1710/10/news040.html

元の記事(上野達弘さんの投稿)
http://rclip.jp/2017/09/09/201708column/

IJCAIでの発表の記事
https://www.businessinsider.jp/post-102878


法律を引用すると
-----------------
http://www.houko.com/00/01/S45/048.HTM
(情報解析のための複製等)
第四七条の七 著作物は、電子計算機による情報解析(多数の著作物その他の大量の情報から、当該情報を構成する言語、音、影像その他の要素に係る情報を抽出し、比較、分類その他の統計的な解析を行うことをいう。以下この条において同じ。)を行うことを目的とする場合には、必要と認められる限度において、記録媒体への記録又は翻案(これにより創作した二次的著作物の記録を含む。)を行うことができる。ただし、情報解析を行う者の用に供するために作成されたデータベースの著作物については、この限りでない。
-----------------
とのことなので、別の法律による制限がないのであれば、確かに機械学習用のデータ利用では著作権の制限を受けないように読める。この規定の例外はデータベースで、おそらく、販売されているデータベース(新聞のコーパスやLDC(https://www.ldc.upenn.edu/language-resources))については、きちんと購入してライセンスを受けてください、ということなのだろう。

2017/10/09

加速度センサーで行動分類(CNN編)

前回は加速度センサーデータをLSTMで分類しましたが、今回はCNNで試してみます。

センサーデータは、LSTMのときと同様にランダムに窓で区切ります。ただし、そのまま入力はせず、3次元のヒストグラムを作り、それをCNNの入力とします。センサーデータは3軸あり、振幅は整数でその範囲は[0,64)です。単純に、64x64x64の3次元配列でヒストグラムを作っても良いのですが、大部分が0になるスカスカのヒストグラムになる点と、CNNの規模が大きくなる点を考慮して、振幅を4分の1にしてヒストグラムを作ります。つまり、16x16x16の3次元配列(テンソル)をCNNの入力にします。これでも要素の数は4096個あり、窓幅を64とするとスカスカです。

使用するCNNのモデルは、3Dの畳み込み層3つと全結合層2つです。具体的には、

BIN_SIZE = 4
si = int(64/BIN_SIZE)
model = Sequential()
model.add(Conv3D(32, kernel_size=(5, 5, 5),
                 input_shape=(si, si, si, 1), 
                 data_format='channels_last'))
model.add(LeakyReLU(0.2))
for _ in range(2):
    model.add(Conv3D(32, kernel_size=(5, 5, 5), data_format='channels_last'))
    model.add(LeakyReLU(0.2))
model.add(Flatten())
model.add(Dense(256))
model.add(LeakyReLU(0.2))
model.add(Dropout(0.5))
model.add(Dense(num_classes))
model.add(Activation('softmax'))
の通りです。Conv3Dのフィルタサイズは32としました。Conv3Dの出力は順に、12x12x12x32、8x8x8x32、4x4x4x32です。したがって、1つ目のDenseへの入力サイズは2048となります。Conv3Dで適度なサイズに絞らないと、全結合層のパラメータ数が増えすぎて、GPUで処理できなくなります。

3次元のヒストグラムはデータを増やすためにランダムに回転させます。回転範囲の制限は、回転行列の対角要素の全てが一定値θより大きいという条件を満たす回転行列のみを使用することにより実現します。

前回と同様、バッチサイズ32、1エポックあたりのステップ数は100で、100エポック学習させました。結果は下表の通りです。100エポック終了後のモデルで訓練、検証ともに32000サンプル生成し、評価しています。32000サンプルはランダムに選択しているので毎回基準は揺れます。

Train accuracy





Window width No rotation Random rotation (θ)
0.95 0.9 0.5 0 Full
32 83.7% 57.5% 52.9% 40.2% 34.8% 29.7%
64 93.1% 70.4% 64.4% 49.8% 41.6% 35.9%
128 98.6% 85.5% 80.1% 63.4% 54.6% 47.6%
256 99.9% 95.0% 91.7% 79.5% 70.0% 61.8%
512 100.0% 98.7% 96.9% 89.6% 83.1% 75.6%

Validation accuracy
Window width No rotation Random rotation (θ)
0.95 0.9 0.5 0 Full
32 51.7% 51.2% 50.2% 38.5% 33.5% 25.3%
64 55.9% 59.9% 57.2% 46.1% 36.2% 31.9%
128 65.0% 70.8% 70.3% 58.6% 50.5% 37.3%
256 74.8% 80.6% 82.2% 77.1% 63.1% 51.9%
512 78.4% 80.7% 82.3% 83.1% 71.4% 65.0%

窓幅(Windows width)が32の場合、1秒のデータだけを見て分類していることになります。表から分かるように、入力データのランダム回転の範囲にかかわらず、窓幅が広いほど精度が高くなります。これは、窓幅が広くなるほど入力データのバリエーションが減り、分類に使えるデータが増えるためです。

入力データのランダムな回転の有効性は、窓幅によって変化します。窓幅が狭いと回転範囲が狭い(θが大きい)ほうが検証用データでの精度が良く、窓幅が広くなるにつれて精度が最も良くなる回転範囲が広めになります。ただし、どの窓幅でも回転範囲を広げすぎると精度が下がります。

今回の実験に使ったコードは次の通りです。util.random_rotateとutil.historyはここにあります。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
# This script trains a neural network model using accelerometer data
import sys,os
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + '/../keras-examples/src')

import keras, datetime
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation, Conv3D, LeakyReLU, Flatten, Dropout
from keras.layers import LSTM
from util.random_rotate import uniform_random_rotation_matrix_3D
from util.history import ExperimentHistory

TRAIN_EPOCH = 100
BIN_SIZE = 4

def make_model(num_classes):
    si = int(64/BIN_SIZE) # Size of input
    model = Sequential()
    model.add(Conv3D(32, kernel_size=(5, 5, 5),
                     input_shape=(si, si, si, 1), 
                     data_format='channels_last'))
    model.add(LeakyReLU(0.2))
    for _ in range(2):
        model.add(Conv3D(32, kernel_size=(5, 5, 5), data_format='channels_last'))
        model.add(LeakyReLU(0.2))
    model.add(Flatten())
    model.add(Dense(256))
    model.add(LeakyReLU(0.2))
    model.add(Dropout(0.5))
    model.add(Dense(num_classes))
    model.add(Activation('softmax'))
    model.summary()
    return model

def read_accel(filename):
    with open(filename) as f:
        data = []
        for line in f:
            for v in line.rstrip().split(" "):
                if not v.isdigit():
                    print(v, "is not digit")
            data.append([int(v) for v in line.rstrip().split(" ")])
        data = np.array(data)
        return data

# Data generator
def data_gen(filelist, batch_size, window_width, random_rotate=False, rotation_range=0):
    # Make a dict for converting a class ID to files
    num_class = 0
    with open(filelist) as f:
        id2f = dict()
        for line in f:
            id, filename = line.split(" ", maxsplit=2)
            id = int(id)
            id2f.setdefault(id, []).append(read_accel(filename.rstrip()))
            num_class = max(num_class, id)
        num_class += 1
    yield num_class # Return the number of classes at first

    # Generate data eternally
    shift = np.array([32, 32, 32]) # Specify the center of rotation
    while True:
        data = []
        labels = []
        for i in range(batch_size):
            while True:
                c = np.random.randint(0, num_class) # Select a class
                if c in id2f:
                    break
            # Select a file in the class
            f = np.random.randint(0, len(id2f[c]))
            # Select a window position
            p = np.random.randint(0, max(0, len(id2f[c][f])-window_width)+1)
            # Make a label for the selected class
            labels.append(keras.utils.to_categorical(c, num_classes=num_class))
            # Add selected data
            if(random_rotate):
                while True:
                    m = uniform_random_rotation_matrix_3D()
                    if m[0,0] > rotation_range and m[1,1] > rotation_range and m[2,2] > rotation_range:
                        break
                data.append(np.dot(m, (id2f[c][f][p:p+window_width,:]-shift).T).T+shift)
            else:
                data.append(id2f[c][f][p:p+window_width,:])
                
            if(len(data[-1]) < window_width):
                padding = window_width - len(data[-1])
                data[-1] = np.append(np.zeros((padding, data[-1].shape[1]), dtype=np.int),
                                              data[-1], axis=0)
        # Make 3D-histogram (x, y and z range are 0-63)
        hist_data = []
        for i in range(batch_size):
            h = np.zeros(shape=(int(64/BIN_SIZE), int(64/BIN_SIZE), int(64/BIN_SIZE), 1))
            for t in range(data[i].shape[0]):
                if data[i][t,0] > 63 or data[i][t,0] < 0 or \
                   data[i][t,1] > 63 or data[i][t,1] < 0 or \
                   data[i][t,2] > 63 or data[i][t,2] < 0:
                    continue
                x = int(data[i][t,0]/BIN_SIZE)
                y = int(data[i][t,1]/BIN_SIZE)
                z = int(data[i][t,2]/BIN_SIZE)
                h[x, y, z, 0] += 1
            hist_data.append(h)
        yield np.array(hist_data), np.array(labels).reshape((batch_size, num_class))

def get_data(data_generator, steps):
    dgr = [next(data_generator) for _ in range(steps)]
    data = np.array([w for v in dgr for w in v[0]])
    labels = np.array([w for v in dgr for w in v[1]])
    return data, labels

def run(eh):
    dg = data_gen("train.list", eh.batch_size, eh.window_width, eh.random_rotate, eh.rotation_range)
    num_classes = next(dg)
    dgv = data_gen("val.list", eh.batch_size, eh.window_width)
    assert num_classes == next(dgv)

    # Make a model
    opt = keras.optimizers.Adamax(lr=0.002, beta_1=0.9, beta_2=0.999,
                                  epsilon=1e-08, decay=1e-4)
    model = make_model(num_classes)
    model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy'])

    # Train the model
    hist = model.fit_generator(dg, steps_per_epoch=100, 
                               validation_data=dgv, validation_steps=100,
                               epochs=TRAIN_EPOCH).history
    # Get data for evaluation
    score_train = model.evaluate(*get_data(dg, 1000))
    score_val = model.evaluate(*get_data(dgv, 1000))
    comments = "#Final model loss_train={0:10.6f} acc_train={1:10.6f} loss_val={2:10.6f} acc_val={3:10.6f}"\
               .format(score_train[0], score_train[1], score_val[0], score_val[1])
    eh.write("history.log", model, opt, hist, comments)

    # Evaluate the last model by train and validation data
    model.save("models/cnn_b{0}_ww{1}_rr{2}.hdf5"\
               .format(BIN_SIZE, eh.window_width, str(eh.rotation_range) if eh.random_rotate else "no"))

if __name__ == '__main__':
    if not os.path.exists("models"):
        os.mkdir("models")

    eh = ExperimentHistory()
    eh.BIN_SIZE = BIN_SIZE
    eh.batch_size = 32
    for (rr, rrange) in [(False, 0), (True, 0.95), (True, 0.9), (True, 0.5), (True, 0), (True, -2)]:
        eh.random_rotate = rr
        eh.rotation_range = rrange
        for ww in [32, 64, 128, 256, 512]:
            eh.window_width = ww
            run(eh)