TSP DP（その２） bit DP / 巡回セールスマン問題 / 動的計画法

前回のTSP DP（動的計画法による巡回セールスマン問題）では、訪問先を表現するためにset()関数を用いましたが、今回はbit DPで試してみました。
bit演算させたほうが高速になるのと、メモ化ライブラリであるfunctools.lru_cacheも使えるようになります。

ビット演算でどこを訪れていて、どこを訪れていないか（訪問フラグ）は以下のように表すようです。
まず訪問先を0, 1, 2, 3とすれば4桁の二進数になります。Pythonでは左側に0bをつければ二進数となります。

S = 0b0001

0ビット目（一番右側）を1にしておいて、ここが出発点となります。
つぎに2に訪れたら（右から0番目、1番目、2番目）、

S = 0b0101

となり、すでに訪れた場所の状態がわかります。
次に1に訪れて、

S = 0b0111

つぎに3に訪れて、

S = 0b1111

最終的にすべてが1になれば完了というわけです。
ビット演算用の記号を使いますが、実際は整数としても認識するので、

S = 1 | 1 << 2 = 0b0101

こんな感じになります。

ということから二進数を用いて前回のTSPを書き直してみました。


functools.lru_cacheを使った場合：（Python 3.6、Jupyter Notebook使用）
メモ化ライブラリであるfunctools.lru_cacheを使うと以下のようになります。
ちなみにこれは合計最短経路のみの出力です。尚、TSP/DPのアルゴリズムについては前回の方が詳しく説明してあるかと。

from functools import lru_cache

@lru_cache(maxsize=None)
def TSP_lru(S, pos):        
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_lru(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# @lru_cache(maxsize=None)
# def TSP_lru(S, pos):  # same as the above
#     if S == all_visited:
#         return dist[pos][0]
#
#     d_min = float('inf')
#     for i in range(num):
#         if S & (1 << i) == 0:     # if not visited 
#             d = dist[pos][i] + TSP_lru(S | (1 << i), i)
#             if d < d_min:
#                 d_min = d
#   return d_min

# Solve by TSP with functools.lru_cache
%time print('Total Distance:', TSP_lru(1, 0))
print(TSP_lru.cache_info())  # show cache info
TSP_lru.cache_clear()        # clear cache otherwise same outputs as before

メモ化されていないファンクションの上に@マークのデコレータをつけるだけでメモ化されます。
maxsize=Noneでメモリ制限なしにしています。デコレータなしだとメモ化しないためかなり遅くなります。
ただし、前回のようなset関数やリストを引数にしている場合はこのlru_cacheが使えません。タプルを引数にするなら大丈夫だと思います。今回は訪問先の集合をビットで表したので使うことができるというわけです。おそらくビット演算にするだけでも、少しは高速になっているのではないでしょうか。
それから、一度runするとcacheに内容が保持されたままになるので、再度TSR_lru(1, 0)をrunすると、前回と同じ答えを出してしまいます。ノード数を変えて試してみたい場合は、毎回cache_clear()する必要があります。何度も同じ答えを呼び出すには高速ですが、キャッシュクリアしないかぎりは同じ答えを出し続けるので、他のプログラムで使うときには注意したほうがいいと思います。


memoize関数を使った場合：
自前で書いたmemoize関数をデコレータにすることもできるようです（こちらを参考）。
合計最短経路のみの出力。

def memoize(f):
    cache = {}
    def func(*args):
        if not args in cache:
            cache[args] = f(*args)
        return cache[args]
    return func

@memoize
def TSP_memo(S, pos):  # re-run this function to clear cache otherwise same outputs as before
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_memo(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# Solve by TSP with memoization decorator
%time print('Total Distance:', TSP_memo(1, 0))

これも先ほどのlru_cacheと同じようなものです。
この場合もキャッシュクリアしない限りは同じ答えを出すので、TSP_memo(1, 0)をrunさせる場合は、def TSP_memo():も一緒に（jupyter上のこのセルごと）runさせたほうがいいです（このキャッシュクリアの方法がわからない）。


リストを使ったDPメモ化：
ビット化した以外はほぼ前回と同じです。
最後にメモ化した内容から経路順を割り出します。

def TSP_BDP(S, pos):
    # if visited before
    if dp[S][pos] >= 0:
        return dp[S][pos]
    
    # if all visited  
    if S == all_visited:         
        dp[S][pos] = dist[pos][0]
        return dp[S][pos]
    
    # if not visited yet
    dp[S][pos] = min([dist[pos][i] + TSP_BDP(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])
    
    return dp[S][pos]

# Solve TSP by bit DP (15 nodes: 0.75 sec, 18 nodes: 8.42 sec, 20 nodes: 44 sec)
dp = [[-1] * num for _ in range(1 << num)]    # initialize memoization
%time print('Total Distance:', TSP_BDP(1, 0))

基本的に、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + TSP_BDP(次の訪問状態, 移動先)

を毎回記録/計算しています。forループのi（移動先）は、未訪問の地点があるかどうか確かめてから、あればそこへ移動する場合の距離を計算しmin()によって最短経路を選択します。
具体的な値を入れてみると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1011 | (1 << 2), 2)

Sについて2ビット目が0となっているので（右から数えて0,1,2番目）、まだ訪れていない地点は2ビット目だけということになり、forループでi=2のときこの処理がなされます。
そうすると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

右辺においてS=0b1111（全て訪問済み）となり、あとは現在地3から2までの距離と残りTSP_BDP(0b1111, 2)による距離を足して左辺のdp[0b1011][3] へ記録することになります。
しかし、dp[0b1011][3] へ記録する前に、TSP_BDP(0b1111, 2)を計算しなければいけないので、TSP_BDP(0b1111, 2)は次の再帰ループに突入してif S == all_visited:の条件分岐により、

dp[0b1111][2] = dist[2][0]

になり距離テーブルdist[2][0]から値を得てdp[0b1111][2]に記録します。このとき再帰ループで一つ先に移動しているので現在地pos=2、移動先は0になっています。距離が記録されたのちdp[0b1111][2]（= dist[2][0]）がリターンされたのでTSP_BDP(0b1111, 2)に戻ります。
この結果、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

この式の右辺側にあるTSP_BDP(0b1111, 2)は、

TSP_BDP(0b1111, 2) = dist[2][0]

とわかったので、それを代入すると

dp[0b1011][3] = dist[3][2] + dist[2][0]

つまり、

未訪問2で現在地3 = 3から2までの距離 + 2から0（ゴール）までの距離

となってゴールから二歩手前までの距離が求まります。

訪問フラグと距離の関係（4ビットの場合、posはその時の現在地）：
dp[ob1111][pos]（4ビット1の場合）ゴールから一歩手前までの距離
dp[0b1011][pos]（3ビット1の場合）ゴールから二歩手前までの距離
dp[0b0011][pos]（2ビット1の場合）ゴールから三歩手前までの距離
dp[0b0001][pos]（1ビット1の場合）ゴールからスタートまでの距離

ということから、この仕組みを理解して経路の生成を以下で行いました。


経路の生成：
上で記録した二次元配列dpを元に経路の順番を生成します。

def gen_dp_tour():
    tour = [0]
    S = 1
    while len(tour) < num:
        d_min = np.inf
        pos = None
        for i in range(num):
            if S & (1 << i) == 0:
                d = dist[tour[-1]][i] + dp[S | (1 << i)][i]
                if d < d_min:
                    d_min = d
                    pos = i
        S |= (1 << pos)
        tour.append(pos)
    return tour + [0]

# Generate Tour
Tour = gen_dp_tour()
print('Tour:', Tour)
plot_path(Tour)

基本的には最短経路を求めるコードと似ています。前のコードのforループ内のTSP_BDP()の部分をdp[][]に置き換えて距離を計算させるかわりに直接メモリから呼び出します。再帰ループの代わりにwhileループにしています。最初出発点にいる状態から全てを訪問するまで、min関数で得た最短距離と同時に現在地も記録しておき、それをtour配列に入れていきます。
先ほどの訪問フラグと距離の関係から、dp[0b0001]にはスタートからゴールまでの合計最短距離が記録されています。これを一つずつ分解していく感じです。

dp[0b0001][pos] = dist[pos][i] + dp[S | (s << i)][i]

というのは、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + dp[次の訪問状態][移動先]

なので、初期状態は0ビット目にいるため0b0001になり現在地pos=0になります。次の移動先iは複数の可能性があるのでforループで各パターンを確かめます。min()関数を使ってdist[][] + dp[][]の合計距離で最短の地点iを次の移動先にしますが、この時最短のiを経路用配列tourにappendしていきます。次のループでは再度訪問フラグをチェックして、複数ある未訪問の地点iについてまたdist[][] + dp[][]の合計距離を調べます。あとはそれの繰り返しとなります。尚、現在地posはその都度追加していった経路用配列tourの最後の値tour[-1]を参照すればいいことになります。初期状態では0から出発することになるので、tour=[0]にしています。

次へ：TSP LP (Linear Programming)：線形計画法による巡回セールスマン問題 / Subtour Elimination 部分巡回路除去

関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）