2018年11月11日日曜日

場所によって異なる興奮性ニューロン

ニューロンは色々種類があるけれど、満遍なく存在するわけではないらしい。 というわけで、以下ざっくりとメモ。

ざっくり


Tasicらによると、マウスの脳の異なる領域にある興奮性ニューロンはタイプが異なるらしい[1, 2]。

Tasicらは、マウスの脳のAnterior lateral motor cortex (ALM, 運動野のどこか)とPrimary visual (VISp, 視覚野)を対象として、遺伝子発現プロファイルに基づいて細胞タイプを調べた。 そうしたところ、133個の細胞タイプが存在し、101個のタイプがALMに、111個のタイプがVISpにあることが分かった。そのうち、79個は両方の領域で存在することもわかった。

また、おおきく3つに分類すると、抑制(GABAergic)タイプが61個、興奮(glutamatergic)タイプが56個、ニューロンでないタイプが16個、となった。

これらのうち、興奮タイプのニューロンの大部分が上記2つの領域間で共有されていなかった。 抑制タイプの中には領域間で遺伝子発現が異なるものがあったが、統計的に異なるというには不十分であった([2]のFig.2 cを参照)。

ちょこっと計算


計算してみると、領域間で共有されていないタイプはALMでは22個、VISpでは32個あり、合計で54個となります。 共有の有無の基準は良く読まないとわかりませんが、この54個がすべて興奮タイプのニューロンだということのようです。[2]のFig.2 cを見ると、確かに興奮性のものばかり遺伝子発現に差異があるようです(図のDE genesはdifferentially expresed genesの略かな?)。

参考


[1] https://www.nature.com/articles/d41586-018-07027-3
[2] https://doi.org/10.1038/s41586-018-0654-5

Firefoxのタブの設定

Firefoxで新しいタブを開く条件を設定する方法です。
動作確認は63.0.1で実施しています。

まず、about:configを開きます。

ポップアップウィンドウをタブで開くには、

browser.link.open_newwindow.restriction 0
とします。

ロケーションバーから開くときに新しいタブで開くには、

browser.urlbar.openintab true
とします。

ブックマークから開くときに新しいタブで開くには、

browser.tabs.loadBookmarksInTabs true
とします。

検索ボックスから開くときに新しいタブで開くには、

browser.search.openintab true
とします。

これで、何か新しいURLを開こうとしたときに、常に新しいタブで開くようになります。

2018年10月6日土曜日

LVMでパーティション作成

はじめに


LVMでディスクを設定します。細かい説明は適当にぐぐってください。

環境


VirtualBox環境で、2GBのHDDを2個目のHDDとして追加しています(/dev/sdb)。 この追加したHDDをLVMにします。OSはUbuntu 16.04.4 LTSです。

まとめ


$ sudo parted /dev/sdb
$ (parted) mklabel gpt
$ (parted) mkpart base-partition 0% 100%
$ (parted) set 1 lvm on
$ (parted) quit
$ sudo pvcreate /dev/sdb1
$ sudo vgcreate vg0 /dev/sdb1
$ sudo lvcreate -l 100%FREE -n lv0 vg0
$ sudo mkfs.ext4 /dev/vg0/lv0
$ sudo mkdir /disk1
$ sudo mount /dev/vg0/lv0 /disk1

詳細


  1. まずは元にするパーティションが必要なので、partedで作成します。
    $ sudo parted --list 
    
    で一覧を取得して、対象となるデバイスを選びます。 不明な点があればman partedを読みましょう。 parted --helpは内容を知っている人向けです。
  2. 選択したデバイスを指定してpartedを起動します。
    $ sudo parted /dev/sdb
    
  3. パーティションの情報を表示します。
    (parted) print
    ...
    パーティションテーブル: unknown
    ...
  4. パーティションの管理方式をGPTに変更します。
    (parted) mklabel gpt
    (parted) print
    ...
    パーティションテーブル: gpt
    ...
  5. パーティションを作成します。
    (parted) mkpart base-partition 0% 100%
    (parted) print
    モデル: ATA VBOX HARDDISK (scsi)
    ディスク /dev/sdb: 2147MB
    セクタサイズ (論理/物理): 512B/512B
    パーティションテーブル: gpt
    ディスクフラグ: 
    
    番号  開始    終了    サイズ  ファイルシステム  名前            フラグ
     1    1049kB  2146MB  2145MB                    base-partition
    
    man partedで表示されるマニュアルでは、
    mkpart part-type [fs-type] start end
                         Make  a  part-type  partition  for filesystem fs-type (if
                         specified), beginning at start  and  ending  at  end  (by
                         default  in megabytes).  part-type should be one of "pri‐
                         mary", "logical", or "extended".
    
    と書かれていますが、第1引数に指定した文字列は、パーティションの名前になります。
  6. LVM用のフラグを立てます。
    (parted) set 1 lvm on
    (parted) p                                                                
    モデル: ATA VBOX HARDDISK (scsi)
    ディスク /dev/sdb: 2147MB
    セクタサイズ (論理/物理): 512B/512B
    パーティションテーブル: gpt
    ディスクフラグ: 
    
    番号  開始    終了    サイズ  ファイルシステム  名前            フラグ
     1    1049kB  2146MB  2145MB                    base-partition  lvm
    

    parted上での作業は完了ですので、quitで抜けます。

  7. Physical Volumeを作ります。
    $ sudo pvcreate /dev/sdb1
    sdbではなく、sdb1を指定します。 結果は
    $ sudo pvdisplay
    で確認できます。
  8. Volume Groupを作ります。
    $ sudo vgcreate vg0 /dev/sdb1
    結果は
    $ sudo vgdisplay
    で確認できます。
  9. Logical Volumeを作ります。サイズ指定なら
    sudo lvcreate -L 1.996GiB -n lv0 vg0
    全て使うなら
    sudo lvcreate -l 100%FREE -n lv0 vg0
    とします。 出来上がったサイズは
    $ sudo lvdisplay
      --- Logical volume ---
      LV Path                /dev/vg0/lv0
      LV Name                lv0
      VG Name                vg0
      ...
      LV Size                2.00 GiB
      ...
    で確認できます。
  10. ext4でフォーマットします。
    $ sudo mkfs.ext4 /dev/vg0/lv0
  11. 出来上がったディスクをマウントします。
    $ sudo mkdir /disk1
    $ sudo mount /dev/vg0/lv0 /disk1
    $ df
    Filesystem                  1K-blocks     Used Available Use% Mounted on
    /dev/mapper/vg0-lv0           2027408     3072   1903300   1% /disk1
    
    マウントするときは、UUIDでも可能です。
    $ blkid
    /dev/mapper/vg0-lv0: UUID="0805787f-b7a0-488a-9321-ed7ed77322ab" TYPE="ext4"
    
    $ sudo mount UUID=3be08a7c-449a-4d74-b576-1064ef72a7ef /disk1
    

    fstabに追加したいときは、

    $ man fstab
    を参照するかぐぐりましょう。

2018年9月15日土曜日

tmuxの画面が固まったら

Windows Subsystem for Linux (WSL)のターミナルでtmuxを使っているときに、
ないしは、byobuのバックエンドでtmuxを使っているときに、
うっかり、SJISの文字コードを表示すると画面が固まります。

復旧するすべは分からないのですが、固まったウィンドウのみを削除することはできます。

表示が乱れるかもしれませんが、固まっているウィンドウを表示し、

 [prefix-key]→ & → y → [prefix-key] → d
を押します。ウィンドウをkillして、tmuxから抜けています。そして、再度、tmux (またはbyobu) に接続すると、固まっていたウィンドウがなくなっています。

または、固まっていないウィンドウに移動し、画面を正常にするために一旦 tmux から抜けて、再度接続し、そこで
$ tmux kill-window -t :[ウィンドウの番号]
を実行すると、固まったウィンドウを削除できます。tmuxバックエンドのbyobuでも同じです。

2018年9月4日火曜日

非可逆圧縮音声ファイルと非圧縮音声ファイルの自動判定

はじめに


人間の耳では感知しにくい、mp3による圧縮音声と非圧縮音声の違いをニューラルネットワークを使って判定できるか試してみます。

mp3とwavのスペクトログラムを下図のように比較すると、明らかに周波数の高い部分がmp3では消えているので、それを特徴として捉えられれば判定できそうです。

方法


Hennequinらの方法(Codec Independent Lossy Audio Compression Detection, 2017)を参考にしました。

今回の実験に使用したニューラルネットワークを以下に示します。 入力は、フレーム長(窓幅)を512サンプルとした対数パワースペクトルです。Nフレームつなげたものを画像とみなしてCNNで処理します。フレームシフトは256サンプルです。活性化関数は全てLeakyReLUです。

Conv2D     Nx256 32-channel
MaxPooling (N/2)x128 32-channel
Conv2D     (N/2)x128 16-channel
MaxPooling (N/4)x64 16-channel
Conv2D     (N/4)x64 16-channel
MaxPooling (N/8)x32 16-channel
Conv2D     (N/8)x32 16-channel
MaxPooling (N/16)x16 16-channel
Conv2D     (N/16)x16 16-channel
Flatten    (N/16)x16x16
Dense      64
Dropout
Dense      64
Dropout
Dense      2
Softmax(2-classes)

データ


訓練データは12曲で、先頭から曲の終端まで1000フレーム(5.8秒)ごとに1回、特徴ベクトルを抽出しました。 合計506個です。検証用データは15曲で特徴ベクトルは792個です。

結果


1エポック当り100ステップ、バッチサイズ4で50エポック訓練した結果、以下のようになりました。 N=1のみ3回実験しました。
NTrain-accValidation-acc
102410.968250
5120.9992500.955000
2560.9982500.902750
1280.9952500.933000
640.9935000.864750
320.9787500.925500
160.9815000.976750
80.9855000.972250
40.9800000.899500
20.9815000.889000
10.9750000.931500
10.9635000.956000
10.9547500.936250
Nが小さくても判定できています。僅か1フレームでも9割前後の精度がでています。 なお、MaxPooling層がない場合は訓練が進まず、精度がまったくでませんでした(0.5前後で運任せ)。

コード


https://bitbucket.org/bluewidz/mp3detect/src/default/ を参照ください。

2018年8月18日土曜日

Linuxbrew installer

gccくらいはインストールされている環境にLinuxbrewをインストールするスクリプトを作りました。
https://bitbucket.org/bluewidz/linuxbrew-bootstrap/
動作確認は、LXDのDebian 9のコンテナ内で実施しています。

Linuxbrewって何?という方は、
https://qiita.com/thermes/items/926b478ff6e3758ecfea
をどうぞ。

簡単にいうと、Linuxbrew
http://linuxbrew.sh/
を使うと、ルート権限のないLinuxマシンで色々なソフトウェアを簡単にインストールできます。

pythonなら、例えば、

$ brew install python
のようにします。

2018年8月5日日曜日

Raspberry Pi で録音と再生

Raspberry Pi で録音と再生をしてみます。

準備


Linuxでオーディオデバイスを扱う場合、ALSAとPulseAudioが使えます。
[Audio device]-[ALSA]-[PulseAudio]-[Application]
の関係があるので、色々と例外はあるようですが、ここではPulseAudioを使うことにします。そのあたりの詳細は
http://mickey-happygolucky.hatenablog.com/entry/2015/04/04/105512
https://cpplover.blogspot.com/2012/03/pulseaudio.html
を参照ください。

Raspberry Pi 3にはマイク端子がないので、USBマイクか、USBオーディオインターフェイスを使ってマイクを接続します。

録音と再生


マイクを接続した後、デバイスの番号を探します。
$ pactl list sources
今回試した環境では、
Source #0
        State: SUSPENDED
        Name: alsa_output.usb-0d8c_C-Media_USB_Headphone_Set-00-Set.analog-stereo.monitor
        Description: Monitor of Audio Adapter アナログステレオ
        <以下略>

Source #1
        State: SUSPENDED
        Name: alsa_input.usb-0d8c_C-Media_USB_Headphone_Set-00-Set.analog-mono
        Description: Audio Adapter アナログモノ
        <以下略>

Source #2
        State: SUSPENDED
        Name: alsa_output.1.analog-stereo.monitor
        Description: Monitor of bcm2835 ALSA アナログステレオ
        <以下略>
と出力されました。USB接続のマイクは、Source #1が該当するので、マイクから録音するには
$ parecord -d 1 test.wav
とします。

録音できていることを確認するには、
$ paplay -d 0 test.wav
とします。paplayで指定するデバイス番号は、
$ pactl list sinks
で調べます。読み方はsourcesの場合と同様です。

デバイス番号を毎回指定したくないときは、スピーカの場合は、
$ pactl set-default-sink [デバイス番号]
で、マイクの場合は、
$ pactl set-default-source [デバイス番号]
でデフォルトで使うデバイスを指定できます。

別のPCのオーディオデバイスの利用


PulseAudioはネットワーク越しのオーディオデバイスも扱えます。
そこで今回は、Windows PCにつながったスピーカやマイクを使う方法を試してみます。

まず、pulseaudio-1.1のWindows用バイナリをどこかから入手して、各種設定をします。バイナリの入手と設定は
https://vogel.at.webry.info/201711/article_9.html
が参考になります。

今回試すWindows側の default.pa は
load-module module-waveout
load-module module-native-protocol-tcp auth-ip-acl=192.168.0.0/16
としました。

Windowsのコマンドプロンプトから
pulseaudio -F ..\default.pa --exit-idle-time=60000
で起動しておいて、Raspberry Piから
$ paplay -s 192.168.0.1 test.wav
とすると、Windows側から再生されます。WindowsのIPアドレスが192.168.0.1だとしています。

同様に録音も出来ます。
$ parecord -s 192.168.0.1 test.wav

GStreamer


PulseAudioだけだと録音できるだけですが、GStreamer (https://gstreamer.freedesktop.org/) を経由すると色々遊べます。

まずは、マイク入力を直接スピーカに流してみます。
$ gst-launch-1.0 pulsesrc device=1 ! pulsesink device=0
デバイス番号は、さきほどpactlで調べた番号と同じです。Ctrl+Cで停止できます。 ちょっと遅延しますが、マイクで録音した音がスピーカから流れます。

deviceの番号については、デバイスが存在していれば間違えてもエラーは起きません。代わりに、全くの無音(振幅0)が録音されます。何も音が出ない場合はデバイス番号を間違っていないか、まずはチェックしましょう。

さて、ただ流すだけでは面白くないので、エコーをかけてみます。

$ gst-launch-1.0 pulsesrc device=1 \
    ! audioecho delay=300000000 intensity=0.5 feedback=0.5 \
    ! pulsesink device=0
delayの単位はナノ秒なので、この設定では300ミリ秒遅延のエコーになります。詳細は、
$ gst-inspect-1.0 audioecho
を実行して、Element Properties 部分の説明を読んでみてください。

今度は、ローパスフィルタを試してみます。

gst-launch-1.0 pulsesrc device=1 ! audiowsinclimit mode=low-pass cutoff=100 ! pulsesink device=0
cutoffの単位はHzなので、100Hz以下を通過させるフィルタになっています。 声がくぐもった感じになります。

せっかくなので、ハイパスフィルタも試してみます。

$ gst-launch-1.0 pulsesrc device=1 ! audiowsinclimit mode=high-pass cutoff=2000 ! pulsesink device=0
低い音が弱くなります。

wav形式で録音するには、

$ gst-launch-1.0 pulsesrc device=1 ! wavenc ! filesink location="test.wav"
とします。これを使うと、加工した音声も簡単に記録できます。

同様に、録音したwavを再生するには、

$ gst-launch-1.0 filesrc location=test.wav ! decodebin ! pulsesink device=0
とします。

もちろん、GStreamerを使うときもPulseAudioのサーバーに接続できます。

$ gst-launch-1.0 pulsesrc server=192.168.0.1 ! pulsesink server=192.168.0.1

GStreamerの不具合


GStreamerのpulsesrc Ver.1.8.3には不具合があり、どう設定しても無音しか記録できません。 古いバージョン(例えば、1.4.4)または新しいバージョン(1.14.2)では動作します。

不具合修正の場所とコミットメッセージは
https://cgit.freedesktop.org/gstreamer/gst-plugins-good/commit/ext/pulse/pulsesrc.c?id=4833f02e47d5b7fc4284517e92d5976b20188a4a
にあります。2016年9月25日にコミットされています。

ブランチ1.8のログ
https://cgit.freedesktop.org/gstreamer/gst-plugins-good/log/?h=1.8
を見ると、1.8.3のリリースの直後に修正されていることがわかります。 tagが付いている範囲だと、Release 1.9.2とRelease 1.9.90の間で修正されています。

Ubuntu 16.04 LTSを使っている人は気をつけましょう。 デフォルトでインストールされるGStreamerがこの不具合を踏んでいます。

2018年7月21日土曜日

OneNote 2016

OneNote 2016がなぜかCPUを1コア消費し続けるようになった。
再起動しても、OneNoteのキャッシュを削除しても、起動してしばらくすると再発する。

何かのアップデートでバグが入ったようだ。

SysinternalsのProcess Explorerを使ってどのスレッドが消費しているのかを調べたところ、

ntdll.dll!RtlReleaseSRWLockExclusive

がCPUを消費し続けていることが分かった。名前とCPUを消費し続けていることから考えると、排他制御をスピンロックを使って実装したものの、ロックを抜ける条件を書くところで何かをミスったのだろう。

何を処理しているスレッドなのかは分からないが、Process Explorerで当該スレッドをSuspendすることで、とりあえずCPUを消費し続ける状況は回避できた。今のところ、OneNoteの操作には影響はないようにみえる。この方法で不具合修正が入るまで回避かな。

2018年5月27日日曜日

Linuxターミナルの色設定

はじめに


Linuxのターミナルの色設定が変わってしまった場合の対処方法についてメモしておきます。 Debian 8.10環境を前提にしています。

環境変数TERM


環境変数TERMに設定することで色設定を変更できます。例えば、
$ TERM=xterm-256color
とします。

なお、xterm-256colorの代わりにscreen-256colorを使うと、lessの検索時にヒットした文字列が斜体(イタリック)になり、ひじょーに見にくくなります。

TERMに設定できる値一覧


TERMに設定できる値は、
$ infocmp -D
で表示されるディレクトリ以下にあるファイル名です。

$ infocmp xterm-256color
のように、TERMに指定できる値を引数に指定することで、内部の設定を表示できます。読み方がわかりませんが…。2個指定すると2つの設定の差分が表示されます。

byobu経由


byobu経由でTERMがおかしくなるときは、手動でTERMを設定するか、byobuの設定を変更します。 tmuxバックエンドの場合は、~/.config/byobu/.tmux.conf に
set-option -g default-terminal xterm-256color
set -g terminal-overrides 'xterm:colors=256'
を追記します。もしかすると2行目は不要かもしれません。 ~/.byobu 以下にも設定ファイルがあるかもしれませんが、有効なほうを変更します。

参考


https://unix.stackexchange.com/questions/179173/make-less-highlight-search-patterns-instead-of-italicizing-them
https://www.tldp.org/HOWTO/Text-Terminal-HOWTO-16.html

追記(2018/9/8)


byobuの起動前後でターミナルの色が異なる場合は、
byobu-disable-prompt
を実行しましょう。byobuの色設定が使われなくなります。逆に、byobuの設定を使いたい場合は、
byobu-enable-prompt
を実行します。

2018年5月4日金曜日

PythonとC++間でのnumpy配列の受け渡し

はじめに


Pythonで作ったnumpy.ndarrayをC++側で加工して、Pythonに返す方法を試します。

PythonからCの関数を呼び出す方法は色々あるようですが、 今回は、numpyを操作したいので、Boost C++ Librariesを使ってみます。

準備


Boost C++ Librariesをインストールします。
https://www.boost.org/doc/libs/1_67_0/more/getting_started/unix-variants.html
の第1節に従ってダウンロードし、第5節に従ってビルドします。

方法(動的リンク版)


まず、C++のコードです。ファイル名は test1.cpp とします。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#include <boost/python/numpy.hpp>
namespace p = boost::python;
namespace np = boost::python::numpy;

void plus_scalar(np::ndarray &a, float b)
{
  auto N = a.shape(0);
  auto *p = reinterpret_cast<float*>(a.get_data());
  for(auto v=p; v!=p+N; ++v){*v=*v+b;}
}

BOOST_PYTHON_MODULE(test1)
{
  Py_Initialize();
  np::initialize();
  p::def("plus_scalar", plus_scalar);
}
Python側から第1引数で渡される1次元NumPy配列に第2引数で渡される値を足すだけの関数plus_scalarを定義しています。

このtest1.cppをコンパイルします($はシェルのプロンプトです)。

$ g++ --std=c++14 --shared -fPIC -I$BOOST_PATH/include/ -I/usr/include/python2.7/ test1.cpp -L$BOOST_PATH/lib -lboost_numpy27 -lboost_python27 -o test1.so
$BOOST_PATHはBoost C++ Librariesをインストールしたパスです。

ライブラリのパスが通っていないので、

$ export LD_LIBRARY_PATH=$BOOST_PATH/lib
を実行し、.soをロードできるようにしておきます。

これをPythonから呼び出します。次のようにPythonスクリプト t1.py を作成します。

1
2
3
4
5
6
7
8
9
import sys
sys.path.append('.')
import test1
import numpy as np

a = np.array([1,2,3], dtype=np.float32)
print(a)
test1.plus_scalar(a, 3.0)
print(a)

そして、実行します。

python t1.py
すると、
[ 1.  2.  3.]
[ 4.  5.  6.]
が出力されます。確かに、3.0が足されています。

方法(静的リンク版)


動的リンク版では LD_LIBRARY_PATH を設定しましたが、使いにくいのでBoost関連のライブラリだけでも静的リンクしておきましょう。次のようにすることで.soを作ることができます。
g++ --std=c++14 --shared -fPIC -I$BOOST_PATH/include/ -I/usr/include/python2.7/ test2.cpp $BOOST_PATH/lib/libboost_numpy27.a $BOOST_PATH/lib/libboost_python27.a -o test2.so

もし、次のエラー

/usr/bin/ld: $BOOST_PATH/lib/libboost_numpy27.a(ndarray.o): 再配置 R_X86_64_32 (`.rodata.str1.8' に対する) は共有オブジェクト作成時には使用できません。-fPIC を付けて再コンパイルしてください。
$BOOST_PATH/lib/libboost_numpy27.a: error adding symbols: 不正な値です
collect2: error: ld returned 1 exit status
が出て、リンクできないときは、エラーメッセージの通りにBoostをビルドし直します。 具体的には、[4]を参考に、
$ ./b2 --clear
$ ./b2 cxxflags=-fPIC cflags=-fPIC install
のように、cxxflagsを追加してビルドし直します。cflagsは不要かもしれません。なお、ここで追加したオプション名は、--cxxflagsでも、-cxxflagsでもなく、ハイフン無しの単なるcxxflagsです。

test2.cppは

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
#include <boost/python/numpy.hpp>

namespace p = boost::python;
namespace np = boost::python::numpy;

void plus_scalar(np::ndarray &a, float b)
{
  auto N = a.shape(0);
  auto *p = reinterpret_cast<float*>(a.get_data());
  for(auto v=p; v!=p+N; ++v){*v=*v+b;}
}

BOOST_PYTHON_MODULE(test2)
{
  Py_Initialize();
  np::initialize();
  p::def("plus_scalar", plus_scalar);
}
t2.pyは
1
2
3
4
5
6
7
8
9
import sys
sys.path.append('.')
import test2
import numpy as np

a = np.array([1,2,3], dtype=np.float32)
print(a)
test2.plus_scalar(a, 3.0)
print(a)
です。test1がtest2に変わっただけです。
$ python t2.py
を実行すると、動的リンク版のときと同様に
[ 1.  2.  3.]
[ 4.  5.  6.]
と出力されます。

参考


[1] http://d.hatena.ne.jp/nihohi/20120306/1331002942
[2] http://tadaoyamaoka.hatenablog.com/entry/2017/05/25/234934
[3] https://stackoverflow.com/questions/10968309/how-to-import-python-module-from-so-file
[4] https://stackoverflow.com/questions/27848105/how-to-force-compilation-of-boost-to-use-fpic

2018年4月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年4月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年4月8日日曜日

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年4月7日土曜日

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年4月1日日曜日

量子機械学習

半年前の記事ですが、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年3月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年2月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年1月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)