Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい!

DeepLearning(3): そして逆伝播(でも全結合層まで)

第1、2回で順伝播の処理は実現した。今回やっと「学習」処理に入る。 (それにしても時間が取れず、またしてもだいぶあいだが空いてしまった。記憶が飛んでいる)

前回作成したプログラムでは、多段のニューラルネットワークを構成したが、 その終盤には2つの全結合層(隠れ層と出力層)が含まれている。本稿では、 逆伝播による学習で全結合層のフィルタ(ウェイト)を更新する部分について 解説し、学習によりフィルタの分類能力が向上することを確認したい。

今回のソースはこちら

1. 全結合層の逆伝播処理

1-1. 処理の流れ

逆伝播の理論や数式の詳細は各種書籍やWebサイトを見ていただくとして、 コード実装の説明に必要な事項をざっと説明しよう。 畳み込み層の逆伝播はちょっと複雑っぽいので次回にし、今回は全結合層だけに しておく。

損失関数

「学習」とは、例えば画像認識では入力画像から「犬」や「猫」を正しく判断 できるように各層のフィルタを修正してくことだ。 (だから本プログラムでは全結合層と畳み込み層がその対象なのだが、 今回は全結合層だけ。)

修正するにはまず「どれぐらい正解とずれているか」を確認しないといけない。 これを損失関数 $E$ で測る。今回は出力層の活性化関数にsoftmaxを使っており、 出力値と正解値の誤差を簡単にするため $E$ には交差エントロピー関数を 採用した。

{ \displaystyle
E = -\sum_k t_k \ln y_k \qquad (1)
}

ここで、$t$ は正解値(ベクトル)、 $y$ は出力値(ベクトル)だ。例えば 本プログラムでは画像を3種類に分類するのでベクトルの要素数は3、 仮に教師データを分類1とするなら、

$t = [1.0, 0.0, 0.0], y = [0.7, 0.1, 0.2]$

といった感じだ。

フィルタの更新

さて損失関数で表される「ずれ」を小さくするため、フィルタ $\boldsymbol{W}$ を更新する。大雑把には次の通り。

{ \displaystyle
 \boldsymbol{W} \leftarrow \boldsymbol{W} - \eta \Delta \boldsymbol{W} \qquad (2)
}

ここで、 \eta は学習率(learning rate)と呼ばれるパラメータだ。 さて問題は $\Delta \boldsymbol{W}$ をどう求めるかだが、損失関数 $E$ の勾配で表す ことができるようだ。

全結合層のフィルタ $\boldsymbol{W}$ は入力数 $\times$ 出力数の行列と見ることができる。 その個々の要素を $w_{ij}$ とすると、

 { \displaystyle
w_{ij} \leftarrow w_{ij} - \eta \frac{\partial E}{\partial w_{ij}} \qquad (3)
}

と書ける(らしい)。だから、最後の偏微分を解けばめでたしめでたしとなる(が、 そこが大変)。直接この偏微分を各層で解くと本当に大変だそうだが、 実は一つ前の層(出力に近い層)の計算結果をうまく使うとこれが簡単にできる、 というのが「誤差逆伝播法」のミソらしい。つまり出力層側から 誤差を計算して次の層に渡す。それを元に次の誤差を計算して さらに次の層に渡す。これを繰り返す。

出力層(Lとする)では、上記の偏微分は以下のようになる。

{ \displaystyle
\frac{\partial E}{\partial w_{ij}^{(L)}} = \frac{\partial E}{\partial y_i^{(L)}}
 \frac{\partial y_i^{(L)}}{\partial a_i^{(L)}}
 \frac{\partial a_i^{(L)}}{\partial w_{ij}^{(L)}}
= \delta_i y_j^{(L-1)} \qquad (4)
}

{ \displaystyle
\left( \delta_i = \frac{\partial E}{\partial y_i^{(L)}}
 \frac{\partial y_i^{(L)}}{\partial a_i^{(L)}} \right)
}

この $\delta_i$ が次の層へ伝える「誤差」、これに一つ手前の層の結果 $y_j$ を掛ける。

順伝播/逆伝播の流れと数式

これらを踏まえ、各層の順伝播・逆伝播を表にしてみた (逆伝播は、今回は全結合層まで)。 右上のカッコ付き数字は層の番号。 太字はベクトル/行列で、 $\boldsymbol{W}_{[3,2]}$ は 3x2行列、のように右下に要素数を記載した。 なお、 $\boldsymbol{c}$ は教師データである。 (層番号と要素数をつけたら数式が非常に醜くなってしまった。。。)

No. 順伝播 逆伝播 勾配($\Delta \boldsymbol{W}$)
0 入力 $\boldsymbol{y}_{[12,12,1]}^{(0)}$
1 畳み込み層1 $\boldsymbol{a}_{[10,10,10]}^{(1)}, a_{cij}^{(1)} = \sum_c \sum_s \sum_t w_{st}^{(1)} y_{(s+i)(t+j)}^{(0)}$
2 活性化層1 $\boldsymbol{y}_{[10,10,10]}^{(2)} = \max (\boldsymbol{a}^{(1)}, 0)$
3 プーリング層1 $\boldsymbol{y}_{[5,5,10]}^{(3)} = \max (y_{(li+s)(lj+t)}^{(2)}), \; {\scriptsize \mathrm{where} : s,t\in : \mid 0,l\mid}$
4 畳み込み層2 $\boldsymbol{a}_{[4,4,20]}^{(4)}, a_{cij}^{(4)} = \sum_c \sum_s \sum_t w_{st}^{(4)} y_{(s+i)(t+j)}^{(3)}$
5 活性化層2 $\boldsymbol{y}_{[4,4,20]}^{(5)} = \max (\boldsymbol{a}^{(4)}, 0)$
6 プーリング層2 $\boldsymbol{y}_{[2,2,20]}^{(6)} = \max (y_{(li+s)(lj+t)}^{(5)}), \; {\scriptsize \mathrm{where} : s,t\in : \mid 0,l\mid}$
7 平坦化層 $\boldsymbol{y}_{[80]}^{(7)} \leftarrow \boldsymbol{y}^{(6)}$ $\boldsymbol{\delta}_{[80]}^{(7)} = {}^t\boldsymbol{W}_{[20,80]}^{(8)} \cdot \boldsymbol{\delta}^{(8)}$ -
8 全結合層1 $\boldsymbol{a}_{[20]}^{(8)} = \boldsymbol{W}_{[80,20]}^{(8)} \cdot \boldsymbol{y}^{(7)} + \boldsymbol{b}_{[20]}^{(8)}$ $\boldsymbol{\delta}_{[20]}^{(8)} = f'(\boldsymbol{a}^{(8)}) \odot \boldsymbol{\delta}^{(9)}, \; {\scriptsize where : f = \max}$ $\Delta w_{kl}^{(8)} = \delta_k^{(8)} \cdot y_l^{(7)},\; b_k^{(8)} = \delta_k^{(8)}$
9 活性化層3 $\boldsymbol{y}_{[20]}^{(9)} = \max (\boldsymbol{a}^{(8)}, 0)$ $\boldsymbol{\delta}_{[20]}^{(9)} = {}^t\boldsymbol{W}_{[3,20]}^{(10)} \cdot \boldsymbol{\delta}^{(10)}$ -
10 全結合層2 $\boldsymbol{a}_{[3]}^{(10)} = \boldsymbol{W}_{[20,3]}^{(10)} \cdot \boldsymbol{y}^{(9)} + \boldsymbol{b}_{[3]}^{(10)}$ $\boldsymbol{\delta}_{[3]}^{(10)} = \boldsymbol{y}^{(11)} - \boldsymbol{c}_{[3]}$ $\Delta w_{ij}^{(10)} = \delta_i^{(10)} \cdot y_j^{(9)},\; b_i^{(10)} = \delta_i^{(10)}$
11 活性化層4(出力) $\boldsymbol{y}_{[3]}^{(11)} = \mathrm{softmax} (\boldsymbol{a}^{(10)})$ - -

※ なぜNo.10:逆伝播で $\boldsymbol{\delta} = \boldsymbol{y} - \boldsymbol{c}$ となるのかはそっち系のサイト(こことか)や書籍を読んでほしい。

誤差 $\delta$ が各層で計算され、次の層へ伝わっていくイメージが伝わるだろうか? これをすべての層に対して遡っていけば、めでたく一個の教師データに対する 学習が完了するわけだ。

実際には複数個(N)の教師データについて $\Delta \boldsymbol{W}$ を求め、それを平均して $\boldsymbol{W}$ の更新に使う。数式(2)は実際は次のような感じになる。 $\boldsymbol{W}$ の右下の $i$ はi番目の教師データから得られたものという意味。

{ \displaystyle
\boldsymbol{W} \leftarrow \boldsymbol{W} - \eta \frac{1}{N} \sum_{i=1}^{N} \Delta \boldsymbol{W}_i \qquad (5)
}

1-2. 誤差逆伝播の実装

面倒な理論・数式の説明はこれぐらいにして、どう実装したか説明しよう。

処理の全体像

下記のmain関数内で呼ばれているtrainLoopが中心部だ。

(Main-train.hs)

main = do
     :
  layers' <- trainLoop getTeachers sampleE putF (learnR st) (layers st) is
     :

学習を繰り返すためにepoch番号を入れた配列isを渡していて、学習済みの状態 (層リスト)を次に引き渡すためちょっと面倒だが再帰になっている。

(Main-train.hs)

trainLoop :: (Int -> IO [Trainer]) -> [Trainer]
  -> (Int -> [Layer] -> [Trainer] -> IO ()) -> Double -> [Layer] -> [Int]
  -> IO [Layer]
trainLoop _ _ _ _ ls []             = return ls
trainLoop getT se putF lr ls (i:is) = do
  ls' <- trainLoop getT se putF lr ls is
  teachers <- getT i
  let
    rls = tail $ map reverseLayer $ reverse ls'    -- fist element isn't used
    dls = map (train ls' rls) teachers
    ls'' = update lr (transpose dls) ls'           -- dls = diff of layers
  putF i ls'' se
  return ls''

trainLoop関数の最初の行でisを一つ減らしてtrainLoopを呼んでいるところ。 このため、isの要素がなくなるところまでtrainLoopの呼び出しが続き、 最後に「学習前の層リスト」が返され、以後少しずつ層リストを学習で修正する。 最後に、学習が完了した層リストがmain中のlayers'に入るという寸法だ。

trainLoopの中では、上記の通り一つ前のepoch番号の学習結果を得てから 今回のepochで使う教師データを生成している。yusugomori実装では教師データ数が 150個に対しバッチ数が150のため毎回全件で学習しているのだ。

次のrlsはちょっとややこしい。学習では誤差逆伝播法を使うが、 そのため層リストを逆順に辿りたい。ただ単純な順序の反転ではダメなので、 まず層リストを反転した上で各層を「逆変換」している(reverseLayer)。 全結合層ではフィルタ行列の転置が逆変換に相当する。

次の行に出てくるtrainがメイン処理だ。これは後述する。 この行で教師データ全部に対し勾配($\Delta \boldsymbol{W}_i$)を計算し、 次の行で層リストを更新して終わりだ。なおputFは学習結果を使って テストデータの認識精度を計算し、表示する処理だ。5 epoch毎に結果を出力する。

次にtrainについて。

(Trainer.hs)

train :: [Layer] -> [Layer] -> Trainer -> [Layer]
train [] _ (i, c) = []
train ls rls (i, c) = dls
  where
    (y, op') = splitAt 1 $ forwardProp ls [i]
    d = (head $ head $ head y) `vsub` c
    (_, dls) = backwardProp (zip (tail op') rls) (d, [])

where句の次の行は前回までに作成した順伝播処理だ。得られた 結果リスト(各層の処理結果)を最後の出力=yとそれ以外に分離している。 次の行、yがImage型のため余計な次元を除去して一次元ベクトルにした後で $\boldsymbol{y} - \boldsymbol{c}$ を計算している (にしても3個もheadがあるのは格好悪い)。 これが最初の誤差 $\boldsymbol{\delta}^{(10)}$ だ。

3行目が逆伝播の本体。一つ目の引数zip (tail op') rlsは 各層の出力値と逆変換用の各層を層毎に組みにしたもののリストだ。 逆伝播はこのリストに対して繰り返す。 二つ目の引数で前層で計算した誤差と各層の勾配、それを更新して次へ渡す。

(Trainer.hs)

backwardProp :: [(Image, Layer)] -> (Delta, [Layer]) -> (Delta, [Layer])
backwardProp [] (_, ls) = ([], ls)
backwardProp ((im,l):ols) (d, ls) = backwardProp ols (d', l':ls)
  where
    (d', l') = backwardLayer l im d

ここまでが、先述の表と図で説明したものを処理している部分である。 うーん、今回も話が長くなりそうだ。

誤差と勾配の計算

逆伝播処理は層の種類で異なるため、backwardLayer関数のパターンマッチで 仕分けている。今のところ活性化層(ActLayer)と全結合層(FullConnLayer)のみ 実装している。他は次回になんとかしようと思う。

(Layer.hs)

backwardLayer :: Layer -> Image -> Delta -> (Delta, Layer)
backwardLayer (NopLayer)           _  d = (d, NopLayer)
backwardLayer (ActLayer f)         im d = deactivate f    im d
backwardLayer (MaxPoolLayer s)     im d = depoolMax  s    im d
backwardLayer (ConvLayer s fs)     im d = deconvolve s fs im d
backwardLayer (FullConnLayer fs)   im d = deconnect    fs im d
backwardLayer l@(FlattenLayer x y) im d = (head $ head $ unflatten x y [[d]], l)

中身を示そう。まず簡単な活性化層の方から。

(ActLayer.hs)

deactivate :: ActFunc -> Image -> Delta -> (Delta, Layer)
deactivate f im delta
  | c == [0.0]  = (zipWith (*) delta f', ActLayer relu')
  | c == [1.0]  = ([], ActLayer f)
  | otherwise = ([], ActLayer f)
  where
    c = f [0.0]
    f' = relu' (head $ head im)

出力に近い方、全結合層の後の活性化層のうち最後の層=出力層では softmax関数を使っているが、誤差を計算する必要がないのでほったらかしている。 一方その二つ前の活性化層3(No.9)はReLUを採用しているので、逆関数 $f'$ は ステップ関数(relu'として別に定義)になり、前層からきた誤差と 掛け合わせている(zipWith (*) delta f'の部分)。

ところで、活性化関数の種類によって処理を変えたいのだが、普通の引数の 型(DoubleとかStringとか)や値でのパターンマッチが使えないみたいだ。 引数fの値(reluとかsoftmaxとか)で分けられたら楽なんだが。。。 それができないので、関数fにダミー値([0.0])を与えた結果で場合分け している。reluなら0.0に、softmaxなら1.0になるからだ(たぶん)。

次は全結合層の説明。

(FullConnLayer.hs)

deconnect :: [FilterF] -> Image -> Delta -> (Delta, Layer)
deconnect fs im delta = (mmul delta fs', FullConnLayer $ calcDiff delta im')
  where
    fs' = tail fs
    im' = head $ head im

calcDiff :: Delta -> [Double] -> [FilterF]
calcDiff delta im = map (mulImage im') delta
  where
    im' = 1.0:im
    mulImage :: [Double] -> Double -> [Double]
    mulImage im d = map (*d) im'

誤差は、フィルタfsと前層からの誤差deltaを掛けている(mmul delta fs')。 fsはすでに転置行列にしてあるが、その最初の要素はバイアスbに相当するため 削る(tail fs)。全結合層は学習によりフィルタを更新するから、誤差から 学習に必要な勾配 $\Delta W$を計算している。calcDiffがそれ。 誤差値deltaと順伝播での入力値imのそれぞれの要素を掛け合わせて 作り出している。本実装ではバイアスbはフィルタに統合されているので、 入力値に1.0を加えているが(im' = 1.0:im)。

全結合層のフィルタ更新

さあ最後に、各教師データから得られた勾配を平均してフィルタを更新する 部分を見てみよう。

(Trainer.hs)

update :: Double -> [[Layer]] -> [Layer] -> [Layer]
update lr _ [] = []
update lr [] ls = ls
update lr (dl:dls) (l:ls) = (updateLayer lr l dl):(update lr dls ls)

dlは処理対象の層lに対する各教師データからの勾配が入っている。 これをupdateLayerに渡して計算している。lrは学習率だ。

次が、全結合層での勾配を集計する処理である。

(FullConnLayer.hs)

updateFullConnFilter :: [FilterF] -> Double -> [Layer] -> Layer
updateFullConnFilter fs lr dl = FullConnLayer fs'
  where
    ms = strip dl
    delta = mscale (lr / (fromIntegral $ length ms))  $ msum ms
    fs' = msub fs delta

ほぼ数式(5)そのままなので、ほぼ説明不要だろう。ms = strip dlに ついてだけ。各教師データから得た勾配のリストから、この層での対象外の データがあれば除外するようにしている。普通は問題ないはずだが。

2. 「学習」の評価

プログラムも出来上がったことだし、どれだけ学習するか評価してみよう。 コンパイルと実行は次の通り。

% cabal build
   :
% ./dist/build/train/train
Initializing...
Training the model...
iter =     0/500 accuracy = 0.3333333333 time = 0.149965s
iter =     5/500 accuracy = 0.3333413554 time = 5.132365s
   :

途中経過を 5 epoch 毎に画面に出している。前回から少しフォーマットを 変更し、正答率を「認識精度=accuracy」に改めている。

2-1. 学習結果(ノーマル版)

上記のプログラムを10回実行し、それぞれ認識精度と実行時間を計測した。 結果は次のグラフの通り。

f:id:eijian:20161108015850p:plain

点線が各試行時の実測値、太い実線が10回の平均値だ。これを見ると、認識精度は 各試行でかなりばらついているのがわかる。また、認識精度の向上が 急峻なものと緩慢なものがあり、かなり差がある。

プログラムでは、学習用データ、テスト用データ、および各フィルタに乱数を 使っているが、それらの初期状態によりこれだけの差が生じてしまっていると 思われる。(学習効果が薄い=うまく特徴が出ない画像群が生成されたとか、 フィルタの初期状態がイケてないとか)

ただ、総じて 500 epoch ぐらい繰り返せば認識精度がほぼ9割を超えることも 分かった。 まだ畳み込み層の学習をしていない段階でここまで精度が上がるとは思っておらず、 嬉しい誤算だった。(こんなものなのかな?)

2-2. 学習結果(学習率の差)

先のグラフで、認識精度の向上は最初からではなく50 epochを超えたあたりから 急に上がりだしている。私は素人なのでよくわからないが、いろいろな情報を見ると 「学習率」はどう設定するかとか、学習進行に応じて変化させるとか書かれて いて、重要なパラメータであるらしい。

このプログラムでは学習率を0.1に固定しているが、初期段階ではこれが小さすぎる のかもしれない。ということで、学習率を0.2と0.05に変えて試行した結果も出して みた。(0.2と0.05の試行は5回、その平均をとった)

f:id:eijian:20161108015855p:plain

予想通りかというとちょっと微妙だが、学習率を大きくした方がより速く認識精度が 向上していることがわかる。ただ、初っ端から急上昇しているわけではないので もっと大きくしないといけないのだろう。ただ、大きいままで学習を続けると 収束しないという話もある(?)ので、やはり動的に変化させるのがいいのだろう。 あとでそれにも取り組みたい。

今回のプログラムはyusugomori実装を「写経」することが第一課題なので 学習率は0.1に戻しておこう。

2-3. yusugomori実装との比較

ここまで、学習が進んだ!と喜んできたわけだが、「写経」元のyusugomori実装と 同じことができているのだろうか?

ということでyusugomori実装との比較をしてみたのが次のグラフだ。

f:id:eijian:20161108015900p:plain

avg-0.1が本プログラム(学習率0.1、試行10回の平均値、青線)、py-cnnが yusugomori実装(学習率0.1、試行1回、赤線)、py-cnnfcがyusugomori実装で 畳み込み層の学習を抑制した改造版(学習率0.1、試行1回、緑線)だ。 yusugomori実装では乱数のseedが固定されているため毎回同じ結果になるのと、 後で述べるようにめちゃくちゃ処理時間が掛かるのとで試行は1回だけにした。

先述の通り乱数の生成値により試行結果はかなりぶれるのだが、青と緑の線は 350 epoch近辺からほとんど同じような値になっていて、(楽観的に見れば) 本プログラムはyusugomori実装と同等な処理ができているのではないかと思う。 一方本来の学習(=畳み込み層も学習)であれば認識精度が急上昇し、かつ 200 epoch ぐらいでも認識精度が98%を超えてくる。それに比べ本プログラムでは 畳み込み層のフィルタがいつまでたっても初期状態(ランダム)のままなので、 この差は致し方ないところだろう。まあ、ここまで性能差があると早く 畳み込み層の学習部分も完成させたいと思うからモチベーションの向上になる。

次に処理時間について見てみる。200 epoch 実行時の結果は次表のとおりである。 (iMac early 2009, 3.0GHz Core 2 Duo, 4GB memory)

avg-0.1 py-cnn py-cnnfc
処理時間(s) 157 24838 6152
処理性能(epoch/s) 1.273 0.008 0.033
認識精度(%) 74.2 98.4 62.9

yusugomori実装の処理時間が半端ではないことがわかる。本プログラムも コンパイルされたネイティブコードであることを考えると決して十分な処理速度 ではないが、それでもまったく比べものになっていない。

NumPyライブラリのインストールに何か抜けている手順があるのだろうか? インストール時に何か特別なことをした記憶もないので、以下の記事にあるように 「ちゃんとした」インストールができていないのかもしれない。

(参考1) Kesin's diary: NumPy, SciPyをちゃんとインストールする

もしくは次の記事のように多重ループと配列要素への直接読み書きが 影響しているのかもしれない。

(参考2) numpyで明示的にループを書くと極端に遅くなる

いずれにせよ、当方は今の所Pythonを使う気はないので深入りはしない...

3. まとめ

今回は全結合層部分の逆伝播処理を実装した。これにより学習ができるように なったので、精度はともかく画像認識ができるようになった。ただ、これはまだ 半分。次回は畳み込み層の逆伝播を実装し、プログラムを完成させようと思う。

DeepLearning(2): まずは順伝播(下)

前回に続き、CNNの順伝播処理について。前回は下準備という感じだったが、 今回は各層の処理を考よう。

なお今回の記事に限らないが、私の記事は自分自身の忘備録であることが 主目的なため、各種技術・理論については正しくないかもしれないし プログラミングに必要なことに絞って書いている。よってCNNやDeep Learningの 正しいところは前回紹介したような、各種書籍・Webサイトにあたってほしい。

CNNの組み立て

まずはyusugomoriさんの Python実装 (yusugomori実装)を写経するので、そこで定義されているCNNを作りたいと思う。

CNNの構造

前回出した図をもう一度使うが、作成するCNNの構造は次のようになっている。

f:id:eijian:20160830205123p:plain

  • 第1畳み込み層
    • 入力画像はモノクロ12x12ドット(1チャネル)
    • 畳み込みフィルタは3x3ドット、10枚
    • 畳み込みにより10x10ドット、10チャネルの「画像」に変換される
  • 第1活性化層
    • 活性化関数はランプ関数(ReLU)、全ドットにReLUを適用
    • 画像サイズ、チャネル数に変化なし
  • 第1プーリング層
    • プーリングサイズは2x2ドット
    • 最大プーリングを適用
    • 出力画像の大きさは10÷2で5x5ドット、10チャネルとなる
  • 第2畳み込み層
    • 入力画像はモノクロ10x10ドット(10チャネル)
    • 畳み込みフィルタは2x2ドット、20枚
    • 畳み込みにより4x4ドット、20チャネルの「画像」に変換される
  • 第2活性化層
    • 活性化関数はReLU、全ドットにReLUを適用
    • 画像サイズ、チャネル数に変化なし
  • 第2プーリング層
    • プーリングサイズは2x2ドット
    • 最大プーリングを適用
    • 出力画像の大きさは4÷2で2x2ドット、20チャネルとなる
  • 一次元化
    • 2x2ドット、20チャネルを一次元化して2x2x20=80個のリストに変換
  • 隠れ層(全結合層)
    • 入力要素は80個、出力要素は20個
  • 第3活性化層
    • 全要素にReLUを適用
  • 出力層(全結合層)
    • 入力要素は20個、出力要素は3個
  • 第4活性化層
    • 全要素にSoftmaxを適用

各層は前回説明したようにmain関数内にLayerのリストとして 定義している。

トレーニング

各教師データ(画像)に対して上記のCNNを適用するのだが、一枚毎に 次のtrain関数で変換している。前回示したloop関数内で呼び出されている。 中身は各層の順伝播処理を呼んでいるだけなので簡単だ。

(Trainer.hs)

train :: [Layer] -> Trainer -> IO (Image, [Layer])
train [] (i, c) = return (i, [])
train ls (i, c) = do
  let op  = forwardProp ls [i]
      ls' = backwordProp ls op
  return (head op, ls')

forwardPropが順伝播、backwordPropが逆伝播だが、まだ順伝播のみの実装だ。 forwardProp再帰により層を順に適用し、それぞれの順伝播処理を呼び出す ようにした。

(Layer.hs)

forwardProp :: [Layer] -> [Image] -> [Image]
forwardProp [] is = is
forwardProp (l:ls) is = forwardProp ls (forward l is)

forward :: Layer -> [Image] -> [Image]
forward _ [] = []
forward NopLayer i = i
forward l@(ActLayer f) im@(i:is)     = (activate l i):im
forward l@(MaxPoolLayer s) im@(i:is) = (poolMax l i):im
forward l@(ConvLayer s fs) im@(i:is) = (convolve s fs i):im
forward l@(HiddenLayer fs) im@(i:is) = (connect fs i):im
forward l@(FlattenLayer f) im@(i:is) = (f i):im

各層の処理

何もしない層(NopLayer)

入力をそのまま返すダミーだ。(プログラムの動作を確認するため)

畳み込み層(ConvLayer)

この層はCNNのキモになるので自分の忘備録のために少し丁寧に記述しておく。 畳み込み層では、画像の各チャネルについてフィルタを適用する。 フィルタの大きさは3×3などの比較的小さいものだ。一般的でないかもしれないが 簡便のため正方形とする。

各チャネルの二次元データ上の各画素について、そこを起点にフィルタと 同じ大きさの正方形を取り出してフィルタの対応する画素と掛け合わせて 合計するのだ。下図では入力を$x$、出力を$u$、フィルタを$w$とする。

f:id:eijian:20160910202832p:plain

この処理、手続き型言語であれば多重ループで難なく書ける処理なのだが、 Haskellでどう書くか? 画像をリストでなくArrayで表現して、各画素を何度も順序ばらばらにアクセス すればすれば多重ループで書けるが、今回は再帰で処理するために次のようにした。

f:id:eijian:20160910202840p:plain

フィルタの大きさをsとすると、画像データの最初のs行をtakeで取り出して 各行の処理(convolveLine)を行う。そのあと、1行目を取り除き(tail)、 また同じことを繰り返すのだ。残りの行数がフィルタの大きさ(s)より小さく なったら終了する。

(ConvLayer.hs)

-- s: filter size

convolvePlain :: Int -> [Double] -> Plain -> Plain
convolvePlain s k ps
  | len < s   = []
  | otherwise = (convolveLine s k p:convolvePlain s k ps')
  where
    len = length ps
    p   = take s ps
    ps' = tail ps

各行についてtakeで行の最初からs個取り出すと$s \times s$個の 画素データが得られる。これにフィルタを適用する。 そのあと、各行の最初の要素を取り除き、それを 使って同じ処理を行の要素がs個未満になるまで繰り返す。 sum $ zipWith (*)の部分がフィルタを適用しているところである。

(ConvLayer.hs)

convolveLine :: Int -> [Double] -> [[Double]] -> [Double]
convolveLine s k ps
  | len < s   = []
  | otherwise = ((sum $ zipWith (*) vs k):convolveLine s k ps')
  where
    len = minimum $ map (length) ps
    vs = concat $ map (take s) ps
    ps' = map (tail) ps

この場合、入力画像のサイズに対し「フィルタのサイズ-1」だけ 小さな画像が出力される。今回の例では、元画像が12x12、フィルタが3x3のため、 12-(3-1)=10となり10x10の画像が得られる。

もう一つ、フィルタサイズより小さくなった時は黒や白のドットを 補うことで入力画像のサイズを維持する手法もあるようだ。 が、一般的な画像認識で使われるであろう数百×数百程度の画像では誤差の ようなものなので前者の手法を使う。

こうすれば、画像の各画素に対して再帰的にフィルタを適用できる。 入力画像の1チャネル毎にこの処理を繰り返すと、全体としては下図のように 複数チャネルの画像に変換されて出てくるわけだ。

f:id:eijian:20160910202845p:plain

その後、入力画像の各チャネルの処理結果を足し合わせる必要がある。 これは各チャネルの同じ位置にある画素データを足すだけなので簡単な ことのはずだが、頭がこんがらがって以外と手間取った。

(ConvLayer.hs)

sumPlains :: Double -> [Plain] -> Plain
sumPlains b [] = repeat $ repeat b  -- 2 Dimensions (X*Y)
sumPlains b (p:ps) = zipWith f p (sumPlains b ps)
  where
    f :: [Double] -> [Double] -> [Double]
    f [] _          = []
    f (a:[]) (b:_) = [a + b]
    f (a:as) (b:bs) = (a + b):(f as bs)

特に全チャネル分を足し合わせ、さらにバイアス値を加える必要があるので 多少複雑になってしまった。全画素に同じバイアス値を加えるために、 repeat $ repeat バイアス値により画像サイズがわからなくても 処理できた。普通ならちゃんと画像サイズ分の二次元データを用意する ところだろう。細かいところだが、遅延評価の恩恵か。

活性化層/活性化(ActLayer)

先述の通り、ニューラルネットについての書籍やサイトでは活性化層という 独立した層として説明されておらず、畳み込み層や隠れ層の後処理的な扱いだ。 が、今回は独立した層としたほうが「私にとって」わかり易そうだったので 分けることにした。 ただ、後でチューニングのために他の層に統合するかもしれない・・・。

活性化関数として、ランプ関数(Rectified Linear Unit, ReLU)とsoftmax関数を 用意した。それぞれ下記で表される。 これらはActLayer層の作成時にどちらか選ぶのだ。

(ReLU) $f(x_i) = \max(x_i, 0)$

(Softmax) $f(\boldsymbol{x}) = \frac{e^{x_i}}{\sum_{k=1}^K e^{x_k}}$

ここで、$x_i$は$i$番目の出力値、$K$は出力値の個数(次元)である。 実装はこう。

(ActLayer.hs)

relu :: ActFunc
relu as = map f as
  where
    f :: Double -> Double
    f a = if a > 0.0 then a else 0.0

softmax :: ActFunc
softmax as = map (\x -> x / sume) es
  where
    amax = maximum as
    es   = map (\x -> exp (x - amax)) as
    sume = sum es

記事を書いていて、ReLUのほうは素直にmaxを使えばいいと思った。 後で直そう。Softmaxの方も、ほぼ定義通りだ。ただし入力xに対し、 xの最大値を引いているがこれは発散しないための工夫だそうだ。 activate本体はすべての画素に対して活性化関数を適用する。 まあ、難しいものではない。

(ActLayer.hs)

activate :: Layer -> Image -> Image
activate (ActLayer f) im = map (map f) im

プーリング層(MaxPoolLayer)

本プログラムでは最大プーリングのみ対応する。他も必要になったらその時考えよう。 プーリングの処理では画像を小ブロックに分け、その中の各値から出力値を 計算する。最大プーリングでは、最大値を選択して出力するわけだ。 幾つかの画素から一つだけ値を出力するのだから、当然出力値の数は少なくなる。 小ブロックが2x2ドットだとすると、元の画像が縦横半分のサイズになるのだ。

f:id:eijian:20160910202851p:plain

これをプログラムにしてみよう。まず、1チャネルの画像の処理について。

(PoolLayer.hs)

poolMaxPlain :: Int -> Plain -> [[Pix]]
poolMaxPlain s p = map (poolMaxLine s) ls
  where
    ls = splitPlain s p

splitPlain :: Int -> Plain -> [[[Double]]]
splitPlain s [] = []
splitPlain s p  = (p':splitPlain s ps)
  where
    (p', ps)  = splitAt s p

poolMaxPlainが1チャネル分の処理の主体だ。入力データ(Plain)を splitPlainによってブロックサイズで縦方向に分割している。 縦sドットの横に細長いデータをそれぞれpoolMaxLineに引き渡して 処理させているのだ。

(PoolLayer.hs)

poolMaxLine :: Int -> [[Double]] -> [Pix]
poolMaxLine _ [] = []
poolMaxLine s ls
  | len == 0  = []
  | otherwise = (pixs:poolMaxLine s ls')
  where
    len  = length $ head ls
    pixs = max' $ zip (concat $ map (take s) ls) [1..]
    ls'  = map (drop s) ls

max' :: [Pix] -> Pix
max' [] = error "empty list!"
max' [x] = x
max' (x:xs) = maximum' x (max' xs)

maximum' :: Pix -> Pix -> Pix
maximum' a@(v1, i1) b@(v2, i2) = if v1 < v2 then b else a

ここは(今は)余計な処理が入っているためちょっと複雑だ。説明しよう。 大まかには、与えられた細長いデータからsドット分を取り出し、max'関数で 最大値を選択している。残りのデータは再帰により同じことを繰り返す。 ここで処理対象がDoubleではなくPixという型に変わっているが、 その定義はこう。

(PoolLayer.hs)

type Pix = (Double, Int)

これは小ブロックの各要素に番号をつけるために用意した型だ。 poolMaxLineでは、$s \times s$ドット分取り出したら一次元化して それに番号を振っている。以下の部分。

    zip (concat $ map (take s) ls) [1..]

これは最大プーリングでどの画素が最大だったかを覚えておいて、 後々逆伝播処理で伝播先(?)とするために必要だとか (とどこかで読んだはず、だがまだ勉強中で本当かどうかわかっていない。 いつのまにかしれっと消えているかもしれない。。。)。

二次元データについてpoolMaxPlainで処理できれば、それを画像データ全部に 適用すればよいだけだ(poolMaxpoolMax')。

一次元化(FlattenLayer)

これは層というよりデータの整形だ。畳み込み層やプーリング層の 出力は二次元平面、複数チャネルという、いわゆる「画像」の形をしているが、 隠れ層では「一次元配列」を入力に取る。ということでその間に変換が必要なのだ。

f:id:eijian:20160910202856p:plain

(Image.hs)

flatten :: Image -> Image
flatten im = [[concat $ concat im]]

単に全チャネルの画素値を一列に並べているだけである。 ただし、Layer型なので結果は「画像=3次元配列」にする必要がある。 なのでわざわざ前後に2つの角カッコをつけて3元配列にみせかけている。

隠れ層(HiddenLayer)

これはいわゆる全結合層なのだが、他の記事やプログラムを見ると隠れ層と いう名で作られているのでそれに倣った。隠れ層としながら、最終の出力層も このHiddenLayerを使っているので名前は変えるべきかもしれない。。。

この層では、入力(channel)c個を出力(kernel)k個に変換する。 c個の要素全てがk個の各出力に関連付けられているため全結合というらしい。 その関連付けをフィルタとして与える必要があるのでそのサイズは$c \times k$に なる。初期はランダム値を設定し、学習によってその値を少しずつ補正することで 画像認識などができるらしい。フィルタは$c \times k$の二次元リストで 与えることにしよう。

(HiddenLayer.hs)

initFilterH :: Int -> Int -> IO [FilterH]
initFilterH k c = do
  f <- forM [1..k] $ \i -> do
    f' <- initKernel c
    return f'
  return f
  where
    a = 1.0 / fromIntegral c
    initKernel :: Int -> IO FilterH
    initKernel c = do
      w <- forM [1..c] $ \i -> do
        r <- MT.randomIO :: IO Double
        return ((r * 2.0 - 1.0) * a)
      return (0.0:w)

各要素は入力要素数cに対し$1/c$ ~ $-1/c$ の間の乱数としている。 ただし、最初の要素はバイアス値で、初期値は0.0とするらしい。

この層では入力値にフィルタ値を掛けてすべて足しこむ(k個分)。 式で書くとこうなる。なおフィルタの要素は $w_{ij}$ としている。

[u_i = \sum{j=0}^c w{ij} x_j + b_i]

ここで$x_j$は$j$番目の入力値、$u_i$は$i$番目の出力値、$b_i$は$i$番目の バイアス値である。

実装はこう。

(HiddenLayer.hs)

connect :: [FilterH] -> Image -> Image
connect [] _ = error "invalid FilterH"
connect _ [] = error "invalid Image"
connect fs [[im]] = [[map (dot im) fs]]
  where
    dot :: [Double] -> [Double] -> Double
    dot im f = sum $ zipWith (*) (1.0:im) f

dot関数のzipWith第一引数に1.0が追加されているのは先の$b_i$を 足し込むためである(FilterHには$b_i$が含まれている)。

実行してみた

以上でひとまず順伝播については実装できた(と思う)。まずはコンパイル

$ cd deeplearning
$ cabal configure
  :
$ cabal build
  :

では実行してみよう。

$ dist/build/deeplearning/deeplearning
Building the model...
Training the model...
iter = 0/200 time = 0s ratio = 0.0
iter = 5/200 time = 0.0030001s ratio = 0.33333713817332733
iter = 10/200 time = 0.333019s ratio = 0.33333713817332733
iter = 15/200 time = 0.682039s ratio = 0.33333713817332733
     :
iter = 200/200 time = 13.0807482s ratio = 0.33333713817332733
Finished!

反復5回毎に状態を出力しており、左から

  • epoch(反復回数)
  • 経過時間(秒)
  • 正解率
    • 正解率は、全テストデータを順伝播にかけ得られた出力値から正解にあたる 要素の値(その要素が正解である確率)を取り出して平均したもの。
    • 3種類に分類されるデータを使ったので何も学習していない状況だと 出力は均等割合となり約1/3になる。

である。 正解率は期待通り約1/3になっているものの、学習してないのでこの結果だけで はなんとも言えない。ただ、何となくそれらしい出力が出てくると楽しい。

まとめ

前回から取り組んでいたCNN画像認識プログラムの「順伝播」のみ実装完了した。 出力結果にはほとんど意味はないが、各層の処理を経てなんとなく結果が出力 されたので雰囲気はつかめたかなと。

次回以降で本当のメインである学習処理を実装し「意味のある結果」が得られる ことを期待しよう。

(ここまでのソースはこちら

DeepLearning(1): まずは順伝播(上)

さて、今回からは新たなネタとしてDeepLearningに取り組もう。 DeepLearningといえばすでにブームは去って、今は現実の課題に どんどん活用されている状況だ。使いやすく性能のいいツールや ライブラリがいろいろあり、ちょっと学べば簡単にその恩恵に 与ることができる(かもしれない)。

しかし、仕組みがわからないのをただ使うのは面白くないし 個人的に写真の自動分類をしたいというニーズもあるし、遅まきながら 自分で作ってみよう、というわけだ。

目標と参考資料

最終的には「一般物体認識」まで行きたい(というかそうでないと 写真の分類にならない)が、目標が高すぎると頓挫するのでここでは 以下の2つのステップまでとする。

  • ステップ1: yusugomoriさんのpython実装を写経する
  • ステップ2: 上を使ってお決まりのMNISTデータを使った性能検証をする

「自分で作ってみよう」と大口を叩いていながら他の方の実装を写経する とは詐欺的だが、そこは当方素人なので。。。以後、yusugomoriさんのpython実装を 「yusugomori実装」と表記する。

作るのは、畳み込みニューラルネットワーク(Convolutional Neural Network, CNN)による画像認識プログラムとする。以下、参考にした書籍や Webサイトを列挙。

ひとつの本やWebサイトだけでは理解が進まないので、多数に目を通して同じ概念を いろいろな説明で読むほうがよいと思う。

まずは型を考える

(ちゃんと理解できているとは言い難いが)CNNによる画像認識は 次の図のように「画像」を各層で変換していって最終的に分類結果 (どれに該当するかを確率で示す)を得るようだ。

f:id:eijian:20160830205123p:plain

そこでまずは「画像(Image)」型と「層(Layer)」型を定義しよう。 カラー画像だと、$X \times Y$二次元のドットで構成された平面(Plain)が 赤緑青の3色(チャネルとする)集まってできている。一般の画像ファイル では各ドットを0-255の整数で表すことが多いが、CNNでは強さを実数で表すようだ。 畳み込み層の変換後のデータも多チャンネルの「画像」とみなせそうなので、 Imageは次のように定義した。実体はDoubleの3次元リストだ。

(Image.hs)

type Plain = [[Double]]    -- 2D: X x Y pixels
type Image = [Plain]       -- n channel

わざわざImagePlainのリストにしたのは、Plainごとに処理することが 多そうだからだ。 あとついでに、「学習」の際に必要となる教師データも定義しておく。

(Image.hs)

type Class = [Double]    -- trainer vector
type Trainer = (Image, Class)

実数のリストだが、教師データは要素のどれか一つが1.0、他が0.0となる。

次にレイヤを考える。CNNでは次のようなレイヤが使用される。

  • 畳み込み層 (convolution layer)
  • プーリング層 (pooling layer)
  • 活性化層/活性化関数 (activation layer)
  • 全結合層 (fully connected layer)

このうち、活性化層というのは一般的な表現かどうかよくわからない。 最初は型クラスとしてLayerを定義し各層を型としてやればいいと考えた。 が、各層をつなげて一連の変換処理を表現するのに「リスト」を使いたかったので (私の知識範囲では)うまい方法が思いつかなかった。 ということで、Layer型を次のように代数データ型で定義することにした。

(LayerType.hs)

data Layer = NopLayer
           | ConvLayer Int [FilterC]
           | MaxPoolLayer Int
           | ActLayer ActFunc
           | HiddenLayer [FilterH]
           | FlattenLayer (Image -> Image)

NopLayerは何も変換しないダミーの層(テスト用)、 プーリング層は実用上Max Poolingだけでよいと考えてMaxPoolLayer、 全結合層は隠れ層ということでHiddenLayerとしている。 各行に出てくる引数の型はそれぞれ以下のように定義してある。 ActFuncは活性化関数、FilterCは畳み込み層のためのフィルタ、 FilterHは全結合層(隠れ層)のためのフィルタ、をそれぞれ表す型だ。

(LayerType.hs)

-- for Activation Layer
type ActFunc = [Double] -> [Double]

-- filter for Convolution Layer
type Kernel = [[Double]]
type Bias   = Double
type FilterC = (Kernel, Bias)

-- filter for Hidden Layer
type FilterH = [Double]

学習用プログラム

deep learningで画像認識させるためのステップは大雑把には以下の手順を踏む 必要があると思っている。

  1. 大量のサンプル画像を用意し、それらを分類する。
  2. サンプル画像から「教師データ」を作成する。
  3. 教師データを使って「学習」を行う。(学習とは、各層で使われるフィルタを 教師データを元に調整していくこと、と理解)
  4. 画像を調整済みフィルタを使って分類する(各分類の確率を計算する)。

1と2はプログラム以前の準備段階だ。が、yusugomori実装では教師データと 学習状況の確認に使うテストデータをプログラム内で生成している。 最終的には学習データとテストデータは外から与える形にしたいが、 まずは同じようにプログラム内で生成するようにしよう。

教師データの生成

ゆくゆく教師データを外から与えることも考え、「画像データの集合」 の型??Poolを作ろう。これらの集合から幾つかの画像を取り出して 学習に使いたい。そこで、データに追番をつけてMapに登録 することにした。とりあえずはオンメモリで処理できれば良いので 以下のようにプール(オンメモリ版、MemPool)を定義した。

(Pool.hs)

class Pool p where
  getImages :: p -> Int -> Int -> IO [Trainer]
  nSample :: p -> Int

newtype MemPool = MemPool { m :: Map.Map Int Trainer }

yusugomori実装に倣って教師データを生成する関数も用意した。

(Pool.hs)

initSamplePool :: Int -> (Int, Int) -> Int -> Double -> Int -> IO MemPool
initSamplePool c (sx, sy) o p n = do
  s0 <- forM [0..(n-1)] $ \i -> do
    let cl = i `mod` o  -- class of this image

    -- Image data
    s1 <- forM [1..c] $ \j -> do
      s2 <- forM [0..(sy-1)] $ \y -> do
        let p' = if y `div` st == cl then p else (1-p)
        s3 <- forM [1..sx] $ \x -> do
          a <- pixel p'
          return a
        return s3
      return s2

    -- Trainer data
    e1 <- forM [0..(o-1)] $ \j -> do
      return $ if j == cl then 1.0 else 0.0
    return (s1, e1)

  return (MemPool (Map.fromList $ zip [0..] s0))
  where
    st = sy `div` o
    pixel :: Double -> IO Double
    pixel p = do
      v <- MT.randomIO :: IO Double
      let v' = if v < p then 0.5 else 0.0
      return v'

結構ごちゃごちゃしているが、以下のような3種類の画像データを生成してプール として返しているのだ。

f:id:eijian:20160830011222p:plain

実際は上記のようなキッチリした画像ではなく、灰色部分に一定割合で白が 混ざっていて、白色部分は逆に灰色が混ざっている。教師データではその 割合を5%、テストデータは10%としている。

つぎに、プールから指定枚数の画像を取り出す関数getImagesとプール内の 画像総数を返す関数nSampleを示そう。

(Pool.hs)

{-
getImages
  IN : pool
       epoch number
       batch size
-}

instance Pool MemPool where
  getImages p@(MemPool m) e b = do
    let s = nSample p
        o = (e-1) * b `mod` s
        mx0 = o + b - 1
        mx2 = mx0 - s
        mx1 = if mx2 < 0 then mx0 else s - 1
        im0 = mapMaybe (\x -> Map.lookup x m) [o..mx1]
        im1 = mapMaybe (\x -> Map.lookup x m) [0..mx2]
    return (im0 ++ im1)

  nSample (MemPool m) = Map.size m

getImagesは若干ややこしいが、要は教師データを満遍なく学習させたいので epoch毎に違う画像を取り出すようにしている(最初のepochで1番目から10個 取り出したら次のepochでは11番目から10個取り出す)。 ゆくゆくは順番に取り出すのではなくランダムに取り出すオプションも実装したい。

main関数(全体の流れ)

ここで、全体の流れを説明しておきたいのでmain関数の話をしよう。 大雑把には次のような流れだ。

  • パラメータ定義(main以前):フィルタの大きさとかバッチサイズとか (将来的にパラメータを設定ファイルから読み込むようにしたい)
  • 教師データ、テストデータの準備
  • 各レイヤの初期フィルタの生成
  • レイヤの組み立て
  • 学習の繰り返し

まずパラメータの説明から。画像は3種類に分類する(k)。教師データが 50x3=150個、テストデータが10x3=30個だ。画像サイズ、フィルタは 先に示した各レイヤの仕様どおり。なおバッチサイズは教師データが少ないので 毎回全教師データを学習するために150個としている。

(Main-cnn.hs)

-- PARAMETERS

-- training/test dataset
k = 3   -- number of class
n = 50  -- number of teacher data for each class
m = 10  -- number of test data for each class

train_N    = n * k
test_N     = m * k

-- image size
image_size = [12, 12]
channel    = 1

-- filter specs
n_kernels    = [10, 20]
kernel_sizes = [3, 2]
pool_sizes   = [2, 2]
n_hidden     = 20
n_out = k

-- loop parameters
epochs = 200
opStep = 5
epoch0 = 1
batch  = 150

-- others
learning_rate = 0.1

次にメインルーチンだ。教師データとテストデータは次のように生成 している。説明はいらないだろう。

(Main-cnn.hs)

main :: IO ()
main = do
  putStrLn "Building the model..."
  poolT <- initSamplePool 1 (12, 12) 3 0.95 train_N   -- trainer data pool
  poolE <- initSamplePool 1 (12, 12) 3 0.90 test_N    -- test data pool
  sampleE <- getImages poolE 1 test_N

次はレイヤの定義だ。まずフィルタを初期化してからレイヤを複数並べる。

(Main-cnn.hs)

  fc1 <- initFilterC 10 1 12 12 3 2
  fc2 <- initFilterC 20 10 5 5 2 2
  fh1 <- initFilterH n_hidden (2*2*20)
  fh2 <- initFilterH n_out n_hidden

  let layers = [ConvLayer 3 fc1, ActLayer relu, MaxPoolLayer 2,
                ConvLayer 2 fc2, ActLayer relu, MaxPoolLayer 2,
                FlattenLayer flatten,
                HiddenLayer fh1, ActLayer relu,
                HiddenLayer fh2, ActLayer softmax]

ここでは、畳み込み層1(活性化関数ReLU)→プーリング層1→ 畳み込み層2(活性化関数ReLU)→プーリング層2→隠れ層(活性化関数ReLU)→ 出力層(活性化関数softmax) という多層構造にした。 なお、FlattenLayerは畳み込み+プーリングの処理から全結合へデータ形式を 変換する処理である。

最後が学習の繰り返しだ。

(Main-cnn.hs)

  tm0 <- getCurrentTime
  putStatus 0 tm0 [([0.0], 0.0)] layers
  loop is tm0 layers batch poolT sampleE

時間の記録と学習前の状態をダミーで表示したあと、繰り返し学習する loopを呼んでいる。loopの詳細は次の通り。

(Main-cnn.hs)

loop :: Pool p => [Int] -> UTCTime -> [Layer] -> Int -> p -> [Trainer]
     -> IO ()
loop [] _ _ _ _ _ = putStr ""
loop (i:is) tm0 ls b pt se = do
  teachers <- getImages pt i b
  ops <- mapM (train ls) teachers
  let (_, dls) = unzip ops
      ls' = updateLayer ls dls   -- dls = diff of layers
  if i `mod` opStep == 0 then putStatus i tm0 (evaluate ls' se) ls'
                         else putStr ""
  loop is tm0 ls' b pt se

大したことはやってなくて、epoch毎に繰り返しloopを呼び、 epoch番号がなくなったら(空リストになったら)、繰り返しを終了する。 中身は、学習する教師データをプールから取り出して、それぞれ学習をし、 各層を学習結果を使って更新する。とはいえ、まだ「学習」は実装しておらず、 この部分はダミーである。

ここまでがmain関数だ。

長くなってしまったので、各層(レイヤ)の処理については次に 回すことにする。

まとめ

今回はCNNの画像認識プログラムの初回として以下に取り組んだ。

  • CNNで使われるいろいろな型を定義した。
  • 入力値を変換するいろいろな「層」を定義した。
  • これらを使って順伝播の処理を実装した。

次回は各層についての説明。

(ここまでのソースはこちら

CPUの創りかた(10): おまけ、アセンブラ

CPU自体は前回までで完成したので、次のネタに行ってもよかったのだが、

  • CPU(の実行プログラム)に与えるのがマシンコードだと、いちいち ハンドアセンブルするのが面倒くさい。
  • 別件でパーサを書く必要がありParsecライブラリに興味を持っていたので そのうち勉強したいと思っていた。

という理由により、今回はおまけとして簡易アセンブラを作ってみよう。

まずは文法の定義から

このCPU(TD4)は命令数が非常に少ないので、何となれば、単純な文字列 変換でもよかったのだが、Parsecの勉強も兼ねているからそこはちゃんと 文法の定義から入らなければなるまい。

文法定義とくればBNFだろう。大学時代に、わからないなりに適当な文法を 定義して遊んでいたのを思い出す。その頃の記憶を頼りにBNFで書き始めたのだが、 どうやら世間には EBNF というのがあるようなのでどうせならこちらにしよう、BNFの問題を 解決しているということだし。

ところでアセンブリ言語(ニーモニック?)の文法はどのように定義したら いいのだろう? と調べたら、 こういうの が見つかった。作りたいのはこれよりだいぶ簡素なものなので、ありがたく 参考にさせていただいた。EBNFはこれが初めてなので書き方が間違っているかも しれないが、とりあえず次のように定義してみた。

# EBNF for TD4 assembler

  program = instcode , { instcode } ;

  instcode = ( code2 | code1 ) , linefeed ;

  code2 = inst2, space, operand2 ;
  code1 = inst1, space, operand1 ;

  inst2 = 'add' | 'mov' ;
  inst1 = 'in' | 'out' | 'jnc' | 'jmp' ;

  operand2 = register , "," , operand1 ;
  operand1 = register | imdata ;

  register = 'a' | 'b' ;

  imdata = 4 * digit2 ;
  digit2 = "0" | "1" ;

  space = white space , { white space } ;
  white space = '\x20' | '\x09' ;

  linefeed = '\x0a' | ( '\x0d' , '\x0a' ) ;

2オペランド命令のANDとMOV、それ以外のは1オペランド命令なので定義が異なる。 あとはレジスタとか直値とかを定義していけばよさそう。ひとまず形になった 気がするので先に進もう。

パーサの前に

どうもParsecを使えば、EBNF(BNF)の定義からある程度容易にパーサの プログラムに落とし込める、という話があるらしい。ならEBNFができている のだから簡単に作れるよね、と思ったのは甘すぎだった。

まず参考となるWebサイトをいろいろ読んでみて、こう作ればよさそうという 感触を得たかったのだが、ちっとも理解できなかった。今なら何を理解できて いなかったのか分かるが、当初は本当にわからなかったのだ。

要は、ソースをパーサが字句解析したあと、それをどのようにして次の処理 (意味解析、コード生成)に結びつけたらいいのかイメージが湧かなかったため。 Web上でありがちな例は、"add 1,2"みたいなのが入力されたら、"add"を 解析して「"add"という文字列を返す」みたいなやつだ。文字列を解析したい のに解析結果が文字列ならそれをさらに解釈する処理が必要で堂々巡りに 思えてしまった。転機はこの記事。 この記事では明快に「返す値を型で定義している」。

つまりパーサからは「解析木」かそれに類する構造化されたデータが出力されるのだ。 それをまず考えて(定義して)おかないとパーサが書けるわけがない。 ではどのような構造のデータがあればマシンコードに変換できるかを考えてみる。

TD4の命令は、命令コードと1つまたは2つのオペランドからなる(オペランドなしの 命令は無い)。そこで、これらを組で表すことにする。最終的にはこう。

data Inst = Add | Mov | In | Out | Jnc | Jmp deriving (Enum, Show)
data Operand = RegA | RegB | Imdata String deriving Show

type Mnemonic = (Inst, (Operand, Maybe Operand))

Mnemonicは二重の組だ。命令コードとオペランドの組からなる。 オペランドの組は、一つ目は必ず存在するので普通にOperandで、 二つ目は無いかもしれないのでMaybe型とした。オペランドはA,Bレジスタか 直値の三種類しかない。パーサはソースを解析してMnemonicを命令の数だけ リストにして返してくれればよいのだ(以後、解析リスト、とする)。

パーサを書く

ではパーサを書いてみようと思う。何となくだが末端から作って積み重ねて 行った方がわかりやすそうなので、個々のオペランドから始める。まずは 直値(immediate data)から。

-- EBNF
--   imdata = 4 * digit2 ;
--   digit2 = "0" | "1" ;

imdata :: Parser Operand
imdata = do
  im <- count 4 (oneOf "01")
  return $ Imdata im

先に定義したように、直値もOperand型の一つなのでパーサの型が Parser Operandになっている。直値は二進数だけを扱うことにし、 かならず4桁と決めた。なので、oneOfで0か1に限定し、それを4つ 連続して取り出したら返すようにした。EBNFの定義と見比べると、 決めた通りにプログラムを書けばよいのがわかる。

同様に、A,Bレジスタはこうなる。

-- EBNF
--   register = 'a' | 'b' ;

register :: Parser Operand
register = do
  rg <- (regA <|> regB)
  return rg

regA :: Parser Operand
regA = do
  rg <- string "a"
  return $ RegA

regB :: Parser Operand
regB = do
  rg <- string "b"
  return $ RegB

AかBかの区別をつけるため、それぞれ別にパーサを定義し、それを<|>で 合わせてレジスタのパーサとした。

あとは同様に、EBNFをもとに以下のようなパーサを作った。

program :: Parser [Mnemonic]
program = do
  pg <- many1 $ instcode
  return pg

instcode :: Parser Mnemonic
instcode = do
  cd <- code2 <|> code1
  many1 $ oneOf "\r\n"
  return cd

code2 :: Parser Mnemonic
code2 = do
  in2 <- inst2
  many1 space
  op2 <- operand2
  return (in2, op2)

code1 :: Parser Mnemonic
code1 = do
  in1 <- inst1
  many1 space
  op1 <- operand1
  return (in1, (op1, Nothing))

inst2 :: Parser Inst
inst2 = do
  i2 <- (string "add" <|> string "mov")
  let i = if i2 == "add" then Add else Mov
  return i

inst1 :: Parser Inst
inst1 = do
  i1 <- (string "in" <|> string "out" <|>
         try (string "jnc") <|> (string "jmp"))
  let i = case i1 of
            "in"  -> In
            "out" -> Out
            "jnc" -> Jnc
            "jmp" -> Jmp
  return i

operand2 :: Parser (Operand, Maybe Operand)
operand2 = do
  op2 <- register
  char ','
  op1 <- operand1
  return $ (op2, Just op1)

operand1 :: Parser Operand
operand1 = do
  op1 <- (register <|> imdata)
  return op1

EBNFの定義と見比べれば、それぞれ何をしているかは分かると思う。 ではテストしてみよう。テスト用プログラムはざっと以下のような感じ。

module Main where

import Control.Applicative hiding ((<|>))
import Data.Char
import Text.Parsec
import Text.Parsec.String
import Text.Parsec.Char

(型、パーサ定義は省略)

main :: IO ()
main = do
  parseTest inst2 "add"
  parseTest inst2 "mov"
  parseTest inst2 "abc"
  parseTest register "a"
  parseTest register "b"
  parseTest register "c"
  parseTest imdata "0100"
  parseTest imdata "10100"
  parseTest imdata "1012"
  parseTest operand1 "aa"
  parseTest operand2 "a 1100"
  parseTest code2 "add a,b"
  parseTest code2 "mov a,0011"
  parseTest code1 "jmp 1011"
  parseTest code1 "in  b"
  parseTest code1 "in  a"
  parseTest code1 "in  0110"
  parseTest instcode "add a,b\n"
  parseTest instcode "add a,b\r\n"
  parseTest instcode "add a,a"
  parseTest instcode "jmp 1100\r\n"

このように、パーサに食わせてみたい文字列とそれを処理するパーサ関数を 組みにして与えればいいようだ。パーサの解析に失敗した(=文法にそぐわない) 場合はエラーが返るのですぐわかる。ざっとみたところうまく動いている ように・・・見えなかった!

文法定義のやり直し

うまくいったと思っていたのに実はいけていなかった。add a,aとか in aという命令はTD4には存在しないがパーサを通ってしまうのだ。 文法定義が甘かったわけだ。2オペランド命令と1オペランド命令を 分けるだけでは個々の命令の細かい制限(Bレジスタしか指定できない等)が 表現できていなかった。

ということで、EBNFの定義を見直す。

EBNF for TD4 assembler (revision 1)

  program  = line , { line } ;
  line      = instcode , linefeed ;

  instcode = inst_add | inst_mov | inst_in | inst_out | inst_jump ;
  inst_add  = 'add' , space , register , "," , imdata ;
  inst_mov  = 'mov' , space , ( op_mov1 | op_mov2 | op_mov3 ) ;
  inst_in   = 'in'  , space , register ;
  inst_out  = 'out' , space , ( reg_b | imdata) ;
  inst_jump = ( 'jnc' | 'jmp' ) , space , imdata ;

  op_mov1   = register , "," , imdata ;
  op_mov2   = reg_a    , "," , reg_b ;
  op_mov3   = reg_b    , "," , reg_a ;

  register  = reg_a | reg_b ;
  reg_a     = 'a' ;
  reg_b     = 'b' ;

  (imdata以降は同じなので割愛)

ポイントは、個々の命令ごとにオペランドのパターンを記述したこと。 結局オペランドが共通なのはジャンプ命令だけだった。これに沿って パーサも書き直す。

instcode :: Parser Mnemonic
instcode = do
  cd <- inst_add <|> inst_mov <|> inst_in <|> inst_out <|> inst_jump
  return cd

inst_add :: Parser Mnemonic
inst_add = do
  in1 <- string "add"
  many1 space
  rg1 <- register
  char ','
  im1 <- imdata
  return (toInst in1, (rg1, Just im1))

inst_mov :: Parser Mnemonic
inst_mov = do
  in1 <- string "mov"
  many1 space
  op  <- try (op_mov1) <|> (try op_mov2 <|> op_mov3)
  return (toInst in1, op)

inst_in :: Parser Mnemonic
inst_in = do
  in1 <- string "in"
  many1 space
  rg <- register
  return (toInst in1, (rg, Nothing))

inst_out :: Parser Mnemonic
inst_out = do
  in1 <- string "out"
  many1 space
  op <- regB <|> imdata
  return (toInst in1, (op, Nothing))

inst_jump :: Parser Mnemonic
inst_jump = do
  in1 <- try (string "jnc") <|> (string "jmp")
  many1 space
  im <- imdata
  return (toInst in1, (im, Nothing))

op_mov1 :: Parser (Operand, Maybe Operand)
op_mov1 = do
  rg <- register
  char ','
  im <- imdata
  return (rg, Just im)

op_mov2 :: Parser (Operand, Maybe Operand)
op_mov2 = do
  op1 <- regA
  char ','
  op2 <- regB
  return (op1, Just op2)

op_mov3 :: Parser (Operand, Maybe Operand)
op_mov3 = do
  op1 <- regB
  char ','
  op2 <- regA
  return (op1, Just op2)

toInst :: String -> Inst
toInst s = case s of
             "add" -> Add
             "mov" -> Mov
             "in"  -> In
             "out" -> Out
             "jnc" -> Jnc
             "jmp" -> Jmp

先ほどと同じテストを流してみると・・・ちゃんとmov a,ain aが エラーになっている! パーサが出来上がった! (Applicativeスタイルへの対応はまた今度。上記の書き方だけでも まだ全然理屈を咀嚼できていない)

コード生成

さて、パーサはソースプログラムの解析をするだけなので、そこから 目的のマシンコードを作り出す処理が必要だ。次はコード生成部分を作ろう。 と言っても、パーサがきちんと解析リストを生成してくれれば あとは簡単だ。AならBというふうに対応するマシンコードを返せば良い。 TD4は命令数が非常に少ないので、全部列挙することにした。

generate :: Either ParseError [Mnemonic] -> [String]
generate (Left s)       = [show s]
generate (Right [])     = []
generate (Right (x:xs)) = (translateOne x):(generate $ Right xs)

translateOne :: Mnemonic -> String
translateOne (Add, (RegA, Just (Imdata s))) = "0000" ++ s
translateOne (Mov, (RegA, Just RegB))       = "00010000"
translateOne (In , (RegA, Nothing))         = "00100000"
translateOne (Mov, (RegA, Just (Imdata s))) = "0011" ++ s
translateOne (Mov, (RegB, Just RegA))       = "01000000"
translateOne (Add, (RegB, Just (Imdata s))) = "0101" ++ s
translateOne (In , (RegB, Nothing))         = "01100000"
translateOne (Mov, (RegB, Just (Imdata s))) = "0111" ++ s
translateOne (Out, (RegB, Nothing))         = "10010000"
translateOne (Out, (Imdata s, Nothing))     = "1011" ++ s
translateOne (Jnc, (Imdata s, Nothing))     = "1110" ++ s
translateOne (Jmp, (Imdata s, Nothing))     = "1111" ++ s
translateOne _                              = error "no such mnemonic"

最後の行は保険だ。パーサがきちんと解析できていれば定義されていない 命令の解析リストが入力されることは無いだろうから。コンパイルして実行 してみる。

$ cabal build td4asm
  :
$ echo "mov a,b" | dist/build/td4asm/td4asm
00010000

おお!正しく変換されている!ちなみに、入力の最後に改行が無いとエラーになる。

$ echo -n "mov a,b" | dist/build/td4asm/td4asm
"TD4 asm" (line 1, column 8):
unexpected end of input
expecting new-line

構文解析で、命令行の最後は改行で終わるように定義しているからだ・・・。 これぐらいの制約は多めに見てもらおう。

再びラーメンタイマーの実行

前回はラーメンタイマープログラムのマシンコードを手で作ってCPUに入れていた。 今回はアセンブラアセンブリ言語のソースを入れてマシンコードを生成させて CPUへ入れたいと思う。アセンブラは生成したマシンコードを標準出力に出すので そのままパイプでCPUへつないであげればいいのだ。

* 見やすくするため実際には指定しているディレクトリなどは省いている *

$ cat ramen.a
out 0111
add a,0001
jnc 0001
add a,0001
jnc 0011
out 0110
add a,0001
jnc 0110
add a,0001
jnc 1000
out 0000
out 0100
add a,0001
jnc 1010
out 1000
jmp 1111

$ td4asm < ramen.a | td4
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0111][PC:0001]
step 2; [CF:0][A:0001][B:0000][OP:0111][PC:0010]
^C

ちゃんと動いている! やはりハンドアセンブルより断然楽ちんだ。

まとめ

さて、今回で本当に"CPU回"はおしまい。記事を細切れにしたせいもあって 10回にもなってしまった。長丁場だったが、CPUネタは面白い取り組みだったなあ。 どこかにi4004の回路図落ちてないかな :-)

(ここまでのソース)

CPUの創りかた(9): CPUはじめました

さあ、前回までで必要なモジュールは出揃った。今回はそれらを組み立てて 動くCPUを作ってしまおう! 一気に最終形はしんどいので少しずつピースを埋めていく感じで進めていきたい。

なお、しつこいようだがここで作っているCPUは以下の本で解説されている TD4という名前のオリジナル4bit CPUだ。説明中にtd4と出てくるのは その名前である。

ステップ 0: 電源、クロックジェネレータ(に相当するところ)

これまで論理回路の細かいところやCPU内の各種モジュールを作ることばかり やってきて、実行できるプログラムにする部分には目を瞑っていた。 しかしさすがに今回は「プログラム」を動かしたいのでそうはいかない。

そこでステップ0として完動させるための周辺部分を作っていこう。 電子工作では電源モジュールだとかクロックジェネレータとかその他の アナログ回路部分に相当するだろうか。

まず仕様を列挙しよう。

  • 「プログラム」は標準入力から投入する。
  • 「プログラム」は'0'と'1'の連続した文字列とする。また間にホワイトスペースが いくつ入ってもよい。
  • 「プログラム」におけるビット並びは(慣れているので) MSB...LSB の順とする。
  • プログラムカウンタが4 bit なので「プログラム」は16 bytes = 128文字。 ただしそれより少ない場合は'0'で補填する。多い場合は切り捨てる。
  • コマンドライン引数は順に"クロック間隔"と"入力ポート"の2つ。省略可能だが、 クロック間隔だけを省略することはできない。
  • クロック間隔の単位は秒、小数も使える。デフォルト値は1.0秒。入力ポートは 4桁の二進数でデフォルト値は"0000"。
    • 例) td4 0.5 0101 < program

ではこの仕様に基づいて作っていこう。

main :: IO ()
main = do
  pg <- getContents
  opts <- getArgs
  let (clock, iport) = parseOpts opts
  putStrLn ("clock " ++ (show clock) ++ " sec; I/P " ++ toStr iport)
  -- CLR(1),CF(1),A(4),B(4),OP(4),PC(4)
  let stat = toBits "011000010011000010"
  loop 0 clock lc_td4 stat iport (createRom pg)

getContentsで標準入力を読み、getArgsコマンドライン引数を取り込む。 どちらも標準で用意されている関数だ。parseOptsでオプションを解析している。

parseOptsは次の通り。 コマンドライン引数の数に応じてその値を読み込んだりデフォルト値を 使ったりしている。

defClock :: Double
defClock = 1.0              -- default clock time = 1 sec
defInput :: [Bin]
defInput  = toBits "0000"   -- default value of Input port

parseOpts :: [String] -> (Double, [Bin])
parseOpts [] = (defClock, defInput)
parseOpts (x:[]) = ((read :: String -> Double) x, defInput)
parseOpts (x:y:_) = ((read :: String -> Double) x, toBits y)

次は「プログラム」の整形についてだ。上記仕様ではビットの並びは MSB...LSBだが、これまで作ってきた論理回路モジュールでは、 入力(Binの配列)が全てLSB...MSBの順だ。そこで前もって順序を入れ替えておこう。 それを行っているのがcreateRom

createRom :: String -> [Bin]
createRom rs = concat $ map reverse $ split8 rs'
  where
    rs' = take 128 (toBits rs ++ repeat sLO)   -- 128 bits = 16 bytes

そのまんまだが。。。入力された文字列をtoBitsBinの配列にし、 足りなければsLO(=0)を付け加えて、先頭から16 bytes取り出している。 このような大雑把な記述が可能なのはHaskellの遅延評価のおかげだなぁ。 ちなみにtoBitsは0と1以外の文字は無視するので、間にスペースや 改行があっても問題ない。あとは8 bits単位に切り出してそれぞれを 逆順に並べ替えれば完成だ。

さあ、いよいよCPUモジュールを駆動する(呼び出す)ところだ。これは 「クロックの立ち上がり」のたびに関数を呼び出す無限ループである。 前回までに解説したように状態はCPUの外で管理することにしたので、 入力はCPUの状態+ROMの内容、出力はCPUの最新状態だ。 それをループにしたいのだ。出力値を次の入力値(の一部)に使うので、 再帰呼び出しが良さそうだ。ということで次のようなloop関数を作ってみた。

loop :: Int -> Double -> LogicCircuit -> [Bin] -> [Bin] -> [Bin] -> IO ()
loop s w lc st ip pg = do
  let os = lc (st ++ ip ++ pg)
  putStatus s os
  threadDelay $ floor (w * 1000 * 1000)
  -- set CLR to HI and take status from output
  let st' = [sHI] ++ (take 17 os)
  loop (s+1) w lc st' ip pg

クロックの度に状態を画面に出力したいので、ループの数(=ステップ数s)を 引数の最初に入れている。次はクロック間隔w、3番目(lc)がCPUを表す関数 (以後、CPU関数と呼ぼう)だ。 CPU関数への入力は「状態」「入力ポート値」「ROM」の3つ。

関数を呼び出して得た出力を画面に表示(putStatus)し、クロック間隔だけ待ち (threadDelay)、状態を更新して次のステップを呼び出す。この繰り返し。

次回の入力値を作っている少々奇妙な部分について。

  let st' = [sHI] ++ (take 17 os)

入力の最初の値が「リセット信号」を表しており、これがLOだとリセットが かかるようになっている。だから一番最初の呼び出し以外はHIにしないといけない。 あと、出力値から先頭の17個を取り出しているがこれには以下が含まれている。

  • carryフラグ(1 bit)
  • Aレジスタ(4 bit)
  • Bレジスタ(4 bit)
  • 出力ポート値(4 bit)
  • プログラムカウンタ(4 bit)

これで"評価ボード"(?)ができた。早速ダミーのCPU関数で動かしてみよう。 中身は何もせず入力を出力に回すだけ。

lc_td4_st0 :: LogicCircuit
lc_td4_st0 xs = concat [cf, a, b, op, pc]
  where
    [_, cf, a, b, op, pc, _, _] = splitInput xs

splitInput :: [Bin] -> [[Bin]]
splitInput xs = [cl, cf, a, b, op, pc, ip, rom]
  where
    (cl, xs0) = splitAt 1 xs
    (cf, xs1) = splitAt 1 xs0
    (a , xs2) = splitAt 4 xs1
    (b , xs3) = splitAt 4 xs2
    (op, xs4) = splitAt 4 xs3
    (pc, xs5) = splitAt 4 xs4
    (ip, rom) = splitAt 4 xs5

splitInputで分割して必要なものを取り出して並べているだけ。 さあコンパイルして実行してみよう。

$ cabal configure
Resolving dependencies...
Configuring mkcpu-0.1.0.0...

$ cabal build
Building mkcpu-0.1.0.0...
Preprocessing executable 'td4' for mkcpu-0.1.0.0...
[7 of 7] Compiling Main             ( src/Main-td4.hs, dist/build/td4/td4-tmp/Main.o )
Linking dist/build/td4/td4 ...

$ echo "0000" | dist/build/td4/td4 
clock 1.0 sec; I/P 0000
step 0; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 1; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 2; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 3; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 4; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 5; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
^C

レジスタなどの状態が表示されている。ちなみにその適当な値は、 実はmainの中で指定してある。

  -- CLR(1),CF(1),A(4),B(4),OP(4),PC(4)
  let stat = toBits "011000010011000010"
  loop 0 clock lc_td4_st0 stat iport (createRom pg)

このstatだ。A、B、OP、PCの値はそれぞれ1、2、3、4にセットしてあるのだ。 先述の通り一番最初のビットはリセット(CLR)であり、最初だけは'0'に してある。が、この何もしないダミーCPUではリセット信号が使われないので Aレジスタなどは初期値が入ったまま(のように見えるの)だ。

兎にも角にも、まずはCPUを駆動する周辺回路に相当する部分は一応 動いたようだ。これを使ってCPUを最終形まで組み立てていこう。

ステップ 1: レジスタの使用

CPUは本来状態を保持したり更新したりして処理を進めていくものだ。 状態はレジスタに保持されているのだが、以前の回で書いたように、CPUの 1サイクルの最終段階でレジスタを更新(もしくは保持)している。 この部分だけを作ってみよう。前回示したブロック図では一番右端にある部分だ。

f:id:eijian:20160504150006p:plain

コードは以下。

lc_td4_st1 :: LogicCircuit
lc_td4_st1 xs = concat [cf', a', b', op', pc']
  where
    [cl, cf, a, b, op, pc, _, _] = splitInput xs
    v0  = toBits "0000"
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ cf)
    a'  = lc_register4 (cl ++ [sHI] ++ a  ++ v0)
    b'  = lc_register4 (cl ++ [sHI] ++ b  ++ v0)
    op' = lc_register4 (cl ++ [sHI] ++ op ++ v0)
    pc' = lc_counter4  (cl ++ [sHI] ++ pc ++ v0)

入力を切り出す部分は同じ。v0はダミー値だ。 フラグやレジスタの入力値をそれぞれレジスタモジュールやカウンタモジュールへ 入れているだけだ。またリセット信号(cl)もそれぞれに入れている。 実行してみよう。

$ echo "0000" | dist/build/td4/td4 
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0000][PC:0001]
step 2; [CF:0][A:0000][B:0000][OP:0000][PC:0010]
step 3; [CF:0][A:0000][B:0000][OP:0000][PC:0011]
step 4; [CF:0][A:0000][B:0000][OP:0000][PC:0100]
step 5; [CF:0][A:0000][B:0000][OP:0000][PC:0101]
^C

ステップ0の結果とはだいぶ変わっている。まず、リセット信号が入ったため、 Aレジスタなどの初期入力値は一旦クリアされて0になっているのがわかる。 さらに、ステップが進む毎にプログラムカウンタ(PC)がカウントアップされて いる!カウンタモジュールは前に作ってテストしているから当然こうなるの だが、実行プログラムとしてこの出力になるのはちょっと嬉しい! (CPUが動いているぞ、という感じがする)

ステップ 2: 加算器の追加

次に加算器を取り付けよう。加算器には入力が2つ必要だが、状態を確認する ためにAレジスタの値を使う。もう一方の値はROMから無理やり取り出そう。 ROMにはプログラムカウンタをつないで0番地から順に値を取り出すようにする。 取り出した8bitから下4bitを使ってAレジスタに足し、結果がAレジスタ入るように 配線する。もちろんcarryフラグも更新する。プログラムはこうだ。

lc_td4_st2 :: LogicCircuit
lc_td4_st2 xs = concat [cf', a', b', op', pc']
  where
    [cl, _, a, b, op, pc, _, rom] = splitInput xs
    rdata = lc_rom16 (pc ++ rom) -- get data addressed by PC
    v0  = toBits "0000"
    im  = take 4 rdata
    (s0, c0) = splitAt 4 $ lc_adder (a ++ im)
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ c0)
    a'  = lc_register4 (cl ++ [sLO] ++ a  ++ s0)
    b'  = lc_register4 (cl ++ [sHI] ++ b  ++ v0)
    op' = lc_register4 (cl ++ [sHI] ++ op ++ v0)
    pc' = lc_counter4  (cl ++ [sHI] ++ pc ++ v0)

4行目でROMの現在番地の値から下4bitを取り出し、5行目(lc_adderのある行)で Aレジスタと足しあわせている。それをs0, c0にしてそれぞれAレジスタと carryフラグへ入れている。 Aレジスタの方は引数の2つ目(レジスタのLD入力)をHIではなくLOにしている。 これは保持している値ではなく外から入った値(s0)をセットするためだ。

プログラム(とは言えないが)は下4桁に加算したい数字を記載している。 上から、1,2,3,4,5,1である。

$ cat program
00000001
00000010
00000011
00000100
00000101
00000001

実行してみよう。

$ dist/build/td4/td4 < program
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0001][B:0000][OP:0000][PC:0001]
step 2; [CF:0][A:0011][B:0000][OP:0000][PC:0010]
step 3; [CF:0][A:0110][B:0000][OP:0000][PC:0011]
step 4; [CF:0][A:1010][B:0000][OP:0000][PC:0100]
step 5; [CF:0][A:1111][B:0000][OP:0000][PC:0101]
step 6; [CF:1][A:0000][B:0000][OP:0000][PC:0110]
^C

1から5まで足すとAレジスタが最大値の15になり、 そこに1を足せばcarryフラグが立ってAが0になるという寸法だが、CPUの 出力も確かにそうなっているのがわかる。

ちなみに、ここまでスラスラ進んでいるように書いているが、実際は 入力値の区切り位置を間違っていたり、入力と出力のパラメータの順序を 間違っていたりして、何度も出力が予想外になってバグ取りが大変だった。 実際の電子工作では「配線間違い」に相当するのだろうか。。。

ステップ3の前に(オペランドの選択)

いよいよ全体を組み上げるわけだが、その前に加算器への入力(ブロック図の 左側に並ぶA,Bレジスタと入力ポート値、および0)を切り替える部分を考えよう。 これには以前作ったmultiplexerが使える。入力値は4bitなので、multiplexerを 4つ、入力値の各桁用に並べればよい。

selectInput :: [Bin] -> [Bin] -> [Bin] -> [Bin] -> [Bin] -> [Bin]
selectInput s a b ip z = concat $ map (\x -> lc_multiplexer4ch (s ++ x)) mi
  where
    mi = buildMultiplexerInput [a, b, ip, z]

    buildMultiplexerInput :: [[Bin]] -> [[Bin]]
    buildMultiplexerInput xs = map (\i -> pickBit i xs) [0..3]

    pickBit :: Int -> [[Bin]] -> [Bin]
    pickBit i xs = map (!!i) xs

最初の引数sでどの入力値を使うかを指定する。あとは入力値の各桁を 集めてきてmultiplexerへ入れてやれば、sが選択する入力値を出力してくれる。 テストは以下。うまくいっているようだ。

>>> let a  = toBits "1000"
>>> let b  = toBits "0110"
>>> let ip = toBits "0001"
>>> toStr $ selectInput [sLO, sLO] a b ip zero
"1000"
>>> toStr $ selectInput [sHI, sLO] a b ip zero
"0110"
>>> toStr $ selectInput [sLO, sHI] a b ip zero
"0001"
>>> toStr $ selectInput [sHI, sHI] a b ip zero
"0000"

いよいよ最後の組み立てを残すのみ。

ステップ 3: CPUの組み立て

さあ、最終段階にきた。ROMから読み込んだ命令を命令デコーダへ入れ、 入力値を選択し、加算器を通して結果をどこに書き出すかを 命令デコーダに指示させればよいのだ。以下が最終のCPUのコードだ。

lc_td4 :: LogicCircuit
lc_td4 xs = concat [cf', a', b', op', pc']
  where
    [cl, cf, a, b, op, pc, ip, rom] = splitInput xs
    rdata = lc_rom16 (pc ++ rom)
    (im, inst) = splitAt 4 rdata
    [sa, sb, ld0, ld1, ld2, ld3] = lc_inst_decorder (inst ++ cf)
    (s0, c0) = splitAt 4 $ lc_adder ((selectInput [sa, sb] a b ip zero) ++ im)
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ c0)
    a'  = lc_register4 (cl ++ [ld0] ++ a  ++ s0)
    b'  = lc_register4 (cl ++ [ld1] ++ b  ++ s0)
    op' = lc_register4 (cl ++ [ld2] ++ op ++ s0)
    pc' = lc_counter4  (cl ++ [ld3] ++ pc ++ s0)

ROMから取り出した命令をデコーダに入れ、その結果を加算器と各レジスタへ つないでいる。ステップ2との差はそれぐらいだが、これで完成だ。意外と あっけなく出来上がった。

さすがに各命令の処理をテストしておく必要があるだろう。以下がテスト用 コードの一部だ(ADD A,Im と MOV A,B)。 具体的な命令コードを与え結果を想定と比較する、これまで作ってきた 論理回路モジュールと同じだ。

>>> let rom0 = take ((16-1) * 8) $ repeat '0'

>>> -- ADD A,Im (A=1, Im=4 -> A=5, CF=0)
>>> toStr $ lc_td4 $ toBits ("10 1000 0000 0000 0000 0000 00100000" ++ rom0)
"01010000000001000"

>>> -- ADD A,Im (A=13, Im=4 -> A=1, CF=1)
>>> toStr $ lc_td4 $ toBits ("10 1011 0000 0000 0000 0000 00100000" ++ rom0)
"11000000000001000"

>>> -- MOV A,B (A=13, B=3 -> A=3)
>>> toStr $ lc_td4 $ toBits ("10 1011 1100 0000 0000 0000 00001000" ++ rom0)
"01100110000001000"

さあ、実際にプログラムを走らせてみよう! 本に記載されているラーメンタイマーを実行してみよう。うまくいけば、 下記のように出力ポートが変化するはずだ。

[0111] -> [0110] -> [0100](点滅) -> [1000]

本には"ニーモニック"しか書いていないのでハンドアセンブルした結果が これ。

10110111
00000001
11100001
00000001
11100011
10110110
00000001
11100110
00000001
11101000
10110000
10110100
00000001
11101010
10111000
11111111

これをファイル(program.ramen)に書いてtd4に食わせればよい。

$ dist/build/td4/td4 < program.ramen
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0111][PC:0001]
step 2; [CF:0][A:0001][B:0000][OP:0111][PC:0010]
  :
step 31; [CF:0][A:1111][B:0000][OP:0111][PC:0001]
step 32; [CF:1][A:0000][B:0000][OP:0111][PC:0010]
step 33; [CF:0][A:0000][B:0000][OP:0111][PC:0011]
  :
step 63; [CF:0][A:1111][B:0000][OP:0111][PC:0011]
step 64; [CF:1][A:0000][B:0000][OP:0111][PC:0100]
step 65; [CF:0][A:0000][B:0000][OP:0111][PC:0101]
step 66; [CF:0][A:0000][B:0000][OP:0110][PC:0110]
  :
step 96; [CF:0][A:1111][B:0000][OP:0110][PC:0110]
step 97; [CF:1][A:0000][B:0000][OP:0110][PC:0111]
step 98; [CF:0][A:0000][B:0000][OP:0110][PC:1000]
  :
step 128; [CF:0][A:1111][B:0000][OP:0110][PC:1000]
step 129; [CF:1][A:0000][B:0000][OP:0110][PC:1001]
step 130; [CF:0][A:0000][B:0000][OP:0110][PC:1010]
step 131; [CF:0][A:0000][B:0000][OP:0000][PC:1011]
step 132; [CF:0][A:0000][B:0000][OP:0100][PC:1100]
  :
step 192; [CF:0][A:1111][B:0000][OP:0100][PC:1100]
step 193; [CF:1][A:0000][B:0000][OP:0100][PC:1101]
step 194; [CF:0][A:0000][B:0000][OP:0100][PC:1110]
step 195; [CF:0][A:0000][B:0000][OP:1000][PC:1111]
  :

命令種やビット数に制限があり正確に3分とはいかないが、ちゃんと想定したとおりの 動きをしているようだ。これだけでも意外にうれしいものだ!

まとめ

やっとCPUが完成した!機能的にはかなり低レベルではあるが、本当に 論理回路の組み合わせだけでCPUという複雑な仕組みが成り立っていて 動くのを確認できた。最初に考え出した人は本当にすごい!

ところで、今回のプログラミングは「Haskellならでは」というのがあまり なかったなあと思う。論理的な処理は論理ゲートとその配線で決まるので、 Haskellらしさは「配線」に相当する処理ぐらいだ(複数の配線を mapで一括処理するとか)。それもあってか、プログラミング的には 淡々と並べただけに終わった気もする。Haskellのもっと高度な 機能を使えば、今よりエレガントな記述ができるのかもしれないが、それは もっと"使える"ようになってから考えよう。

さてこの先だが、

  • 8bit化: レジスタや加算器で1bitの部品を4個から8個に増やせばよい。 プログラム的には繰り返し回数を増やすだけなのでかなり簡単なはず。
  • 加算以外の命令: 本にも記載があるが本格的なALUを用意すればもっと いろいろできることが増える。
  • i4004の製作: 回路図がわかれば今回と同様になんでも製作できそうだ (山のように時間があれば)。ならば実在する有名どころを作ってみるのも 楽しそう。

などは手がつけられそうだ。こういった拡張は頭で考えるだけでも楽しいものだ。

次回だが、このネタのおまけで何かCPU関係にするか、新しいネタにするか、、、。

(ここまでのソース)

CPUの創りかた(8): すべては足し算だった

CPUの製作もいよいよゴールが近づいてきた。残る主要な部品は 命令デコーダ(instruction decorder)のみ。その前にCPUの全体構造を 確認しておきたい。そうすれば命令デコーダをどう作ればいいかがよく 分かると思う。次に命令デコーダを考えよう。

CPUのブロック図と命令デコーダ

CPU(TD4)を入力、処理、出力に分けて考えると次のようなブロック図に なると思う。

f:id:eijian:20160504150006p:plain

この図では、入力値(入力ポート(I/P)、A,Bレジスタ、直値)は加算器 (adder)に入り、結果は全て加算器から出てきて出力ポート(O/P)やレジスタの いずれかへ書き込まれるような構成になっている。 前回作った加算器では、入力は二つの4 bit値、出力は一つの4 bit値とcarry フラグだった。ということは、加算器に入れる値をいくつかの候補から選択する 必要がある。それが図のSelector INだ。 一方、加算結果を書き出す先も選択する必要があるが、これはSelector OUTの 仕事だ。

では、どれを足してどこに結果を書くかはどう決まるのかというと、もちろん 命令(operation code)で決まるわけだ。ここまでくると命令デコーダの やっていることが見えてくる。このCPU(TD4)においては、 「命令を解釈して、どの値を足して、どこへ書き出すかを決める」 ということをやっているわけだ。ROMから命令デコーダへ入っている矢印が 「命令」を、命令デコーダからSelector IN/OUTへ伸びている矢印が どれを選択するかを表す情報というわけだ。

すべては足し算

さて、上記の図と説明を見て「あれ?」と思っただろうか? そう、TD4の構造では全ての命令(処理)は加算器を通る。それはつまり「足し算しか できない」ということだ。しかしCPUの処理というものは、代入とかポートの 入出力とか、ジャンプとかいろいろあるわけで、それはどこへ行ったんだ、 となる。

これについては例の本(CPUの創りかた)を読むとよくわかる。一見全く別の 命令に見えたものが、実は「全部足し算で表現できる」のだ。次の表を 参照してほしい。これは命令コード(opcode)とニーモニック、それに 対応する足し算を並べたものだ。

opcode mnemonic equation
0000xxxx ADD A,Im A ← A + Im
00010000 MOV A,B A ← B + 0000
00100000 IN A A ← I/P + 0000
0011xxxx MOV A,Im A ← 0000 + Im
01000000 MOV B,A B ← A + 0000
0101xxxx ADD B,Im B ← B + Im
01100000 IN B B ← I/P + 0000
0111xxxx MOV B,Im B ← 0000 + Im
10010000 OUT B O/P ← B + 0000
1011xxxx OUT Im O/P ← 0000 + Im
1110xxxx JNC Im PC ← 0000 + Im (cf=0)
1111xxxx JMP Im PC ← 0000 + Im

確かに全部足し算で表現できている!恥ずかしながら、今までCPUの内部でこういう処理を しているとは知らなかった。。。代入は足し算を使って(片方を0にして) 実現しており、出力ポートへの書き出しとは書き込む先がポートになっている、 ただの代入だ。ジャンプは飛び先アドレスをPCへ代入しているだけだ。

式中のプラス記号の左と右を比べてほしい。プラス記号の左側はSelector INから 来るもの、右側はROMから直接入ってくるものだ。左側には図に書かれていない '0000'もある。これも実際に作るときにはSelector INの入力の一つと してつくり込む必要がある。一方、右側は直値(Im)か0だ。それをどうやって 使い分けているか、は簡単だ。命令コードの下4 bitがそれで、 0を入れる時は命令コードの下4 bitが0なのだ。

もちろん一般のCPUはもっと複雑な構造をしているので、 こんな単純に全部足し算で処理しているのではなく、ALU内で複雑な処理を しているのだろう。とはいえそれも、加算器とか乗算器とか、論理演算 などを処理する部分にどういったものを入力するか、それらから出てくる 結果をどこへ書き出すかを命令デコーダが選択すればよいのだろうから、 大差無いだろう。

「命令を解釈する」と言われると、なんとなく神秘的でとても高度なことを しているというイメージがあったのだが。 今回の製作でその辺のカラクリが見えたのは収穫だ。 (知っている人にとっては当たり前すぎて、バカバカしいと思われるだろうが)

命令デコーダ

では、最後の大物(?)にとりかかろう。先述の通り、このCPUの命令デコーダは 「命令を解釈してSelector IN/OUTを適切に制御する」ものだ。 ではSelectr IN/OUTは何者かということになる。まあ本に書いてあることだが、

  • Selector IN: 4chマルチプレクサを4 bit分並べたもの。一つのマルチプレクサで 1 bit分が選択できるので4つ必要。
  • Selector OUT: どのレジスタやポート(実体はFlip Flop)のLoadを有効に するかを指定する。実体があるわけではなく命令デコーダの一部と言える。 (わかりやすくするためSelector OUTと表現しているが)

ということだ。結局命令デコーダは、Selector IN/OUTに必要な 「選択するための情報」を作り出せばいいのだ。

さて、命令デコーダの入力は命令コード(上4 bit、OP0-OP3)とcarryフラグ(CF)だ。 Select INへは(SA、SBの)2 bit、4種類の情報を与えてやれば良い。 Selector OUTについては、4つある出力先(Flip Flop)のどれか一つの'Load'を 有効に、それ以外を無効にする必要があるので、4 bitの情報で表す(LD0-LD3)。 まとめると、入力は5 bit(OP0-OP3,CF)、出力が6 bit(SA,SB,LD0-LD3)だ。

SA,SBの値と、選択される入力値の対応は下表のとおり。

SA SB input
0 0 A register
1 0 B register
0 1 input port
1 1 value '0000'

LD0-LD3と出力先の対応は下表のとおり。

LOAD output
LD0 A register
LD1 B register
LD2 output port
LD3 program counter

これらを踏まえ入出力を真理値表にまとめたものがp.242にある。さらに カルノー図などで分析した結果、各出力は次のように表現できる。

{
  S_A = OP_0 + OP_3
}

{
  S_B = OP_1
}

{
  LD_0 = OP_2 + OP_3
}

{
  LD_1 = \overline{OP_2} + OP_3
}

{
  LD_2 = OP_2 + \overline{OP_3}
}

{
  LD_3 = \overline{OP_2} + \overline{OP_3} + \overline{OP_0} \cdot CF
}

非常に簡潔になった、素晴らしい。もちろん偶然ではなく、この本の著者は こういう単純な回路になるように綿密に考えて、命令コードを割り振っているわけだ。 ここまでくればコードに落とすのは簡単だ。

{- |
instruction decorder

  IN : [OP0,OP1,OP2,OP3,CF]
  OUT: [SA,SB,LD0,LD1,LD2,LD3]

    OPn: operation code (upper 4 bit)
    CF : carry flag
    SA : Select A
    SB : Select B
    LDn: LOAD0 - LOAD3
-}

lc_inst_decorder :: LogicCircuit
lc_inst_decorder (op0:op1:op2:op3:c:_) = [sa, sb, l0, l1, l2, l3]
  where
    sa = op0 |> op3
    sb = op1
    nop0 = (!>) op0
    nop2 = (!>) op2
    nop3 = (!>) op3
    l0 = op2  |> op3
    l1 = nop2 |> op3
    l2 = op2  |> nop3
    l3 = nop2 |> nop3 |> (nop0 &> c)

なんと命令デコーダの小さいことか!わずか12行だ!

各opcodeについてテストを書いてみる。

>>> toStr $ lc_inst_decorder $ toBits "00000"  -- ADD A,Im
"000111"
>>> toStr $ lc_inst_decorder $ toBits "10000"  -- MOV A,B
"100111"
>>> toStr $ lc_inst_decorder $ toBits "01000"  -- IN A
"010111"
>>> toStr $ lc_inst_decorder $ toBits "11000"  -- MOV A,Im
"110111"
>>> toStr $ lc_inst_decorder $ toBits "00100"  -- MOV B,A
"001011"
>>> toStr $ lc_inst_decorder $ toBits "10100"  -- ADD B,Im
"101011"
>>> toStr $ lc_inst_decorder $ toBits "01100"  -- IN B
"011011"
>>> toStr $ lc_inst_decorder $ toBits "11100"  -- MOV B,Im
"111011"
>>> toStr $ lc_inst_decorder $ toBits "10010"  -- OUT B
"101101"
>>> toStr $ lc_inst_decorder $ toBits "11010"  -- OUT Im
"111101"
>>> toStr $ lc_inst_decorder $ toBits "01110"  -- JNC(C=0) Im
"111110"
>>> toStr $ lc_inst_decorder $ toBits "01111"  -- JNC(C=1) Im
"111111"
>>> toStr $ lc_inst_decorder $ toBits "11110"  -- JMP Im
"111110"

ちゃんとテストも全部通る。これは嬉しい。 (本の真理値表とはビット順序がいくらか違っているので混乱しないように)

さあ、部品が出揃った!いよいよCPUとして全体を組み上げる頃合いだ! が、続きは次回。

まとめ

今回は最後の重要なモジュールである命令デコーダを作った。 次回はCPUを組み上げて簡単な命令の実行を試そう。うまく動くかな?

(ここまでのソース)

CPUの創りかた(7): 加算器を作る

今回はCPUの中でも「純粋に」計算するところ、加算器を作ろう。ALUを作る!と 行きたいところだが、残念ながらTD4(しつこいようだが下記の本で実装されるCPU) では加算器しかないので今の所は諦めよう。

1 bit 全加算器

ものの本では半加算器から説明されるようだが、どうせ使わないので 最初から全加算器を作る。まずは1 bitから。

f:id:eijian:20160411200537p:plain

入力は、2つの値$A$と$B$、あとくり上がり$C_i$(carry)の3つ、出力は加算の結果 $S$と上の桁へのくり上がり$C$だ。

真理値表は次のようになるだろう。以下の説明では従来のHI、LOの代わりに 1と0を使うことにする。

$A$ $B$ $C_i$ $S$ $C$
0 0 0 0 0
0 1 0 1 0
1 0 0 1 0
1 1 0 0 1
0 0 1 1 0
0 1 1 0 1
1 0 1 0 1
1 1 1 1 1

例によって、これをSとCについてそれぞれカルノー図を描いて考えると、 最終的に以下のように表される(面倒なので途中の計算は省略する)。

{
S= Ci \oplus (A \oplus B)
}

{
C= A \cdot B + Ci \cdot (A \oplus B)
}

コードはこの式をそのまま表してみた。

{- |
  1 bit full adder

  IN : [Ci,A,B]
  OUT: [S,C]

    Ci  : carry in
    A,B : value
    S   : answer
    C   : carry out
-}

lc_adder1 :: LogicCircuit
lc_adder1 (ci:a:b:_) = [s, c]
  where
    a_xor_b = a <+> b
    s = ci <+> a_xor_b
    c = (a &> b) |> (ci &> a_xor_b)

全加算器は大した回路ではないので、テストもすんなりと通った。

>>> lc_adder1 [sLO, sLO, sLO] == [sLO, sLO]
True
>>> lc_adder1 [sLO, sLO, sHI] == [sHI, sLO]
True
>>> lc_adder1 [sLO, sHI, sLO] == [sHI, sLO]
True
>>> lc_adder1 [sLO, sHI, sHI] == [sLO, sHI]
True
>>> lc_adder1 [sHI, sLO, sLO] == [sHI, sLO]
True
>>> lc_adder1 [sHI, sLO, sHI] == [sLO, sHI]
True
>>> lc_adder1 [sHI, sHI, sLO] == [sLO, sHI]
True
>>> lc_adder1 [sHI, sHI, sHI] == [sHI, sHI]
True

4 bit 加算器

1 bit 全加算器ができたので、これを数珠つなぎにすれば$n$ bitの加算器を 作ることができる。本CPU(TD4)はレジスタが4 bitなので4 bit 加算器を作ろう。 言わずもがなだが、全加算器を次のようにつなげれば良い。入力$A_n$、$B_n$ および$C_i$、出力が$S_n$と$C$だ。

f:id:eijian:20160411200544p:plain

よって、先に作った全加算器(lc_adder1)をこの図のとおりにつなげて 計算すればよいのだが、数珠つなぎだからこれは再帰で表現できそうだ。 まず再帰処理の本体から考える。入力は$n$ bitの$A$と$B$だが、それらの同じ桁の 値を組にしたリストを用意しよう。これを順に加算器に入れ、各桁の結果$S_i$を 求めたい。大雑把には次の通り。

[(a0,b0), (a1,b1), ... , (an,bn)] → [s0,s1, ... , sn]

ただし、くり上がり(carry)があるので少々複雑になる。くり上がりも含め 計算するために、その桁に足し込みたいcarryも引数で与える。さらに、 計算結果$S$を次に渡していくので、これも引数にする。これらを踏まえると、 つぎのような関数adder_nができる。

{-
  IN : Ci, [Si], [(Ai,Bi)]
  OUT: [Si,C]

    Ci    : carry in
    Ai, Bi: value
    Si    : Answer
    C     : carry out
-}

adder_n :: Bin -> [Bin] -> [(Bin, Bin)] -> [Bin]
adder_n ci ss [] = [ci] ++ ss
adder_n ci ss ((a, b):ds) = adder_n c (s:ss) ds
  where
    [s, c] = lc_adder1 [ci, a, b]

$A$と$B$の各桁の組のリストを先頭から一組取り出して、 carryと合わせて加算する。残りのリストと計算結果を使ってadder_nを 呼出す(再帰)。処理するリストがなくなったら、$S$の前に最後の くり上がりを追加して終わり。この関数では「最後の桁」が先頭に来るから、 出来上がった$S$のリストは逆順にしないといけない。加算器の出力は、 $S_n$,$C$の順にしたいので、再帰の終端で最後のくり上がりを「先頭に追加」 したわけだ。逆順にしたら最後になるように。結局出力は以下のようになる。

[s0,s1, ... , sn, c]

あとは、このadder_nを呼び出す「外側」の関数を作れば良い。 lc_adderとする。

{- |
  n bit full adder

  IN : [A,B]
  OUT: [S,C]

    A,B : value
    S   : answer
    C   : carry out
-}

lc_adder :: LogicCircuit
lc_adder ds = reverse $ adder_n sLO [] $ zip a b
  where
    l2 = length ds `div` 2
    a = take l2 ds
    b = drop l2 ds

加算器のビット数を明示的に与えていないことに注意。ビット数は、 入力(リスト)の長さの半分とした。$A_n$と$B_n$が繋がったリストで入力される 前提だ(4 bitなら要素は8個のはず)。奇数個だったら余りは無視される。 これで$n$ bit加算器の完成だ。テストも上々だ(長くなるので割愛)。

まとめ

今回は少々軽めの内容だったが、「一番計算機らしい回路」と言えなくもない(?) 加算器を作った。先述の通りこのCPUは加算器しか持たないので乗算器などに 手は出さないが、将来CPUをパワーアップするとしたら要検討かもしれない。

これでCPUで必要となる論理回路のモジュールがほぼ出揃った。あと残すは 「本丸」の命令デコーダのみだ。だが次回はその前にCPUの全体構成を 確認したいと思う。その方が命令デコーダの作り方がよく分かるとおもうから。

(ここまでのソース)