にせねこメモ

はてなダイアリーがUTF-8じゃないので移ってきました。

このブログについて

文字・フォント・プログラム・技術・趣味などについて、Twitterでは書きづらい長い内容などをまとめるためのブログです。基本的には自分用のメモとして書いている部分が多いです。

リンク等

Pixiv http://pixiv.me/nixeneko
Tumblr http://nixeneko.tumblr.com/ 絵。Pixivアカウント持ってなくても見れます
Twitter https://twitter.com/nixeneko  
MediaMarker http://mediamarker.net/u/nixeneko/ 主に文字・言語関係の蔵書
GitHub https://github.com/nixeneko プログラム用、あまり使ってない
Amazon
欲しい物リスト
http://www.amazon.co.jp/registry/wishlist/1C43ZFBA4IL6Z プレゼントを下さい
ナナシスID ZhRYMnA
デレステID 421820148

同人誌(無料公開)

http://nixeneko.hatenablog.com/entry/c88_russian_alphabethttp://nixeneko.hatenablog.com/entry/c90_greek_latin_cyrillic

WindowsでChainerCVのサンプルを動かしFaster R-CNNをトレーニングしてみる

先日のChainer Meetupにて、ChainerCVというライブラリを知った。

Chainerの上で動作する、コンピュータビジョンタスクのためのディープラーニングライブラリとのことである。

このChainerCVで、Faster R-CNNによる物体検出が簡単に利用できるらしく、exampleを動かすことで試してみた。

環境

記事時点でのChainer, ChainerCVのバージョンは以下の通り。

  • Chainer: 2.0.0
  • ChainerCV: 0.5.1

ディープラーニングWindows使うとつまらない所でひっかかって割とつらいのでLinux使うべきだと思う。

インストー

  1. Chainer (Version 2.0.0)とCuPyをインストールする。
  2. 必要なライブラリが入ってなければpipで入れる
    • 必須: Cython, Pillow
    • 任意: Matplotlib, OpenCV
  3. ChainerCV (Version 0.5.1)をインストールする
    • Anacondaプロンプトを開き、次を実行する:
pip install chainercv

これでインストールは完了。


あと、(v0.5.1の) examplesを動かしてる時に、データをダウンロードするところでZeroDivisionErrorで止まったので、
chainercv/utils/download.py (C:\Anaconda3\Lib\site-packages\chainercv\utils\download.py などにある)の25行目

    speed = int(progress_size / (1024 * duration))

がdurationが0での時でも動くように書き換える:

    try:
        speed = int(progress_size / (1024 * duration))
    except ZeroDivisionError:
        speed = float('inf')

これは最新のリポジトリでは直ってるので、新しいバージョンがリリースされれば解決される問題であると思う。

examplesを動かしてみる

ソースコードをとってくる

git clone https://github.com/chainer/chainercv.git
cd chainercv
git checkout v0.5.1
cd examples/faster_rcnn

demo.py

demo.pyは物体検出を行い、結果を表示する。デフォルトでは、トレーニング済みのモデルをネットからダウンロードしてきてそれを利用する。

適当な画像hoge.jpgを用意して

python demo.py hoge.jpg

とすると、検出結果が表示される🐱
f:id:nixeneko:20170623002835p:plain

train.py

これはモデルのトレーニング用プログラムで、次のようにして実行する。0番目のGPUを使用する設定にしている*1

python train.py --gpu 0

しかし、Windowsでは実行すると次のようなエラーが発生する。

AttributeError: Can't pickle local object 'main.<locals>.transform'

このサイトによると、multiprocessingがサブプロセスを生成する方法がUnixWindowsで異なっていることがエラーの原因らしい。

Unixではデフォルトでforkによりサブプロセスを作るが、Windowsではforkは使えず、spawnによってサブプロセスを生成する。
このspawnでサブプロセスを生成する際、targetはpickableである(pickleでシリアライズできる)必要があるとのこと。特に、関数内に定義された関数はpickableでないので、モジュールのトップレベルに持ってくる必要がある。

エラーメッセージを見ると、train.pyのmain関数内に定義されたtransform関数が引っかかっているようである。なので、transform関数をトップレベルまで持ってくれば良さそうだ。

train.pyのfaster_rcnnの定義とtransform関数の定義をmain関数の定義の前に持ってくると、とりあえず動くようになった。faster_rcnnはtransformの中で参照されているため、これもトップレベルにもってくる。

faster_rcnn = FasterRCNNVGG16(n_fg_class=len(voc_detection_label_names),
                              pretrained_model='imagenet')
def transform(in_data):
    img, bbox, label = in_data
    _, H, W = img.shape
    img = faster_rcnn.prepare(img)
    _, o_H, o_W = img.shape
    scale = o_H / H
    bbox = transforms.resize_bbox(bbox, (H, W), (o_H, o_W))

    # horizontally flip
    img, params = transforms.random_flip(
        img, x_random=True, return_param=True)
    bbox = transforms.flip_bbox(
        bbox, (o_H, o_W), x_flip=params['x_flip'])

    return img, bbox, label, scale
    
def main():
    #(略)

この様に変更した後に、

python train.py --gpu 0

を実行し、トレーニングループを回し、終了するまで待つ。

トレーニングが終わると、デフォルトではresult/snapshot_model.npzにトレーニング後のモデルが保存される。
なので、それを指定してdemo.pyを実行する:

python demo.py --pretrained_model "result/snapshot_model.npz" hoge.jpg

すると、先ほどと同様の結果が出力された。
f:id:nixeneko:20170624093754p:plain

*1:GTX 1060で10時間ほど、CPUだと60日程度かかるらしい。

アウトラインがぶれるフォント

TrueType命令で遊ぶシリーズ。

  1. 初めてのTrueType命令: Windowsでは見えないフォントをつくる - にせねこメモ
  2. フォントサイズに合わせて回転するフォントを作る(1) - にせねこメモ
  3. PPEM・ポイントサイズを表示するフォント - にせねこメモ
  4. TrueType命令で三角関数(sin, cos)を計算する - にせねこメモ
  5. フォントサイズに合わせて回転するフォントを作る(2) - にせねこメモ
  6. TrueType命令で擬似乱数: 線形合同法 - にせねこメモ
  7. TrueType命令でビット演算 - にせねこメモ

擬似乱数を使って何か

前に擬似乱数を生成する関数を作ったので、それを使って何かやってみようというのが今回の目的である。
アウトラインをぶれさせてみたら、ランダムさが効果的に使えるのではないかと思うので、それをやってみる。

完成品

f:id:nixeneko:20170606181753p:plain
M+フォント(mplus-1p-regular.ttf)のアウトラインをガタガタにしている。フォントサイズによって文字の形が変わっているところに注目してほしい。

イラレとかで大きさを変えていくのを見ても面白い。


ダウンロード

なお、例によってTrueType命令を利用しているのでMacでは動かない。

実装

やっていることは単純で、各制御点の座標を読み出し、それに対して乱数生成関数によって生成した乱数を(適当にスケールして)足し合わせ、その座標に制御点を動かす。これをX, Y軸、および全ての制御点に対して行う。

さて、実際に実装してみる。

初期化

各種値の初期化は'prep'で行う。

まず、擬似乱数関数に与える初期seedを適当に初期化する。できるだけ乱雑になってほしかったのでMPPEM命令でPPEMを取得してそれを初期seedとした。これにより、フォントサイズによって乱数列が異なり、最終的に得られるアウトラインもフォントサイズ依存になる。
このseedをStorage Area 0番地に保存する。

次に、乱数のスケールに使用するscale factor sを準備する。
XまたはY軸方向に移動する距離の最大値(つまり、スケール後の乱数の絶対値の最大)をCVTテーブルの0番に書いておく。ここでは30 (FUnits)とした。
これをRCVTで読みだすとpixel単位になるので、それを pとする。
乱数の最大値は0x7FFFFFFFであるので、 s = \mathrm{0x7FFFFFFF} / pとする。
このとき、 p < 1の場合に、 sを計算する際にオーバーフローしてしまうため、それを防ぐために p := 1.0 \mbox{ (if } p < 1 \mbox{)}という処理を入れている。

読みだした乱数を sで割ることで、乱数(の絶対値)が0~pの範囲に含まれる様になる。
最後にscale factor  sをStorage Areaの1番地に保存している。

#seed初期化
PUSHB_1
 0
MPPEM   # 初期seed (>0)
WS      # StorageArea[0] = PPEM

#スケール係数の初期化
PUSHW_7 #0xFF|256.0|0xFF|256.0|0xFF|256.0|0x7F|
 255   # 0xFF
 16384 # 256.0
 255   # 0xFF
 16384 # 256.0
 255   # 0xFF
 16384 # 256.0
 127   # 0x7F
MUL     #0xFF|256.0|0xFF|256.0|0xFF|0x7F00|
ADD     #0xFF|256.0|0xFF|256.0|0x7FFF|
MUL     #0xFF|256.0|0xFF|0x7FFF00|
ADD     #0xFF|256.0|0x7FFFFF|
MUL     #0xFF|0x7FFFFF00|
ADD     #0x7FFFFFFF| #2147483647 = 0x7FFFFFFF
PUSHB_1
 0
RCVT    #0x7FFFFFFF|CVT[0]|
DUP     #0x7FFFFFFF|CVT[0]|CVT[0]|
PUSHB_1 #0x7FFFFFFF|CVT[0]|CVT[0]|1.0|
 64
LT      #0x7FFFFFFF|CVT[0]|CVT[0]<1.0|
IF      #0x7FFFFFFF|CVT[0]|  #if CVT[0]<1.0:
 POP     #0x7FFFFFFF|
 PUSHB_1 #0x7FFFFFFF|1.0|
  64
EIF
DIV     #0x7FFFFFFF/CVT[0]|
PUSHB_1 #0x7FFFFFFF/CVT[0]|1|
 1
SWAP    #1|0x7FFFFFFF/CVT[0]|
WS    #StorageArea[1] = 0x7FFFFFFF/CVT[0]

Storage Areaへの割り付け

番号 value
0 seed (PPEMで初期化)
1 scale factor

関数群

関数0: 乱数の生成

前に作った乱数関数(改良版の関数1)をもってきて、関数0とする。
関数番号以外は同一なので詳しくは前の記事を参照。

実行する度に1~0x7FFFFFFFの乱数列が返ってくる。

関数1: スケールした乱数の取得

関数0を呼び出して乱数を取得し、'prep'でStorage Areaの1番地に保存したscale factor  sで割ることでスケールする。
この際、取得した乱数を1ビット左シフトし、取得した乱数の上から2ビット目を符号ビットとして扱うことで、負の数も得られる様にしている。
最終的に得られる乱数は、最大動き幅pについて-p~pの範囲となるはず。

/* Function 1: returns scaled random value */
PUSHB_1
 1
FDEF          /* ..|      *//* ← initial stack */
PUSHB_1
 0
CALL          /* ..|rand| */
DUP
ADD           /* ..|rand<<1| */ /* to make it signed */
PUSHB_1
 1
RS            /* ..|rand<<1|s| */ /*s = StorageArea[1], scaling factor */
DIV           /* ..|(rand<<1)/s| */
ENDF

関数2: 指定された制御点をランダムな大きさだけ動かす

指定された制御点番号kに対応する制御点のX座標、Y座標を取得し、それに対して関数1で取得したスケール済み乱数を足しあわせ、計算された座標にSCFSで移動している。

/* Function 2: moves the control point k         */
/* initial stack ..| k |       k: 編集する制御点番号                           */
/* final stack   ..|                                                      */
PUSHB_1
 2
FDEF     /* ..|k| */
DUP      /* ..|k|k| */
DUP      /* ..|k|k|k| */
DUP      /* ..|k|k|k|k| */
SVTCA[x-axis]                                  /* X座標に設定 */
GC[cur]  /* ..|k|k|k|x_k| */
PUSHB_1
 1
CALL     /* ..|k|k|k|x_k|rand| */ /* get scaled random value*/
ADD      /* ..|k|k|k|x_k+rand| */
SCFS     /* ..|k|k| */        /* 制御点kのX座標を x_k + rand に */

SVTCA[y-axis]                                  /* Y座標に設定 */
GC[cur]  /* ..|k|y_k| */
PUSHB_1
 1
CALL     /* ..|k|y_k|rand| */
ADD      /* ..|k|y_k+rand| */
SCFS     /* ..| */            /* 制御点kのY座標を y_k + rand に */

ENDF

関数3: カウントアップしつつ関数2を呼び出す

LOOPCALLで呼び出される用の関数。スタックトップの番号(カウンター)を引数として関数2を呼び出し、最後にカウンターを1増やす。

/* function 3: call func2 with value n and increment n */
/* initial stack: ..|n|   */
/* final stack:   ..|n+1| */
PUSHB_1
 3
FDEF      /* ..|n| */ /* repeat n times */
DUP       /* ..|n|n| */
PUSHB_1
 2
CALL      /* ..|n| */ /* moves the point n */
PUSHB_1   /* ..|n|1| */
 1
ADD       /* ..|n+1| */
ENDF

関数4: LOOPCALLで関数3を指定回数実行する

グリフから呼び出され、指定回数(=制御点の個数回)関数3を実行する。関数3では呼び出す度にスタックトップの値(カウンター)を1ずつ増やしていくので、すべての制御点に対して操作が行われることになる。

/* function 4: call func3 n times */
/* initial stack: ..|n|  n: number of repetition */
/* final stack:   ..|  */
PUSHB_1
 4
FDEF      /* ..|n| */ /* n: num of points in the glyph */
PUSHB_1   /* ..|n|0| */ /* initialize the counter by 0 */
 0
SWAP      /* ..|0|n| */
PUSHB_1   /* ..|0|n|3| */
 3
LOOPCALL  /* ..| */ /* call function 3, n times */
ENDF

グリフ

関数4を、グリフから次のように呼び出す。
例えば、グリフに含まれる制御点の個数が10個であれば、

PUSHB_2
 10
 4
CALL

となる。

すべてのグリフへの適用

'cvt ', 'fpgm', 'prep'を設定して書き出したフォントAmovepointsrandom.ttfをttxを使ってXMLファイル(.ttx)にダンプする。

ttx Amovepointsrandom.ttf

その後、次のPythonプログラムを実行し、Amovepointsrandom.ttxのすべてのグリフに対して関数を呼び出すTrueType命令を適用した結果(Amovepointsrandom-out.ttx)を得る。
やっていることは、グリフ毎に制御点の個数を調べ、その値を組み込んだ関数呼び出しの命令を付加するという感じ。

#!/usr/bin/env python3
# conding: utf-8

import xml.etree.ElementTree as ET

INFILE = "Amovepointsrandom.ttx"
OUTFILE = "Amovepointsrandom-out.ttx"
xmltree = ET.parse(INFILE)
xmlroot = xmltree.getroot()

for glyph in xmlroot.find('glyf').findall('TTGlyph'):
    cnt = 0
    for contour in glyph.findall('contour'):
        cnt += len(contour.findall('pt'))
    if cnt > 0:
        prog ="""
          PUSHB[ ]      /* 2 values pushed */
          {} 4
          CALL[ ]       /* CallFunction */
          """.format(cnt)
        glyph.find('instructions').find('assembly').text = prog

with open(OUTFILE, 'w') as w:
    w.write('<?xml version="1.0" encoding="UTF-8"?>\n')
    xmlstr = ET.tostring(xmlroot, method='xml', encoding="unicode")
    #xmltree.write(OUTFILE)
    w.write(xmlstr)

最後に得られた.ttxファイルをttxで.ttfに変換する。

ttx Amovepointsrandom-out.ttx

これによって生成されるAmovepointsrandom-out.ttfが完成品である。

コード

'cvt '

number value
0 30

'fpgm'

PUSHB_1
 0
FDEF
PUSHB_1
 0
RS
DUP
PUSHW_2
 30000
 14488
ADD
DUP
ROLL
SWAP
PUSHW_1
 4096
MUL
DIV
ROLL
ROLL
PUSHB_1
 3
CINDEX
PUSHW_1
 4096
MUL
MUL
SUB
PUSHW_3
 4096
 30000
 18271
ADD
MUL
MUL
SWAP
PUSHW_2
 3399
 4096
MUL
MUL
SUB
DUP
PUSHB_1
 0
GT
IF
ELSE
PUSHW_7
 255
 16384
 255
 16384
 255
 16384
 127
MUL
ADD
MUL
ADD
MUL
ADD
ADD
EIF
DUP
PUSHB_1
 0
SWAP
WS
ENDF
PUSHB_1
 1
FDEF
PUSHB_1
 0
CALL
DUP
ADD
PUSHB_1
 1
RS
DIV
ENDF
PUSHB_1
 2
FDEF
DUP
DUP
DUP
SVTCA[x-axis]
GC[cur]
PUSHB_1
 1
CALL
ADD
SCFS
SVTCA[y-axis]
GC[cur]
PUSHB_1
 1
CALL
ADD
SCFS
ENDF
PUSHB_1
 3
FDEF
DUP
PUSHB_1
 2
CALL
PUSHB_1
 1
ADD
ENDF
PUSHB_1
 4
FDEF
PUSHB_1
 0
SWAP
PUSHB_1
 3
LOOPCALL
ENDF

'prep'

PUSHB_1
 0
MPPEM
WS
PUSHW_7
 255
 16384
 255
 16384
 255
 16384
 127
MUL
ADD
MUL
ADD
MUL
ADD
PUSHB_1
 0
RCVT
DUP
PUSHB_1
 64
LT
IF
POP
PUSHB_1
 64
EIF
DIV
PUSHB_1
 1
SWAP
WS

三相3線式の交流電源の1線を接地した場合の電圧

三相3線式の交流電源について考えてみる。

三相交流とは、3種類の位相の単相交流電源を組み合わせたもので、どの2つをとっても位相差の大きさが \displaystyle 120\mbox{°} = \frac{2}{3}\piとなるというものらしい。

試しに、3つの相R, S, Tの時間 tにおける電圧 V_R(t), V_S(t), V_T(t)を次のようにしてみる。ここで、 V_0は最大電圧、 \omega = 2 \pi fは角周波数( fは電源の周波数)を示す。
 \displaystyle \begin{align}
V_R(t) &= V_0 \sin \left(\omega t - \frac{2}{3}\pi\right) \\
V_S(t) &= V_0 \sin \omega t \\
V_T(t) &= V_0 \sin \left(\omega t + \frac{2}{3}\pi\right)
\end{align}
グラフにかくと次のようになる。
f:id:nixeneko:20170601042038p:plain


さて、三相3線式ではこのうち一線を接地する。ここでは、Sの線を接地してみることにする。
電圧は接地したところを基準に測られるので、接地後のSの電圧は V_S'(t) = 0である。
同様に、接地後の各電圧は、
 \displaystyle \begin{align}
V_R'(t) &= V_0 \sin \left(\omega t - \frac{2}{3}\pi\right) - V_0 \sin \omega t\\
V_S'(t) &= 0 \\
V_T'(t) &= V_0 \sin \left(\omega t + \frac{2}{3}\pi\right) - V_0 \sin \omega t
\end{align}
として測定されるはずである。


ここで、加法定理から
 \displaystyle \begin{align}
\sin \left(\omega t - \frac{2}{3}\pi\right) &= \sin\omega t \cos\frac{2}{3}\pi - \cos\omega t \sin\frac{2}{3}\pi \\
&= -\frac{1}{2}\sin\omega t -\frac{\sqrt{3}}{2} \cos\omega t \\
\sin \left(\omega t + \frac{2}{3}\pi\right) &= \sin\omega t \cos\frac{2}{3}\pi + \cos\omega t \sin\frac{2}{3}\pi \\
&= -\frac{1}{2}\sin\omega t +\frac{\sqrt{3}}{2} \cos\omega t 
\end{align}
であるから、よって
 \displaystyle \begin{align}
V_R'(t) &= V_0 \sin \left(\omega t - \frac{2}{3}\pi\right) - V_0 \sin \omega t\\
&= V_0 \left( -\frac{1}{2}\sin\omega t -\frac{\sqrt{3}}{2} \cos\omega t - \sin \omega t \right) \\
&= V_0 \left( -\frac{3}{2}\sin\omega t -\frac{\sqrt{3}}{2} \cos\omega t \right)
\end{align}
となる。
三角関数の合成より
 \displaystyle \begin{align}
 -\frac{3}{2}\sin\omega t -\frac{\sqrt{3}}{2} \cos\omega t = \sqrt{3}\sin\left(\omega t - \frac{5}{6}\pi\right)
\end{align}
したがって、
 \displaystyle \begin{align}
V_R'(t) = \sqrt{3} V_0 \sin\left(\omega t + \frac{5}{6}\pi\right)
\end{align}
である。
同様に、
 \displaystyle \begin{align}
V_T'(t) = \sqrt{3} V_0 \sin\left(\omega t - \frac{5}{6}\pi\right)
\end{align}
も計算できる。

まとめると、
 \displaystyle \begin{align}
V_R'(t) &= \sqrt{3} V_0 \sin\left(\omega t + \frac{5}{6}\pi\right) \\
V_S'(t) &= 0 \\
V_T'(t) &= \sqrt{3} V_0 \sin\left(\omega t - \frac{5}{6}\pi\right)
\end{align}
で、これらは対地電圧になっている。

グラフにすると次のような感じになる。
f:id:nixeneko:20170601042205p:plain

また、 V_R'(t) - V_T'(t)についても計算してみると、
 \displaystyle \begin{align}
&V_R'(t) - V_T'(t) \\
&= V_0 \sin \left(\omega t - \frac{2}{3}\pi\right) - V_0 \sin \omega t - \left\{V_0 \sin \left(\omega t + \frac{2}{3}\pi\right) - V_0 \sin \omega t \right\} \\
&= V_0 \left\{ \sin \left(\omega t - \frac{2}{3}\pi\right) - \sin \left(\omega t + \frac{2}{3}\pi\right) \right\} \\
&= V_0 \left\{ -\frac{1}{2}\sin\omega t -\frac{\sqrt{3}}{2} \cos\omega t - \left(-\frac{1}{2}\sin\omega t +\frac{\sqrt{3}}{2} \cos\omega t \right) \right\} \\
&= -\sqrt{3}V_0\cos\omega t \\
&= \sqrt{3}V_0\sin\left(\omega t - \frac{\pi}{2}\right)
\end{align}
となり、振幅は V_R'(t) V_T'(t)と等しい \sqrt{3}V_0となっている。

ってなわけで、有効電圧は最大電圧の \displaystyle \frac{1}{\sqrt{2}}なので、三相3線式でどの2線をとっても線間電圧が \displaystyle \sqrt{\frac{3}{2}}V_0になるようである。


これらのことから、電源電圧が200Vの場合は、
 \displaystyle \begin{align}
\sqrt{\frac{3}{2}}V_0 = 200 \mbox{ [V]} \\
\therefore V_0 = 200 \times \sqrt{\frac{2}{3}} \approx 163 \mbox{ [V]}
\end{align}
となることがわかる。

Chainerを触ってみる: XOR関数を訓練する

※この記事ではPython 3.5.3とChainer v1.20.0.1を使っています。

調べながら書いています。間違い等あればご指摘願います。

はじめに

『ゼロから作るDeep Learning』を一通り読んで、実際にライブラリ使って機械学習してみようという段、Chainerを触ってみたものの、deep learning界のHello worldとも呼ばれるMNISTのexampleを見ても、なにをしているのかよくわからない。

なので、もっと単純なネットワークを作って訓練してみることで、chainerの使い方を勉強して、MNISTの例で何をやってるのかが理解できるようにすることを目指した。

Chainerの基本的なオブジェクト

chainer.Variable

  • 変数に相当する。
  • 後述するFunctionやLinkは入力がVariableで、出力もVariableである。
    • 入力に関してはnumpy.ndarrayを与えると内部で変換してくれるのでどちらでも良い感じかも。
  • .dataでデータ配列numpy.ndarray (あるいはcupy.ndarray)がとれる。

chainer.Function

  • 関数に相当する。訓練すべきパラメータをもたない。
  • chainer.functions以下に(chainer.Functionを継承したものが)色々定義されている。
  • activation関数, loss関数, accuracy関数などなど

chainer.Link

  • 訓練すべきパラメータをもつ関数を表す。
  • chainer.links以下にいろいろ定義されている。例:
    • Linear (全結合)
    • Convolution2D (2次元CNN)
  • ミニバッチを入力とするので入力の次元に注意。

パーセプトロン

f:id:nixeneko:20170531234148p:plain
いくつかの入力を重みづけして足し合わせ、さらに定数(バイアス)を足し、その結果に適当な非線形関数(activation関数という)を適用したものを(単純)パーセプトロンという。

AND関数を作ってみる

ANDは単純パーセプトロンで実装できる。
f:id:nixeneko:20170601000053p:plain

  • AND関数は入力2ノード、出力1ノードの関数であり、全結合のLinearを使って実装できる。
  • initialW, initial_biasを設定することで、ネットワークのパラメータを決め打ちしている。
  • Linearはミニバッチを入力とすることに注意する。
  • activation関数にはReLUを使った。
import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np

# L.Linear(in_size, out_size, wscale=1, bias=0, nobias=False, 
#          initialW=None, initial_bias=None)
W = np.array([[1,1]]) # weight
b = -1.0 # bias
a = L.Linear(2, 1, initialW=W, initial_bias=b) # full-connected layer
f = F.relu # activation function
def AND(x): # 順伝播を定義
    return f(a(x))

# use network as a function
x = np.array([[1,1],[0,1],[1,0],[0,0]], dtype=np.float32) # input data
print(AND(x).data)
#[[ 1.]
# [ 0.]
# [ 0.]
# [ 0.]]

実行すると、[1,1]の入力に対しては[1]が、それ以外には[0]が帰ってきていることがわかる。これはまさにANDの挙動である。

多層パーセプトロン(MLP)

単純パーセプトロンを複数、層状に組み合わせたものを多層パーセプトロンという。各層はベクトルを入力としベクトルを出力とする関数になっていて、全体で複雑な非線形関数を実現できる。

次は3層の例で、バイアスやactivation関数の記載は省略した。*1
f:id:nixeneko:20170531234740p:plain

XOR関数を作ってみる

同様にXOR関数を作ってみる。
XORは線形分離不可能であるので単純パーセプトロンでは再現できず、全結合層を2つ用意する必要がある。中間層(隠れ層)のノード数は2にした。
f:id:nixeneko:20170601000835p:plain

  • 全結合層(Linear)2層
  • 入力層はノード2つ, 中間層ノード2つ, 出力層ノード1つ
  • 層の重みをinitialW, initial_biasで設定。うまくいく値を決め打ちで使った。
  • activation関数として第一層にはReLUを使った。(第二層はidentity functionを使ったというか、そのままの値)
  • 前層の出力と次層の入力でノードの数を合わせる。
    • Linearの第1引数(in_size)をNoneにすると、最初に入力が与えられたときに自動的にサイズを推定してくれる。
import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np

# L.Linear(in_size, out_size, wscale=1, bias=0, nobias=False, 
#          initialW=None, initial_bias=None)
# 第1層
W1 = np.array([[1,-1],[-1,1]]) # weight parameters
b1 = 0.0 # bias
a1 = L.Linear(2, 2, initialW=W1, initial_bias=b1) #入力2, 出力2
f1 = F.relu # activation function

# 第2層
W2 = np.array([[1,1]]) # weight parameters
b2 = 0.0 # bias
a2 = L.Linear(2, 1, initialW=W2, initial_bias=b2) #入力2, 出力1

# 順伝播を定義
def XOR(x):
    h1 = f1(a1(x))
    return a2(h1)

# 試してみる
x = np.array([[1,1],[1,0],[0,1],[0,0]], dtype=np.float32)
print(XOR(x).data)
#[[ 0.]
# [ 1.]
# [ 1.]
# [ 0.]]

これも実際に動かしてみると、入力の[1,1]と[0,0]には[0]を、[0,1]と[1,0]には[1]を返していて、XORを表現した関数となっている。

chainer.Chain

ChainでLinkオブジェクト(=パラメータ付き関数)をまとめることができる。
実際にはChainを継承したクラスをつくって、順伝播をその中の__call__メソッドに定義するという使い方をするようだ。
先ほどのXORを書き換えると次のようになる。

import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np

# parameters
W1 = np.array([[1,-1],[-1,1]])
b1 = 0.0
W2 = np.array([[1,1]])
b2 = 0.0

class Xor(chainer.Chain):  #chainer.Chainを継承したクラスで全体の関数を表現
    def __init__(self): # constructor
        #親クラスchainer.Chainのコンストラクタでlinkをまとめる
        super(Xor, self).__init__( 
            a1 = L.Linear(2, 2, initialW=W1, initial_bias=b1),
            a2 = L.Linear(2, 1, initialW=W2, initial_bias=b2)
        )

    def __call__(self, x): # 順伝播の定義
        h1 = F.relu(self.a1(x))
        return self.a2(h1)

# 定義したものを使ってみる
XOR = Xor()
x = np.array([[1,1],[1,0],[0,1],[0,0]], dtype=np.float32)
print(XOR(x).data)
#[[ 0.]
# [ 1.]
# [ 1.]
# [ 0.]]

chainer.Chainを継承したXorクラスのインスタンスに入力を与えて呼び出すと計算結果が出力される。

さて、ここではXORネットワークのパラメータを手で予め求めた値に設定したが、これをデータから学習させてみたい。

Chainerにおける訓練の流れ

多分こんな感じ。

  1. 学習用データセットの用意
  2. lossを計算するネットワークモデルを用意(これを訓練する)
  3. Optimizerを用意してモデルをセット
  4. Updaterを用意
  5. Trainerを準備
  6. 訓練ループの実行

データセットはモデルを定義するのと同時に考えないといけないと思うので最初に入れているが、MNISTのサンプルなどではOptimizerとUpdaterの間に入っている。

これより後のコードでは

import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np
from chainer.training import extensions

を前提とする。
最後の方に訓練用のコード全体を載せておいた。

1. 学習用データセットの準備

学習用データセットは、入力データと正解データの組である。
今回はXORの2入力がそれぞれ0か1のどちらかの値をとるので、入力は4種類ある。それぞれに対して対応する出力を用意する。
なお、全結合Linearが入力・出力としてnumpy.float32を要求するのでデータ形式をそのように設定している。

indata = np.array([[0,0],[0,1],[1,0],[1,1]], dtype=np.float32)
labels = np.array([ [0],  [1],  [1],  [0] ], dtype=np.float32)

dataset = chainer.datasets.TupleDataset(indata, labels)
train_iter = chainer.iterators.SerialIterator(dataset, 4) # 訓練用. batch size = 4
test_iter = chainer.iterators.SerialIterator(dataset, 4, repeat=False, shuffle=False)
                                                          # 評価用. batch size = 4
  • TupleDatasetは内部で単に2つの入力をタプルにするくらいのことしかしてなくて、あまり本質ではないっぽい。
  • 重要なのはその後のSerialIteratorなどのイテレータで、これがあるとミニバッチ単位でデータを取り出すのが簡単になる。
  • 訓練用と評価用のイテレータを用意する。
    • 今回は訓練セットと評価セットを全く同じにしているが、普通は適当に分ける。
    • どちらもミニバッチサイズを4とした。
  • 訓練用のイテレータは後のUpdaterに渡される。
  • 評価用のイテレータは特になくても訓練はできるはずだけれど、訓練中のモデルの途中経過のaccuracyとかを見るのに使われる。
    • repeat=False, shuffle=False を設定しておく(デフォルトではTrue)。特に repeat=False としないと評価が終わらなくなる。

2. lossを返すモデルの用意

訓練するモデルは、呼び出すと、つまり、 .__call__ メソッドが:

  • 用意したデータセット(のミニバッチ)を入力として受け付ける
  • lossを返す

ようなネットワークである必要がある*2

#L.Linear(in_size, out_size, wscale=1, bias=0, nobias=False,
#         initialW=None, initial_bias=None)
class Xor(chainer.Chain):  #chainer.Chainを継承したクラスで全体の関数を表現
    def __init__(self): # constructor
        super(Xor, self).__init__( #chainer.Chainでlinkをまとめる
            l1 = L.Linear(2, 10),
            l2 = L.Linear(10, 1)
        )

    def __call__(self, x): #forward propagation definition
        h1 = F.relu(self.l1(x))
        return self.l2(h1)

my_xor = Xor() # 最終的に求めたい関数を示すモデル
accfunc = lambda x, t: F.sum(1 - abs(x-t))/x.size #ないと動かないのでとりあえず用意
model = L.Classifier(my_xor, lossfun=F.mean_squared_error, accfun=accfunc) 
                                      # datasetのデータを入力するとlossを返すモデル
  • ここで、Xorは上で定義したクラスとほとんど同じだけれど、中間層のノードを10に増やしてある*3
  • また、initialWやinitial_biasの設定を行わず、パラメータはランダムな値で初期化される様になっている。
  • Xorのインスタンスmy_xorを作成する。
    • これは用意したデータセットのデータをそのまま入力とすることができず、さらに呼び出すとXORの計算結果を返すため、そのままではOptimizerに渡すことはできない*4
    • そのため、用意したデータセットを入力として受け付けることができ、lossを返すようなネットワーク(ここではClassifier)でmy_xorをくるむ。
    • Classifierのデフォルトのaccuracy評価関数はnumpy.int32を入力としてとるので、代わりにnumpy.float32を入力としてとれる関数accfuncを用意して設定している。訓練には必須ではないのでダミーでもよい。

ここでは、最終的に得たい関数を表すネットワークをlossを返すネットワークでくるんでいるが(MNISTのサンプルも同様である)、他方PaintsChainerなどでは、.callメソッドで最終的に得たい関数を計算するようにして、.__call__ではlossを計算して返すようなモデルを使う実装になっていた。

3. Optimizerの用意

lossを見てモデルのパラメータを更新するのがOptimizerである。lossの計算もOptimizerの中で行うことができる*5
これはUpdaterから呼び出される。

  • 色々な最適化手法がchainer.optimizers以下に用意されていて、適当なものを使う(ここではAdam)。
  • インスタンスを生成して.setupメソッドに2.で用意したlossを返すモデルをセットする。
  • また、必要があればOptimizerにフック関数を設定して、例えば重み減衰による正則化などが設定できる。
optimizer = chainer.optimizers.Adam()
optimizer.setup(model) #lossを返すモデルをセット
#optimizer.add_hook(chainer.optimizer.WeightDecay(0.0005))

4. Updaterの用意

Updaterの役割は、訓練用データセットイテレータからデータのミニバッチを読み込んで、lossを計算してoptimizerでモデルのパラメータを更新することである。

今回はStandardUpdaterをそのまま使うが、複雑なモデルとかだと、StandardUpdaterを継承したUpdaterを作成してパラメータを更新するということもできる。

updater = chainer.training.StandardUpdater(train_iter, optimizer)

5. Trainerの準備

Trainerは全体の訓練ループを管理し、進捗状況表示やログ出力などを行う。

trainer = chainer.training.Trainer(updater, (3000, 'epoch'), out="test_result")

trainer.extend(extensions.Evaluator(test_iter, model))
trainer.extend(extensions.LogReport())
trainer.extend(extensions.PrintReport(['epoch', 'main/loss', 'main/accuracy']))
trainer.extend(extensions.ProgressBar())

まずTrainerを作成する。(3000, 'epoch')は終了条件で、3000エポックで終了するという意味。
その後.extendメソッドにchainer.training.extensions以下にあるようなExtensionを渡して実行させることができ、ここでは

  • Evaluatorで現在のモデルの評価
  • LogReportでログファイルへの書き出し
  • PrintReportで指定したデータの画面表示
    • 'epoch'がエポック, 'main/loss'はloss, 'main/accuracy'はaccuracyに対応する。
      • main/はおまじないみたいな感じで、Optimizerが訓練対象とするモデルを表すらしい。
  • ProgressBarで進捗状況を示すプログレスバーの表示

などを設定している。他にも、訓練中のモデルやOptimizerの状態を定期的にファイルに保存することもでき、時間のかかる訓練では、途中でやめても保存した状態から訓練を再開できるようにしておくと良さそう。

6. 訓練ループの実行

Trainerを動かして訓練を行う。

trainer.run()

lossが下がっていくのを眺めながら、訓練が終わるまで待つ。
初期値によっては局所解に囚われてしまってlossが減少しなくなったりするので、その場合は何回かやり直してみるといいかもしれない。

訓練のソース全体

import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np
from chainer.training import extensions

# データセットの準備
indata = np.array([[0,0],[0,1],[1,0],[1,1]], dtype=np.float32)
labels = np.array([ [0],  [1],  [1],  [0] ], dtype=np.float32)

dataset = chainer.datasets.TupleDataset(indata, labels)
train_iter = chainer.iterators.SerialIterator(dataset, 4) # 訓練用 batch size = 4
test_iter = chainer.iterators.SerialIterator(dataset, 4, repeat=False, shuffle=False)
                                                          # 評価用 batch size = 4

#L.Linear(in_size, out_size, wscale=1, bias=0, nobias=False,
#         initialW=None, initial_bias=None)
class Xor(chainer.Chain):  #chainer.Chainを継承したクラスで全体の関数を表現
    def __init__(self): # constructor
        super(Xor, self).__init__( #chainer.Chainでlinkをまとめる
            l1 = L.Linear(2, 10),
            l2 = L.Linear(10, 1)
        )

    def __call__(self, x): #forward propagation definition
        h1 = F.relu(self.l1(x))
        return self.l2(h1)


my_xor = Xor() # 訓練したい関数を示すモデル
accfun = lambda x, t: F.sum(1 - abs(x-t))/x.size #ないと動かないのでとりあえず用意
model = L.Classifier(my_xor, lossfun=F.mean_squared_error, accfun=accfun) 
                                     # datasetのデータを入力するとlossを返すモデル

optimizer = chainer.optimizers.Adam()
optimizer.setup(model) #lossを返すモデルをセット
#optimizer.add_hook(chainer.optimizer.WeightDecay(0.0005))

# 訓練の準備
updater = chainer.training.StandardUpdater(train_iter, optimizer)
trainer = chainer.training.Trainer(updater, (3000, 'epoch'), out="test_result")
trainer.extend(extensions.Evaluator(test_iter, model))
trainer.extend(extensions.PrintReport(['epoch', 'main/loss', 'main/accuracy']))
trainer.extend(extensions.LogReport())
trainer.extend(extensions.ProgressBar())
#print('start')
trainer.run()
#print('done')

訓練したモデルを使う

my_xorが最終的に得たい関数であるXORを表すものになっているはずである。

x = np.array([[1,1],[1,0],[0,1],[0,0]], dtype=np.float32)
result = my_xor(x)
print(result.data)
#出力例:
#[[  4.03262675e-07]
# [  9.99999642e-01]
# [  9.99999642e-01]
# [  3.80910933e-07]]

訓練する度に値は変わるが、うまいこといっていれば出力に0か1に近い値がでてくる。

また、パラメータは次のようにして見ることができる。

print("l1.W:")
print(my_xor.l1.W.data)
print("l1.b:")
print(my_xor.l1.b.data)
print("l2.W:")
print(my_xor.l2.W.data)
print("l2.b:")
print(my_xor.l2.b.data)
#出力例:
#l1.W:
#[[-0.7456618   1.03302336]
# [ 1.10809541 -0.91810882]
# [ 0.47287568  0.47287571]
# [-0.73910308 -0.74221992]
# [ 0.88031632 -0.73166984]
# [-0.73261809 -1.85475707]
# [ 0.70534587  0.70534694]
# [-0.15760551 -0.42765161]
# [ 0.07897222  0.93817586]
# [ 1.17522132  0.32095835]]
#l1.b:
#[ -1.48954052e-06  -2.97201019e-09  -4.72875684e-01   0.00000000e+00
#  -5.84823523e-09   0.00000000e+00  -7.05345809e-01   0.00000000e+00
#   9.12837006e-09   6.72058091e-02]
#l2.W:
#[[ 0.49185807  0.19722475 -0.84104687  0.05196132  0.75117517 -0.09614659
#  -0.7157371   0.72252452  0.50084352  0.06860919]]
#l2.b:
#[-0.00461056]

訓練したモデルのパラメータを保存する

chainer.serializers以下に、オブジェクトをファイルに保存・読み込みする関数が用意されている。保存形式にはNPZ (numpyの形式)とHDF5の2種類があり、それぞれ別々の関数によって読み書きを行う。

#NPZ用
chainer.serializers.save_npz(filename, obj, compression=True)
chainer.serializers.load_npz(filename, obj)

#HDF5用
chainer.serializers.save_hdf5(filename, obj, compression=4)
chainer.serializers.load_hdf5(filename, obj)

訓練の途中経過のスナップショットをとってそこから再開するといった場合にはtrainerを保存・読み込みすると良いようだ。
また、さらなる訓練が要らなければ、今回はmy_xorを保存しておけばよさそうである。
つまり、保存は

chainer.serializers.save_npz('my_xor.npz', my_xor)

で、これを別のファイルで読み出すのであれば、

import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np

class Xor(chainer.Chain):  #chainer.Chainを継承したクラスで全体の関数を表現
    def __init__(self): # constructor
        super(Xor, self).__init__( #chainer.Chainでlinkをまとめる
            l1 = L.Linear(2, 10),
            l2 = L.Linear(10, 1)
        )

    def __call__(self, x): #forward propagation definition
        h1 = F.relu(self.l1(x))
        return self.l2(h1)

my_xor = Xor()

chainer.serializers.load_npz('my_xor.npz', my_xor)

などとすればよい。

*1:outputは第3層の出力と等しい。色々図示されたものを見るに、最後のoutputは第3層と一つにまとめて書かれるのが普通な感じがある。

*2:必ずしもそうでなくても問題ないのだが、Trainerを使った訓練で扱うには色々面倒っぽい

*3:これは訓練したら極小解に陥ってしまって、何を入力しても0.5を返すようになってしまうことが多かったため。

*4:Optimizerに渡すことはできるが、独自のUpdaterを用意してlossを計算するなどしないといけない。

*5:loss関数と入力を与えて呼び出すと、内部でlossを計算し、パラメータを更新する

TrueType命令でビット演算

TrueType命令で遊ぶシリーズ。

  1. 初めてのTrueType命令: Windowsでは見えないフォントをつくる - にせねこメモ
  2. フォントサイズに合わせて回転するフォントを作る(1) - にせねこメモ
  3. PPEM・ポイントサイズを表示するフォント - にせねこメモ
  4. TrueType命令で三角関数(sin, cos)を計算する - にせねこメモ
  5. フォントサイズに合わせて回転するフォントを作る(2) - にせねこメモ
  6. TrueType命令で擬似乱数: 線形合同法 - にせねこメモ

ビット演算

TrueType命令にはビット演算がない。ヒンティングには重要でないからだろう。
それでも、変なことしようとするとあると便利なのがビット演算である。

それっぽい名前の命令としてNEG, AND, ORなどがあるが、それぞれ数値の符号反転 (-a), 論理AND (a&&b), 論理OR (a||b)である。
仕方がないので自分で実装してしまおうというのがこの記事の目的である。

bitwise NOT

符号反転して-1を足すだけ。(別に1を引くのでもいい)

NEG
PUSHW_1 -1
ADD


左シフト

1ビット左シフト

要するに2を掛ければよいのだけれど、掛け算命令使うとF26DOT6との折り合いを付けないといけなくて面倒なのとオーバーフロー時の動作が心配なので、同じものを複製して足し合わせて計算することにする。

DUP
ADD

多ビット左シフト

1ビット左シフトをシフトしたいだけ繰り返せばよい。

右1ビット算術シフト

2.0 (つまり、F26DOT6の6ビット分だけ下駄を履かせて128)で割る。ただ、任意の数について適切に扱えるかはわからない。

PUSHB_1 128
DIV

多ビット(nビットとする)の場合は 2^n \times 64で割るといいと思う。

右1ビット論理シフト

最上位ビットに気を付けて2で割る。
まず数値が負かどうかを見て最上位ビット(MSB)によって場合分けする。

  • 非負(MSB==0)ならそのまま2で割る。
  • 負(MSB==1)なら0x80000000を足してMSBを0にした後、2で割り、もともとのMSBに相当する部分のビットを立てる(つまり、0x40000000を足す)。
             #x
DUP          #x|x|
PUSHB_1 0    #x|x|0|
LT           #x|x<0| #MSB
IF           #x| #when MSB==1
 PUSHW_3
  16384  #0x100*64
  16384  #0x100*64
  -32768 #0x8000
 MUL
 MUL         #x|0x80000000|
 ADD         #x+0x80000000|
 PUSHB_1 128 #x+0x80000000|2.0|
 DIV         #(x+0x80000000)/2|
 PUSHW_3
  16384 #0x100*64
  16384 #0x100*64
  16384 #0x4000
 MUL
 MUL         #(x+0x80000000)/2|0x40000000|
 ADD         #(x+0x80000000)/2+0x40000000|
ELSE         #x| #when MSB==0
 PUSHB_1 128 #x|2.0|
 DIV         #x/2|
EIF

これを一度計算すると最上位ビットが0であることが保証されるので、多ビット(nビット)シフトする場合はこれを実行した後に 2^{n-1} \times 64で割るといいと思う。

bitwise AND

AND, ORって命令はあるけれど、0をfalse, 0以外をtrueとしたlogical演算(演算結果がtrueの場合は1, falseの場合は0が返される)であって、bit演算ではない。

MSBは0だと非負整数、1だと負正数となる。そのため、値が0より小さければそのMSBが1だとわかる。また、左シフトもできるので、1ビットごと計算をしていけばよさそうだ。

演算対象の2つの数値についてMSBがわかるので、それらのMSBに対して(logical) ANDをとる。返す値を1ビット左シフトしながらその結果を加算する(LSBとなる)。
演算対象の2つの数値について1ビット左シフトして、この全体のプロセスを32回繰り返すと、bitwise ANDの結果が得られる。

#func0 - calculates bitwise AND
#initial stack: ...|x|y|
#final stack:   ...|x & b| (bitwise and)
PUSHB_1 0
FDEF        #x|y|
PUSHB_1 0   #x|y|r=0| #initialize return value r by 0
ROLL        #y|r|x|
PUSHB_1 32  #y|r|x|c=32| #initialize counter c by 32
ROLL        #y|x|c|r|
ROLL        #y|c|r|x|
PUSHB_1 4
MINDEX      #c|r|x|y|
#ここからループする(32回実行) #c|r|x|y|   #基準の形*1
# jump_to_here:
PUSHB_1 4
MINDEX      #r|x|y|c|
PUSHB_1 1   #r|x|y|c|1|
SUB         #r|x|y|c-1|
PUSHB_1 4
MINDEX      #x|y|c-1|r|
PUSHB_1 4
MINDEX      #y|c-1|r|x|
PUSHB_1 4
MINDEX      #c-1|r|x|y|
DUP         #c-1|r|x|y|y| 
PUSHB_1 0   #c-1|r|x|y|y|0|
LT          #c-1|r|x|y|y<0|
PUSHB_1 3
CINDEX      #c-1|r|x|y|y<0|x|
PUSHB_1 0   #c-1|r|x|y|y<0|x|0|
LT          #c-1|r|x|y|y<0|x<0|
AND         #c-1|r|x|y|y<0 && x<0|
PUSHB_1 4
MINDEX      #c-1|x|y|y<0 && x<0|r|
DUP         #c-1|x|y|y<0 && x<0|r|r|
ADD         #c-1|x|y|y<0 && x<0|r<<1|
ADD         #c-1|x|y|r<<1 + (y<0 && x<0)|
ROLL        #c-1|y|r<<1 + (y<0 && x<0)|x|
DUP         #c-1|y|r<<1 + (y<0 && x<0)|x|x|
ADD         #c-1|y|r<<1 + (y<0 && x<0)|x<<1|
ROLL        #c-1|r<<1 + (y<0 && x<0)|x<<1|y|
DUP         #c-1|r<<1 + (y<0 && x<0)|x<<1|y|y|
ADD         #c-1|r<<1 + (y<0 && x<0)|x<<1|y<<1| #*1と同じ形
#ここまでループする
PUSHW_1 -44 #c-1|r<<1 + (y<0 && x<0)|x<<1|y<<1|jump_delta| #ジャンプ先
PUSHB_1 5
CINDEX      #c-1|r<<1 + (y<0 && x<0)|x<<1|y<<1|jump_delta|c-1|
JROT        #c-1|r<<1 + (y<0 && x<0)|x<<1|y<<1| #jump if c-1 != 0
POP         #c|r|x|
POP         #c|r|
SWAP        #r|c|
POP         #r|
ENDF

bitwise OR

bitwise ANDの場合と同様だが、ANDの代わりにORを使えばbitwise ORになる。一命令入れ替えただけ。

#func1 - calculates bitwise OR
#initial stack: ...|x|y|
#final stack:   ...|x | b| (bitwise or)
PUSHB_1 1
FDEF        #x|y|
PUSHB_1 0   #x|y|r=0| #initialize return value r by 0
ROLL        #y|r|x|
PUSHB_1 32  #y|r|x|c=32| #initialize counter c by 32
ROLL        #y|x|c|r|
ROLL        #y|c|r|x|
PUSHB_1 4
MINDEX      #c|r|x|y|
#ここからループする(32回実行) #c|r|x|y|  #基準の形*1
# jump_to_here:
PUSHB_1 4
MINDEX      #r|x|y|c|
PUSHB_1 1   #r|x|y|c|1|
SUB         #r|x|y|c-1|
PUSHB_1 4
MINDEX      #x|y|c-1|r|
PUSHB_1 4
MINDEX      #y|c-1|r|x|
PUSHB_1 4
MINDEX      #c-1|r|x|y|
DUP         #c-1|r|x|y|y| 
PUSHB_1 0   #c-1|r|x|y|y|0|
LT          #c-1|r|x|y|y<0|
PUSHB_1 3
CINDEX      #c-1|r|x|y|y<0|x|
PUSHB_1 0   #c-1|r|x|y|y<0|x|0|
LT          #c-1|r|x|y|y<0|x<0|
OR          #c-1|r|x|y|y<0 || x<0|
PUSHB_1 4
MINDEX      #c-1|x|y|y<0 || x<0|r|
DUP         #c-1|x|y|y<0 || x<0|r|r|
ADD         #c-1|x|y|y<0 || x<0|r<<1|
ADD         #c-1|x|y|r<<1 + (y<0 || x<0)|
ROLL        #c-1|y|r<<1 + (y<0 || x<0)|x|
DUP         #c-1|y|r<<1 + (y<0 || x<0)|x|x|
ADD         #c-1|y|r<<1 + (y<0 || x<0)|x<<1|
ROLL        #c-1|r<<1 + (y<0 || x<0)|x<<1|y|
DUP         #c-1|r<<1 + (y<0 || x<0)|x<<1|y|y|
ADD         #c-1|r<<1 + (y<0 || x<0)|x<<1|y<<1| #*1と同じ形
#ここまでループする
PUSHW_1 -44 #c-1|r<<1 + (y<0 || x<0)|x<<1|y<<1|jump_delta| #ジャンプ先
PUSHB_1 5
CINDEX      #c-1|r<<1 + (y<0 || x<0)|x<<1|y<<1|jump_delta|c-1|
JROT        #c-1|r<<1 + (y<0 || x<0)|x<<1|y<<1| #jump if c-1 != 0
POP         #c|r|x|
POP         #c|r|
SWAP        #r|c|
POP         #r|
ENDF

bitwise XOR

 a,b\in \{0,1\}でXORしたければ、ANDやORの代わりに a-b != 0 をチェックすればよさそう。
なので、bitwise ANDの場合と同様であるが、ANDの代わりに a-b != 0 を計算する命令が入っている。

#func2 - calculates bitwise XOR
#initial stack: ...|x|y|
#final stack:   ...|x ^ b| (bitwise xor)
PUSHB_1 2
FDEF        #x|y|
PUSHB_1 0   #x|y|r=0| #initialize return value r by 0
ROLL        #y|r|x|
PUSHB_1 32  #y|r|x|c=32| #initialize counter c by 32
ROLL        #y|x|c|r|
ROLL        #y|c|r|x|
PUSHB_1 4
MINDEX      #c|r|x|y|
#ここからループする(32回実行) #c|r|x|y|  #基準の形*1
# jump_to_here:
PUSHB_1 4
MINDEX      #r|x|y|c|
PUSHB_1 1   #r|x|y|c|1|
SUB         #r|x|y|c-1|
PUSHB_1 4
MINDEX      #x|y|c-1|r|
PUSHB_1 4
MINDEX      #y|c-1|r|x|
PUSHB_1 4
MINDEX      #c-1|r|x|y| 
DUP         #c-1|r|x|y|y| 
PUSHB_1 0   #c-1|r|x|y|y|0|
LT          #c-1|r|x|y|y<0|
PUSHB_1 3
CINDEX      #c-1|r|x|y|y<0|x|
PUSHB_1 0   #c-1|r|x|y|y<0|x|0|
LT          #c-1|r|x|y|y<0|x<0|
SUB         #c-1|r|x|y|(y<0)-(x<0)|
PUSHB_1 0   #c-1|r|x|y|(y<0)-(x<0)|0|
NEQ         #c-1|r|x|y|(y<0)-(x<0)==0|
PUSHB_1 4
MINDEX      #c-1|x|y|(y<0)-(x<0)==0|r|
DUP         #c-1|x|y|(y<0)-(x<0)==0|r|r|
ADD         #c-1|x|y|(y<0)-(x<0)==0|r<<1|
ADD         #c-1|x|y|r<<1 + ((y<0)-(x<0)==0)|
ROLL        #c-1|y|r<<1 + ((y<0)-(x<0)==0)|x|
DUP         #c-1|y|r<<1 + ((y<0)-(x<0)==0)|x|x|
ADD         #c-1|y|r<<1 + ((y<0)-(x<0)==0)|x<<1|
ROLL        #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y|
DUP         #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y|y|
ADD         #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y<<1| #*1と同じ形
#ここまでループする
PUSHW_1 -47 #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y<<1|jump_delta| #ジャンプ先
PUSHB_1 5
CINDEX      #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y<<1|jump_delta|c-1|
JROT        #c-1|r<<1 + ((y<0)-(x<0)==0)|x<<1|y<<1| #jump if c-1 != 0
POP         #c|r|x|
POP         #c|r|
SWAP        #r|c|
POP         #r|
ENDF

感想

  • もっと効率化する方法はあるんだろうか。
  • あと、任意の数に対して正しい結果を返すかも自信ないし心配。
    • 検証どうすればいいのだろう。
    • 実装によって異なる動作をするみたいなのがあったら怖いなあ。

TrueType命令で擬似乱数: 線形合同法

TrueType命令で遊ぶシリーズ。

  1. 初めてのTrueType命令: Windowsでは見えないフォントをつくる - にせねこメモ
  2. フォントサイズに合わせて回転するフォントを作る(1) - にせねこメモ
  3. PPEM・ポイントサイズを表示するフォント - にせねこメモ
  4. TrueType命令で三角関数(sin, cos)を計算する - にせねこメモ
  5. フォントサイズに合わせて回転するフォントを作る(2) - にせねこメモ

擬似乱数の生成

TrueType命令で擬似乱数が計算できれば夢が広がるのではないかなーと思ってやってみた。

擬似乱数線形合同法を使うのが楽っぽい。実際に全然乱数ではないのだが、シミュレーションとかに使うのではないので、そこそこ乱数っぽく見えれば問題ないのではと思う。

線形合同法

線形合同法 - Wikipedia
Wikipediaに書いてあった、Park & Millerによる「最低基準」とされる次の式を計算する。

 \displaystyle
X_{n+1} = (48,271 \times X_n ) \mod (2^{31}-1)

http://c-faq.com/lib/rand.htmlにこれを実装したCコードがあり、次のようになっている。(16807→48271に変更してある)

#define a 48271
#define m 2147483647
#define q (m / a)
#define r (m % a)

static long int seed = 1;

long int PMrand()
{
	long int hi = seed / q;
	long int lo = seed % q;
	long int test = a * lo - r * hi;
	if(test > 0)
		seed = test;
	else	seed = test + m;
	return seed;
}

ここで、

2^{31} = 44,488 \times 48,271 + 3,399
なので、

hi = seed / 44488;
lo = seed % 44488;
test = 48271 * lo - 3399 * hi;
if(test>0) seed = test;
else       seed = test+2147483647;

みたいな計算を実際にすればいいことがわかる。

実装

seedは正数とする。
seedをスタックトップにおいて呼び出すとseedから計算される次の(擬似)乱数が返ってくる関数0は次のようになった:

# func0: calculate random number (linear congruential generator)
# initial stack: ...|seed(>0)| 
# final stack:   ...|pseudo_random_number|

PUSHB_1 # func0
 0
FDEF     #seed|
DUP      #seed|seed|
PUSHW_2
 30000
 14488
ADD      #seed|seed|44488|
DUP      #seed|seed|44488|44488|
ROLL     #seed|44488|44488|seed|
SWAP     #seed|44488|seed|44488|
PUSHW_1  #seed|44488|seed|44488|64.0|
 4096
MUL      #seed|44488|seed|44488*64|
DIV      #seed|44488|seed/44488|
ROLL     #44488|seed/44488|seed|
ROLL     #seed/44488|seed|44488|
PUSHB_1
 3
CINDEX   #seed/44488|seed|44488|seed/44488|
PUSHW_1  #seed/44488|seed|44488|seed/44488|64.0|
 4096
MUL      #seed/44488|seed|44488|(seed/44488)*64|
MUL      #seed/44488|seed|44488*(seed/44488)|
SUB      #seed/44488|seed%44488| #seed%44488 = seed - 44488*(seed/44488)
PUSHW_3  #seed/44488|seed%44488|64.0|30000|18271|
 4096
 30000
 18271
ADD      #seed/44488|seed%44488|64.0|48271|
MUL      #seed/44488|seed%44488|48271*64|
MUL      #seed/44488|(seed%44488)*48271|
SWAP     #(seed%44488)*48271|seed/44488|
PUSHW_2  #(seed%44488)*48271|seed/44488|3399|64.0|
 3399
 4096
MUL      #(seed%44488)*48271|seed/44488|3399*64|
MUL      #(seed%44488)*48271|(seed/44488)*3399|
SUB      #test| # test = (seed%44488)*48271-(seed/44488)*3399
DUP      #test|test|
PUSHB_1
 0       #test|test|0|
GT       #test|test>0|
IF       #test|  #if test>0
  # do nothing here
ELSE
 PUSHW_7 #test|0xFF|256.0|0xFF|256.0|0xFF|256.0|0x7F|
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  127   # 0x7F
 MUL     #test|0xFF|256.0|0xFF|256.0|0xFF|0x7F00|
 ADD     #test|0xFF|256.0|0xFF|256.0|0x7FFF|
 MUL     #test|0xFF|256.0|0xFF|0x7FFF00|
 ADD     #test|0xFF|256.0|0x7FFFFF|
 MUL     #test|0xFF|0x7FFFFF00|
 ADD     #test|0x7FFFFFFF| #2147483647 = 0x7FFFFFFF
 ADD     #test+2147483647|
EIF
ENDF

これを、例えば

PUSHB_1
 1
0 CALL

などとしてseed=1として呼び出すとスタックトップに48271が積まれる。

改良

実際に乱数として使うのであれば、seedはスタックに積まずにstorage areaに格納することにするのが良いのではないかと思う。
prepで初期seedをstorage areaに書きこみ、乱数関数ではstorage areaからseedを読み込み、次の乱数を計算して値を返す(この乱数は次のseedとなるためstorage areaに書きこまれる)。
ここでは、storage areaの0番をseedの記録に使用することにする。
このような関数1を先ほどの関数を弄って作成した。

# func0: calculate random number (linear congruential generator)
# initial stack: ...| 
# final stack:   ...|pseudo_random_number|

PUSHB_1 # func1
 1
FDEF     #|
PUSHB_1
 0   # storage area number
RS       #seed|  # seed = StorageArea[0]
DUP      #seed|seed|
PUSHW_2
 30000
 14488
ADD      #seed|seed|44488|
DUP      #seed|seed|44488|44488|
ROLL     #seed|44488|44488|seed|
SWAP     #seed|44488|seed|44488|
PUSHW_1  #seed|44488|seed|44488|64.0|
 4096
MUL      #seed|44488|seed|44488*64|
DIV      #seed|44488|seed/44488|
ROLL     #44488|seed/44488|seed|
ROLL     #seed/44488|seed|44488|
PUSHB_1
 3
CINDEX   #seed/44488|seed|44488|seed/44488|
PUSHW_1  #seed/44488|seed|44488|seed/44488|64.0|
 4096
MUL      #seed/44488|seed|44488|(seed/44488)*64|
MUL      #seed/44488|seed|44488*(seed/44488)|
SUB      #seed/44488|seed%44488| #seed%44488 = seed - 44488*(seed/44488)
PUSHW_3  #seed/44488|seed%44488|64.0|30000|18271|
 4096
 30000
 18271
ADD      #seed/44488|seed%44488|64.0|48271|
MUL      #seed/44488|seed%44488|48271*64|
MUL      #seed/44488|(seed%44488)*48271|
SWAP     #(seed%44488)*48271|seed/44488|
PUSHW_2  #(seed%44488)*48271|seed/44488|3399|64.0|
 3399
 4096
MUL      #(seed%44488)*48271|seed/44488|3399*64|
MUL      #(seed%44488)*48271|(seed/44488)*3399|
SUB      #test| # test = (seed%44488)*48271-(seed/44488)*3399
DUP      #test|test|
PUSHB_1
 0       #test|test|0|
GT       #test|test>0|
IF       #test|  #if test>0
  #do nothing here
ELSE
 PUSHW_7 #test|0xFF|256.0|0xFF|256.0|0xFF|256.0|0x7F|
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  127   # 0x7F
 MUL     #test|0xFF|256.0|0xFF|256.0|0xFF|0x7F00|
 ADD     #test|0xFF|256.0|0xFF|256.0|0x7FFF|
 MUL     #test|0xFF|256.0|0xFF|0x7FFF00|
 ADD     #test|0xFF|256.0|0x7FFFFF|
 MUL     #test|0xFF|0x7FFFFF00|
 ADD     #test|0x7FFFFFFF| #2147483647 = 0x7FFFFFFF
 ADD     #test+2147483647|
EIF
DUP      #seed|seed|
PUSHB_1  #seed|seed|0|
 0
SWAP     #seed|0|seed|
WS
ENDF

これを'prep'で、適当な値で初期化する。

PUSHB_1
 0
PUSHB_1
 1  #初期seed (>0)
WS      # StorageArea[0] = 1

すると、関数1を呼び出すと、スタックに乱数が積まれていくようになるはずである。
試しに、グリフから

10 1 LOOPCALL

と呼び出してみると、スタックに乱数が積まれていく。
f:id:nixeneko:20170528031648p:plain
右側中段のスタック表示ウィンドウを見ると、上で記載したC言語のコードの出力と同じになっていることが確かめられる。

実際に使う段には、最大値が大きいので、適当な数で割ってスケールして使う必要があると思う。

コード(コピペ用)

'fpgm'
PUSHB_1
 1
FDEF
PUSHB_1
 0
RS
DUP
PUSHW_2
 30000
 14488
ADD
DUP
ROLL
SWAP
PUSHW_1
 4096
MUL
DIV
ROLL
ROLL
PUSHB_1
 3
CINDEX
PUSHW_1
 4096
MUL
MUL
SUB
PUSHW_3
 4096
 30000
 18271
ADD
MUL
MUL
SWAP
PUSHW_2
 3399
 4096
MUL
MUL
SUB
DUP
PUSHB_1
 0
GT
IF
ELSE
PUSHW_7
 255
 16384
 255
 16384
 255
 16384
 127
MUL
ADD
MUL
ADD
MUL
ADD
ADD
EIF
DUP
PUSHB_1
 0
SWAP
WS
ENDF
'prep'
PUSHB_1
 0
PUSHB_1
 1
WS