2018/04/30

CppUTest を試してみる

CppUTestを試してみます。CppUTestの本家はここ。
http://cpputest.github.io/

現時点の最新をダウンロードします。

$ wget https://github.com/cpputest/cpputest/releases/download/v3.8/cpputest-3.8.zip
展開して、
$ unzip cpputest-3.8.zip
$ cd cpputest-3.8/cpputest_build/
ビルドするだけです。なお、cmake版もあります。
$ ../configure
$ make
一応、CppUTest自身のテストもしておきます。
$ make tdd
以下のようなコードを書いて、
#include "CppUTest/CommandLineTestRunner.h"
#include "CppUTest/TestHarness.h"

TEST_GROUP(TestSuite){};
TEST(TestSuite, Case1)
{
  LONGS_EQUAL(3, 4);
}

int main(int argc, char *argv[])
{
  return RUN_ALL_TESTS(argc, argv);
}
ビルドします。
$ g++ -Icpputest-3.8/include test.cpp cpputest-3.8/cpputest_build/lib/libCppUTest.a
実行すると、
$ ./a.out
test.cpp:7: error: Failure in TEST(TestSuite, Case1)
        expected <3 0x3>
        but was  <4 0x4>

.
Errors (1 failures, 1 tests, 1 ran, 1 checks, 0 ignored, 0 filtered out, 3 ms)
となり、無事、エラーが検出できました。


なお、次のようにファイルの指定順序をうっかり間違うと、リンカエラーがでます。
$ g++ -Icpputest-3.8/include cpputest-3.8/cpputest_build/lib/libCppUTest.a test.cpp
理由を知りたい方は、
http://www.yunabe.jp/docs/static_library.html
をどうぞ。

2018/04/16

映像を利用した音声分離

複数の音声が混じった音から個別の音声を抜き出す音声分離の性能を映像を使って向上させる方法をgoogleが発表。

評価には、signal-to-distortion ratio (SDR)を利用。
arXivに投稿されている論文のTable 3によると、音声のみを用いた最新手法より分離性能が高い。


arXiv
https://arxiv.org/abs/1804.03619

GIGAZINE
https://gigazine.net/news/20180412-looking-to-listen-google/

Google Research Blog
https://research.googleblog.com/2018/04/looking-to-listen-audio-visual-speech.html


SDRの解説は
https://hal.inria.fr/inria-00564760/document
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.412.3918&rep=rep1&type=pdf
に記載があるが、きちんとは読んでいない。

黙読?の認識

MITが声を出さずに頭の中でしゃべった声を認識する装置を開発したとのこと。

口を動かしたつもりがないのに、実際には微弱ながらも筋肉が反応しており
それを捕らえて認識している模様。CNNベースの孤立単語認識で実現している。

深層学習のおかげで、脳内でしゃべった単語に対応する電圧(電流?)が
筋電計から取得できることさえ分かってしまえば、あとは装置を作ってデータを集めるだけ。
EEGやNIRSに比べれば装置は作りやすく、リアルタイム性も良さそう。



論文
http://dam-prod.media.mit.edu/x/2018/03/23/p43-kapur_BRjFwE6.pdf

スラド
https://science.srad.jp/story/18/04/10/0326229/

GIGAZINE
https://gigazine.net/news/20180405-mit-internal-verbalization-transcribe-computer-system/

2018/04/08

PulseAudioで音量表示

はじめに


PulseAudioを使って録音しながらリアルタイムで音量を表示します。

方法とコード


PulseAudioのSimple API [1] を使って録音し、音量を表示します。音量の計算はざっくりです。20秒間録音して終了します。
 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
#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <pulse/simple.h>

std::string getbar(double v)
{
  v = 10*(v-4); // 10 and -4 are determined by heuristic
  if(v < 0){v = 0;}
  return std::string(v, '*');
}

int main(void)
{
  pa_sample_spec ss;
  ss.format = PA_SAMPLE_S16LE; // Signed, 16bits, Little endian
  ss.channels = 1;
  ss.rate = 8000; // Sampling frequency is 8kHz
  pa_simple *s = pa_simple_new(NULL, "Volume meter", PA_STREAM_RECORD, NULL, "Microphone vm",
          &ss, NULL, NULL, NULL);
  std::vector<short> buffer(800); // For 100ms
  for(int i=0; i<200; ++i) // Record for 20s
  {
    // Read audio for 100ms
    pa_simple_read(s, &buffer[0], buffer.size()*sizeof(short), nullptr);

    // Calculate power for the range of 100ms
    double power = 0;
    for(auto &v : buffer)
    {
      power += v*v;
    }
    // Output a result with a bar
    double logp = log10(power/buffer.size());
    std::cout << logp << "\t" << getbar(logp) << std::endl;
  }
  pa_simple_free(s);
  return 0;
}
このコードのファイル名をrec.cppとすると、コンパイルは、
g++ --std=c++14 -lpulse-simple rec.cpp
で行えます。g++はバージョン4.9.2を利用しています。libpulse-dev が必要ですので、コンパイルする前にapt-getやaptitudeを使ってインストールしておきましょう。

マイクがリモートにある場合は、リモート側でPulseAudioのサーバの設定を行い、

PULSE_SERVER=192.168.1.2 ./a.out
のようにして実行します。ここでは、リモートマシンのIPアドレスが192.168.1.2であるとしています。 リモートマシンの設定は前の記事を参照ください。

結果


実行すると録音が始まり、コンソールに音量がバーで表示されます。
5.09851 **********
4.93314 *********
4.76426 *******
4.96767 *********
6.26456 **********************
6.17523 *********************
5.34311 *************
6.9875 *****************************
4.98664 *********
7.23169 ********************************
5.60815 ****************
6.16099 *********************
5.55161 ***************
5.16174 ***********
5.62978 ****************
5.19768 ***********
5.73633 *****************
5.85854 ******************
6.54446 *************************
5.53201 ***************
5.42636 **************
5.40482 **************
5.77586 *****************
4.9568 *********
4.82601 ********
7.56319 ***********************************
8.19005 *****************************************

参考


[1] https://freedesktop.org/software/pulseaudio/doxygen/simple.html

2018/04/07

PulseAudioを使った音の転送

はじめに


ネットワーク経由で音の再生・録音をしてみます。OSはLinuxです。遠いところに置いた省電力のデバイスから簡単に音を出せます。無線でつなげば、どこでも音楽が聴けるようになりますね。

構成


ここでは、音の再生・録音の操作をする計算機をクライアント、実際にスピーカやマイクがつながっている計算機をサーバとします。

クライアントは、VirualBox内で稼動しているDebian 8.1です。
サーバは、Raspberry Pi 3 Model Bです。こちらのOSはDebian 8.0です。

音の転送には、PulseAudio 5.0を用います。ALSAのひとつ上の層で動くサウンドサーバです。音をファイルから再生する場合は、

[クライアント]→[サーバ]→[スピーカ]→音
というふうにデータが流れ、録音する場合は、
音→[マイク]→[サーバ]→[クライアント]
のように流れます。

サーバの設定


まず、設定ファイルをコピーします。
cp /etc/pulse/default.pa ~/.config/pulse/
この中から、コメントアウトされているmodule-native-protocol-tcpの項目を見つけてきて、コメントをはずし、次のようにオプションを追記します。
load-module module-native-protocol-tcp auth-ip-acl=127.0.0.1;192.168.1.0/24
auth-ip-acl=の右側はアクセスを許可するクライアントのIPアドレスです。複数指定する場合は、;で区切ります。マスクも使えます。

設定を反映するために、PulseAudioを再起動します。

$ pulseaudio --kill
$ pulseaudio --start
これで、サーバ側の準備は終わりました。

クライアントから操作


動作テストに使うwavファイルを適当に入手します。例えば、
$ wget https://freewavesamples.com/files/1980s-Casio-Piano-C5.wav

毎回サーバを指定


環境変数PULSE_SERVERを指定するだけで、サーバのスピーカから音を出力することができます。 サーバ(Raspberry Pi)のアドレスが192.168.1.2であれば、
PULSE_SERVER=192.168.1.2 aplay 1980s-Casio-Piano-C5.wav
とすると、サーバにつなげたスピーカから音が出ます。

同じように、録音もできます。

PULSE_SERVER=192.168.1.2 arecord test.wav

デフォルトサーバの指定


環境変数を毎回指定するのが手間な場合は、
export PULSE_SERVER=192.168.1.2
とします。または、client.confを設定します。具体的には、
cp /etc/pulse/client.conf ~/.config/pulse/
のように、クライアント用の設定ファイルをコピーし、コメントアウトされているdefault-serverを
default-server = 192.168.1.2
のように有効化します。

なお、これらの設定をするとクライアントにしている計算機でPulseAudioのサーバを起動できなくなります。例えば、

$ pulseaudio --start
N: [pulseaudio] main.c: ユーザーが設定したサーバー 192.168.1.2 は start/autospawn を拒否しています。
というエラーが出るようになります。

参考


[1] http://nishio.hateblo.jp/entry/20111215/1323962652
[2] https://wiki.archlinux.jp/index.php/PulseAudio/%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB
[3] https://qiita.com/tukiyo3/items/d39c90f91d782001e7d7
[4] http://d.hatena.ne.jp/kakurasan/20081212/p2

2018/04/01

量子機械学習

半年前の記事ですが、Quantum machine learning についてのレビュー記事があったので、メモ。

https://www.nature.com/articles/nature23474
(doi:10.1038/nature23474)

普通の機械学習の量子版が色々と考えられているようです。例えば、
- Quantum principal component analysis (PCAの量子版)
- Quantum support vector machines (SVMの量子版)
- Deep quantum learning (深層学習の量子版)
- Quantum reinforcement learning (強化学習の量子版)
があるようです。

2018/03/31

VirtualBox+cygwinで画面をWindows側に

目的


LinuxのウィンドウをWindowsの画面として表示させます。

環境


仮想マシンにはVirtualBoxを使います。ここではホストにWindows、ゲストにDebianを用います。

X Window Systemには、VcXsrv (https://sourceforge.net/projects/vcxsrv/) を使います。また、ssh接続のためにcygwinも用います。

方法


環境が準備できたら、以下のコマンドで簡単にWindows側でウィンドウを表示できます。 以下の例は、Visual Studio Codeを起動する例です。

まず、ゲストOSにsshでログインします。

ssh -R 6000:localhost:6000 -p [port] username@localhost
次に、ディスプレイの表示先を変更し、ウィンドウを使うアプリケーションを起動します。
export DISPLAY=localhost:0.0
code --disable-gpu
[port]はVirtualBoxのsshのポートフォワーディングルールとして設定しているホストポートの番号を入れます。具体的には、VirutalBoxの仮想マシンごとの設定の[ネットワーク]→[高度]→[ポートフォワーディング]で、
設定項目
名前ssh
プロトコルTCP
ホストIP(空欄)
ホストポート3222
ゲストIP(空欄)
ゲストポート22
としていれば、[port]は3222になります。6000はX Window Systemが使うデフォルトのポートです。

GPUによるアクセラレーションは使えないので、Visual Studio Codeはgpuを使用しないモードで起動する必要があります。ちょっと描画が遅いのが気になりますが、一応、表示や操作はできます。

2018/02/10

Fréchet Inception Distance のサンプル数依存性

はじめに

少し前に試したFréchet Inception Distance (FID)ですが、[1]によると、サンプル数依存性があるとのことです。今回は、本当にサンプル数依存性、つまり、距離計算に使う画像の数にFIDが影響を受けるのか確認してみます。

実験

例によって、手書き数字のMNISTで試してみます。訓練用の画像とテスト用の画像の距離を計算します。結果は、下表のようになりました。訓練用の画像の数とテスト用の画像の数は同じで、それぞれ表の通りの数を先頭から順に選んでいます。

# of imagesFID
30007.05
40005.51
50004.73
60003.76
70002.87
80002.26
90001.96
100001.63
減ってる!

グラフにしてみると、

となりました。

まとめ

というわけで、[1]のFigure 1 (b)と同じような結果が得られました。画像集合間の近さをFIDで比較するときは、数をあわせましょう。

参考文献

[1] Mikołaj Bińkowski, Dougal J. Sutherland, Michael Arbel, Arthur Gretton, "Demystifying MMD GANs," 2018, https://arxiv.org/abs/1801.01401

2018/01/27

Coulomb GANでMNIST

はじめに


Coulomb GANでMNISTの手書き数字を生成してみます。生成した画像の良さはFréchet Inception Distance (FID)を使ってとりあえずは測定できるので、普通のGANと比較してみます。

生成器と識別器


生成器と識別器は、以前の利用したものとほぼ同じものを利用します。 Kerasのコードで、生成器は、
g = Sequential()
g.add(Dense(1024, input_dim=100))
g.add(BatchNormalization())
g.add(Activation('relu'))
g.add(Dense(128*7*7))
g.add(BatchNormalization())
g.add(Activation('relu'))
g.add(Reshape((128, 7, 7), input_shape=(128*7*7,)))
g.add(UpSampling2D((2, 2), data_format='channels_first'))
g.add(Conv2D(64, (5, 5), padding='same', data_format='channels_first'))
g.add(BatchNormalization())
g.add(Activation('relu'))
g.add(UpSampling2D((2, 2), data_format='channels_first'))
g.add(Conv2D(1, (5, 5), padding='same', data_format='channels_first'))
g.add(Activation('tanh'))
識別器は、
d = Sequential()
d.add(Conv2D(64, (5, 5), strides=(2, 2), padding='same', input_shape=(1, 28, 28), 
                 data_format='channels_first'))
d.add(LeakyReLU(0.2))
d.add(Conv2D(128, (5, 5), strides=(2, 2), data_format='channels_first'))
d.add(LeakyReLU(0.2))
d.add(Flatten())
d.add(Dense(256))
d.add(LeakyReLU(0.2))
d.add(Dropout(0.5))
d.add(Dense(1))
を使います。最後のsigmoidは省略しています。Coulomb GANはポテンシャルを模擬するので、sigmoidがあると模擬ができません。一方、普通のGANでは真偽を1と0で表すのでsigmoidがあったほうが良さそうですが、試してみたところ問題なく学習できました。

実験結果


さっそく、Coulomb GANと普通のGANを比較してみます。結果は下表の通りです。 FIDの比較対象は全てMNISTの訓練用のデータです。Coulomb GANは自身のポテンシャルを無視する設定で動作させました。また、Plummer kernelの\(\varepsilon\)の半減期は5000、次元は3に固定しました。FIDは3000サンプルで計算しました。
TypeParametersFIDGenerated images
Standard GAN Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-4
\(D\), Adam lr =1e-5
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-4
\(G\), Adam lr =1e-4
14.66 Epoch=99, Batch=400
Coulomb GAN
No. 1
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-3
\(D\), Adam lr =1e-3
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-3
\(G\), Adam lr =1e-3
Plummer kernel \(\varepsilon\)=1.0
52.32 Epoch=99, Batch=400
Coulomb GAN
No. 2
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-3
\(D\), Adam lr =1e-4
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-3
\(G\), Adam lr =1e-3
Plummer kernel \(\varepsilon\)=1.0
90.82 Epoch=99, Batch=400
Coulomb GAN
No. 3
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-4
\(D\), Adam lr =1e-3
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-3
\(G\), Adam lr =1e-4
Plummer kernel \(\varepsilon\)=1.0
106.5 Epoch=99, Batch=400
Coulomb GAN
No. 4
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-3
\(D\), Adam lr =1e-3
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-3
\(G\), Adam lr =1e-3
Plummer kernel \(\varepsilon\)=0.1
108.0 Epoch=99, Batch=400
Coulomb GAN
No. 5
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-2
\(D\), Adam lr =1e-3
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-2
\(G\), Adam lr =1e-3
Plummer kernel \(\varepsilon\)=10.0
118.4 Epoch=70, Batch=0
Coulomb GAN
No. 6
Batch size=128
\(D\), Adam \(\beta_1\) =0.5
\(D\), Adam decay =1e-4
\(D\), Adam lr =1e-3
\(G\), Adam \(\beta_1\) =0.5
\(G\), Adam decay =1e-4
\(G\), Adam lr =1e-3
Plummer kernel \(\varepsilon\)=1.0
149.9 Epoch=99, Batch=400
FIDで比較してみると、普通のGANはCoulomb GANよりFIDが小さく、元の分布を良く再現できているようです。実際、画像をよく見てみると、Coulomb GANのうち最もFIDが小さい実験1の場合でも、崩れ気味かつ1ピクセルだけ輝度が高い部分が残っています。Coulomb GANの実験4はFIDが108.0と大きめですが、数字の背景の点々(ノイズ)も少なく、線もはっきりとしています。ただ、なんというか、数字を一定のルールで崩したような感じに見えます。他の実験結果は、ぜんぜんダメです。

\(G\)の学習速度については、普通のGANでは、15エポックもすると数字っぽく主観で綺麗な画像が出てくるのに対し、Coulomb GANは実験1のケースであっても95エポックくらいまで到達しないと崩れたりノイズが多い画像が出力されます。パラメータはここに挙げたもの以外も試してみましたが、普通のGANと同じ速さで\(G\)を学習できるパラメータを見つけることはできませんでした。

まとめ


点の分布の学習では安定して良さげな結果を出していたCoulomb GANでしたが、今回のMNISTの画像生成ではいまいちな結果となりました。真の画像よりバリエーションが豊富な画像を生成したいときはFIDは適度に大きくなる必要がありますが、どのくらいのFIDになると良いのかは今回の実験では分かりませんでした。

コード


今回の実験で使ったコードです。最後のパラメータ部分は実験ごとに書き換えています。
  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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# -*- 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 util.history import ExperimentHistory
from gan.coulomb import CoulombPotentials
from gan.gaussian_mixture.datagen import RandomSampler
from keras.models import Sequential
from keras.optimizers import Adam
import keras.backend as K
from keras.models import Sequential
from keras.layers import Dense, Activation, Reshape
from keras.layers.normalization import BatchNormalization
from keras.layers.convolutional import UpSampling2D, Conv2D
from keras.layers.advanced_activations import LeakyReLU
from keras.layers import Flatten, Dropout
from keras.datasets import mnist
from PIL import Image

GENERATED_IMAGE_PATH = 'coulomb_images/'
BATCH_SIZE = 32*4
NUM_EPOCH = 100

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

def combine_images(generated_images):
    total = generated_images.shape[0]
    cols = int(math.sqrt(total))
    rows = math.ceil(float(total)/cols)
    width, height = generated_images.shape[2:]
    combined_image = np.zeros((height*rows, width*cols), dtype=generated_images.dtype)

    for index, image in enumerate(generated_images):
        i = int(index/cols)
        j = index % cols
        combined_image[width*i:width*(i+1), height*j:height*(j+1)] = image[0, :, :]
    return combined_image

def gan_model_mnist(initializer='glorot_uniform'):
    """
    return: (generator, discriminator)
    """
    g = Sequential()
    g.add(Dense(1024, input_dim=100))
    g.add(BatchNormalization())
    g.add(Activation('relu'))
    g.add(Dense(128*7*7))
    g.add(BatchNormalization())
    g.add(Activation('relu'))
    g.add(Reshape((128, 7, 7), input_shape=(128*7*7,)))
    g.add(UpSampling2D((2, 2), data_format='channels_first'))
    g.add(Conv2D(64, (5, 5), padding='same', data_format='channels_first'))
    g.add(BatchNormalization())
    g.add(Activation('relu'))
    g.add(UpSampling2D((2, 2), data_format='channels_first'))
    g.add(Conv2D(1, (5, 5), padding='same', data_format='channels_first'))
    g.add(Activation('tanh'))
    print(g.summary())

    d = Sequential()
    d.add(Conv2D(64, (5, 5), strides=(2, 2), padding='same', input_shape=(1, 28, 28), 
                     data_format='channels_first'))
    d.add(LeakyReLU(0.2))
    d.add(Conv2D(128, (5, 5), strides=(2, 2), data_format='channels_first'))
    d.add(LeakyReLU(0.2))
    d.add(Flatten())
    d.add(Dense(256))
    d.add(LeakyReLU(0.2))
    d.add(Dropout(0.5))
    d.add(Dense(1))
    #d.add(Activation('sigmoid'))
    print(d.summary())
    return g, d

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 DiscriminatorLabelCoulomb:
    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):
    (X_train, y_train), (_, _) = mnist.load_data()
    X_train = (X_train.astype(np.float32) - 127.5)/127.5
    X_train = X_train.reshape(X_train.shape[0], 1, X_train.shape[1], X_train.shape[2])

    # Random sampler for each batch
    rs = RandomSampler(100, "normal" if eh.random_with_normal_dist else "uniform")

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

    # Make a generator G and a discriminator D
    generator, discriminator = gan_model_mnist()
    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)

            # Output generated images
            if index % 200 == 0:
                image = combine_images(points_fake)
                image = image*127.5 + 127.5
                if not os.path.exists(GENERATED_IMAGE_PATH):
                    os.mkdir(GENERATED_IMAGE_PATH)
                Image.fromarray(image.astype(np.uint8))\
                    .save(GENERATED_IMAGE_PATH+"%04d_%04d.png" % (epoch, index))

                generator.save(GENERATED_IMAGE_PATH+"generator.model")
                discriminator.save(GENERATED_IMAGE_PATH+"discriminator.model")

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

            # Update a 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()))

if __name__ == '__main__':
    GENERATED_IMAGE_PATH = 'coulomb_images/'
    eh = ExperimentHistory()
    eh.batch_size = BATCH_SIZE
    eh.random_with_normal_dist = False
    eh.X_train_is_shuffled = True
    eh.plummer_kernel_dim = 3.0
    eh.plummer_kernel_eps = 0.1
    eh.plummer_kernel_eps_half_life = 5000.0
    eh.plummer_kernel_ignore_self_potential = True
    eh.disc_Adam_decay = 1e-3
    eh.disc_Adam_lr = 1e-3
    eh.disc_Adam_beta_1 = 0.5
    eh.gen_Adam_decay = 1e-3
    eh.gen_Adam_lr = 1e-3
    eh.gen_Adam_beta_1 = 0.5
    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)

2017/12/24

Fréchet Inception Distance

はじめに

Fréchet Inception Distance (FID)と呼ばれる、Generative Adversarial Network (GAN)が生成する画像の品質を評価する指標を試してみます[1]。この指標は、画像の集合間の距離を表します。前回試したInception Scoreは画像の集合そのものの良さを表すスコアでしたので1つ画像の集合を与えるだけで計算できましたが、FIDはそのようには計算できません。GANで再現したい真の分布から生成された画像の集合と、GANで再現した分布から生成した画像の集合との距離を計算することになります。距離が近ければ近いほど良い画像であると判断します。FIDは、Google Brainが実施したGANの大規模評価の評価指標にも用いられています[2]

計算方法

FIDは、Inceptionモデルの途中の層の出力から得られるベクトル\(h\)を使ってFréchet Distance [3, 4]を計算することで求められます。Fréchet Distanceは曲線同士の距離のため、\(h\)のままでは距離を計算できません。そこで、画像から得られるベクトル\(h\)の分布が多変量正規分布(Multivariate normal distribution)に従うと仮定します。多変量正規分布は曲線なので、2つの多変量正規分布を求めると、その分布間のFréchet Distanceを計算できます。平均ベクトルと共分散行列が分かっている多変量正規分布間のFréchet Distanceは[5]で計算できます。

具体的な計算方法は[1]の著者らが[6]で公開しています。TensorFlowベースです。画像の集合を\(A\)、その要素を\(a \in A\)とし、Inceptionモデルの途中の層まで計算する関数を\(f_{\rm inception} : A \rightarrow H\)とします。\(H\)は\(h\)の集合です。まず、平均ベクトル\(\mu\)と共分散行列\(\Sigma\)を計算します。\(H\)の各要素は\(f_{\rm inception}\)で計算済みであるとしています。 \[ \mu = \frac{1}{|A|} \sum_{h \in H} h \] \[ \Sigma = \frac{1}{|A|-1} \sum_{h \in H} (h-\mu)(h-\mu)^{T} \] \(H\)から推定した分布間の距離を計算するので、[6]に倣って、\(\Sigma\)を不偏共分散行列として計算しています。少々距離が短くなりますが、標本共分散行列で計算してもGANで生成した画像の評価指標として使う分には特に問題ないでしょう。

2つの画像集合\(A_1\)と\(A_2\)の距離を計算したいので、それぞれの平均ベクトルを\(\mu_1, \mu_2\)、共分散行列を\(\Sigma_1, \Sigma_2\)とすると、Fréchet Distanceは \[ d^2 = |\mu_1-\mu_2|^2 + {\rm tr}\left (\Sigma_1 + \Sigma_2 - 2(\Sigma_1 \Sigma_2)^{\frac{1}{2}}\right )\] で計算できます。

[1]のp.7 L.13-14によると、\(h\)には最後のプーリング層を使っているとのことです。これは、Inception-v3の論文[7]のTable 1の下から3行目のpoolのことを指しています。この層の出力は1x1x2048なので、\(h\)は2048次元のベクトルということになります。[8]にも書いてありますが、画像のサンプル数が2048個より多くないと\(d^2\)が計算できないので、注意が必要です(ちょっと試すだけでも2000枚強の画像が必要とは、なんて面倒な指標なんだ!)。

実験結果

MNISTとImageNetの一部の画像を使ってFIDを計算します。MNISTは、訓練用と検証用ともに先頭3000枚を利用します。ImageNetは、10クラスからランダムに2956枚選んだ画像と、6クラスから順に2956枚に達するまで選んだ画像を利用します。

10クラスは、具体的には{n02066245, n02096294, n02100735, n02119789, n02123394, n02124075, n02125311, n02417914, n02423022, n02509815}です。6クラスは{n02066245, n02096294, n02100735, n02119789, n02123394, n02124075}です。

計算結果は下表の通りです。

\(A_1\)\(A_2\)FID (\(d^2\))
MNIST trainMNIST train6.8959054997e-11
MNIST trainMNIST val7.05136659744
MNIST valMNIST train7.0513665973
MNIST trainImageNet 10 classes338.586062646
MNIST trainImageNet 6 classes346.218188602
ImageNet 10 classesImageNet 6 classes67.1025959331
同じ画像集合間の距離は計算誤差を無視すると0になっています(1行目)。MNISTのtrainとvalの距離は7.05です(2, 3行目)。MNISTとImageNetの距離はMNIST間の距離より遠く、約340です(4, 5行目)。クラス数が異なるImageNetの画像集合間の距離は67.1で、MNISTのtrainとvalの間の距離より遠くなっています(6行目)。

距離の大きさはともかく、その大小関係は期待通りになっています。具体的には、

  • 同じ集合間の距離は0。
  • クラスが同じ画像集合間の距離は、クラスが多少異なる画像集合間より近い。
  • 全く異なる画像集合間の距離は、クラスが多少異なる画像集合間の距離より遠い。
ということです。

コード

今回の実験に使ったコードは以下です。Kerasを使っています。実験の都合上、計算した\(H\)は一旦ファイルに保存するようにしています。

 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
# -*- coding: utf-8 -*-
import os, glob
import glob
import numpy as np
from keras.applications.inception_v3 import InceptionV3, preprocess_input
from keras.applications.imagenet_utils import decode_predictions
from keras.preprocessing import image
from keras.datasets import mnist
from keras.models import Model
from PIL import Image as pil_image
from scipy.linalg import sqrtm

model = InceptionV3() # Load a model and its weights
model4fid = Model(inputs=model.input, outputs=model.get_layer("avg_pool").output)
def resize_mnist(x):
    x_list = []
    for i in range(x.shape[0]):
        img = image.array_to_img(x[i, :, :, :].reshape(28, 28, -1))
        #img.save("mnist-{0:03d}.png".format(i))
        img = img.resize(size=(299, 299), resample=pil_image.LANCZOS)
        x_list.append(image.img_to_array(img))
    return np.array(x_list)

def resize_do_nothing(x):
    return x

def frechet_distance(m1, c1, m2, c2):
    return np.sum((m1 - m2)**2) + np.trace(c1 + c2 - 2*(sqrtm(np.dot(c1, c2))))

def mean_cov(x):
    mean = np.mean(x, axis=0)
    sigma = np.cov(x, rowvar=False)
    return mean, sigma

def fid(h1, h2):
    m1, c1 = mean_cov(h1)
    m2, c2 = mean_cov(h2)
    return frechet_distance(m1, c1, m2, c2)

def calc_h(x, resizer, batch_size=8):
    r = None
    n_batch = (x.shape[0]+batch_size-1) // batch_size
    for j in range(n_batch):
        x_batch = resizer(x[j*batch_size:(j+1)*batch_size, :, :, :])
        r_batch = model4fid.predict(preprocess_input(x_batch))
        r = r_batch if r is None else np.concatenate([r, r_batch], axis=0)
        if j % 10 == 0:
            print("i =", j)
    return r

def mnist_h(n_train, n_val):
    x = [0, 0]; h = [0, 0]; n = [n_train, n_val]
    (x[0], _), (x[1], _) = mnist.load_data()
    for i in range(2):
        x[i] = np.expand_dims(x[i], axis=3) # shape=(60000, 28, 28) --> (60000, 28, 28, 1)
        x[i] = np.tile(x[i], (1, 1, 1, 3)) # shape=(60000, 28, 28, 1) --> (60000, 28, 28, 3)
        h[i] = calc_h(x[i][0:n[i], :, :, :], resize_mnist)
    return h[0], h[1]

def imagenet_h(files, batch_size=8):
    xs = []; hs = []
    for f in files:
        img = image.load_img(f, target_size=(299, 299))
        x = image.img_to_array(img) # x.shape=(299, 299, 3)
        xs.append(x)
        if len(xs) == batch_size:
            hs.append(calc_h(np.array(xs), resize_do_nothing))
            xs = []
    if len(xs) > 0:
        hs.append(calc_h(np.array(xs), resize_do_nothing))
    return np.concatenate(hs, axis=0)

# Calculate and save H of MNIST
h_train, h_val = mnist_h(3000, 3000)
np.save("mnist_h_train.npy", h_train)
np.save("mnist_h_val.npy", h_val)

# Calculate and save H of the part of Imagenet 
h_imagenet = imagenet_h(glob.glob("from_imagenet/*.jpg")) # 10 classes
h_imagenet_seq = imagenet_h(sorted(glob.glob("from_imagenet_seq/*.jpg"))[0:2956]) # 6 classes
np.save("imagenet_h.npy", h_imagenet)
np.save("imagenet_h_seq.npy", h_imagenet_seq)

# Load H and calculate FID
h_train = np.load("mnist_h_train.npy")
h_val = np.load("mnist_h_val.npy")
h_imagenet = np.load("imagenet_h.npy")
h_imagenet_seq = np.load("imagenet_h_seq.npy")
print("FID between MNIST train and val :", fid(h_train, h_val))
print("FID between MNIST val and train :", fid(h_val, h_train))
print("FID between MNIST train and train :", fid(h_train, h_train))
print("FID between MNIST train and imagenet :", fid(h_train, h_imagenet))
print("FID between MNIST train and imagenet_seq :", fid(h_train, h_imagenet_seq))
print("FID between imagenet and imagenet_seq :", fid(h_imagenet, h_imagenet_seq))

参考

[1] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, Sepp Hochreiter, "GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium," arXiv:1706.08500, 2017, https://arxiv.org/abs/1706.08500
[2] Mario Lucic, Karol Kurach, Marcin Michalski, Sylvain Gelly, Olivier Bousquet, "Are GANs Created Equal? A Large-Scale Study," arXiv:1711.10337, 2017, https://arxiv.org/abs/1711.10337
[3] Fréchet, M. "Sur la distance de deux lois de probabilité," C. R. Acad. Sci. Paris, 244, 689-692, 1957 (内容を確認したわけではない)
[4] http://www.thothchildren.com/chapter/59b4f81975704408bd430061 (Fréchet Distanceの解説記事)
[5] D. C. Dowson and B. V. Landau, "The Fréchet Distance between Multivariate Normal Distributions," Journal of multivariate analysis, 12, 450-455, 1982, http://www.sciencedirect.com/science/article/pii/0047259X8290077X
[6] https://github.com/bioinf-jku/TTUR/blob/master/fid.py ([1]の著者らが作成したコード)
[7] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jonathon Shlens, Zbigniew Wojna, "Rethinking the Inception Architecture for Computer Vision," arXiv:1512.00567, 2015, https://arxiv.org/abs/1512.00567

2017/12/17

XGBoost

はじめに

2年以上前に流行っていた[1]らしいXGBoost[2]をいまさら試してみました。勾配ブースティング決定木(Gradient boosting decision tree)と呼ばれている手法です。手法自体の説明は探すと色々出てきますので、そちらを参照ください。今回は、回帰問題と識別問題を簡単なデータで試してみます。

インストール

基本的には[2] を参考にインストールします。今回はWindows+VC2015のコマンドラインツールがインストールされている環境下で試しました。ソースコードの入手までは出来ている状態で、まず、VS2015 x64 Native Tools Command Prompt を立ち上げ、次のコマンドを実行すると、xgboostをビルドできます。

vcvarsall amd64
cd (path to xgboost)
mkdir build64
cd build64
cmake .. -G "Visual Studio 14 2015 Win64"
cmake --build . --config Release
Python環境下で利用するために、
cd (path to xgboost)
cd python-package
python setup.py install
を実行します。これで、import xgboostが使えるようになります。

回帰

回帰には、XGBRegressorを利用します。実験データは正弦波にノイズを加えたものです。

コード

コードは次の通りです。

%matplotlib inline
import xgboost as xgb
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
from IPython.display import HTML
x = []
y = []
for i in range(2):
    xi = np.linspace(0,1,1000) + np.random.normal(0,0.1,1000)
    yi = np.cos(xi*10) + np.random.normal(0,0.1,1000)
    xi = xi.reshape(-1, 1)
    yi = yi.reshape(-1, 1)
    x.append(xi)
    y.append(yi)

mod = xgb.XGBRegressor(learning_rate=0.1, max_depth=2, n_estimators=100)
mod.fit(x[0], y[0])
y_train_pred = mod.predict(x[0])
y_test_pred = mod.predict(x[1])
from sklearn.metrics import mean_squared_error
print('MSE train : {0:.3f}, test : {1:.3f}'.
      format(mean_squared_error(y[0], y_train_pred),
             mean_squared_error(y[1], y_test_pred)))

x_pred = np.linspace(-0.5,1.5,10000).reshape(-1, 1)
y_pred = mod.predict(x_pred)
#plt.scatter(x[0], y[0], s=1)
plt.scatter(x[1], y[1], s=1)
plt.plot(x_pred, y_pred, 'C3')
plt.show()
plt.close()

結果

Jupyter notebookで実行すると、以下のような結果が得られます。

MSE train : 0.009, test : 0.013
テストデータに対する回帰折れ線(?)は次のようになりました。ギザギザではありますが、中心付近を通過するように学習できていることが分かります。

識別

識別にはXGBClassifierを使います。データは、1個のガウス分布から1000個サンプリングした点をクラス-1、2個のガウス分布から500個ずつサンプリングした点をクラス1として作成しています。

コード

コードは次の通りです。

%matplotlib inline
import xgboost as xgb
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from IPython import display

# Gaussian distribution for class -1
mu = [0, 0]
sigma = [[30, 0], [0, 50]]
x1 = np.random.multivariate_normal(mu, sigma, 1000) # Train
x1t = np.random.multivariate_normal(mu, sigma, 1000) # Test

# 2 Gaussian distribution for class 1
mu = [5, 5]
sigma = [[5, 0], [0, 3]]
x2 = np.random.multivariate_normal(mu, sigma, 500) # Train
x2t = np.random.multivariate_normal(mu, sigma, 500) # Test

mu = [-10, 15]
sigma = [[3, 2], [2, 3]]
x2 = np.concatenate([x2, np.random.multivariate_normal(mu, sigma, 500)]) # Train
x2t = np.concatenate([x2t, np.random.multivariate_normal(mu, sigma, 500)]) # Test

# Make labels of train and test
y = np.concatenate([np.ones(1000)*-1, np.ones(1000)*1])

# classify by xgboost
mod = xgb.XGBClassifier(
        learning_rate=0.1,
        max_depth=2, # Maximum depth of each tree
        n_estimators=100) # The number of trees
x = np.concatenate([x1, x2])
xt = np.concatenate([x1t, x2t])
mod.fit(x, y)
y_pred = mod.predict(x)
yt_pred = mod.predict(xt)

from sklearn.metrics import accuracy_score
print('MSE train : {0:.3f}, test : {1:.3f}'.
      format(accuracy_score(y, y_pred), accuracy_score(y, yt_pred)))

# Show boundaries between class -1 and 1
gx1 = np.linspace(-20, 20, 1000)
gx2 = np.linspace(-20, 20, 1000)
gx1, gx2 = np.meshgrid(gx1, gx2)
y = mod.predict(np.array([gx1, gx2]).T.reshape(-1, 2))
plt.scatter(x1t[:,0], x1t[:,1], s=0.2)
plt.scatter(x2t[:,0], x2t[:,1], s=0.2)
cm = LinearSegmentedColormap.from_list('cm', [(0, 'blue'),(1, 'blue')])
plt.contour(gx1, gx2, y.reshape(1000, -1).T, [0], alpha=1, cmap=cm)
plt.show()

結果

Jupyter notebookで実行すると、以下のような結果が得られます。

MSE train : 0.926, test : 0.913
クラス間の識別境界(青線)と、テストデータのクラス-1(青点)と1(橙点)の分布は次の図のようになります。決定木ベースなのでギザギザですが、見た目では分割できているように見えます。

参考

[1] http://smrmkt.hatenablog.jp/entry/2015/04/28/210039
[2] http://xgboost.readthedocs.io/en/latest/build.html