動的計画法４ 連鎖行列積
問題定義
連鎖行列積問題 n個の行列の連鎖が与えられたとき，スカラーの乗算回数を最小化するように積 A1, A2, .... An を完全に括弧付けする問題です． （与えられた行列計算を実行するために，掛け算の回数が最小となる計算順序を決定することが目的．）
なかなか上の定義だと難解なので，自分解釈します．
０．行列の掛け算
二つの行列の掛け算をする時に必要な回数は，下記のように求まります．
行列演算は結合法則が成り立つので，どの順番で掛け算をしても結果は同じになりますが，結果を求めるための掛け算回数は下記のようにどの順番で計算するかに依存します．
上の例では，A1, A2を計算して，そのあとA3との掛け算をした場合，7500回，A2, A3の掛け算をして，A1の掛け算をした場合，75000回となり，演算回数が10倍も異なります．この問題は，行列の掛け算問題が与えられたときに掛け算回数が最小となる順序を求めるものです．
ということで，動的計画法で問題を解決する下記３ステップを通して，問題を解いていきたいと思います．
１．再帰的に求める方法を考える． ２．再帰関数の分岐から漸化式を作る． ３．漸化式を使って，ボトムアップで表を埋めていく
１．再帰的に求める方法を考える
行列の連鎖積が下記のように与えられたとき，
最小掛け算回数を返してくれる下記の関数があったとします．
get_min_mult_calc(i, j, mat_list) i : mat_list 内の行列開始端インデックス． j : mat_list 内の行列終了端インデックス． mat_list : 行列リスト．
で，この魔法の関数がうまく動くとします．であるとすると，(i, j) の連鎖積の最適値を求めるときには (i, k), (k+1, j) に対するこの関数のアウトプットを使えばよく，この k を変化させて最適値を調べます．ここで，(i, j) の連鎖積を最小にする k の位置は i + 1 ~ j - 1 のどこかにあります．この関数のコールイメージは下記のようになります．
i = j の時， 0 i ≠ j の時， min( get_min_mult_calc(i, i+1, mat_list) + get_min_mult_calc(i+2, j, mat_list) + pi-1pi+1pj, .... get_min_mult_calc(i, k, mat_list) + get_min_mult_calc(k+1, j, mat_list) + pi-1pkpj, .... get_min_mult_calc(i, j-2, mat_list) + get_min_mult_calc(j-1, j, mat_list) + pi-1pj-2pj, )
という形で，i ≠ j の時は k をひとつづつ動かして最小なものを選ぶループが必要になります．上記をコードに落とすと，下記のようになります．
import sys
class Mat(object):
def __init__(self, row, col):
self.row = row
self.col = col
def get_min_mult_cnt(i, j, mat_list):
val = 0
if (i == j):
val = 0
else:
val = sys.maxsize
for k in range(i, j):
p = mat_list[i].row
n = mat_list[k].col
q = mat_list[j].col
tmp = get_min_mult_cnt(i, k, mat_list) + get_min_mult_cnt(k+1, j, mat_list) + p*n*q
if (tmp < val):
val = tmp
return val
if __name__ == '__main__':
print("Matrix Multiplication Problem")
mat_list = [Mat(30, 35), Mat(35, 15), Mat(15, 5), Mat(5, 10), Mat(10, 20), Mat(20, 25)]
print(get_min_mult_cnt(0, 5, mat_list))
２．再帰関数の分岐から漸化式を作る
ということで，ここまでくればあとは機械的に行けそうです．漸化式は，下記のようになります．m[i, j] は，i から j までの行列連鎖積の最小掛け算回数です．
m[i, j] = 0 : （i = j の時．） m[i, j] = : （i ≠ j の時．） mini≦k≦j(m[i, k] + m[k+1, j] + pi-1pkpj)
漸化式の中にループが入っている点が今までと違うところだと思います．この漸化式を使って，実際に表を作っていきます．
３．漸化式を使って，ボトムアップで表を埋めていく
今回の場合，漸化式の中にループが入っているので，テーブルの計算方法が少し複雑になります．ここで，テーブルの構造を下記のように定義します．
m[i, j] : 連鎖行列積列の i 番目から j 番目までの最小掛け算回数を収めたテーブル．
三重のループになるので，テーブルの計算方法がちょっと複雑です．参考に4つの行列連鎖積の場合について，テーブルの計算順序を書き出してみました．３重になっている一番内側のループは，kを動かすためのループになります．kを動かすことによって，m[i, j]がそれ以前に計算されている要素からどのように求められているかは，下記の図の最後に書いてあります．
で，これをコードに落とすと，，，
import sys
import numpy as np
class Mat(object):
def __init__(self, row, col):
self.row = row
self.col = col
def get_min_mult_cnt_rec(i, j, mat_list):
val = 0
if (i == j):
val = 0
else:
val = sys.maxsize
for k in range(i, j):
p = mat_list[i].row
n = mat_list[k].col
q = mat_list[j].col
tmp = get_min_mult_cnt_rec(i, k, mat_list) + get_min_mult_cnt_rec(k+1, j, mat_list) + p*n*q
if (tmp < val):
val = tmp
return val
def get_min_mult_cnt_dp(i, j, mat_list):
multp = np.zeros((len(mat_list), len(mat_list)))
multp.fill(sys.maxsize)
# 長さ0連鎖の演算量
for i in range(len(mat_list)):
multp[i, i] = 0
for le in range(1, len(mat_list)):
for i in range(len(mat_list)-le):
j = i + le
tmp = sys.maxsize
for k in range(i, j):
p = mat_list[i].row
n = mat_list[k].col
q = mat_list[j].col
tmp = multp[i, k] + multp[k+1, j] + p*n*q
if (tmp < multp[i, j]):
multp[i, j] = tmp
print(multp)
return multp[i, j]
if __name__ == '__main__':
print("Matrix Multiplication Problem")
mat_list = [Mat(30, 35), Mat(35, 15), Mat(15, 5), Mat(5, 10), Mat(10, 20), Mat(20, 25)]
print("Solve Recursively")
print(get_min_mult_cnt_rec(0, 5, mat_list))
print("Solve Dynamic Programing")
print(get_min_mult_cnt_dp(0, 5, mat_list))
動的計画法３ 最長共通部分列
動的計画法シリーズ，３個目の例題です．前回のエントリの最後で説明したステップ通りに問題を解いていきたいと思います．今回の問題も有名な問題で，”最長共通部分列問題”と呼ばれるものです．問題の定義は，，，
問題定義
二つの与えられた文字列に対して，最長共通部分列の長さを出力しなさい． ※部分列とは，与えられた文字列を Z = {z1, z2, ... zn} とすると， zi1, zi2, ... zim (i1 < i2 < ... < im) を表すことのできる文字列であり， 二つの文字列の最長共通部分列とは，例えば "abcd", "becd" という文字列を考えた場合， "bcd"が最長共通部分列となる．
ということで，動的計画法で問題を解決する下記３ステップを通して，問題を解いていきたいと思います．
１．再帰的に求める方法を考える． ２．再帰関数の分岐から漸化式を作る． ３．漸化式を使って，ボトムアップで表を埋めていく
まず”１．再帰的に求める方法を考える”で問題を解いていきたいと思います．
１．再帰的に求める方法を考える
参考書を何冊か読んでみたんですが，なかなかピンとくるまで時間がかかりました．どの本を読んでも
”より小さい文字列の最長共通部分列を求める問題が，大きい文字列の問題の部分問題になっているので，再帰的に解ける”
とかって書いてあるんですが，ここがどうしてもピンとこず．．．自分なりに解釈してみました．
下記のように二つの文字列があったとします．この時，比較している文字の添え字をそれぞれ i, j とします．
比較文字列A：abcd
比較文字列B：becd
今，比較している文字列が (i, j) = (2, 2) だとすると，A[2] = 'c', B[2] = 'c' となるのでこのペアは共通文字列の一部になれることがわかります．よって次は (i-1, j-1) = (1, 1) を調べます．今度は A[1] = 'b', B[1] = 'e' となるため，このペアは共通文字列の一部にはなれません．ここの再帰関数の分岐のさせ方がキモで，次の調査候補は (i-1, j), (i, j-1) の二つのパターンになります．一致しなかった場合，片方の文字列を一文字捨てて，もう一度他方の文字列と比較します．文字列が二つあるので，結局分岐は二つになります．
関数を次のように定義すると
get_common_longest_subseq(a_idx, b_idx, seq_a, seq_b) a_idx : 関数コールで比較対象となっている，文字列Aの文字のインデックス． b_idx : 関数コールで比較対象となっている，文字列Bの文字のインデックス． seq_a : 文字列A seq_b : 文字列B
で，関数のコールイメージは下記のようになります．
比較対象文字列が一致した場合：get_common_longest_subseq(a_idx - 1, b_idx - 1, seq_a, seq_b) + 1 比較対象文字列が一致しなかった場合：max(get_common_longest_subseq(a_idx, b_idx - 1, seq_a, seq_b), get_common_longest_subseq(a_idx - 1, b_idx, seq_a, seq_b))
上記をコードに落とすと，下記のようになります．
def longest_common_subsequences(seq1, seq2):
return longest_common_subsequences_int(len(seq1)-1, len(seq2)-1, seq1, seq2)
def longest_common_subsequences_int(i, j, seq1, seq2):
if (i == -1 or j == -1):
return 0
if (seq1[i] == seq2[j]):
val = longest_common_subsequences_int(i-1, j-1, seq1, seq2) + 1
else:
val1 = longest_common_subsequences_int(i-1, j, seq1, seq2)
val2 = longest_common_subsequences_int(i, j-1, seq1, seq2)
val = max(val1, val2)
return val
if __name__ == '__main__':
seq1 = "abcdbab"
seq2 = "bdcaba"
print(longest_common_subsequences(seq1, seq2))
２．再帰関数の分岐から漸化式を作る
ということで，ここまでくればあとは機械的に行けそうです．漸化式は，下記のようになります．
length_tbl[0][:] = 0, length_tbl[:][0] = 0 length_tbl[i][j] = length_tbl[i-1][j-1] + 1 (比較文字列が一致したとき) length_tbl[i][j] = max(length_tbl[i-1][j], length_tbl[i][j-1]) (比較文字列が不一致の時)
３．漸化式を使って，ボトムアップで表を埋めていく
あとは２の漸化式を使って，(i, j) に関する二次元表を末端から埋めていくだけです．
def longest_common_subsequences_dp(seq1, seq2):
tbl = np.zeros((len(seq1)+1, len(seq2)+1), dtype=int)
tbl.fill(-sys.maxsize)
tbl[:, 0] = 0
tbl[0, :] = 0
for i in range(1, len(seq1)+1):
for j in range(1, len(seq2)+1):
if (seq1[i-1] == seq2[j-1]):
tbl[i,j] = tbl[i-1, j-1]+1
else:
tbl[i,j] = max(tbl[i-1, j], tbl[i, j-1])
return tbl[len(seq1), len(seq2)]
if __name__ == '__main__':
seq1 = "abcdbab"
seq2 = "bdcaba"
print(longest_common_subsequences_dp(seq1, seq2))
完成したテーブルは下記のようになります．
[[0 0 0 0 0 0 0]
[0 0 0 0 1 1 1]
[0 1 1 1 1 2 2]
[0 1 1 2 2 2 2]
[0 1 2 2 2 2 2]
[0 1 2 2 2 3 3]
[0 1 2 2 3 3 4]
[0 1 2 2 3 4 4]]
動的計画法２ コイン問題
動的計画法シリーズ，２個目の例題です．この問題も結構有名な問題らしく，ググるといろいろ出てきます．今回は再帰を使ってこの問題を解く方法を考えて，そこから動的計画法の表を作ってみたいと思います．
問題定義
額面が C1, C2, ... , Cm のコインを使って，n 円を支払うときの最小のコインの枚数をもとめよ． ただし，各額面のコインは何回使ってもよいものとする．
問題のイメージは下記のとおりです．
１．問題の解き方（統治分割法＆貪欲法）
まずはいつものごとく，貪欲に解いてみます．ナップザック問題とこの問題が違うところは，同一のコインを何枚使ってもいいというところです．この違いは再帰関数内の分岐条件の違いとして現れます．一つの額面を一回しか使えないという制限がついている場合，分岐は次のようになります．
関数を次のように定義します．
get_coin_num(cur_idx, amt, coin_list) cur_idx : 0 ~ cur_idx-1 が支払いに使用できるコインの種類 amt : 支払い金額 coin_list : コインの額面リスト
つまり，get_coin_num(3, 10, coin_list) は，coin_list の 0 ~ 2 までのコインを使って amt の金額を支払う際のコインの最小枚数です．
自分を使った場合の枚数：get_coin_num(cur_idx-1, amt-coin_list[cur_idx], coin_list) + 1 自分を使わなかった場合の枚数：get_coin_num(cur_idx-1, amt, coin_list)
このコインを使おうが，使うまいが，同一 idx のコインを考慮するのは一回だけなので，再帰関数を呼び出すときには cur_idx + 1 としてそれ以降の再帰関数ではこの額面のコインが扱われることはありません．一方で，同一のコインを何回も使っていい場合，分岐を次のように書き替える必要があります．
自分を使った場合の枚数：get_coin_num(cur_idx, amt-coin_list[cur_idx], coin_list) + 1 自分を使わなかった場合の枚数：get_coin_num(cur_idx-1, amt, coin_list)
これをコードに落とすと，下記のようになります．
def get_num_coin_rec(amt, coin_list):
return get_num_coin_rec_int(len(coin_list)-1, amt, coin_list)
def get_num_coin_rec_int(cur_idx, amt, coin_list):
if (amt == 0):
return 0
elif (amt < 0 or cur_idx == -1):
return sys.maxsize
val1 = get_num_coin_rec_int(cur_idx, amt-coin_list[cur_idx], coin_list) + 1
val2 = get_num_coin_rec_int(cur_idx-1, amt, coin_list)
return min(val1, val2)
if __name__ == '__main__':
print("Coin Changing Problem")
amt = 15
coin_list = [1, 2, 7, 8, 12, 50]
print(get_num_coin_rec(amt, coin_list))
２．問題の解き方（統治分割法＆メモ化再帰）
関数を見てわかる通り，コールされるたびに変わる変数は cur_idx と amt です．よって，この 2 変数を添え字としたメモを作成すれば，上記関数をメモ化再帰に買えることができます．
def get_num_coin_rec_memo(amt, coin_list):
memo = np.zeros((len(coin_list)+1, amt+1))
memo.fill(sys.maxsize)
return get_num_coin_rec_memo_int(len(coin_list)-1, amt, coin_list, memo)
def get_num_coin_rec_memo_int(cur_idx, amt, coin_list, memo):
if (amt == 0):
val = 0
elif (cur_idx == -1 or amt < 0):
val = sys.maxsize
elif (memo[cur_idx, amt]!=sys.maxsize):
val = memo[cur_idx, amt]
else:
val1 = get_num_coin_rec_memo_int(cur_idx, amt-coin_list[cur_idx], coin_list, memo) + 1
val2 = get_num_coin_rec_memo_int(cur_idx-1, amt, coin_list, memo)
memo[cur_idx, amt] = min(val1, val2)
val = memo[cur_idx, amt]
return val
if __name__ == '__main__':
print("Coin Changing Problem")
amt = 15
coin_list = [1, 2, 7, 8, 12, 50]
print(get_num_coin_rec_memo(amt, coin_list))
３．問題の解き方（動的計画法）
１と２のコードの分岐のさせ方から，下記のように漸化式が求まります．
table[0][j] = INFINITY (０円の額面のコインを使うと，無限枚必要．) table[i][j] = table[i-1][j] (if j < amt[i]) table[i][j] = min(table[i][j-amt[i]]+1, table[i-1][j]) (if j < amt[i])
漸化式が求まっていったので，結局は i, j が小さいところから表を埋めていけば，DPでこの問題を解けます．実装は，下記のようになります．
def get_num_coin_dp(amt, coin_list):
cnt_tbl = np.zeros((len(coin_list)+1, amt+1), dtype=int)
cnt_tbl.fill(255)
cnt_tbl[:, 0] = 0
for i in range(1, len(coin_list)+1):
for j in range(1, amt+1):
if (j < coin_list[i-1]):
cnt_tbl[i, j] = cnt_tbl[i-1, j]
else:
cnt_tbl[i, j] = min(cnt_tbl[i-1, j], cnt_tbl[i, j-coin_list[i-1]] + 1)
return cnt_tbl[len(coin_list), amt]
if __name__ == '__main__':
print("Coin Changing Problem")
amt = 15
coin_list = [1, 2, 7, 8, 12, 50]
print(get_num_coin_dp(amt, coin_list))上記の計算でも止まるテーブルを可視化すると，下記のようになります．
[[ 0 INF INF INF INF INF INF INF INF INF INF INF INF INF INF INF] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15] [ 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8] [ 0 1 1 2 2 3 3 1 2 2 3 3 4 4 2 3] [ 0 1 1 2 2 3 3 1 1 2 2 3 3 4 2 2] [ 0 1 1 2 2 3 3 1 1 2 2 3 1 2 2 2] [ 0 1 1 2 2 3 3 1 1 2 2 3 1 2 2 2]]
だいぶ動的計画法の考え方に慣れてきましたが，問題を解くコツは
１．再帰的に求める方法を考える．
２．再帰関数の分岐から漸化式を作る．
３．漸化式を使って，ボトムアップで表を埋めていく
という感じですかね．
動的計画法１ ナップザック問題
引き続き，動的計画法です．今回はフィボナッチ数の問題よりも少し難しくなりますが，動的計画法のありがたさがよくわかる例かと思います．取り上げる問題は，有名な "ナップザック問題" です．
問題は，
重さが wi, 価値が vi であるような n 個の品物がある． これらの品物から，重さの総和が W を超えないように選んだ時の，価値の総和の最大値を求めなさい．
わかりやすく言うと，手持ちのナップザックが重さWまで入れてOKだとしたときに，入れたものの価値が最大になるような品物の組み合わせを見つけることです．
この問題を図示すると，下記のようになります．（とっても汚い図で恐縮です(ﾟДﾟ)ﾉ）
上記の図では，ナップザックに入れていい重さが10の時に，与えられた4つの品物（□，〇，△，☆）を選択して価値を高めるとすると，（□，〇，△）の組み合わせを選んでやれば重さが10，価値が14となって価値の総和が最大になります．
１．問題の解き方（統治分割法＆貪欲法）
結局はそれぞれの品物をナップザックに”入れる”か”入れない”かの選択になるので，品物がn個あった場合はそれぞれの品物について”入れる”か”入れない”を考えていくと，組み合わせ数は2nになります．
上図では，ナップザックに"入れる"を1に，"入れない"を0として記載しましたが，全部で24=16通りのパターンがあることがわかると思います．この方法で解く場合，下記のような手順になります．
１．組み合わせをもれなく洗い出す．
２．１で出てきた組み合わせのそれぞれに対して，重さと価値を計算する．
３．計算した組み合わせの中から，重さが１０以下で価値が最大となるものを選択する．
この手順を実装すると，下記のようになります．
class Product(object):
def __init__(self, val, weight):
self.value = val
self.weight = weight
comb_cnt = 0
def find_maximum_knapsack_greedy(max_weight, product_list):
return find_maximum_knapsack_greedy_int(0, 0, max_weight, product_list)
def find_maximum_knapsack_greedy_int(cur_idx, cur_weight, max_weight, plist):
if (cur_idx >= len(product_list)):
global comb_cnt
comb_cnt = comb_cnt + 1
if (cur_weight > max_weight):
return -sys.maxsize
return 0
return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)
if __name__ == '__main__':
print("Knapsack Problem")
product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
max_weight = 5
print(find_maximum_knapsack_greedy(max_weight, product_list))
print("Combination Cnt : " + str(comb_cnt))
上の実装では，全探索した時の組み合わせ数を計算したかったので，重さの条件チェックを最後までやっていません．この関数を実行すると組み合わせ数がたしかに16となっていることがわかります．
ただ，重さがナップサックの許容量を超えている時点で計算をやめてもいいはずなので，この条件を付け加えるとコードは下記のようになります．
class Product(object):
def __init__(self, val, weight):
self.value = val
self.weight = weight
def find_maximum_knapsack_greedy_int(cur_idx, cur_weight, max_weight, plist):
if (cur_weight > max_weight):
return -sys.maxsize
elif (cur_idx >= len(product_list)):
global comb_cnt
comb_cnt = comb_cnt + 1
return 0
return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)重さを計算の打ち切り条件にしたことで，組み合わせ数が12となりました．この問題を図示すると，下記のような二分木になります．
結局は，「4つある品物の組み合わせを全部洗い出して，重さ条件を満たす組み合わせのものの中で価値が最大のものを探す」ということをやっていることがわかります．ここで，「品物nを選ぶ，選ばない」という選択が上記のコードのどこで行われているかというと，下記の記述になります．
return max(find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist), \
find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)今考えている品物を選ばない：find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight, max_weight, plist)
今考えている品物を選らぶ ：find_maximum_knapsack_greedy_int(cur_idx+1, cur_weight+plist[cur_idx].weight, max_weight, plist) + plist[cur_idx].value)
となります．この分岐を書くこと自体は動的計画法というよりは統治分割法の話になるのかもしれませんが，この分岐がはっきりすれば動的計画法での実装に一歩近づくものと思われます．
２．問題の解き方（統治分割法＆メモ化再帰）
１の統治分割法＆貪欲法で書いた二分木ですが，もう少し詳しく書いて見ました．それぞれのノードが（i, j）という形で書かれていますが，ここでは下記のように解釈します．（※ここをどう解釈するかが動的計画法のキモです．）
検討すべき品物数がnであるとしたときに，
（i, j）： 許容重量が j [kg] のときの， i+1 ~ n 番目までの品物を組み合わせて作った価値の総和の最大値．
です．下記に具体的な図を描いてみましたが，二分木の上の決定（例えば，品物1を入れる）が許容重量の形でツリーの下側に伝播していく様子がわかります．また，注意深く見ていくと，同じ構造をしている木が複数あることに気づくと思います．
例えば，(2, 3) 「許容重量が 3[kg] の時の，3 ~ 4 番目までの品物を組み合わせて作った価値の総和の最大値」が2回出てきています．この計算自体は何度やっても同じ結果になるので，この結果を”メモ”することで不要な再帰計算を省略することができます．
ここからメモ化の出番です．計算したパスの結果を保存して再利用することで，再帰関数の呼び出しを抑えることが目的になります．
class Product(object):
def __init__(self, val, weight):
self.value = val
self.weight = weight
def find_maximum_kanpsack_rec_memo(max_weight, plist):
memo = np.zeros((len(plist), max_weight+1))
memo.fill(-1)
memo[:, 0] = 0
return find_maximum_kanpsack_rec_memo_int(0, max_weight, plist, memo)
def find_maximum_kanpsack_rec_memo_int(cur_idx, cur_weight, plist, memo):
if (cur_weight < 0):
return -sys.maxsize
elif (cur_idx == len(plist)):
return 0
elif (memo[cur_idx, cur_weight] != -1):
return memo[cur_idx, cur_weight]
inc_value = find_maximum_kanpsack_rec_memo_int(cur_idx+1, cur_weight-plist[cur_idx].weight, plist, memo) + plist[cur_idx].value
exc_value = find_maximum_kanpsack_rec_memo_int(cur_idx+1, cur_weight, plist, memo)
memo[cur_idx, cur_weight] = max(inc_value, exc_value)
return memo[cur_idx, cur_weight]
if __name__ == '__main__':
print("Knapsack Problem")
product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
max_weight = 5
print(find_maximum_kanpsack_rec_memo(max_weight, product_list))
３．問題の解き方（動的計画法）
二分木による可視化とメモ化ができれば，あとはこれを再帰を使わずに for ループで実現するための表を作ることは簡単です．2で可視化したツリーを見るとわかるように，末端以外のノードは自分の下に2つの子供を持っています．当該ノードが持つ値は，下記のどちらか大きいほうの値になります．
１．自分を含めた時の価値の最大値
２．自分を含めなかった時の価値の最大値
(1, 3) のノードを例にして書くと，下記のようになります．
(1, 3) = max((2, 1) + 品物2の価値, (2, 3))
つまり，漸化式は下記のようになります．
val[n+1][j] = 0 val[i][j] = max(val[i+1][j]) (j < w[i]) : 品物の重量が許容重量より大きく，考慮に入れられない時． val[i][j] = max(val[i+1][j-w[i]]+value[i], val[i+1][j]) : それ以外
といった具合になります．また，ノードの値を求めるには当該ノードの子供の値が確定していないと求まらないので，for ループでテーブルを作成するには末端のノード，つまり4番目のノードの値から求めてあげる必要があります．この場合，テーブルは下記のようになります．
※ちなみに，5番目の品物というものは存在しないので，ダミーみたいなもんです．
これをコードに落とすと．．．
class Product(object):
def __init__(self, val, weight):
self.value = val
self.weight = weight
def find_maximum_knapsack_dp(max_weight, plist):
val_table = np.zeros((len(plist)+1, max_weight+1))
val_table[4, :] = 0
for i in range(len(plist)-1, -1, -1):
for j in range(max_weight + 1):
if (plist[i].weight > j):
val_table[i, j] = val_table[i+1, j]
else:
val1 = val_table[i+1, j-plist[i].weight] + plist[i].value
val2 = val_table[i+1, j]
val_table[i, j] = max(val1, val2)
return val_table[0, max_weight]
if __name__ == '__main__':
print("Knapsack Problem")
product_list = [Product(4, 2), Product(5, 2), Product(2, 1), Product(8, 3)]
max_weight = 5
print(find_maximum_knapsack_dp(max_weight, product_list))ちなみに，上記のコードだとDPの表を求める時の一階層目のループが逆方向になってます．順方向のほうがすっきり感じる場合は（自分は少なくともそうでした），問題の見方を少し変えれば順方向のループでテーブルが作れます．
この二分木では，(i, j) を下記のように解釈しています．
（i, j）： 許容重量が j [kg] のときの， i 番目までの品物を組み合わせて作った価値の総和の最大値．
こっちのほうが自然ですよね？こう解釈すると，表が順方向に計算できます．
これをコードに落とすと，，，
def find_maximum_knapsack_dp_low2high(max_weight, plist):
val_table = np.zeros((len(plist)+1, max_weight+1))
val_table[0, :] = 0
for i in range(1, len(plist)+1):
for j in range(max_weight + 1):
if (plist[i-1].weight > j):
val_table[i, j] = val_table[i-1, j]
else:
val1 = val_table[i-1, j-plist[i-1].weight] + plist[i-1].value
val2 = val_table[i-1, j]
val_table[i, j] = max(val1, val2)
return val_table[len(plist), max_weight]
となります．引き続き，もう少し動的計画法をやっていきます．
動的計画法０ フィボナッチ数
つくばチャレンジ本番まであとわずかですが，，，カメラとLidarのフュージョンをするにあたってせっかくなのでSGMとStixelをちゃんと理解しておこうと思い論文を読んでます．
SGMの論文
https://core.ac.uk/download/pdf/11134866.pdf
Stixelの論文
http://www.lelaps.de/papers/badino_dagm09.pdf
が，，，SGMもStixelも動的計画法がキーなんですが，これがいまいち理解できず．．．このあたり，やっぱり情報系の基礎・基本がなってないんだろうなあと思います．ということで，遠回りではあるんですが，一念発起して今週は動的計画法の勉強をすることにしました．
動的計画法の定義
で，積読してあったアルゴリズムの本を開いて定義から確認です．動的計画法は英語では Dynamic Programming というそうですが，この Programming という単語はプログラミングの意ではなく，【表】を意味するんだそうです．ということで，厳密な意味で言うと動的計画法は
"計算結果を表に保存しておき，必要になったらその結果を利用することで計算の効率化を図る方法．"
といえるのかなと思います．動的計画法のキモは，”【表】をどうやって作るか？”というところかなと理解しました．この定義だけ書いても意味不明だと思うので，最初の例題としてフィボナッチ数列をやってみます．
フィボナッチ数列
フィボナッチ数列はググればいろんな情報が出てきます．ひとまず，Wikipediaのリンクを載せておきます．
フィボナッチ数 - Wikipedia
で，実際のフィボナッチ数列の定義ですが，下記のようになります．
こいつをプログラムで求めるとすると，下記のようになります．
def fibonacci_recursive(n):
print("fibonacci(" + str(n) + ")")
if (n == 0):
return 0
elif (n == 1):
return 1
return fibonacci_recursive(n-2) + fibonacci_recursive(n-1) この再帰関数は，nが0か1になるまでスタックに関数が積まれていきます．試しにn=6として上記関数を実行すると，下記のように関数がコールされます．
fibonacci(6), fibonacci(4), fibonacci(2), fibonacci(0), fibonacci(1), fibonacci(3), fibonacci(1), fibonacci(2), fibonacci(0), fibonacci(1), fibonacci(5), fibonacci(3), fibonacci(1), fibonacci(2), fibonacci(0), fibonacci(1), fibonacci(4), fibonacci(2), fibonacci(0), fibonacci(1),
fibonacci(3), fibonacci(1), fibonacci(2), fibonacci(0), fibonacci(1)
たかだかフィボナッチ数列の6項目を求めるのために，かなりの回数関数が呼ばれていることがわかります．ひとつの関数が2つの関数を再帰的に呼び出すので，図にすると下記のように2分木になります．
で，上図を見てもらうとわかるんですが，一度計算している関数が何回も呼ばれていることがわかります．例えば，fib(4)とか．一度計算した値に関しては保存しておいて再利用すれば，関数を呼び出す回数はかなり削減できそうです．
メモ化再帰を用いたフィボナッチ数列の計算
ということで，一度計算済のフィボナッチ数をメモしながらフィボナッチ数を求めます．
def fibonacci_recurcive_memo(n):
memo = np.zeros(n+1)
memo.fill(-1)
return fibonacci_recurcive_memo_int(n, memo)
def fibonacci_recurcive_memo_int(n, memo):
print("fibonacci(" + str(n) + ")")
if (memo[n] != -1):
return memo[n]
elif (n == 0):
memo[n] = 0
return memo[n]
elif (n == 1):
memo[n] = 1
return memo[n]
else:
memo[n-2] = fibonacci_recurcive_memo_int(n-2, memo)
memo[n-1] = fibonacci_recurcive_memo_int(n-1, memo)
return memo[n-2] + memo[n-1]同じようにn=6として実行すると，下記のように関数が呼ばれていきます．
fibonacci(6), fibonacci(4), fibonacci(2), fibonacci(0), fibonacci(1), fibonacci(3), fibonacci(1), fibonacci(2), fibonacci(5), fibonacci(3), fibonacci(4)
一度計算した値を保存して再利用することで，関数の呼び出し回数がかなり減りました．再帰呼出しの中で途中結果を保存しておく方法を"メモ化再帰"と呼ぶようです．おそらく厳密な意味ではこれは動的計画法ではなく，”分割統治法・再帰法+メモ化再帰”なのだと思います．ただ，メモ化再帰ができるということは，表（＝メモ）を作ることができます．
動的計画法を用いたフィボナッチ数列の計算
やっとこさ本題の動的計画法です．最初に説明したように，動的計画法で問題を解くには，”表”を作る必要があります．この問題の場合は，fib(0), fib(1), fib(2), ...というようにそれぞれのフィボナッチ数を表として作成する必要があります．で，ここでポイントなんですが，フィボナッチ数を求める際に，値が定義されているのは n=0とn=1だけです．あとは漸化式を解いて求める必要があります．つまり，表を計算するときにも，nが小さい方から計算していく必要があります．動的計画法（＝表作成）は，この例のようにボトムから解ける問題に対して応用可能な技法になります．基本的には，問題を漸化式の形に落とし込めることができれば，動的計画法で問題が解ける可能性が大ということかと思います．
初自動走行
本日9/15日，第三回試走会の日でしたが，，，，あいにくの雨だったので参加をあきらめました．ロボットの防水がしっかりできておらず，壊れてしまうとつらいので．．．結局次回の試走会が10/13ということで，だいぶ後になってしまうのですが，確認走行のトライはそこまでおあずけです．試走会はあと五回．．．だんだん追い込まれてきました．ただ，当初の目標だった下記6項目は完了しました．
1. 確認走行コースのデータ取得：完了
2. 確認走行コースのマップ作成：完了
3. AMCLを使っての自己位置推定：完了
4. Yp-Spurを使ってロボット制御実施．制御パラメータ調整：完了
5. Teb-Local Planner と Yp-Spur を使ってロボット制御実施：完了
6. move_base をつかって，センシングしながらのロボット制御を実施：完了
コース（リビング？）の写真です．段ボールの真ん中を目指して，障害物をよけながら走ります．
ちょっとまだ Local Planner が使いこなせてない感じでロボットの軌跡がふらふらしてますが，ここはもう少しチューニングが必要です．ただ，Lidar + ROS の威力はすさまじいですね．ハードを作ればあれよあれよというまに自動走行まで来てしまいました．ここからは，カメラとLidarのFusionに取り組んでいきます！
ROS Teb Local Planner を使った位置制御
ということで第三回試走会前最後の週末，割と頑張りましたが ROS Navigation Stack の理解で思ったより時間がかかってしまい，，，自宅の地図を作ってセンシングしながらの自律走行まで完了させたかったんですが，Teb Local Planner を使ってオドメトリだけを見ながら移動するとこまでで終わってしまいました．金曜までに何とか6まで行ければいいのですが．．．
1. 確認走行コースのデータ取得：完了
2. 確認走行コースのマップ作成：完了
3. AMCLを使っての自己位置推定：完了
4. Yp-Spurを使ってロボット制御実施．制御パラメータ調整：完了
5. Teb-Local Planner と Yp-Spur を使ってロボット制御実施：完了
6. move_base をつかって，センシングしながらのロボット制御を実施：これから
下記，PS3のコントローラでロボットが動くようになりました．
www.youtube.com
下記，Teb Local Planner で指定した４点を順々に移動していきます．差動式のロボットだと位置を変えずに回転できるので，車型のロボットよりだいぶ有利ですね．
www.youtube.com
Teb Local Planner と Yp-Spurの制御がうまくかみ合わず，パラメータ調整にも結構時間がかかりました．部屋の中だけで動かしている状態なので，また外に持っていくといろいろと調整をしなくてはいけないんだろうなと思います．
