蟻本をpythonで: Sec. 2-3, 動的計画法(DP)

セクションごとにして編集で追加していくスタイルに変更。

page 52: ナップサック問題

# -*- coding: utf-8 -*-

n = int(raw_input())
w = map(int,raw_input().split())
v =  map(int,raw_input().split())
W =  int(raw_input())
################ naiveな全探索 ########################

def val_tot(i, j):
    if i == n:
        res = 0
    elif j < w[i]:
        res = val_tot(i+1,j)
    else:
        res = max(val_tot(i+1,j), val_tot(i+1, j-w[i]) + v[i])

    return res

print val_tot(0, W)

################ メモ化 ########################
import numpy as np
val_tot_list = -1*np.ones((n+1,W+1))

def val_tot_memo(i, j):
    if val_tot_list[i,j] !=-1:
        return val_tot_list[i,j]
    else:
        if i == n:
            res = 0
        elif j < w[i]:
            res = val_tot_memo(i+1,j)
        else:
            res = max(val_tot_memo(i+1,j), val_tot_memo(i+1, j-w[i]) + v[i])

        val_tot_list[i,j] = res

        return res

print int(val_tot_memo(0, W))

################ 漸化式 (DP) ########################

## dp[i][j] = dp[i+1][j] ( w[i] > j)
##          = max(dp[i+1][j], dp[i+1][j-w[i]] + v[i])
## val_tot_list[n,W] = 0

# import numpy as np
val_tot_list = -1*np.ones((n+1,W+1))

val_tot_list[n,:] = 0*val_tot_list[n,:]

for i in range(n-1,-1,-1):
    for j in range(W+1):
        if i ==n:
            val_tot_list[n,j] = 0
        else:
            if w[i]  > j:
                val_tot_list[i][j] = val_tot_list[i+1][j]
            else:
                val_tot_list[i][j] = max(val_tot_list[i+1][j], val_tot_list[i+1][j - w[i]] + v[i])

print int(val_tot_list[0][W])

計算量はそれぞれ O(2^n), O(nW), O(nW).
いきなり漸化式を立ててコードに持っていくのはまだ慣れないので(添字間違えそう)、ひとまずnaive全探索を書いてからメモ化するのを身に着けたい。

#page 56: 最長共通部分問題(LCS)

# -*- coding: utf-8 -*-
s = raw_input()
t = raw_input()

#### dp[i][j]: sのi番目以前、tのj番目以前の文字列についての最大共通部分列の長さ


########################## naive #########################
def lcs(i,j):
    if i == 0 or j == 0:
        return 0
    else:
        if s[i-1] == t[j-1]:
            return lcs(i-1,j-1) + 1
        else:
            return max(lcs(i-1,j-1+1), lcs(i-1+1,j-1))


print lcs(len(s),len(t))


##########################メモ化##################################
import numpy as np
lcs_list = -1*np.ones((len(s)+1, len(t)+1),dtype = int)

def lcs_memo(i,j):
    if lcs_list[i][j] != -1:
        return lcs_list[i][j]
    else:
        if i == 0 or j == 0:
            lcs_list[i][j] = 0
            return 0
        else:
            if s[i-1] == t[j-1]:
                lcs_list[i][j] =  lcs_memo(i-1,j-1)+1
            else:
                lcs_list[i][j] = max([lcs_memo(i-1,j-1+1), lcs_memo(i-1+1,j-1)])
            return lcs_list[i][j]


print lcs_memo(len(s),len(t))



########################## DP ##################################
import numpy as np
lcs_list = -1*np.ones((len(s)+1, len(t)+1),dtype = int)

lcs_list[0,:] = 0*lcs_list[0,:]
lcs_list[:,0] = 0*lcs_list[:,0]   # initialize

for i in range(len(s)):
    for j in range(len(t)):
        if s[i] == t[j]:
            lcs_list[i+1,j+1] = lcs_list[i,j]+1
        else:
            lcs_list[i+1,j+1] = np.max([lcs_list[i+1,j], lcs_list[i,j+1] ])

print lcs_list[-1,-1]

「ひとまずnaive全探索を書いてからメモ化するのを身に着けたい。」と言ったものの、この問題に関してはDPが一番すんなり書けた。ナップサック問題の場合と違ってi, j のループが順方向だったので初期化が分かりやすかったのが一因っぽい。

他二つはlcsの引数と配列のindexingの対応づけでちょっと混乱した。

page58: 個数制限なしナップサック問題

# -*- coding: utf-8 -*-
n = int(raw_input())
w = map(int,raw_input().split())
v = map(int,raw_input().split())
W = int(raw_input())

 #### dp[i][j] : i番目までの品物の組み合わせで重さj以下にしたときの最大価値, 最終的な計算結果は d[n][W]
 #### 初期値 d[0][:] = 0, dp[:][0] = 0 (w_i >= 1なので)

 #### 漸化式は dp[i+1][j] = max(dp[i][j-k*w[i]] + k*v[i]| k ... ith　品物の個数)


########################## 全探索   ############################
def dp(i,j):
    if i <= 0 or j <= 0:
        return 0
    else:
        return max([dp(i-1,j-k*w[i-1])+ k*v[i-1] for k in xrange(j//w[i-1] + 1)])

print dp(n, W)

##########################  メモ化   ############################
import numpy as np

dplist = -1*np.ones((n+1, W+1),dtype=int)


def dp_memo(i,j):

    if dplist[i,j] != -1:
        return dplist[i,j]
    else:
        if i <= 0 or j <= 0:
            dplist[i,j] =  0
            return 0
        else:
            dplist[i,j] =  max([dp_memo(i-1,j-k*w[i-1])+ k*v[i-1] for k in xrange(j//w[i-1] + 1)])
            return dplist[i,j]

print dp_memo(n,W)

############################  DP_1   ##########################

import numpy as np

dplist = -1*np.ones((n+1, W+1),dtype=int)

dplist[0,:] = 0
dplist[:,0] = 0

for i in xrange(n):
    for j in xrange(W+1):
        dplist[i+1,j] = max([dplist[i,j-k*w[i]]+k*v[i] for k in xrange(j//w[i]+1)])

print dplist[n][W]

############################  DP_2   ##########################

# import numpy as np

dplist = -1*np.ones((n+1, W+1),dtype=int)

dplist[0,:] = 0
dplist[:,0] = 0

for i in xrange(n):
    for j in xrange(W+1):
        if j < w[i]:   #### i+1番目の品物はもう入らない
            dplist[i+1,j] = dplist[i,j]
        else:
            dplist[i+1,j] = max(dplist[i,j], dplist[i+1,j-w[i]] + v[i])  #### 漸化式の変形(page 59). ループが減らせる.

print dplist[n][W]

蟻本にもあるように、漸化式からそのままDPを書くとループが三重になってしまう。漸化式の式変形で重複した計算を回避すれば2重ループにできる。

メモ化してもループしていることに変わりはないのでいずれにせよ処理は必要っぽい。
配列の扱いが楽だからという理由だけでnumpyを使ってしまっているがビルトインのもろもろだけで書けるようになっておいたほうが良い気がする。

#page60: ナップサック問題その2

# -*- coding: utf-8 -*-

n = int(raw_input())
w = map(int,raw_input().split())
v =  map(int,raw_input().split())
W =  int(raw_input())

######### 重さの制約がきつい場合に対応して、DPの対象を価値から重さにした場合 ##########

########## dp[i + 1][j]: i番目の品物までで価格jにした時の最小の重さ
###########　最終的な出力は、 dp[n][j] <= W を満たす最大のj
########### 漸化式は dp[i+1][j] = min(dp[i][j], dp[i][j-v[i]]+w[i])
###########  初期値は dp[0][0] = 0, dp[0][j] = INF =float("INF")


###################  全探索 ################
def dp(i,j):
    if i == 0 and j == 0:
        return 0
    elif j != 0 and i == 0:
        return float("INF")
    else:
        if j < v[i-1]:
            return dp(i-1,j)
        else:
            return min(dp(i-1,j), dp(i-1,j-v[i-1])+w[i-1])

res = 0
for j_val in xrange(sum(v)):
    if dp(n,j_val)  <= W:
        res = j_val
print res


################### メモ化 ################

import numpy as np

dplist = np.ones((n+1, sum(v)+1)) * (-1)

def dp_memo(i,j):

    if i == 0 and j == 0:
        dplist[i,j] = 0
        return 0
    elif j != 0 and i == 0:
        return float("INF")
    elif dplist[i,j] != -1:
        return dplist[i,j]
    else:
        if j < v[i-1]:
            dplist[i,j] = dp_memo(i-1,j)
            return dplist[i,j]
        else:
            dplist[i,j] =  min(dp_memo(i-1,j), dp_memo(i-1,j-v[i-1])+w[i-1])
            return dplist[i,j]

res = 0
for j_val in xrange(sum(v)):
    if dp_memo(n,j_val)  <= W:
        res = j_val
print res



################### DP ##################

import numpy as np

dplist = np.ones((n+1, sum(v)+1)) * (-1)

dplist[0,:] = float("INF")
dplist[0,0] = 0

for i in xrange(n):
    for j in xrange(sum(v) +1):
        if j < v[i]:
            dplist[i+1,j] =dplist[i,j]
        else:
            dplist[i+1,j] = min(dplist[i,j], dplist[i,j-v[i]]+w[i])

res = 0
for j_val in xrange(sum(v)):
    if dplist[n,j_val]  <= W:
        res = j_val
print res

これもDPが一番書きやすい。

漸化式が dp[i] = hogehoge(dp[i+1])かdp[i+1] = hogehoge(dp[i])のどちらか。

再帰を使って書く場合は前者の方が楽で、DPの場合は後者が楽 (配列や関数の引数で混乱が少ない)。どちらで書くかに応じて立てる漸化式を変えたほうが良さそう。
今の問題にしても、dp[i,j]を「i番目以降の品物で〜」、とすれば
漸化式がdp[i,j] = min(dp[i+1,j], dp[i+1,j-v[i]]+w[i])となるので再帰では書きやすい。

#page62: 個数制限付き部分和問題

# -*- coding: utf-8 -*-


n = int(raw_input())
a = map(int,raw_input().split())
m = map(int,raw_input().split())
K = int(raw_input())

import numpy as np


################ dp[i+1][j]: i番目まででjを作るときに余る最大のi番目の個数
################ 漸化式は dp[i+1][j] = dp[i+1][j-a[i]] -1
################                    = -1
################                    = m[i]      (dp[i][j] >= 0)

################ 作れない場合は-1を入れておく

dplist = np.ones((len(a)+1,K+1))*(-1)

dplist[0,0] = 0

for i in range(len(a)):
    for j in range(K+1):
        if dplist[i,j] >= 0:
            dplist[i+1,j] = m[i]
        elif dplist[i+1,j-a[i]]<= 0:
            dplist[i+1,j] = -1
        elif j - a[i] < 0:
            dplist[i+1,j] = -1         #######  j  - a[i] <0のとき、i番目を入れてもjが作れるかは変わらない。(i-1)番目までで作れる場合は処理されているので、ここでは作れない場合のみ考えれば良い。

        else:
            dplist[i+1][j] = dplist[i+1,j-a[i]] -1

print dplist[n,K] >= 0

boolを求めるDPは無駄が多いらしい。

page 62のbool値DPでは、a[i]をいくつ入れるかを処理するループがあるため三重ループになってしまっている(O(K \sum_i m[i])と書いてあるけどO(nK \sum_i m[i])が正しそう)。

個数制限なしナップサック問題のときもそうだったが、三重ループが出てくるようならどこかに無駄があると考えたほうが良さそう。a[i]の個数についてのループを消すためにはどうすればよいかと考える。

dp[i+1][j]を、i番目まででjが作れるか、としたboolのDPでは、
「dp[i+1, j-k*a[i]] == Trueなるkが存在するか」という処理が三重ループの原因になっている。なので、dp[i+1, j-a[i]]の結果を使ってdp[i+1][j]が判定できるようになっていれば三重ループは必要無くなる。

しかし、bool値だけを持っていると、dp[i+1, j-a[i]] = Trueが与えられても「a[i]がまだ残っているのか」が分からないのでdp[i+1, j]について判断できずにループを回すハメになる。だったら配列として残りのa[i]の個数を持っておけばよい、というのが考え方っぽい。

流石に問題文から一足飛びにここまでたどり着くにはまだ修行が足りない。上のように、素直な(愚直な)方法で書いてみて(考えてみて)改善点を探していくのが筋だろうか(基本的な問題だからパターンを覚えろというだけかもしれない)。

#page 63: 最長部分増加列問題

# -*- coding: utf-8 -*-
import numpy as np

a = map(int,raw_input().split())


#####################page 63#########################

### dp[i]: a[i]が最後尾にいる最長部分列の長さ

dp = np.ones(len(a),dtype = int) * (-1)

dp[0] = 0
res = 0
for i in range(len(a)):
    dp[i] = 1
    for j in range(i):
        if a[j] < a[i]:
            dp[i] = max(dp[i],dp[j] + 1)

    res = max(res,dp[i])

print res


#####################page 64: O(n^2)#########################

### dp[i]: 長さがi+1の部分列の最終要素の最小値

dp = np.ones(len(a),dtype = int) * float("INF")
n = len(a)
for j in xrange(n):
    for i in xrange(n):
        if i == 0 or dp[i-1] <a[j]:
            dp[i] = min(dp[i],a[j])
print len(dp[dp<float("INF")])


#####################page 65: O(n log n)#########################

import bisect   #### C++でのlower_boundと同等の動作を実現するために。
dp = np.ones(len(a),dtype = int) * float("INF")
n = len(a)

for j in xrange(n):
    dp[bisect.bisect_left(dp,a[j])] = a[j]
print len(dp[dp<float("INF")])

下二つが難しい。

とりあえずpage 64の方針の場合、jについてループを回すと以下のように配列が埋まっていく。

a = {2,1,3,2}

1. dp[0] = 2, dp[1] = INF
2. dp[0] = min(2,1) = 1, dp[1] = INF
3. dp[0] = min(1,3) = 1, dp[1] = min(INF,3) = 3, dp[2] = INF
4. dp[0] = min(1,2) = 1, dp[1] = min(3,2) = 2, dp[2] = INF

aの末尾に新しく値を追加したときにdpテーブルがどう更新されるべきか考える。

まず、dp[0]は長さ1の増加部分列なので、新しいaの最小要素を持ってくればよい。

dp[1]は、長さ2の増加部分列の最後尾の最小要素。新しい要素の追加で値が更新されるかどうかは、追加前のaでのdp[1]と新しい要素の小さい方で決めれば良い。ただし、dp[0]が新しい要素で更新されてしまった場合、新しい要素は増加部分列の最後尾になりえないので排除する (if i ==0 or dp[i-1] < a[j]の部分)。

dp[2]について、dp[1]が新しい要素で更新されたとする。このとき、dp[2]もまた更新されてしまうとすると、それより小さい値で長さ2の部分列の最後尾に来るものがいることになるので矛盾。なので、dp[1]が更新されていない場合に限りもとのdp[1]と新しい要素を比較する。

というわけで、一般に新しい要素をaに追加した場合、dp[i]は、dp[i-1]がその値で更新されていないときにだけmin(dp[i], a_new)で更新すればよい。

初期値をINFで与えておけば、aの要素についてループを回し、dpの頭から上記のように判定していけば良い。これがO(n^2)での解答。

page 65の二分探索を用いた方法は、dpの頭から見て判定していく部分をO(log n)にしたもの。最終的なdp配列は値が単調増加になるので、a[j]未満のdp[i]については見る必要がない。よって、bisect.bisect_left (STLではlower_bound)によってdp[i]>=a[j]となるようなiを探せばよい。探し方から、dp[i-1] < a[j], dp[i] = min(dp[i],a[j]) = a[j]が自動的に満たされているので、条件判定などせずにdp[i] = a[j]で更新することができる。

#page66: 分割数

# -*- coding: utf-8 -*-
import numpy as np


n = int(raw_input())
m = int(raw_input())
M = int(raw_input())


####### dp[i][j]を、i個をj個に分割するときの方法の数とする。

####### 漸化式はdp[i][j] = dp[i-j][j] + dp[i][j-1]

dp = np.zeros((n+1,m+1),dtype = int)

dp[:,1] = 1

for i in xrange(n+1):
    for j in xrange(2,m+1):
        if i-j >=0:
            dp[i][j] = dp[i][j-1] + dp[i-j][j]
        else:
            dp[i][j] = dp[i][j-1]

print dp[-1][-1]

一番個数の多い箱を抜き出した残りについて考えて、という感じの漸化式を最初考えたものの、重複があってダメだった。上手いことやれば重複を排除できたかもしれないが、どのみちO(mn^2)になってたのでよろしくない。(プログラミング以前に数学力が足りない)

漸化式からコードに落とし込むのは慣れてきた。

#page 67: 重複組合せ

# -*- coding: utf-8 -*-
import numpy as np

n = int(raw_input())
m = int(raw_input())
a = map(int,raw_input().split())
M = int(raw_input())

dp = np.zeros((n+1,m+1),dtype = int)

dp[0,:a[0]] = 1
dp[:,0] = 1

for i in xrange(n):
    for j in xrange(1,m+1):
        if j -1 - a[i] >=0:
            dp[i+1][j] = dp[i+1][j-1] + dp[i][j] -dp[i][j-1-a[i]]
        else:
            dp[i+1][j] = dp[i+1][j-1] + dp[i][j]

print dp[-1][-1]

漸化式は単純だが、そこからどうやって計算量を落とすか。
三重ループを消すために、すでに計算してある量で和を書き直す作業が必要。

#Sec.2-3を終えて

実装力というか、漸化式からDPでコードを書く部分はスムーズにできるようになってきた。
添字の割り当て間違いに対してのカンが働くようになってきた感じ。

適切で計算量の小さい漸化式を立てる部分はまだ練習が必要。思えば受験のころから組合せ論とかはちょっと苦手だった。

Chapter 2の最後に練習問題としてPOJの問題番号が載ってるのでそちらもそのうちやっていきたい。